- Email: [email protected]
Estimating the direction of an unknown air pollution source using a digital elevation model and a sample of deposition
Ecological Modelling 124 (1999) 85 – 95 www.elsevier.com/locate/ecolmodel
Estimating the direction of an unknown air pollution source using a digital elevation model and a sample of deposition Oleg Antonic´ *, Tarzan Legovic´ Rudjer Bosˇko6ic´ Institute, P.O. Box 1016, Bijenicˇka 54, HR-10001 Zagreb, Croatia Received 7 August 1998; received in revised form 5 May 1999; accepted 21 June 1999
Abstract A new method is described for estimating the direction of an unknown air pollution source. A heterogeneous topographic exposure to the wind on the polluted surface is assumed. Hence, the method is applicable to mountainous or hilly areas. Estimators of the topographic exposure to wind are modelled using a digital elevation model for different hypothetical directions of wind using an azimuthal step. These estimators are correlated to soil pollution but any other relevant field sampling can be used. If significant, maximum positive correlation indicates the direction towards pollution source(s). When spatial sampling intensity is low, estimators of topographic exposure can be used to significantly improve spatial interpolation of measured values. The first application of the method were performed on a sample of heavy metals soil concentrations in the Risnjak National Park in Croatia. Maximum correlations between soil pollution and topographic exposure are obtained toward the west. This is the main direction of incoming cyclones and also the direction toward major industrial pollution sources. © 1999 Elsevier Science B.V. All rights reserved. Keywords: Exposure; GIS; Heavy metals; Landscape ecotoxicology; Soil; Spatial interpolation; Topography; Wind flux
1. Introduction A study of environmental risk assessment, ecotoxicological or ecosystem health problem, is usually based on sampling of response variables. They have to be biologically relevant, sensitive, easily measured and appropriately scaled (Cairns et al., 1994). When such an environmental problem exists on a landscape, spatially explicit approach and modelling procedures are used (e.g. Stein et al. * Corresponding author. Fax: +385-1-4680-084. E-mail address: [email protected] (O. Antonic´)
1995; Cairns and Niederlehner, 1996). There are two basic reasons to include models in the landscape ecotoxicology research: (1) to explain spatial variability of measured responses (using appropriate sample); and (2) to predict responses at locations or times when measurements have not be made, using spatial or temporal interpolation, extrapolation or a process model.
1.1. Aim In this paper the following case is examined: the pollution source is not a priori known
0304-3800/99/$ - see front matter © 1999 Elsevier Science B.V. All rights reserved. PII: S 0 3 0 4 - 3 8 0 0 ( 9 9 ) 0 0 1 4 9 - 0
86
O. Antonic´, T. Lego6ic´ / Ecological Modelling 124 (1999) 85–95
and the field sample of responses is too small for the usual geostatistical approaches (e.g. interpola-
Fig. 2. Estimators of topographic exposure to wind. N indicates north. Terrain in the point T has the maximum slope m and terrain aspect g. The first estimator s is relative terrain aspect at point T for a given azimuth of the wind flux (d). The third estimator b is the horizon angle of the point T for a given d, determined by the point M under the chosen search distance d. Note that a change of d may change b. a is the maximum angle between a plane orthogonal to the wind and regression plane trough the terrain at the point T (point T’ is a projection of point T on the orthogonal plane). Cosine of a is the second estimator if b or search distance are set to zero (terrain exposure to the horizontal component of wind flux), otherwise it is the fourth estimator (terrain exposure to the sloped wind flux).
tion by kriging). This situation occurs often in pollution and ecological monitoring. The question is: how can we explain spatial variability of measured responses and interpolate measured values over the research area under the already mentioned constraints?
1.2. Basic hypotheses Fig. 1. Upper part: location of the research area. Lower part: elevation map derived from DEM used in the case study. Values of altitude above sea level range from black (400 m) to white (1526 m). The boundary of the Risnjak National park and the measurement stations for heavy metal soil concentration are superimposed. Labeling of stations follows original data set (Vrbek and Gasˇparac, 1992).
When the research area is uneven, its spatial units are differently exposed to the air pollution that is coming from some distant pollution source. Direction of the source can be estimated by maximizing correlation between the sampled
O. Antonic´, T. Lego6ic´ / Ecological Modelling 124 (1999) 85–95
pollution variables and topographic exposure to a particular direction. Suitable estimators of the topographic exposure to the given direction can be derived from the digital elevation model (DEM). This could enable the explanation of the pollution spatial variability and also spatial prediction of contamination values, even in cases
87
when pollution source, atmospheric conditions and/or real wind flux are not known.
2. Model description
2.1. Research area and pollution data The research area is located in the Risnjak National Park which is situated in Western Croatia (Fig. 1), near the city of Delnice (f= 45.400° N, l= 14.804° E). Regarding geological conditions, this area belongs to the karstic region. The macroclimate is perhumid: moderately cold, with high precipitation (over 3500 mm annually), high air moisture, high and enduring snow cover and frequent frost over the entire, relatively short growing season. Spatial topoclimatic contrasts are very strong because the relief is extremely rugged, due to the geological porosity and relatively high precipitation (Antonic´, 1996a). The vegetation complex of the Risnjak National Park belongs to the North Dinaric inland type (Antonic´ and Lovric´, 1996). The largest part of the area is covered by beech, fir and spruce monodominant forests and by beech-fir mixed forests and mountain pine woods. The vegetation pattern is mosaical due to the topoclimatic contrasts (Antonic´, 1996b; Antonic´ et al., 1998). The area of Risnjak National Park belongs to the region of Gorski Kotar (Fig. 1). This region is one of the major forestry resources in Croatia.
Fig. 3.
Fig. 3. Spatial distribution of estimators derived from the digital elevation model: 1st row — relative terrain aspect (°); 2nd row — terrain exposure to the horizontal wind flux (following Eq. (1), using b= 0°); 3rd row — horizon angle (°) for the search distance of 300 m; 4th row — terrain exposure to the sloped wind flux (following Eq. (1), using the previous variable as b). Estimators represent topographic exposure to the wind flux for different wind azimuths: 1st column toward the south and 2nd column toward the west. Spatial distributions of topographic exposures for other azimuths are omitted. To improve visualization, values of variables are stretched in a grey scale from the minimum (white) to the maximum (black), with the exception of the horizon angle (3rd variable) where a reverse grey scale is used. The boundary of the Risnjak National Park and the sample points of heavy metal soil concentration are superimposed.
O. Antonic´, T. Lego6ic´ / Ecological Modelling 124 (1999) 85–95
88
During the last three decades the forests of this region have been declining (Prpic´, 1987; Prpic´ et al., 1991). This phenomenon is especially marked in the fir and beech-fir forests while fir is the most endangered of the tree species. A hypothesis of possible causes of the forest decline emphasizes the influence of air pollution, which results in acid rain and soil contamination by heavy metals (Prpic´, 1987; Glavacˇ et al., 1987a; Prpic´ et al., 1991). Several studies of soil acidity and heavy metal soil concentrations (Glavacˇ et al., 1987a; Vrbek and Gasˇparac, 1992; Vrbek et al., 1994) demonstrated that the region of Gorski Kotar is subject to air pollution. Possible pollution sources must be relatively distant because industrial activities in the Gorski Kotar are negligible. Prpic´ et al., (1991) hypothesize that air pollution is coming from industry located in Rijeka Bay (in Croatia) along with the more distant industry of northern Italy. Both are located approximately west of the research area.
In this study, soil pollution data by Vrbek and Gasˇparac (1992) is used. This data set consists of concentrations of lead, copper and zinc (mg/g) in the lower part of the organic soil layer on the 31 point (using a regular square grid 1× 1 km) in Risnjak National Park (Fig. 1). The concentration on each point is the average of ten random measurements in a 10 m radius around the point. Vrbek and Gasˇparac (1992) and Vrbek et al. (1994) interpreted these data qualitatively and concluded that heavy-metal concentrations are lower at a lower altitude and in the leeward side. Similar results are reported by Glavacˇ et al. (1987a) for a wider area of Gorski Kotar, and by Reiners et al. (1975) and Friedland et al. (1983) for areas in the USA. Glavacˇ et al. (1987b) also showed that correlations between altitude and lead-soil concentration in the forests of the North Hessen in Germany are significantly different on the western than on the eastern slopes.
Table 1 Univariate correlation coefficient (r) between four independent variables and soil concentrations of lead, copper and zinca
Pb vs. Asp Cu vs. Asp Zn vs. Asp Pb vs. Exp0 Cu vs. Exp0 Zn vs. Exp0 Pb vs. Hor300 Cu vs. Hor300 Zn vs. Hor300 Pb vs. Exp300 Cu vs. Exp300 Zn vs. Exp300 a
0°
45°
90°
135°
180°
225°
270°
315°
0.0358 0.0076 −0.0368 −0.0586 −0.0023 0.0956 0.0433
−0.3692* −0.3986* −0.212 −0.4202* −0.4368* −0.1198 −0.3508
−0.4535* −0.4985** −0.155 −0.3605* −0.4223* −0.2138 −0.5842**
−0.3313 −0.3164 −0.0816 −0.3233 −0.347 −0.1537 −0.3994*
−0.0358 −0.0076 0.0368 −0.0556 0.108 0.1753 0.009
0.3692* 0.3986* 0.212 0.2817 0.4617** 0.3775* 0.3645*
0.4534* 0.4985** 0.155 0.4994** 0.6078*** 0.3695* 0.3578*
0.3313 0.3164 0.0816 0.311 0.3263 0.2081 0.3097
0.997 0.998 0.981 0.984 0.994 0.979 0.987
−0.1042
−0.4072*
−0.6281***
−0.4631**
−0.0244
0.3595*
0.2814
0.0661
0.975
−0.2164
−0.3585*
−0.3419
−0.2612
−0.0471
0.112
0.0577
−0.0348
0.994
−0.0705
−0.2807
−0.2596
−0.2334
−0.0257
0.1224
0.3816*
0.209
0.971
−0.0333
−0.4041*
−0.3675*
−0.215
0.1383
0.3494
0.6606***
0.4626**
0.981
0.132
−0.1064
−0.2236
−0.1223
0.1523
0.3915*
0.5344**
0.3256
0.995
R
Independent variables estimate topographic exposure to the wind from a given direction (0° means north): ‘Asp’ is relative terrain aspect, ‘Exp0’ is terrain exposure to the horizontal wind flux (following Eq. (1), under b= 0°), ‘Hor300’ is horizon angle (°) for the search distance of 300 m, ‘Exp300’ is terrain exposure to the sloped wind flux (following Eq. (1), using ‘Hor300’ as b). R is proportion of variance of obtained correlation coefficients accounted for by Eq. (2). * Significant at P = 0.05, ** significant at P= 0.01, *** significant at P =0.001, N =31).
O. Antonic´, T. Lego6ic´ / Ecological Modelling 124 (1999) 85–95
89
2.2. Topographic exposure to the wind flux A DEM for the test area is created by digitizing the topographic map (1:5000), generating a triangular network and rasterizing in a spatial resolution of 10×10 m (Fig. 1). Four different estimators of the topographic exposure to the hypothetical wind flux are selected and modelled. Each of them is examined for eight different directions using the azimuthal step of 45° (N, NE, E, SE, S, SW, W and NW). The first estimator is the relati6e terrain aspect. Consider the angle, v, between azimuth g of the maximum terrain slope (terrain aspect) and the azimuth of the hypothetical wind flux d. If v5 180°, then the first estimator (which we denote by s) is chosen to be s= v, otherwise s= 360− v (Fig. 2). The terrain aspect is calculated from the DEM in the usual way, using a regression plane through the 3× 3 neighbourhood of each point. The second estimator is the terrain exposure toward the horizontal component of the wind flux. It is defined as a cosine of the angle a, between the regression plane through the 3× 3 neighbourhood of a terrain point and the plane orthogonal to the wind flux (Fig. 2). The second estimator can be computed from: cos a= cos(m) sin(b)+ sin(m) cos(b) cos(d − g) (1) where m is the terrain slope, g is the terrain aspect, d is the azimuth of wind flux and b is the horizon angle of wind flux (it is zero for all points because horizontal wind flux is hypothesized). The estimator given by Eq. (1) is commonly used in topographic solar radiation modelling (see e.g. Dubayah and Rich 1995; Antonic´ 1998), where a is the solar illumination angle for a given surface
Fig. 4.
Fig. 4. Correlation coefficients (r) between estimators and heavy metal soil concentration versus azimuth: circles — relative terrain aspect; squares — terrain exposure to the horizontal wind flux; rhombs — horizon angle for the search distance of 300 m; triangles — terrain exposure to the sloped wind flux. Horizontally dotted lines indicate significance of correlation coefficient at P= 0.05. Vertical dotted lines indicate direction of two probable pollution sources. Obtained correlations are fitted by the general sinusoid (rEST — solid line) as a function of azimuth of the hypothetical wind blow (Eq. (2)).
90
O. Antonic´, T. Lego6ic´ / Ecological Modelling 124 (1999) 85–95
(defined by m and g) and for a given sun position on the sky (Keith and Kreider, 1978). The third estimator is the horizon angle b. For a given point T, b is maximum vertical angle (determined by the point M) in the azimuth of the hypothetical wind flux, for the chosen search distance a (Fig. 2). Finding the horizon angle for each point requires a search of the DEM in the azimuth of the wind flux. A similar procedure is used in topographic solar radiation modelling (the calculation of the so-called ‘sky-view factor’, see also Dozier and Frew, 1990; Dubayah and Rich, 1995; Antonic´, 1998). The search distance of 300 m is chosen because a further increase in the study area does not change the horizon angle significantly (Antonic´, 1996b). For points where the horizon angle is negative (below the horizontal), it is set to zero. A horizon angle below zero would be dominantly influenced by altitude. The fourth estimator becomes terrain exposure toward the sloped wind flux, following Eq. (1), where b is the horizon angle of wind flux for the search distance of 300 m (the third estimator). Note, that when the search distance is set to zero,
horizon angle is zero and consequently the fourth estimator becomes terrain exposure toward the horizontal component of the wind flux (the second estimator). Before testing the basic hypotheses by the correlation analysis, the relative terrain aspect (the first estimator) and horizon angle (the third estimator) were multiplied by − 1. Consequently all estimators increased by the increase of the topographic exposure to the hypothetical wind flux. Spatial distributions of the four estimators on the research area are shown in Fig. 3, using the two azimuths of the hypothetical wind flux (south and west) as examples.
2.3. Estimation of critical direction Spatial distributions of the four estimators and locations of sample points have been overlayed into a simple geographic information system. A point database (heavy metal soil concentrations) was completed with the corresponding values of the estimators, using local linear interpolation by the inverse distance weight method from four adjacent pixels.
Fig. 5. Correlation between heavy metal soil concentration and estimators of topographic exposure as a function of search distance: lines with circles, lead; squares, copper; triangles, zinc. Filled symbols represent correlation of horizon angle toward the east. Empty symbols represent the correlation of terrain exposure to the wind flux from west (when search distance is zero, the wind flux is horizontal, otherwise it is sloped using horizon angle toward the west). Hatched lines represent significance of correlation coefficient at P =0.05 and P =0.01.
O. Antonic´, T. Lego6ic´ / Ecological Modelling 124 (1999) 85–95
91
Fig. 6. Selection of the best general linear model (Eq. (3)) for the spatial prediction of soil lead concentration as a function of altitude (a), terrain exposure to the horizontal wind flux from the west (b) and horizon angle for the search distance of 300 m toward the east (c). Used independent variables (including interactions terms) are denoted on the x-axis for each particular model (the best model is bold). The line with circles represents regression coefficients and the line with squares represents differences between Cp statistics and the number of independent variables in the particular model.
Soil concentration for each heavy metal is correlated to each estimator of topographic exposure for each of the eight chosen wind directions using linear regression analysis. Results clearly suggest that correlations significantly vary with direction (Table 1). Obtained correlation coefficients (r) are estimated by the following expression: rEST =b0 +b1sin(di + b2)
(2)
where rEST is the estimated r, di is the azimuth of ith hypothesized wind flux, and bj ( j =0, 1, 2) are parameters. The proportion of variance (R) of obtained correlation coefficients accounted for by Eq. (2) is very high (the last column of Table 1). Correlations increase toward the west for each heavy metal and for each estimator of the topographic exposure (Fig. 4). Positive correlations are maxi-
mized toward the west which means that increasing topographic exposure toward the west increases soil concentration of heavy metals. This agrees with the direction of probable pollution sources hypothesized by Prpic´ et al. (1991). Moreover, the major Van Bebber’s cyclone path for this region is the Genova cyclone. It crosses the research area from the west to the east (Pandzˇic´, 1989). Following Robinson (1984), this cyclone could be a major transporter of air pollutants that are being deposited on the geomorphological barrier of the Gorski Kotar region. The absolute value of negative correlations increases toward the direction opposite to the hypothesized pollution sources, i.e. toward the east (Fig. 4). For the first estimator this is a consequence of the linear dependence between relative aspects obtained in opposite directions. For the
92
O. Antonic´, T. Lego6ic´ / Ecological Modelling 124 (1999) 85–95
Fig. 7. Relation between observed values of soil lead concentration (mg/g) and respective values predicted by the selected general linear model in terms of altitude, terrain exposure to the horizontal wind flux from the west and horizon angle using the search distance of 300 m toward the east (Eq. (3)). Points are labelled using the station numbers shown in Fig. 1. The solid line represents equality between observed and predicted values and hatched line represents respective confidence bands at P =0.95.
second and the fourth estimator (terrain exposure to the horizontal and sloped wind flux) a similar opposite effect exists, although correlations toward the west are dominant (see Table 1). However, correlations between soil heavy metals and the third estimator (horizon angle) are higher toward the east than toward the west and this cannot be explained by the opposite effect. In fact, the horizon angle in the azimuth opposite to the wind azimuth quantifies a topographic barrier to the flux of air masses behind a particular point. This turns out to be a better estimator than the horizon angle toward the azimuth of wind that quantifies the barrier in front of the point. Predictive power of a particular estimator is not the same for different metals. Although to some extent this may be explained by the specific behaviour of a metal in the atmosphere, it may also point to limits in the presented method which is based only on the geomorphology. Important parameters such as total pollution output of a particular source and its temporal distribution, wind speed and frequency, atmospheric conditions and processes, as well as physical and chemical attributes of pollutants are ignored. On the other hand, the method is simple, fast and able to
produce useful results when data for more exact models are not available. Moreover, the method is suitable for testing where estimators based on geomorphology alone can determine direction of incoming pollution.
Table 2 Parameters and statistics of the finally selected modela Variable
kib
t(24)
p(t)
Intercept a b c ab ac bc
−0.9615 0.2070 7740.3853 −86.2442 −7.5451 0.0799 39.5343
−0.01 1.90 2.17 −1.73 −2.33 1.84 2.14
0.9936 0.0692 0.0405 0.0964 0.0284 0.0774 0.0424
a Following Eq. (3) for soil lead concentration (mg/g) as a function of altitude above sea level (a, in meters), terrain exposure to the horizontal flux from the west (b, following Eq. (1), using b= 0°) and horizon angle toward the east using the search distance of 300 m (c, in degrees). b Parameters are denoted by ki. The t-value and resulting P-value are used to test the hypothesis that ki is equal to zero. The regression coefficient is R =0.7866. The ratio between regression mean square and residual mean square is F (6,24) = 6.4902 with P(F) =0.00036.
O. Antonic´, T. Lego6ic´ / Ecological Modelling 124 (1999) 85–95
93
distance used for calculation of the respective horizon angle (Fig. 5). Correlation by the soil lead concentration is maximized when the search distance is zero (horizontal wind flux) and this is the reason for selection of the second independent variable. Horizon angle toward the opposite direction (east) is included according to the above discussion. Because the predictive power of the horizon angle for lead is maximized under the search distance of 300 m (Fig. 5), it is selected as a third independent variable. It must be emphasized that a selection of the second and the third independent variable (regarding the search distance) would differ for copper and zinc (Fig. 5). To test the hypothesis that the correlation between pollution and one estimator is not the same for different values of the other estimator, interaction terms are also included. The full tested model is: Fig. 8. Spatial distribution of the predicted soil lead concentration (mg/g) using the selected general linear model in terms of altitude, terrain exposure to the horizontal wind flux from the west and horizon angle for the search distance of 300 m toward the east (Eq. (3)). Predicted values are limited to the range of observed values and stretched in a grey scale from white (100 mg/g or less) to the black (600 mg/g or more). Boundary of the Risnjak National Park and sampling stations are superimposed.
2.4. Predicting the spatial distribution of soil lead concentration Let us assume that estimators of topographic exposure to the wind flux could also be used for spatial prediction of the examined environmental pollution. This is tested for soil lead concentration. The general linear modelling method is used (Ott, 1993). The following independent variables are selected: altitude, terrain exposure to the horizontal wind flux from the west and horizon angle toward the east (without multiplication by the − 1). The altitude is included as the first independent variable because greater air mass flux and consequently greater amount of deposited air pollutants could be supposed at higher sites (Glavacˇ et al., 1987b). Predictive power of terrain exposure to the wind flux from the west changes with the search
y= k0 + k1a+ k2b+ k3c+ k4ab+k5ac+ k6bc (3) where y is soil lead concentration (mg/g), a is altitude (m), b is terrain exposure to the horizontal wind flux from the west (following Eq. (1), using b= 0°), c is horizon angle (°) toward the east for the search distance of 300 m and ki (i= 0, 1, 2, 3, 4, 5, 6) are parameters. Interaction terms are treated as separate independent variables. Cp statistics (Mallows, 1973) is used to find a regression equation which best fits the data. Note that the best model should have the value of Cp statistics closest to the number of independent variables. It has been found that the difference between Cp and a number of independent variables is minimized for the full tested model (Fig. 6). This finally chosen model (Fig. 7 and Table 2) is used to estimate the spatial distribution of lead in the lower part of the organic soil layer (Fig. 8). Obtained spatial distribution is very different from that which has been obtained by Antonic´ et al. (1994) using ordinary spatial interpolation of the same field sample. It is clear that using ordinary spatial interpolation in this case yields a greater error. This is because the used sampling intensity is very low relative to the existing spatial variability conditioned by the extremely rugged
94
O. Antonic´, T. Lego6ic´ / Ecological Modelling 124 (1999) 85–95
karst relief. Note that to obtain acceptable results (involving only interpolation errors; e.g. Phillips and Marks, 1996) while neglecting topographic effects, a much larger sampling intensity is needed.
3. Discussion The presented results suggest that spatially complex problems of air pollution deposition, including explanation and interpolation, can be successfully solved using digital terrain modelling and field pollution sampling, even with a relatively low spatial sampling intensity. The proposed method assumes heterogeneous topographic exposures to the wind on the polluted surface and it is applicable to mountainous or hilly areas. Pollution source or sources, atmospheric conditions and/or real wind flux do not need to be known. An analogous problem arises with two or more pollution sources. When sources are located in different directions the problem would be more difficult to solve, especially in the case of the same pollutant. In general, the approach would be the same: modelling topographic exposure to the wind flux toward different directions, correlating these models and response variables sampled in the field and finding directions where these correlations are maximized. An increase in the number of pollution sources requires an increase in field sampling intensity. But in the case of more sources of the same pollutant, the influence of some particular pollution source could be unclear and blurred by the influences of other sources, regardless of the sample intensity. This could be solved by decreasing the azimuthal step between directions for which topographic exposures are being modelled. If directions toward sources are hypothesized a priori, weights for different sources proportional to the total pollution output of particular source would be helpful. The most simple estimator of the topographic exposure to the wind flux is the relative terrain aspect. It ignores the influence of terrain slope and also the influence of the horizon angle. Terrain exposure toward the horizontal component
of the wind flux ignores only the influence of the horizon angle. Both estimators are available as standard functions in practically every raster–GIS software (e.g. the aspect and hillshade function in the ARC-Info software package). These estimators can be used for the preliminary search of the direction of an unknown pollution source and an explanation of the spatial variability of pollution. More detailed analysis (including horizon angle modelling and optimizing of the search distance for the specific pollutant) is needed when the method is used for the spatial interpolation of the measured pollution values. It is especially recommended when sampling intensity is low relative to the existing spatial variability and when conventional methods of spatial interpolation yield unrealistic results. The same estimators of the topographic exposure to the wind flux derived from DEM could also be used in other ecological studies where wind is an important independent variable, for instance in explaining the spatial distribution of alpine vegetation types (e.g. Brown, 1994) or species (e.g. Guisan et al., 1998).
Acknowledgements This research has been supported by the Croatian Ministry for Research and Technology and by OIKON Ltd., Zagreb. The authors are grateful to two anonymous reviewers for their comments.
References Antonic´, O., Belusˇic´, R., Ananic´, M., Kusˇan V., 1994 GRID 1.0: Development of Landscape Ecology GIS Functions. Proceedings of the 16th International Conference on Information Technology Interfaces, Pula, Croatia, June 1994, pp. 467 – 470 Antonic´, O., 1996a. Application of the spatial modelling in the karst bioclimatology. Cro. Met. J. 31, 95 – 102. Antonic´, O., 1996b. Models of the topoclimatic impact on the karst vegetation, Ph.D. dissertation, University of Zagreb, Faculty of Natural Sciences, Croatia. Antonic´, O., Lovric´, A.Z& ., 1996. Numerical analysis of vegetation complexes and community diversity of major coastal Dinaric mountains. J. Veg. Sci. 7, 73 – 80.
O. Antonic´, T. Lego6ic´ / Ecological Modelling 124 (1999) 85–95 Antonic´, O., 1998. Modelling daily topographic solar radiation without site-specific hourly radiation data. Ecol. Model. 113, 31 – 40. Antonic´, O., Kusˇan, V., Hrasˇovec, B., 1998. Microclimatic and topoclimatic differences between the phytocoenoses in the ‘Viljska ponikva’ sinkhole, Mt. Risnjak, Croatia. Cro. Met. J., in press. Brown, D.G., 1994. Predicting vegetation types at treeline using topography and biophysical disturbance variables. J. Veg. Sci. 5, 641 – 656. Cairns, J. Jr., McCormick, P.V., Niederlehner, B.R., 1994. A proposed framework for developing indicators of ecosystem health. Hydrobiologia 263, 1–44. Cairns, J. Jr., Niederlehner, B.R., 1996. Developing a field of landscape ecotoxicology. Ecol. Applic. 6, 790–796. Dozier, J., Frew, J., 1990. Rapid calculation of terrain parameters for radiation modeling from digital terrain data. IEEE Trans. Geosci. Remote Sens. 28, 963–969. Dubayah, R., Rich, P.M., 1995. Topographic solar radiation models for GIS. Int. J. Geogr. Inf. Syst. 9, 405–419. Keith, F. and Kreider, J.F., 1978. Principles of Solar Engineering. Hemisphere Publishing Corporation, Washington, DC, 778 pp. Friedland, A.J., Johnson, A.H., Siccama, T.G., 1983. Trace metal content of the forest floor in the Green Mountains of Vermont: Spatial and temporal patterns. Water, Air, Soil Pollut. 21, 161 – 170. Glavacˇ, V., Koenies, H. and Prpic´, B., 1987a. Zur Imissionsbelastung der industriefernen Buchen- und Buchen-Tannenwa¨lder in den Dinarischen Gebirgen Nordwest-juogoslawiens. Verh. Ges. O8 kologie, 15. Jahrestagung in Graz 1985, Band XV, 237–247. Glavac´, V., Koenies, H., Jochheim, H., Heimerich, R., 1987b. Relief effects on the deposition of air pollutants in forest stands-lead deposition as an example. Proceedings of the International Symposium held in Grenoble, France, May, 1987b, pp. 520 – 529
.
95
Guisan, A., Theurillat, J.-P., Kienast, F., 1998. Predicting the potential distribution of plant species in an alpine environment. J. Veg. Sci. 9, 65 – 74. Mallows, C.L., 1973. Some comments on Cp. Technometrics 15, 661 – 675. Ott, R.L., 1993. An Introduction to Statistical Methods and Data Analysis. Duxbury press, Belmont, 1051 pp. Pandzˇic´, K., 1989. Factor analysis of the local climatic fields within global stream, Ph.D. dissertation, University of Zagreb, Faculty of Natural Sciences, Croatia. Phillips, L.D., Marks, D.G., 1996. Spatial uncertainty analysis: propagation of interpolation errors in spatially distributed models. Ecol. Model. 91, 213 – 229. Prpic´, B., 1987. Decline of forest trees in Croatia with the special emphasis on the load of ‘Gorski Kotar’ area with the acid rain and heavy metals (in Croatian with English abstract). S& umarski list 1-2, 53 – 60. Prpic´, B., Seletkovic´, Z., Ivkov, M., 1991. Forest decline in Croatia and its relation to the biotical and abiotical factors today and in the past (in Croatian with English abstract). S& umarski list 3-5, 107 – 129. Reiners, W.A., Marks, R.H., Vitousek, P.M., 1975. Heavy metals in subalpine and alpine soils of New Hampshire. Oikos 26, 264 – 275. Robinson, E., 1984. Dispersion and fate of atmospheric pollutants. In: Treshow, M., (Ed.): Air Pollution and Plant Life. J. Wiley, Chicester, pp. 15 – 37. Stein, A., Staritsky, I., Bouma, J., Van Groenigen, J.W., 1995. Interactive GIS for environmental risk assessment. Int. J. Geogr. Inf. Syst. 9, 509 – 525. Vrbek, B., Gasˇparac, M. 1992.New findings of heavy metals input in National park ‘Risnjak’ (in Croatian with English abstract). Radovi S& umarskog Inst. 27-11, pp. 65 – 75. Vrbek, B., Vrbek, M., Vukelic´, J., Gasˇparac, M., 1994. The results of study of heavy metals input (Pb, Cu, Zn) in National park ‘Risnjak’ (in Croatian with English abstract). Proceeding of the Conference on ‘40 years of NP Risnjak’, Crni Lug, Croatia, September, 1993, pp. 140 – 142.
Estimating the direction of an unknown air pollution source using a digital elevation model and a sample of deposition Oleg Antonic´ *, Tarzan Legovic´ Rudjer Bosˇko6ic´ Institute, P.O. Box 1016, Bijenicˇka 54, HR-10001 Zagreb, Croatia Received 7 August 1998; received in revised form 5 May 1999; accepted 21 June 1999
Abstract A new method is described for estimating the direction of an unknown air pollution source. A heterogeneous topographic exposure to the wind on the polluted surface is assumed. Hence, the method is applicable to mountainous or hilly areas. Estimators of the topographic exposure to wind are modelled using a digital elevation model for different hypothetical directions of wind using an azimuthal step. These estimators are correlated to soil pollution but any other relevant field sampling can be used. If significant, maximum positive correlation indicates the direction towards pollution source(s). When spatial sampling intensity is low, estimators of topographic exposure can be used to significantly improve spatial interpolation of measured values. The first application of the method were performed on a sample of heavy metals soil concentrations in the Risnjak National Park in Croatia. Maximum correlations between soil pollution and topographic exposure are obtained toward the west. This is the main direction of incoming cyclones and also the direction toward major industrial pollution sources. © 1999 Elsevier Science B.V. All rights reserved. Keywords: Exposure; GIS; Heavy metals; Landscape ecotoxicology; Soil; Spatial interpolation; Topography; Wind flux
1. Introduction A study of environmental risk assessment, ecotoxicological or ecosystem health problem, is usually based on sampling of response variables. They have to be biologically relevant, sensitive, easily measured and appropriately scaled (Cairns et al., 1994). When such an environmental problem exists on a landscape, spatially explicit approach and modelling procedures are used (e.g. Stein et al. * Corresponding author. Fax: +385-1-4680-084. E-mail address: [email protected] (O. Antonic´)
1995; Cairns and Niederlehner, 1996). There are two basic reasons to include models in the landscape ecotoxicology research: (1) to explain spatial variability of measured responses (using appropriate sample); and (2) to predict responses at locations or times when measurements have not be made, using spatial or temporal interpolation, extrapolation or a process model.
1.1. Aim In this paper the following case is examined: the pollution source is not a priori known
0304-3800/99/$ - see front matter © 1999 Elsevier Science B.V. All rights reserved. PII: S 0 3 0 4 - 3 8 0 0 ( 9 9 ) 0 0 1 4 9 - 0
86
O. Antonic´, T. Lego6ic´ / Ecological Modelling 124 (1999) 85–95
and the field sample of responses is too small for the usual geostatistical approaches (e.g. interpola-
Fig. 2. Estimators of topographic exposure to wind. N indicates north. Terrain in the point T has the maximum slope m and terrain aspect g. The first estimator s is relative terrain aspect at point T for a given azimuth of the wind flux (d). The third estimator b is the horizon angle of the point T for a given d, determined by the point M under the chosen search distance d. Note that a change of d may change b. a is the maximum angle between a plane orthogonal to the wind and regression plane trough the terrain at the point T (point T’ is a projection of point T on the orthogonal plane). Cosine of a is the second estimator if b or search distance are set to zero (terrain exposure to the horizontal component of wind flux), otherwise it is the fourth estimator (terrain exposure to the sloped wind flux).
tion by kriging). This situation occurs often in pollution and ecological monitoring. The question is: how can we explain spatial variability of measured responses and interpolate measured values over the research area under the already mentioned constraints?
1.2. Basic hypotheses Fig. 1. Upper part: location of the research area. Lower part: elevation map derived from DEM used in the case study. Values of altitude above sea level range from black (400 m) to white (1526 m). The boundary of the Risnjak National park and the measurement stations for heavy metal soil concentration are superimposed. Labeling of stations follows original data set (Vrbek and Gasˇparac, 1992).
When the research area is uneven, its spatial units are differently exposed to the air pollution that is coming from some distant pollution source. Direction of the source can be estimated by maximizing correlation between the sampled
O. Antonic´, T. Lego6ic´ / Ecological Modelling 124 (1999) 85–95
pollution variables and topographic exposure to a particular direction. Suitable estimators of the topographic exposure to the given direction can be derived from the digital elevation model (DEM). This could enable the explanation of the pollution spatial variability and also spatial prediction of contamination values, even in cases
87
when pollution source, atmospheric conditions and/or real wind flux are not known.
2. Model description
2.1. Research area and pollution data The research area is located in the Risnjak National Park which is situated in Western Croatia (Fig. 1), near the city of Delnice (f= 45.400° N, l= 14.804° E). Regarding geological conditions, this area belongs to the karstic region. The macroclimate is perhumid: moderately cold, with high precipitation (over 3500 mm annually), high air moisture, high and enduring snow cover and frequent frost over the entire, relatively short growing season. Spatial topoclimatic contrasts are very strong because the relief is extremely rugged, due to the geological porosity and relatively high precipitation (Antonic´, 1996a). The vegetation complex of the Risnjak National Park belongs to the North Dinaric inland type (Antonic´ and Lovric´, 1996). The largest part of the area is covered by beech, fir and spruce monodominant forests and by beech-fir mixed forests and mountain pine woods. The vegetation pattern is mosaical due to the topoclimatic contrasts (Antonic´, 1996b; Antonic´ et al., 1998). The area of Risnjak National Park belongs to the region of Gorski Kotar (Fig. 1). This region is one of the major forestry resources in Croatia.
Fig. 3.
Fig. 3. Spatial distribution of estimators derived from the digital elevation model: 1st row — relative terrain aspect (°); 2nd row — terrain exposure to the horizontal wind flux (following Eq. (1), using b= 0°); 3rd row — horizon angle (°) for the search distance of 300 m; 4th row — terrain exposure to the sloped wind flux (following Eq. (1), using the previous variable as b). Estimators represent topographic exposure to the wind flux for different wind azimuths: 1st column toward the south and 2nd column toward the west. Spatial distributions of topographic exposures for other azimuths are omitted. To improve visualization, values of variables are stretched in a grey scale from the minimum (white) to the maximum (black), with the exception of the horizon angle (3rd variable) where a reverse grey scale is used. The boundary of the Risnjak National Park and the sample points of heavy metal soil concentration are superimposed.
O. Antonic´, T. Lego6ic´ / Ecological Modelling 124 (1999) 85–95
88
During the last three decades the forests of this region have been declining (Prpic´, 1987; Prpic´ et al., 1991). This phenomenon is especially marked in the fir and beech-fir forests while fir is the most endangered of the tree species. A hypothesis of possible causes of the forest decline emphasizes the influence of air pollution, which results in acid rain and soil contamination by heavy metals (Prpic´, 1987; Glavacˇ et al., 1987a; Prpic´ et al., 1991). Several studies of soil acidity and heavy metal soil concentrations (Glavacˇ et al., 1987a; Vrbek and Gasˇparac, 1992; Vrbek et al., 1994) demonstrated that the region of Gorski Kotar is subject to air pollution. Possible pollution sources must be relatively distant because industrial activities in the Gorski Kotar are negligible. Prpic´ et al., (1991) hypothesize that air pollution is coming from industry located in Rijeka Bay (in Croatia) along with the more distant industry of northern Italy. Both are located approximately west of the research area.
In this study, soil pollution data by Vrbek and Gasˇparac (1992) is used. This data set consists of concentrations of lead, copper and zinc (mg/g) in the lower part of the organic soil layer on the 31 point (using a regular square grid 1× 1 km) in Risnjak National Park (Fig. 1). The concentration on each point is the average of ten random measurements in a 10 m radius around the point. Vrbek and Gasˇparac (1992) and Vrbek et al. (1994) interpreted these data qualitatively and concluded that heavy-metal concentrations are lower at a lower altitude and in the leeward side. Similar results are reported by Glavacˇ et al. (1987a) for a wider area of Gorski Kotar, and by Reiners et al. (1975) and Friedland et al. (1983) for areas in the USA. Glavacˇ et al. (1987b) also showed that correlations between altitude and lead-soil concentration in the forests of the North Hessen in Germany are significantly different on the western than on the eastern slopes.
Table 1 Univariate correlation coefficient (r) between four independent variables and soil concentrations of lead, copper and zinca
Pb vs. Asp Cu vs. Asp Zn vs. Asp Pb vs. Exp0 Cu vs. Exp0 Zn vs. Exp0 Pb vs. Hor300 Cu vs. Hor300 Zn vs. Hor300 Pb vs. Exp300 Cu vs. Exp300 Zn vs. Exp300 a
0°
45°
90°
135°
180°
225°
270°
315°
0.0358 0.0076 −0.0368 −0.0586 −0.0023 0.0956 0.0433
−0.3692* −0.3986* −0.212 −0.4202* −0.4368* −0.1198 −0.3508
−0.4535* −0.4985** −0.155 −0.3605* −0.4223* −0.2138 −0.5842**
−0.3313 −0.3164 −0.0816 −0.3233 −0.347 −0.1537 −0.3994*
−0.0358 −0.0076 0.0368 −0.0556 0.108 0.1753 0.009
0.3692* 0.3986* 0.212 0.2817 0.4617** 0.3775* 0.3645*
0.4534* 0.4985** 0.155 0.4994** 0.6078*** 0.3695* 0.3578*
0.3313 0.3164 0.0816 0.311 0.3263 0.2081 0.3097
0.997 0.998 0.981 0.984 0.994 0.979 0.987
−0.1042
−0.4072*
−0.6281***
−0.4631**
−0.0244
0.3595*
0.2814
0.0661
0.975
−0.2164
−0.3585*
−0.3419
−0.2612
−0.0471
0.112
0.0577
−0.0348
0.994
−0.0705
−0.2807
−0.2596
−0.2334
−0.0257
0.1224
0.3816*
0.209
0.971
−0.0333
−0.4041*
−0.3675*
−0.215
0.1383
0.3494
0.6606***
0.4626**
0.981
0.132
−0.1064
−0.2236
−0.1223
0.1523
0.3915*
0.5344**
0.3256
0.995
R
Independent variables estimate topographic exposure to the wind from a given direction (0° means north): ‘Asp’ is relative terrain aspect, ‘Exp0’ is terrain exposure to the horizontal wind flux (following Eq. (1), under b= 0°), ‘Hor300’ is horizon angle (°) for the search distance of 300 m, ‘Exp300’ is terrain exposure to the sloped wind flux (following Eq. (1), using ‘Hor300’ as b). R is proportion of variance of obtained correlation coefficients accounted for by Eq. (2). * Significant at P = 0.05, ** significant at P= 0.01, *** significant at P =0.001, N =31).
O. Antonic´, T. Lego6ic´ / Ecological Modelling 124 (1999) 85–95
89
2.2. Topographic exposure to the wind flux A DEM for the test area is created by digitizing the topographic map (1:5000), generating a triangular network and rasterizing in a spatial resolution of 10×10 m (Fig. 1). Four different estimators of the topographic exposure to the hypothetical wind flux are selected and modelled. Each of them is examined for eight different directions using the azimuthal step of 45° (N, NE, E, SE, S, SW, W and NW). The first estimator is the relati6e terrain aspect. Consider the angle, v, between azimuth g of the maximum terrain slope (terrain aspect) and the azimuth of the hypothetical wind flux d. If v5 180°, then the first estimator (which we denote by s) is chosen to be s= v, otherwise s= 360− v (Fig. 2). The terrain aspect is calculated from the DEM in the usual way, using a regression plane through the 3× 3 neighbourhood of each point. The second estimator is the terrain exposure toward the horizontal component of the wind flux. It is defined as a cosine of the angle a, between the regression plane through the 3× 3 neighbourhood of a terrain point and the plane orthogonal to the wind flux (Fig. 2). The second estimator can be computed from: cos a= cos(m) sin(b)+ sin(m) cos(b) cos(d − g) (1) where m is the terrain slope, g is the terrain aspect, d is the azimuth of wind flux and b is the horizon angle of wind flux (it is zero for all points because horizontal wind flux is hypothesized). The estimator given by Eq. (1) is commonly used in topographic solar radiation modelling (see e.g. Dubayah and Rich 1995; Antonic´ 1998), where a is the solar illumination angle for a given surface
Fig. 4.
Fig. 4. Correlation coefficients (r) between estimators and heavy metal soil concentration versus azimuth: circles — relative terrain aspect; squares — terrain exposure to the horizontal wind flux; rhombs — horizon angle for the search distance of 300 m; triangles — terrain exposure to the sloped wind flux. Horizontally dotted lines indicate significance of correlation coefficient at P= 0.05. Vertical dotted lines indicate direction of two probable pollution sources. Obtained correlations are fitted by the general sinusoid (rEST — solid line) as a function of azimuth of the hypothetical wind blow (Eq. (2)).
90
O. Antonic´, T. Lego6ic´ / Ecological Modelling 124 (1999) 85–95
(defined by m and g) and for a given sun position on the sky (Keith and Kreider, 1978). The third estimator is the horizon angle b. For a given point T, b is maximum vertical angle (determined by the point M) in the azimuth of the hypothetical wind flux, for the chosen search distance a (Fig. 2). Finding the horizon angle for each point requires a search of the DEM in the azimuth of the wind flux. A similar procedure is used in topographic solar radiation modelling (the calculation of the so-called ‘sky-view factor’, see also Dozier and Frew, 1990; Dubayah and Rich, 1995; Antonic´, 1998). The search distance of 300 m is chosen because a further increase in the study area does not change the horizon angle significantly (Antonic´, 1996b). For points where the horizon angle is negative (below the horizontal), it is set to zero. A horizon angle below zero would be dominantly influenced by altitude. The fourth estimator becomes terrain exposure toward the sloped wind flux, following Eq. (1), where b is the horizon angle of wind flux for the search distance of 300 m (the third estimator). Note, that when the search distance is set to zero,
horizon angle is zero and consequently the fourth estimator becomes terrain exposure toward the horizontal component of the wind flux (the second estimator). Before testing the basic hypotheses by the correlation analysis, the relative terrain aspect (the first estimator) and horizon angle (the third estimator) were multiplied by − 1. Consequently all estimators increased by the increase of the topographic exposure to the hypothetical wind flux. Spatial distributions of the four estimators on the research area are shown in Fig. 3, using the two azimuths of the hypothetical wind flux (south and west) as examples.
2.3. Estimation of critical direction Spatial distributions of the four estimators and locations of sample points have been overlayed into a simple geographic information system. A point database (heavy metal soil concentrations) was completed with the corresponding values of the estimators, using local linear interpolation by the inverse distance weight method from four adjacent pixels.
Fig. 5. Correlation between heavy metal soil concentration and estimators of topographic exposure as a function of search distance: lines with circles, lead; squares, copper; triangles, zinc. Filled symbols represent correlation of horizon angle toward the east. Empty symbols represent the correlation of terrain exposure to the wind flux from west (when search distance is zero, the wind flux is horizontal, otherwise it is sloped using horizon angle toward the west). Hatched lines represent significance of correlation coefficient at P =0.05 and P =0.01.
O. Antonic´, T. Lego6ic´ / Ecological Modelling 124 (1999) 85–95
91
Fig. 6. Selection of the best general linear model (Eq. (3)) for the spatial prediction of soil lead concentration as a function of altitude (a), terrain exposure to the horizontal wind flux from the west (b) and horizon angle for the search distance of 300 m toward the east (c). Used independent variables (including interactions terms) are denoted on the x-axis for each particular model (the best model is bold). The line with circles represents regression coefficients and the line with squares represents differences between Cp statistics and the number of independent variables in the particular model.
Soil concentration for each heavy metal is correlated to each estimator of topographic exposure for each of the eight chosen wind directions using linear regression analysis. Results clearly suggest that correlations significantly vary with direction (Table 1). Obtained correlation coefficients (r) are estimated by the following expression: rEST =b0 +b1sin(di + b2)
(2)
where rEST is the estimated r, di is the azimuth of ith hypothesized wind flux, and bj ( j =0, 1, 2) are parameters. The proportion of variance (R) of obtained correlation coefficients accounted for by Eq. (2) is very high (the last column of Table 1). Correlations increase toward the west for each heavy metal and for each estimator of the topographic exposure (Fig. 4). Positive correlations are maxi-
mized toward the west which means that increasing topographic exposure toward the west increases soil concentration of heavy metals. This agrees with the direction of probable pollution sources hypothesized by Prpic´ et al. (1991). Moreover, the major Van Bebber’s cyclone path for this region is the Genova cyclone. It crosses the research area from the west to the east (Pandzˇic´, 1989). Following Robinson (1984), this cyclone could be a major transporter of air pollutants that are being deposited on the geomorphological barrier of the Gorski Kotar region. The absolute value of negative correlations increases toward the direction opposite to the hypothesized pollution sources, i.e. toward the east (Fig. 4). For the first estimator this is a consequence of the linear dependence between relative aspects obtained in opposite directions. For the
92
O. Antonic´, T. Lego6ic´ / Ecological Modelling 124 (1999) 85–95
Fig. 7. Relation between observed values of soil lead concentration (mg/g) and respective values predicted by the selected general linear model in terms of altitude, terrain exposure to the horizontal wind flux from the west and horizon angle using the search distance of 300 m toward the east (Eq. (3)). Points are labelled using the station numbers shown in Fig. 1. The solid line represents equality between observed and predicted values and hatched line represents respective confidence bands at P =0.95.
second and the fourth estimator (terrain exposure to the horizontal and sloped wind flux) a similar opposite effect exists, although correlations toward the west are dominant (see Table 1). However, correlations between soil heavy metals and the third estimator (horizon angle) are higher toward the east than toward the west and this cannot be explained by the opposite effect. In fact, the horizon angle in the azimuth opposite to the wind azimuth quantifies a topographic barrier to the flux of air masses behind a particular point. This turns out to be a better estimator than the horizon angle toward the azimuth of wind that quantifies the barrier in front of the point. Predictive power of a particular estimator is not the same for different metals. Although to some extent this may be explained by the specific behaviour of a metal in the atmosphere, it may also point to limits in the presented method which is based only on the geomorphology. Important parameters such as total pollution output of a particular source and its temporal distribution, wind speed and frequency, atmospheric conditions and processes, as well as physical and chemical attributes of pollutants are ignored. On the other hand, the method is simple, fast and able to
produce useful results when data for more exact models are not available. Moreover, the method is suitable for testing where estimators based on geomorphology alone can determine direction of incoming pollution.
Table 2 Parameters and statistics of the finally selected modela Variable
kib
t(24)
p(t)
Intercept a b c ab ac bc
−0.9615 0.2070 7740.3853 −86.2442 −7.5451 0.0799 39.5343
−0.01 1.90 2.17 −1.73 −2.33 1.84 2.14
0.9936 0.0692 0.0405 0.0964 0.0284 0.0774 0.0424
a Following Eq. (3) for soil lead concentration (mg/g) as a function of altitude above sea level (a, in meters), terrain exposure to the horizontal flux from the west (b, following Eq. (1), using b= 0°) and horizon angle toward the east using the search distance of 300 m (c, in degrees). b Parameters are denoted by ki. The t-value and resulting P-value are used to test the hypothesis that ki is equal to zero. The regression coefficient is R =0.7866. The ratio between regression mean square and residual mean square is F (6,24) = 6.4902 with P(F) =0.00036.
O. Antonic´, T. Lego6ic´ / Ecological Modelling 124 (1999) 85–95
93
distance used for calculation of the respective horizon angle (Fig. 5). Correlation by the soil lead concentration is maximized when the search distance is zero (horizontal wind flux) and this is the reason for selection of the second independent variable. Horizon angle toward the opposite direction (east) is included according to the above discussion. Because the predictive power of the horizon angle for lead is maximized under the search distance of 300 m (Fig. 5), it is selected as a third independent variable. It must be emphasized that a selection of the second and the third independent variable (regarding the search distance) would differ for copper and zinc (Fig. 5). To test the hypothesis that the correlation between pollution and one estimator is not the same for different values of the other estimator, interaction terms are also included. The full tested model is: Fig. 8. Spatial distribution of the predicted soil lead concentration (mg/g) using the selected general linear model in terms of altitude, terrain exposure to the horizontal wind flux from the west and horizon angle for the search distance of 300 m toward the east (Eq. (3)). Predicted values are limited to the range of observed values and stretched in a grey scale from white (100 mg/g or less) to the black (600 mg/g or more). Boundary of the Risnjak National Park and sampling stations are superimposed.
2.4. Predicting the spatial distribution of soil lead concentration Let us assume that estimators of topographic exposure to the wind flux could also be used for spatial prediction of the examined environmental pollution. This is tested for soil lead concentration. The general linear modelling method is used (Ott, 1993). The following independent variables are selected: altitude, terrain exposure to the horizontal wind flux from the west and horizon angle toward the east (without multiplication by the − 1). The altitude is included as the first independent variable because greater air mass flux and consequently greater amount of deposited air pollutants could be supposed at higher sites (Glavacˇ et al., 1987b). Predictive power of terrain exposure to the wind flux from the west changes with the search
y= k0 + k1a+ k2b+ k3c+ k4ab+k5ac+ k6bc (3) where y is soil lead concentration (mg/g), a is altitude (m), b is terrain exposure to the horizontal wind flux from the west (following Eq. (1), using b= 0°), c is horizon angle (°) toward the east for the search distance of 300 m and ki (i= 0, 1, 2, 3, 4, 5, 6) are parameters. Interaction terms are treated as separate independent variables. Cp statistics (Mallows, 1973) is used to find a regression equation which best fits the data. Note that the best model should have the value of Cp statistics closest to the number of independent variables. It has been found that the difference between Cp and a number of independent variables is minimized for the full tested model (Fig. 6). This finally chosen model (Fig. 7 and Table 2) is used to estimate the spatial distribution of lead in the lower part of the organic soil layer (Fig. 8). Obtained spatial distribution is very different from that which has been obtained by Antonic´ et al. (1994) using ordinary spatial interpolation of the same field sample. It is clear that using ordinary spatial interpolation in this case yields a greater error. This is because the used sampling intensity is very low relative to the existing spatial variability conditioned by the extremely rugged
94
O. Antonic´, T. Lego6ic´ / Ecological Modelling 124 (1999) 85–95
karst relief. Note that to obtain acceptable results (involving only interpolation errors; e.g. Phillips and Marks, 1996) while neglecting topographic effects, a much larger sampling intensity is needed.
3. Discussion The presented results suggest that spatially complex problems of air pollution deposition, including explanation and interpolation, can be successfully solved using digital terrain modelling and field pollution sampling, even with a relatively low spatial sampling intensity. The proposed method assumes heterogeneous topographic exposures to the wind on the polluted surface and it is applicable to mountainous or hilly areas. Pollution source or sources, atmospheric conditions and/or real wind flux do not need to be known. An analogous problem arises with two or more pollution sources. When sources are located in different directions the problem would be more difficult to solve, especially in the case of the same pollutant. In general, the approach would be the same: modelling topographic exposure to the wind flux toward different directions, correlating these models and response variables sampled in the field and finding directions where these correlations are maximized. An increase in the number of pollution sources requires an increase in field sampling intensity. But in the case of more sources of the same pollutant, the influence of some particular pollution source could be unclear and blurred by the influences of other sources, regardless of the sample intensity. This could be solved by decreasing the azimuthal step between directions for which topographic exposures are being modelled. If directions toward sources are hypothesized a priori, weights for different sources proportional to the total pollution output of particular source would be helpful. The most simple estimator of the topographic exposure to the wind flux is the relative terrain aspect. It ignores the influence of terrain slope and also the influence of the horizon angle. Terrain exposure toward the horizontal component
of the wind flux ignores only the influence of the horizon angle. Both estimators are available as standard functions in practically every raster–GIS software (e.g. the aspect and hillshade function in the ARC-Info software package). These estimators can be used for the preliminary search of the direction of an unknown pollution source and an explanation of the spatial variability of pollution. More detailed analysis (including horizon angle modelling and optimizing of the search distance for the specific pollutant) is needed when the method is used for the spatial interpolation of the measured pollution values. It is especially recommended when sampling intensity is low relative to the existing spatial variability and when conventional methods of spatial interpolation yield unrealistic results. The same estimators of the topographic exposure to the wind flux derived from DEM could also be used in other ecological studies where wind is an important independent variable, for instance in explaining the spatial distribution of alpine vegetation types (e.g. Brown, 1994) or species (e.g. Guisan et al., 1998).
Acknowledgements This research has been supported by the Croatian Ministry for Research and Technology and by OIKON Ltd., Zagreb. The authors are grateful to two anonymous reviewers for their comments.
References Antonic´, O., Belusˇic´, R., Ananic´, M., Kusˇan V., 1994 GRID 1.0: Development of Landscape Ecology GIS Functions. Proceedings of the 16th International Conference on Information Technology Interfaces, Pula, Croatia, June 1994, pp. 467 – 470 Antonic´, O., 1996a. Application of the spatial modelling in the karst bioclimatology. Cro. Met. J. 31, 95 – 102. Antonic´, O., 1996b. Models of the topoclimatic impact on the karst vegetation, Ph.D. dissertation, University of Zagreb, Faculty of Natural Sciences, Croatia. Antonic´, O., Lovric´, A.Z& ., 1996. Numerical analysis of vegetation complexes and community diversity of major coastal Dinaric mountains. J. Veg. Sci. 7, 73 – 80.
O. Antonic´, T. Lego6ic´ / Ecological Modelling 124 (1999) 85–95 Antonic´, O., 1998. Modelling daily topographic solar radiation without site-specific hourly radiation data. Ecol. Model. 113, 31 – 40. Antonic´, O., Kusˇan, V., Hrasˇovec, B., 1998. Microclimatic and topoclimatic differences between the phytocoenoses in the ‘Viljska ponikva’ sinkhole, Mt. Risnjak, Croatia. Cro. Met. J., in press. Brown, D.G., 1994. Predicting vegetation types at treeline using topography and biophysical disturbance variables. J. Veg. Sci. 5, 641 – 656. Cairns, J. Jr., McCormick, P.V., Niederlehner, B.R., 1994. A proposed framework for developing indicators of ecosystem health. Hydrobiologia 263, 1–44. Cairns, J. Jr., Niederlehner, B.R., 1996. Developing a field of landscape ecotoxicology. Ecol. Applic. 6, 790–796. Dozier, J., Frew, J., 1990. Rapid calculation of terrain parameters for radiation modeling from digital terrain data. IEEE Trans. Geosci. Remote Sens. 28, 963–969. Dubayah, R., Rich, P.M., 1995. Topographic solar radiation models for GIS. Int. J. Geogr. Inf. Syst. 9, 405–419. Keith, F. and Kreider, J.F., 1978. Principles of Solar Engineering. Hemisphere Publishing Corporation, Washington, DC, 778 pp. Friedland, A.J., Johnson, A.H., Siccama, T.G., 1983. Trace metal content of the forest floor in the Green Mountains of Vermont: Spatial and temporal patterns. Water, Air, Soil Pollut. 21, 161 – 170. Glavacˇ, V., Koenies, H. and Prpic´, B., 1987a. Zur Imissionsbelastung der industriefernen Buchen- und Buchen-Tannenwa¨lder in den Dinarischen Gebirgen Nordwest-juogoslawiens. Verh. Ges. O8 kologie, 15. Jahrestagung in Graz 1985, Band XV, 237–247. Glavac´, V., Koenies, H., Jochheim, H., Heimerich, R., 1987b. Relief effects on the deposition of air pollutants in forest stands-lead deposition as an example. Proceedings of the International Symposium held in Grenoble, France, May, 1987b, pp. 520 – 529
.
95
Guisan, A., Theurillat, J.-P., Kienast, F., 1998. Predicting the potential distribution of plant species in an alpine environment. J. Veg. Sci. 9, 65 – 74. Mallows, C.L., 1973. Some comments on Cp. Technometrics 15, 661 – 675. Ott, R.L., 1993. An Introduction to Statistical Methods and Data Analysis. Duxbury press, Belmont, 1051 pp. Pandzˇic´, K., 1989. Factor analysis of the local climatic fields within global stream, Ph.D. dissertation, University of Zagreb, Faculty of Natural Sciences, Croatia. Phillips, L.D., Marks, D.G., 1996. Spatial uncertainty analysis: propagation of interpolation errors in spatially distributed models. Ecol. Model. 91, 213 – 229. Prpic´, B., 1987. Decline of forest trees in Croatia with the special emphasis on the load of ‘Gorski Kotar’ area with the acid rain and heavy metals (in Croatian with English abstract). S& umarski list 1-2, 53 – 60. Prpic´, B., Seletkovic´, Z., Ivkov, M., 1991. Forest decline in Croatia and its relation to the biotical and abiotical factors today and in the past (in Croatian with English abstract). S& umarski list 3-5, 107 – 129. Reiners, W.A., Marks, R.H., Vitousek, P.M., 1975. Heavy metals in subalpine and alpine soils of New Hampshire. Oikos 26, 264 – 275. Robinson, E., 1984. Dispersion and fate of atmospheric pollutants. In: Treshow, M., (Ed.): Air Pollution and Plant Life. J. Wiley, Chicester, pp. 15 – 37. Stein, A., Staritsky, I., Bouma, J., Van Groenigen, J.W., 1995. Interactive GIS for environmental risk assessment. Int. J. Geogr. Inf. Syst. 9, 509 – 525. Vrbek, B., Gasˇparac, M. 1992.New findings of heavy metals input in National park ‘Risnjak’ (in Croatian with English abstract). Radovi S& umarskog Inst. 27-11, pp. 65 – 75. Vrbek, B., Vrbek, M., Vukelic´, J., Gasˇparac, M., 1994. The results of study of heavy metals input (Pb, Cu, Zn) in National park ‘Risnjak’ (in Croatian with English abstract). Proceeding of the Conference on ‘40 years of NP Risnjak’, Crni Lug, Croatia, September, 1993, pp. 140 – 142.