GIS/EM4 - Integration of spatial analysis and ecological risk assessment in a GIS environment

Integration of spatial analysis and ecological risk assessment in a GIS environment:

The case study of the Venetian lagoon contaminated sediments

GIS/EM4 No. 122

Stefania Bertazzon
Claudio Carlon
Andrea Critto
Antonio Marcomini
Gabriele Zanetto

Abstract

A stepwise, GIS-based procedure was developed for the integration of geostatistic and spatial analytic techniques to assist the screening stage of ecological risk assessment. The procedure was applied to the case study of the risk posed to the ecosystem and to humans by the contaminated sediments of the Venice lagoon. The ecological risk was estimated by the toxic unit method, encompassing contaminant concentration, bioavailability and toxicity. The benthic community was selected as assessment endpoint, because of its bioaccumulation capability. In accordance with the probabilistic nature of risk assessment, the spatial distribution of environmental data was analyzed by applying geostatistical methods, which provide an indication of the uncertainty of the estimate. In the case of mercury (Hg), shown to cause the highest risk, the contaminant bioaccumulation by the clam Tapes philippinarum (species chosen also because of its large human consumption and commercial value) was modeled by a spatial regression method. Despite a limited sample size, a predictive model is proposed for mercury (Hg) bioaccumulation by lagoon clams, depending only on environmental parameters (sediment concentration, pH, total organic carbon), since clam-specific characteristics appeared not significant. The promising potential of spatial analysis of environmental data in a GIS environment to support the assessment, representation and management of the ecological and human risk posed by contaminated sediments, is illustrated and its outlooks are discussed.

Keywords

Ecological risk assessment, spatial analysis, GIS, Venice, Venetian Lagoon, spatial correlation, spatial regression, contaminated sediment bioaccumulation factor, geostatistic.

Background

Both the historical center of Venice and the surrounding lagoon have been declared world heritage sites by UNESCO. The Venice lagoon, (ca. 550 Km² wide) hosts a great variety of natural habitats. The diversity of the natural environment was enhanced and endangered by man-made transformations throughout the centuries. These transformations were, at beginning, mainly of hydraulic type, aimed at preserving the lagoon from getting eroded or filled with riverborne material, as it served as a natural defense and privileged harbor for the Republic of Venice. In the early 20th century new major transformations took place as large chemical industrial plants were established just on the lagoon border (Porto Marghera). As result, a great amount of polluting substances have been discharged in the lagoon and, to a significant extent, accumulated in the sediment. This gave rise to a series of recent decrees for the Venetian lagoon, aiming at setting the water quality objectives, the maximum allowable contaminant loadings and new industrial and municipal discharge concentration limits. Since the Sixties, several analytical campaigns have been conducted to assess the contamination of the Venice lagoon sediments, but the spatial analysis of the obtained data has never been performed. In this study, the utilized data concern concentrations of metals and organic pollutants measured in the top 15-cm sediment and in the tissue of clam Tapes philippinarum. All data were collected within a project funded by the Water Authority of Venice and managed by its concessionaire Consorzio Venezia Nuova. As for the sediment, the data of 95 stations were selected out of the whole set of data after a scrutiny of quality control and quality assurance examination (US-EPA, 1990). Moreover, the selected stations are distributed on a 1x1 km grid defined for the specification of hydrodynamic models of the lagoon. The frequency of sampling stations decreased in the seaward direction from 1 sample every 3 cells to 1 sample every 5 cells of the grid. The clam Tapes philippinarum was collected in 23 stations located in areas used for the mussel capture and, whenever possible, in the same areas where the surface sediment was also collected. In addition to the determination of organic and inorganic pollutants, bulk parameters were also acquired for both sediment (i.e. total organic carbon, fine fraction) and clam (i.e. weight and tissue fat content).

Problem statement and work objectives

Pollutants stored in sediments represent a potential risk for the benthic and fish communities mainly through bioaccumulation processes. Moreover, the human intake of toxic substances by consumption of clams and fish poses human health risks, thus strongly endangering the feasibility of important economic activities, such as aquaculture (mainly clam) and fisheries in the lagoon. In this context, several environmental management options, such as the regulation of the contaminant concentration of industrial and municipal discharge and the definition of lagoon zones for fishing and aquaculture, are still not solved. Ecological risk assessment provides a critical element for environmental decision making by providing an approach for considering available scientific information along with the other factors required (e.g., social, legal, political, or economic) in the selection of a course of action (US-EPA, 1998). Exploratory data analysis (EDA) of the available experimental information, the decision making process, and the communication of results may be enhanced by spatial analysis and visualization. In this direction, Geographic Information Systems (GIS) represent a potentially useful tool for the management of many environmental problems. The integration of specific techniques and spatial analytical tools (statistical spatial analysis, geostatistic analysis) within GIS would greatly increase its potentiality. In the present work some of these spatial analysis techniques and tools were tested, which would be useful, but are not currently available in commercial GIS packages. Specifically these geostatistic and spatial analytic techniques assist the screening stage of the Ecological risk assessment (ERA) for:

selecting index contaminants (Chemical of Potential Ecological Concern, COPEC) to be considered in ERA;
assessing the exposure of lagoon organisms to sediment contaminants;
mapping an ecological risk index based on the organism exposure to the contaminants and the contaminant toxicity;
providing a probabilistic model, i.e. spatial regression, capable of interpreting, and possibly predicting trends of spatial data;
developing the analytical framework for the regulation of contaminant discharge and human activities in the lagoon.

Conceptual approach

A risk-based approach is integrated with spatial analysis in a GIS environment. Ecological Risk Assessment (ERA) is defined as "the process that evaluates the likelihood that adverse ecological effects may occur or are occurring as a result of exposure to one or more stressors" (US-EPA, 1998). The ecological risk associated with each pollutant depends on the pollutant concentration in environmental media, its capability of bioaccumulating throughout the trophic chain and its toxicological effects. With respect to the traditional approach, which is simply based on the hazard of the pollutant, the risk-based approach aims at probabilistically assessing the effect of the pollutant, which always depends on the environmental conditions and the endangered target (assessment endpoint). In this case study, the benthic community was considered as assessment endpoint, because it has a critical role within the bioaccumulation of contaminants through the trophic chain. Subsequently, the analysis was focused on the contaminant bioaccumulation by the clam Tapes philippinarum, a benthic organism largely consummated and commercially highly prized by humans. Spatial analysis of data allows for detection of spatial patterns, as well as understanding and modeling of the processes that might be responsible for the observed patterns. Despite their great potential, these techniques have received only limited attention in ecological applications, often owing to the small size of environmental datasets, even though some ecological processes that have been traditionally approached by linear or multivariate statistics, might be better explored by accounting for their spatial dimension. In the study of bioaccumulation of contaminants by marine organisms, it is fundamental to consider the spatial distribution of the contaminants to which the organisms are exposed. In this work, spatial regression was used to predict the clam bioaccumulation of contaminant based on environmental parameters, such as contaminant concentration, pH and organic carbon in sediment, and clam-specific parameters, such as weight and lipid concentration. Spatial regression is a technique for interpreting, or fitting, a dependent variable, based on a set of independent, or explanatory variables. Unlike other analytical techniques (e.g. spatial econometrics), considered model driven, spatial regression is considered a data driven approach, where a greater importance is given to the characteristics of the observed data rather than to the theoretical nature of the observed phenomenon (Getis, 1999, Anselin, 1999). Conceptually similar to standard regression, spatial regression differs from the former for two important reasons: the first is that it makes use of spatially distributed data, the second is that it attempts at dealing explicitly with the complexities of space (Cressie, 1993). Subsequently, in accordance with the probabilistic nature of risk assessment, sediment contamination and bioaccumulation by clams and ecological risk were analyzed in spatial terms by means of geostatistic techniques. For this reason, in this work GIS was not simply employed as a tool to represent the spatial pattern of the environmental phenomenon, but rather as a framework for data management and analysis, to describe and analyze the spatial relations of data by means of models.

Methods

Geostatistics

Spatial correlation of each sediment contaminant was analyzed and described using geostatistic techniques. Geostatistics is based on Tobler's (1979) First Law of Geography, stating that 'everything is related to everything else, but near things are more related than distant things' (Theory of Regionalised Random Variables, Matheron, 1971). Variography and kriging are the geostatistical tools applied. Variography allows us to quantify and to model the spatial correlation among sample locations. In addition, kriging permitted to obtain interpolations from observed values and their spatial relationships, as inferred from the variography (Isaaks and Srivastava, 1989; De Marsily, 1986). As opposed to deterministic data interpolation methods, such as inverse distance functions or triangulation methods, kriging is based on spatial correlation among variables, as empirically tested and modeled on sample data; it also aims at minimizing the error variance, and provides an indication of the uncertainty of the estimate. In this work, spatial correlation analysis was performed by GS plus^© software (Gamma Design Software, 1998), while kriging and spatial distribution maps of pollutants were obtained by using the GIS software Surfer^© (Golden Software, Inc., 1997).

Ecological Risk Analysis (ERA)

The Toxic Units (TU) Method is an implementation of the well established quotient method, which is widely applied to the characterization of the environmental risk (Chapman et al., 1997). In this work, it was used to obtain a preliminary assessment of the ecological risk for the benthic community. TUs are quotients of the concentration of a given chemical found in a medium divided by the standard test endpoint concentration for that chemical (e.g., Threshold Effect Level for benthic community, TEL, MacDonald, 1994). TUs for all sediment pollutants are summed over the investigated stations (Sum(TUs)) and those chemicals accounting for a fixed percentage, usually 99%, of the Sum(TUs) across sites are retained as Contaminants of Potential Ecological Concern (COPECs). Then, specific longitudinal and transversal sections of the lagoon, i.e. transects, were defined in order to include subgroups of stations (e.g., transects including the river inflows into the lagoon, or transects including the outermost stations), and to obtain plots allowing the matching of Sum(TUs) for each station (Fig. 4). To aid interpretation, the relative scores (percentage on Sum(TUs)) for each chemical in several stations were visually represented by pie charts in a GIS framework. The Sum(TUs) provides a score of cumulative total risk (CRS) for each sampling station. The individual contribution of each pollutant to the total risk (CRS) was calculated by dividing the contaminant relative score (UTij) of each j-pollutant by the CRS of each i-station (CRSi). The GIS software MapInfo^© (MapInfo Corporation, 1996) was used to calculate these values and to obtain the map of Fig. 5. The latter shows, for each sampling station, a pie chart with a diameter proportional to the corresponding CRS. This allows for an easy visualization of those stations posing a major risk for the benthic community and an easy identification, through the composition of each pie, of the pollutants that give the highest contribution to the total risk.

Spatial Regression Analysis

Standard regression methods generally do not consider, for their estimates, the irregularities and discontinuities of spatially distributed variables. Also, spatial data are often characterized by spatial autocorrelation (which was tested, for our sample, using geostatistic analysis) which would affects standard estimates. To address the spatial autocorrelation issue and to account for uneven spatial distributions, a fundamental preliminary step of spatial analysis is the specification of a neighborhood structure, or contiguity matrix. The definition of a contiguity matrix is always a difficult operation involving problematic decisions regarding the nature of the data and its distribution, as well as the nature and properties of the process under scrutiny (Griffith, 1995). After defining the contiguity matrix, which contains information on spatial interaction among observations, a vector of spatial weights is calculated, which is used, in the coefficient estimation, to calibrate the influence of each data-point based on its distance and the underlying spatial autocorrelation model. The spatial autocorrelation model for our data was based on the variography analysis described in the geostatistics section. Spatial regression and related statistical routines were calculated using S-Plus^© spatial statistics software (MathSoft, Inc., 1996). To the best of our knowledge, in fact, the spatial analysis module of S-Plus and Anselin's SpaceStat are the only commercial packages capable of performing such analysis.

Results and Discussion

Ecological exposure characterization based on spatial distribution of concentration data

After analyzing the whole set of data, only As, Cr, Cu, and, to a minor extent, Zn exhibited normal or log-normal distributions. Cd and most organic pollutants, with the exceptions of dioxins (PCDD/F) and dioxin-like PCB, showed highly complex distributions not easily explainable. The distributions of Ni, Pb, and Hg turned out to be normal or log-normal only after the depletion of one sampling station where the concentration values were much higher and appear to be outliers. As for metals, the variography produces well defined spatial variogram of exponential type (Fig. 2), with the exception of Pb, that showed a spherical type variogram. A satisfactory spatial correlation was also found for some organic pollutants: aliphatic hydrocarbons exhibit a spherical variogram while PCDD/F obey an exponential type correlation. Noticeably, all spatial correlations could be defined within a radius (i.e. correlation range) varying between 5,000 and 10,000 meters, except for Hg showing a correlation range of 20,000 meters. This feature encourages further investigation, since it is possible, in principle, to associate the correlation range with the environmental mobility of a given pollutant. In our case study, Hg might be more mobile compared with the other pollutants. Once constructed the variograms, the experimental data were interpolated by the kriging which led to obtain the maps of exposure and of the uncertainty associated with the interpolation, which are showed in Fig. 1 and 3 for Hg. The metals were affected by a low uncertainty level (e.g. a standard deviation in between one half and one quarter of the interpolated value), while the organic pollutants with a good correlation (PCDD/F, dioxin-like PCB and HCB) were affected by a medium-high level of uncertainty (e.g. 1 to 2 times the standard deviation associated with the interpolated value).


	Figure 2. Mercury Omnidirectional Variogram

Figure 1. Kriging interpolation map of mercury (Hg) concentration in surface sediment (i.e. sediment exposure map).	Figure 3. Map of Hg kriging interpolation uncertainty (standard deviation).

Ecological risk characterization

The graphs obtained by the application of the Toxic Units Method (Fig. 4) outline the relevant contribution of Hg to the Sum(TU) in most sampling stations. In several stations, located in the northern and central lagoon, Hg exceeded the benchmark value by a factor up to 10. The examination of the longitudinal and transversal transects shows a twofold to threefold reduction of Sum(TU) by moving from the innermost to the outermost lagoon stations. This means that sediments located along the lagoon littoral pose a lower risk to the benthic community mainly because of the decreased exposure to Hg. Such a conclusion can also be drawn by inspecting the Fig. 5 where a lower cumulative risk in the southern and in the seaward direction is shown.

Figure 4. Example of Toxic Units method application.

Figure 5. Pie chart visualization of the individual contaminant contribution to the cumulative risk (CRS) at each sampling station.

Clam bioaccumulation modeling: preliminary results

Clams live in the surface sediment and filter the overlying water for respiration and nourishment. The concentration of low solubility contaminants in water is mainly the result of partitioning processes between water and particulate matter including resuspended sediment. The contaminant partitioning greatly depends upon sediment characteristics, such as concentration of contaminant, pH, organic carbon content and fine fraction. In the case of metals, an important role is also played by the chemical speciation. As both water and resuspended sediment are constantly moving throughout the lagoon, we regarded the spatial regression as a particularly appropriate technique, in that it allows for explicit consideration of areas within which the sediment properties are expected to impact upon the clam bioaccumulation. For each Contaminant of Potential Ecological Concern (COPEC), the bioaccumulation in clam tissues can be fitted by using a spatial regression function. A simple bioaccumulation model typically contains, as explanatory variables, some characteristics of the observed organism, such as weight and fat percentage, and some environmental variables, primarily the contaminant concentration in the sediment, followed by other environmental indicators, such as the sediment total organic carbon and the fine fraction (i.e. the fraction <0.063 mm). After selecting mercury (Hg) as the contaminant, the resulting model can be expressed as follows:

(Notation 1)
Hg(c) = f { W(c), L(c) ; Hg(s), TOC(s), ff(s), pH(s) }

where the dependent variable is Hg(c), i.e. mercury in the clam tissue (c) (mg/Kg ww); W and L are weight and lipid content of the organism, i.e. in grams and percentage, respectively; and the remaining are the environment variables: sedimentary mercury (Hg(s); mg/Kg), total organic carbon (TOC(s); %), fine fraction (ff(s); %), and pH in the sediment (s). Before undertaking the model estimation, the data needed some refinement, as clam samples did not spatially coincide with sediment samples. Empirical approaches were tested, such as the selection of pairs of nearest neighboring points taken from each set, which were subsequently spatially averaged. Finally, a model approach was preferred, based on the kriging interpolations previously performed. In this way, a unique criterion is provided, which was also consistent with the rest of the analysis. The regression was thus performed on the spatial locations of the clam samples, to which the sediment values were attributed, as read from the kriged grid. The original goal was to perform a 2-step regression, with the first step consisting of a model based on organism characteristics (W and L). The scope of this model is to isolate the share of the total variability which is only attributable to the organisms, and to provide a residual vector that can be used as the dependent variable for the second-step model. Scope of the latter is the modeling of the variability in the bioaccumulation only attributable to environmental factors, expressed in the equation by the sediment variables. After several testings on Equation 1, we found that the first-step equation had generally very poor explanatory capability, indicating, in our view, that little variation is due to individual organism features. This finding may also be related with the paucity of the sample data. Thus the data driven approach led to the selection of a model only based on environment variables for the present sample. The 2-step model retains its conceptual value and might be further considered and possibly tested on larger samples. We finally selected a model based on mercury, total organic carbon, and pH, as shown in notation 2.

(Notation 2)
Hg(c) = -0.63 + 0.02 [Hg(s)] + 0.09 [pH(s)] + 0.003 [TOC(s)]

t[intercept]= - 1.9

t[Hg(s)]= 0.97

t[pH(s)]= 2.07

t[TOC(s)]= 1.11

Res. Std. Error = 0.005 on 18 d.o.f.

Correlation (Hg, HgFit) = 0.40

The final specification was reached after a selection process based on statistical significance (t tests) and data characteristics (i.e. ff, was excluded because of its complex distribution and scarce significance in the model). All the variables in the selected model are statistically significant, with the sign of the coefficients consistent with expectations. The estimation method used by S-Plus is the Maximum Likelihood, which provides no indication of the model's goodness of fit. In these cases Anselin (1993) suggests a pseudo-Rsquare, defined as the squared value of the correlation between observed and predicted variable. The above model was used to obtain the map of Hg bioaccumulation in clam from the kriging-based maps of Hg, pH and TOC in the lagoon sediments (Fig. 6). The resulting distribution of Hg bioaccumulation appears similar to the Hg distribution in sediment (Fig. 1) The highest bioaccumulation values are located in the central lagoon, which is affected by industrial and urban discharges, from the Porto Marghera area and the city of Venice, respectively.

Figure 6. Spatial regression based map of Hg bioaccumulation in clam tissue.

Conclusion and perspectives

The work illustrates the benefit obtained by the integration of GIS procedures within the ecological risk assessment approach, in terms of more efficient data processing capabilities and interpretation of results, as well as in terms of more effective communication of the results obtained. Spatial analysis permitted an optimal characterization of the exposure conditions to which the benthic community is exposed, highlighting the areas of greater risk (i.e. areas near the sources of contamination: the industrial zone of Porto Marghera and the river mouths) and the relative contribution of each single pollutant to the overall risk (the highest risk arose mainly from mercury). This suggests the need for a risk analysis scheme completely performed within a GIS support. The preliminary spatial regression function developed to model the bioaccumulation process of mercury from sediment to clams showed spatial regression analysis to be a promising and innovative tool for the investigation of ecological phenomena on the basis of spatially distributed experimental data. The estimated regressions model was affected by a limited sample size. Further estimations of the models are proposed, based on a refined and increased sample. The analysis should also be extended to include a 2-step regression for the investigation of the specific contributions of individual organisms and the environment, respectively. Finally, the need emerged for developing an analytical framework for the joint consideration of spatial and temporal variability in ERA models. Beyond the significance of the approach to ecological risk analysis, the study poses some conceptual questions for the development of spatial regression. One issue (not uncommon in applications) is the availability of two spatially distinct data sets for the dependent and explanatory variables, which had to be merged before performing the regression. Further questions are the definition of spatial contiguity in the case of point samples, whose area of significance extends well beyond the point itself, and to the relationship of such areas with the spatial correlation among sample points. Moreover, the model of spatial autocorrelation based on a single variable semivariogram appears as a poor tool: some algorithm for the integration of multiple variables either in the calculation of semivariograms or in their subsequent use in the spatial regression might represent a great advance. Such issues stimulate the development of analytical tools, but, more importantly, they arise conceptual questions that tackle fundamental issues such as the conceptualization of space, at the core of Geographical Information Science.

References used

Anselin L. 1999. Interactive Techniques and Exploratory Spatial Data Analysis. In: Longley PA, Goodchild MF, Maguire DJ, Rhind DW, editors. Geographical Information Systems. Second Edition. New York: Wiley. p 253-266.

Anselin L. 1993. SpaceStat Tutorial. Regional Research Institute. West Virginia University. Morgantown, West Virginia.

Cressie NAC. 1993. Statistics for Spatial Data. New York: Wiley.

Getis A. 1999. Spatial Statistics. In: Longley PA, Goodchild MF, Maguire DJ, Rhind DW, editors. Geographical Information Systems. Second Edition. New York: Wiley. p. 239-251.

Griffith DA. 1995. Some Guidelines for specifying the geographic weights matrix contained in Spatial Statistical Models. In: Arlinghaus SL, editor. Practical Handbook of Spatial Statistics. CRC Press.

Tobler WR. 1979. Cellular Geography. In: Gale S, Olsson G, editors. Philosophy in Geography. Dordecht: Reidel. p.379-386.

Isaaks EH, Srivastava RM. 1989. An introduction to applied geostatistics. Oxford University Press.

Chapman PM, Cano M, Fritz AT, Gaudet C, Menzie CA, Sprenger M, Stubblefield WA. 1997. Summary report on contaminated site cleanup decisions. In: Ecological risk assessment of contaminated sediments. Pensacola, FL: SETAC Press. p 83-114.

MacDonald DD. 1994. Approach to the Assessment of Sediment Quality in Florida Coastal Waters. Florida Department of Environmental Protection. Tallahassee, FL.

Matheron G. 1971. The Theory of Regionalised Random Variables and its Applications. Report nr 5. Centre de Morphologie Mathematique de Fontainebleau.

De Marsily G. 1986. Quantitative Hydrogeology. London: Academic Press.

U.S. Environmental Protection Agency [US-EPA]. 1990. Guidance for Data Usability in Risk Assessment. Washington, DC. EPA/540/G-90008.

U.S. Environmental Protection Agency [US-EPA]. 1998. Guidelines for Ecological Risk Assessment. Washington, DC. EPA/630/R-95/002F.

Authors

Stefania Bertazzon, Ph.D., Research Associate, Department of Economics
University of Venice Ca' Foscari, Fondamenta S. Giobbe, Cannaregio 873, Venice, Italy.
Email: sbertazz@unive.it, Tel: +39-0412574160, Fax: +39-0412574176.

Claudio Carlon, Ph.D., Research Associate, Department of Environmental Sciences
University of Venice Ca' Foscari, Calle Larga S. Marta 2137, Venice, Italy.
Email: carlon@unive.it, Tel: +39-0412578690, Fax: +39-0412578584.

Andrea Critto, Ph.D. student, Junior Researcher
Fondazione ENI Enrico Mattei of Venice, Campo S. Maria Formosa, Castello 4778, Venice, Italy.
Email: critto@unive.it, Tel: +39-0412711460, Fax: +39-0412711461.

Antonio Marcomini, Associate Professor, Department of Environmental Sciences
University of Venice Ca' Foscari, Calle Larga S. Marta 2137, Venice, Italy.
Email: marcom@unive.it, Tel: +39-0412578548, Fax: +39-0412578584.

Gabriele Zanetto, Professor, Department of Environmental Sciences
University of Venice Ca' Foscari, Calle Larga S. Marta 2137, Venice, Italy.
Email: gzanetto@unive.it, Tel: +39-0412578672, Fax: +39-0412578584.

(Notation 2)	Hg(c) = -0.63 + 0.02 [Hg(s)] + 0.09 [pH(s)] + 0.003 [TOC(s)]
	t[intercept]= - 1.9	t[Hg(s)]= 0.97
	t[pH(s)]= 2.07	t[TOC(s)]= 1.11
	Res. Std. Error = 0.005 on 18 d.o.f.
	Correlation (Hg, HgFit) = 0.40