A New Spatiotemporal Risk Index for Heavy Metals: Application in Cyprus

The main aim of this research was to improve risk mapping of heavy metals by taking account of erosion effects. A new spatiotemporal index, namely the G2met index, is introduced, with integration of pre-existing methodologies (Hakanson, EPM, and G2). The G2met index is depicted as a series of risk maps for each heavy metal on a month-time step. The southern part of Cyprus Island was selected as a study area. Concentration of major heavy metals was extracted with soil sampling in a grid of 5350 sites. A large number of regional-scale risk maps (with 500-m cell size) were created: one for each heavy metal and totally per month and annually; in addition, choropleth maps in terms of statistics per river basin were produced for every metal. Generally, the G2met maps resulted in different spatial patterns in comparison to those depicted by the Hakanson index alone.


Introduction
Heavy metal contamination of soils is regarded as a potential hazard of food safety and public health.Exposure to polluted soils may take place in different ways, e.g., through the consumption of vegetables grown on contaminated soils, drainage of rich in heavy metals waste disposal, dust inhalation, etc.
The off-site transport of chemicals bounded to sediments may cause pollution and silting of water resources.Sediment is defined as the loose sand, clay, silt or other soil particles that settle at the bottom of a body of water [1] and is considered as a habitat and major nutrient source for aquatic organisms.Sediment analysis is important in evaluating qualities of total ecosystem of a water body in addition to water sample analysis, because it reflects the long term quality situation independent of the current inputs [2].
Contamination of soil and water is indirectly referred to by numerous existing legislative policies in the European Union, such as the Directives regarding the disposal of wastes (Landfill, Mining waste), the Sewage sludge Directives and directives regarding the application of chemicals, the Regulation for Plant protection and Pesticide use Directives, the Industrial emissions Directive, the Water Framework Directive, and the Air quality Framework Directive.
A variety of methods has been developed to estimate heavy metal accumulation into soils and sediments.Among them, pollution risk indices are considered to be a powerful tool for ecological geochemistry assessment [3].In 1969, Müller introduced the Geoaccumulation Index [4], originally used with river bottom sediments and more recently for soil contamination [5]; the same index was used as an evaluation tool for pollution from an abandoned mine by [6], or for comparing differences between current and preindustrial concentrations [7].The Average Daily Intake (ADI) is an index focusing on the impact component of risk assessment, while Monte Carlo simulations have been applied in order to reduce distribution uncertainties, using data from different sources [8].Pollution Index (PI) is an index defined as the ratio of the metal concentration to the background concentration of the corresponding metal [7] and is equivalent to the contamination factor of metal as defined by [9].In his review, Gong et al. [3] report also several indices as different expressions of the fuzzy logic theory applied in risk assessment [10].
A range of risk definitions can be found in the literature [11].A common attribute in all these definitions, however, is the probabilistic nature of the risk, which in mathematical terms is defined usually as a multiplicative equation [12]: According to [13], quite commonly risk assessment is lacking of spatial analysis, whereas is based only on static measurements of pollutant concentrations in some critical sites.However, as Lahr et al. [14] notify, spatial analysis help interpretation of potential impacts of environmental stressors.Especially, considering that remote areas may be connected to the pollution sources in ways that are not immediately apparent [15].Moreover, is important for sustainable territorial planning and prevention of natural disasters [16].
Risk maps can be of different types, such as contamination maps, exposure maps, or hazard maps.Critical parameters of mapping specifications with regard to the detail, minimum mapping unit (MMU), value ranging or classification and symbolization, are discussed by [17].The approach for the estimation of the exposure component of risk is different for point from non-point (diffused) pollution sources and for different scales as well.Risk mapping for soil heavy metal pollution from mines (a diffused case) in China was applied by [8].For a point-source pollution problem (olive mill wastewaters), Karydas et al. [18] has followed a multi-scale approach through pathway analysis.
In all currently used indices for heavy metals, however, there is lack of a component which would convincingly express the effect of soil erosion processes in the potential sediment loads.Defined as the process of detachment and transport of soil material by wind or water [19], the soil erosion processes (soil loss, deposition downwards, and sedimentation in surface waters) [20] could be considered as the exposure component of risk assessment for heavy metals and -at the same time-, the linking mechanism between pollution sources and receptors.According to [21], sedimentation data analysis could also support impact assessment of distant human activities.
The main aim of this research was to improve risk assessment of heavy metals potentially reaching surface water courses and water bodies.Improvement could be achieved by a proposed new risk index which would: i. Account for soil erosion processes, thus linking heavy metal concentrations in soils with downstream potential sediment loads.ii.Be dynamic (temporal), in terms of providing assessments in regular short time steps (and not simply static lumped values).iii.Combine all required risk assessment components, i.e., threat, vulnerability, and impact.iv.Be spatial in terms of providing standardized risk maps.

Index definition
The development of the new risk index for heavy metals was achieved through the conceptual, algorithmic, and spatial integration of three pre-existing risk assessment methods, each dedicated to a very particular role (Figure 1): i.The Potential Ecological Risk Index (PERI, or Hakanson index) for the estimation of the risk associated with the concentration of heavy metals in soils or sediments [9].ii.The G2 model for mapping risk of soil loss on a monthly-time scale [22].iii.The Erosion Potential Model (EPM, or Gavrilovic) and specifically its sedimentation component, for estimating sediment delivery ratio (SDR) on a basin scale [23].The new index measures the risk of a specific heavy metal (or all of them) to reach surface waters.The new index is named "G2met" and is defined mathematically as follows: where 2 : spatial risk index for metal during a specific month (dimensionless); : Hakanson index of metal m or the total index value for all heavy metals found in a site (dimensionless); : soil loss (t•ha −1 ) at the site; and : sediment delivery ratio of the river basin in which the contaminated and eroded site belongs (decimal in the range 0-1); = 1•ha•t −1 (dimensional constant).
Defining G2met as a product of the Hakanson index and erosion effect, G2met values remain comparable with Hakanson values, thus indicating the effect of erosion in terms of a numerical multiplier.As is derived from Equation (2), a G2met value would be equal with the Hm value when one ton of polluted soil per ha is transferred to the stream course as sediment.
Τhe G2met index is expressed in terms of a time series of risk maps for each heavy metal alone and entirely, therefore is by default a spatiotemporal risk index.In their review on pollutant risk mapping [14] conclude that risk maps can substantially assist in environmental risk analysis and communication.
The new index follows the fundamental principle of the environmental risk theory, according to which risk is a combined output of three components: threat (or frequency), exposure, and impact (or vulnerability) [10,14].In the composition of the G2met index, threat is represented by the term i f C of the Hakanson index (called "contamination factor"), while impact is represented by the term i f T of the Hakanson index (called "toxicity factor"); therefore, the Hakanson index provides by default the two out of three risk components for heavy metals.
As a consequence, the only missing component of the specific risk is "exposure" of the hazard to the environment, in terms of probability to connect an existing threat with potential impact on ecosystems.This missing component is provided-in G2met-by the incorporated erosion factors, i.e., soil loss (E) and sediment delivery ratio (SDR).A parameter or process bridging sources and receptors in large scale risk mapping is required for an integrative risk assessment at the basin scale.Moreover, incorporation of erosion processes in pollution risk by heavy metals is in line also with the nature of this risk as a non-point source pollution risk.According to the US Environmental Protection Agency (EPA), nonpoint source pollution generally results from land runoff, precipitation, atmospheric deposition, drainage, seepage or hydrologic modification [24].

Study Area
The Mediterranean island of Cyprus, which is rich in heavy metals and highly risky to erosion, was selected as a study area.Cyprus is the third largest Mediterranean island-and the biggest in the eastern Mediterranean basin-with a surface area of 9251 km 2 .However, the current study was restricted to the southern part of the island due to inaccessibility to required heavy metal data in the northern part.The extent of the study area covers mainly the south-west and the south-central part of the island resulting in an extent of about 5947 km 2 (Figure 2).Considering the volcanic geologic feature of Cyprus, the high slope terrain and the intensification of the human activities during the last decades, examination of environmental risk from heavy metal contamination and high potential of these polluted soils to be eroded and mobilized [26] should be urgently assessed.According to the geochemical Atlas of Cyprus, geochemical values are mostly related to parent lithology and subsequent regolith processes, since there is strong correlation between top soil and sub soil geochemical values.However, sediment quality assessment is more complex than water quality assessment alone, due to several site-specific parameters.Stabilization and other forms of long-term self-containing barriers, could reduce the mobility and availability of critical pollutants [27].

Calculation of the Hakanson Index
Hakanson's ecological risk index (RI) evaluates the potential ecological risks associated with the metal contaminant concentrations found in soils and sediment samples [9].The index reflects not only the single ecological risk of individual contaminants, but also addresses the integrated ecological-toxicological effects of multiple pollutants through dividing the ecological risk levels of soil contamination and is calculated using the following equations: where The heavy metal concentration of major heavy metals (Cd, Cr, Cu, Ni, Pb, Zn, Hg, and As) was extracted by the national geochemical map of Cyprus.The soil and sediment monitoring has been conducted for the period May 2006 to January 2009 at the high sampling density of 1 site per 1 km 2 , taking samples of top soil (0-25 cm depth) and sub soil (50-75 cm depth) from a grid of over 5350 sites across the study area.The results of the soil analysis are presented in detail by [28].Inverse-distance weighting (IDW) was used to generate the geochemical maps.The selected IDW model had a grid-cell size of 0.33 km, a distance weighting exponent of 1.6 and a maximum search radius from the center of each grid-cell of 2 km, in order to preserve the point data characteristics of the sample and limit propagation of anomalies from points with very high or low values [28].
Table 1 shows mean concentrations of heavy metal found in the seven main geologic formations, representing 78% of Cyprus extent.All units showed similar heavy metal range with the exception of Basal group that experiences higher Zn (109.8 mg kg −1 ) and Cu (167.5 mg kg −1 ) content.Industrial zones, urban areas and some other locations display higher contents of Pb, Cu, Hg, Sn, and Zn, especially in topsoil.Reference, or background, or quality reference value refers to the natural concentration of an element or a substance in soils that have not been modified by anthropogenic impacts [29].The heavy metals concentration of the sub soil data was used to establish the reference values in each geologic formation of Cyprus based on the statistical analysis (75th percentile of the frequency distribution of the data series [30]) of the laboratory results of samples.The weighted average background value of all examined metals (Table 2) are lower than the values suggested by [31].Chromium, Copper and Nickel values are found to have higher values than those reported by JRC/EU [31] in some geologic units.For instance higher copper values were observed in Basal, Sheeted Dykes and pillow lavas.High background chromium values were found in Serpentine and sheared Serpentine soils.The background values were estimated lower for Cyprus in comparison to Brazilian soils (Table 2) with the exception of the Ni background value [29].
In order to keep the uncertainty of background metal computations low: (1) Estimated values were compared to those established in various international case studies (Table 2); (2) The ratio of top-soil versus sub-soil data of Ni (a heavy metal mainly originated by pedogenic and in general non-anthropogenic processes) was checked and estimated close to one (average value 0.98), thus indicating subsoils data sufficient to estimate background values; and (3) Strong dependency of the geochemical atlas of Cyprus on geologic formations was evaluated through systematic observations [28].Using Equations ( 3)-( 5), a set risk maps according to the Hakanson index were produced from interpolated concentration layers (surface and reference concentrations) and toxicity constant values for all reported heavy metals in the study area and for their accumulative risk (named "total") (Figure 3).Due to the fact that Cd and Hg have distinctly high biological toxicity factors, their contribution to the total risk is significantly higher than the rest of the metals; together they contribute by about 66%.

Soil Loss Mapping Using the G2 model
The G2 model inherits its main principles from the Universal Soil Loss Equation (USLE).According to [32], G2 is classified as a "Watershed to landscape-Pathway type-Averaged-Empirical" model.Soil loss values are calculated by G2 using the following equation [33]: where : soil loss (t•ha −1 ); : rainfall erosivity (MJ•cm•ha −1 •h −1 ); : vegetation retention (V≥1; dimensionless); : soil erodibility (t•ha•h•MJ −1 •ha −1 •cm −1 ); : topographic influence (T ≥ 0; dimensionless); : slope intercept (1 ≤ I ≤ 2; dimensionless).All units are defined by the empirical equations of USLE for R and S (the latter denoted by K in USLE).G2 is designed to run in a GIS environment provided that one spatial layer is prepared for S, one for T, and one for I; whereas, 12 spatial layers are required for R and another 12 spatial layers are required for V (i.e., one layer per month for each of R and V).Since 2010 when it was first introduced, G2 has evolved into a set of alternative formulas for the calculation of all erosion factors.The G2 model has been made available to decision makers by the Soil Erosion section of the European Soil Data Center (ESDAC) of the Joint Research Center, through provision of guidance, datasets and support [33].
The monthly R-factor layers of Cyprus were calculated from 35 precipitation stations with 30-min rainfall data, which are part of the Rainfall Erosivity Database of Europe (REDES) [34].According to [35], the spatial interpolation of monthly erosivity point data use the precipitation data from the WorldClim database as covariates [36] and in accordance to the Generalized Adaptive Model (GAM) [37].
Estimation of the V-factor by G2 is based on the combination of two parameters: the fraction of the surface covered by vegetation (Fcover) and the land use effect (LU).Fcover layers were downloaded for years 2011-2013 from: catftp.vgt.vito.be(last accessed: 06-03-2015) [38,39].The LU values were calculated by compiling CORINE 2006 Land Cover database with the Gavrilovic (or EPC) model empirical data [26].
The S-factor layer resulted from the interpolation of 89 soil samples using the cubist regression model [40], where S at samples was calculated with the original USLE K-equation [41].
For the estimation of topographic influence on erosion (T), the method of [42] was applied.Required hydrological and topographic parameters were calculated from a mosaic of digital elevation model (derived from set of six tiles) downloaded from the ASTER GDEM/METI-NASA geoportal.
The slope intercept factor expresses the potential of an alternating landscape to intercept rainfall runoff by reducing the slope length.Two Landsat 8 images covering Cyprus and acquired during summer 2013 and 2014, were downloaded from the United States Geological Survey (USGS) portal [43].The images were mosaicked prior to their processing.
The output of the G2 model for the study area comprises a set of 12 monthly soil loss layers at 100-m resolution.The average annual erosion rate in the study area was estimated at 11.75 t•ha −1 with a standard deviation of 10.17 t•ha −1 (Figure 4).Three seasons of different risk level can be discriminated: a high risk season during October-December (with averaged monthly rates larger than 1•t•ha −1 ), three medium risk seasons during May-July, January, and September (with averaged monthly rates between 0.5 and 1•t•ha −1 ), and two low risk seasons during February-April and August (with averaged monthly rates smaller than 0.5•t•ha −1 ).In an erosion study in Yialias basin, located in central Cyprus (111 km 2 , outlet next to Potamia village), a mean value of 20.95 t•ha −1 •yr −1 is reported by [44].In the same basin, G2 resulted in a mean value of 10.58 t•ha −1 •year −1 , which is close to the mean value of the entire study area (10.32 t•ha −1 •year −1 ).However, there are not any systematic experimental data in Cyprus, to which the created maps could be compared; therefore, the results could not be validated against ground truth data.

Sediment Delivery Ratio Using the EPM
EPM comprises a soil loss and a soil sedimentation component.The model was developed in the Morava river basin (Serbia) and later was implemented in other areas of Europe.The equations used for calculating SDR are the following [24]: where SD : sediment delivery ratio (decimal in range 0-1), : basin perimeter (km), : difference of mean altitude from minimum altitude of the basin (km), : total length of the primary stream segments (highest order in the specific basin) (km), : total length of the secondary stream segments (km), and : basin extent (km 2 ).All the hydrologic elements, i.e., river basins and stream networks, required to extract input parameter to Equation ( 7), were mapped using the D8 methodology.The stream network was mapped with a flow accumulation threshold set at 500 after trial and error, while stream order was defined using the Strahler method [45].
The data source for extracting all hydrological and topographic information (e.g., Z) was an ASTER GDEM (global digital elevation model), the same used for extracting the required parameters for T factor calculation in soil loss mapping with G2.The implementation of the EPM equation in the spatial domain resulted in a per basin averaged SDR value of 0.28 with a maximum value of 0.71 (Figure 5).

G2met Risk Maps and Statistics
Using Equation (2), twelve G2met layers and another one for the annual risk assessment were created per heavy metal.In addition to that, another twelve layers (i.e., one per month) plus an annual layer were created for the total risk in the study area.As a result, 117 risk layers at 500-m cell size were produced (Figure 6).This big set of layers allows for a detailed monitoring of contamination risk associated with heavy metals in the study area in the spatiotemporal domain (Table 3).
Due to the linear character of the main G2met equation, the monthly distribution of the total risk by G2met follows the same temporal trend as that of erosion.Also, the two layers are overall spatially equivalent, as the mean ratio of normalized (with the mean) annual G2met total risk layer over the normalized (with the mean) Hakanson index layer, was found to be 0.98.Therefore, the main effect of incorporating erosion processes in the calculation of the G2met index (in comparison to the use of the Hakanson index alone) is the spatial redistribution of risk (Figure 7).A mean value per river basin was calculated from the produced G2met layers of monthly and annual risk levels, either for each heavy metal individually or for all heavy metals together.An exploration of the Hakanson and G2met layers per basin shows a dramatic effect of the erosion component of the G2met model on the mapping result, in comparison to the (static) Hakanson index alone.In an indicative comparison for the total risk using normalized risk values (range 0-1), many basins appearing to have low risks with the Hakanson index, show highly risky with the G2met index; and the inverse (Figures 8 and 9).
Radical changes over months can be understood for every metal, as a result of the temporal component of the G2met model (i.e., the monthly step of the assessments).Indicatively, the risk maps of Cu for January and May, respectively, per river basin (mean values) exhibit very different spatial distribution of risk values per basin; e.g., there are basins which are classified in the high risk range in January, whereas in May the risk decreases substantially; and the inverse (see for example, Kouris basin) (Figures 10 and 11

Evaluation Study
One of the largest reservoirs of Cyprus is the Kouris reservoir, located downstream of the Troodos Ophiolite complex and covering significant drinking water needs.The basin draining to Kouris reservoir covers an extent of about 304 km 2 (Figure 12).Covering a surface of about 1.95 km 2 on average, Kouris reservoir has a storage capacity of 115 Mm 3 and the usable surface and groundwater capacity is 370 Mm 3 .Thus, Kouris reservoir stores the 31% of the usable water mass and is of high importance serving as a drinking water supply.4).
The mean G2met index values for Ni, Cr, and Cu in Kouris basin were calculated and their temporal distribution (on a monthly basis) was graphed (Figure 13).From this graph, a large temporal fluctuation of risk for both metals can be detected, with October to December being distinctly the highest risk months, January, June, July, and September being moderate in risk, and the other months being at lower risk levels.These findings are in accordance with the relatively low concentration values indicated for this season of the year (May to June) in the sampled heavy metals of Kouris reservoir.(Note that Cr values appear lower than those of Ni and Cu due to the fact that Cr has been assigned a lower toxicity factor).

Conclusions and Outlook
In this research, a new risk index for surface water contamination by heavy metal concentrations in soils was developed.The new index, namely the G2met index, is by default spatiotemporal, in terms of producing standardized month-time step risk maps.Considering that the new index incorporates the Hakanson index (as a term in the main equation), future studies (using G2met) will be comparable to previous ones (using Hakanson).
An application of the new index in Cyprus showed that it is a useful tool for supporting analysis of ecosystems threatened by heavy metals.In water shortage areas, such as in Cyprus, construction of reservoirs is a common solution for water storage and management, even though the eroded sediment is trapped into their bed.Therefore, in these cases, special attention should be given in examining the potential of sediment to act as pollution source or sink.A preliminary (though limited) evaluation study of the G2met index for Ni, Cr, and Cu in Kouris basin showed that measured values were in accordance with the estimated ones.It seems that G2met is more realistic than previous attempts.
More extended and detailed experiments should be conducted in order to assist in elaboration of an extra factor in the G2met methodology, accounting for soil fraction in the sediment loads.As it is known, fine soil fraction facilitates adsorption and metal binding to iron and manganese oxides and to organic matter and experiences usually high metal concentrations.

Figure 1 .
Figure 1.An overview of the methodology for calculating the G2met index.

Figure 2 .
Figure 2. Cyprus Island is located in the Eastern Mediterranean Basin (see inset); the study area was addressed in the southern part of the island (basemap from ArcMap©).
i f C : contamination factor of metal i; i s C : measured concentration of metal i in the sample, with the following metals used in this method: Cu, Ni, Pb, As, Zn, Hg, Cd, and Cr; i r C : background (reference) concentration of metal i in the sediments; i r E : specific metal potential ecological risk factor; i f T : biological toxicity factor for a single metal, which is used to reflect the toxic levels of heavy metals and the water sensitive to the metal contamination; and RI : potential ecological risk index, describing a comprehensive value of multiple pollutants.Referring to [9], the following i f T values were employed in the current study: Cd = 30; Cr = 2; Cu = 5; Ni = 5; Pb = 5; Zn = 1; Hg = 40; and As = 10.The parameter RI can be divided into four grades: Low risk (RI ≤ 50), Moderate risk (50 < RI ≤ 100), Considerable risk (100 < RI ≤ 200) and High risk (RI > 200).

Figure 4 .
Figure 4.The annual soil loss map of Cyprus (cell size: 100 m; value ranging: user defined).

Figure 7 .
Figure 7.The layer of ratio between normalized (with the mean) layers of annual G2met total risk and Hakanson index (cell size: 500 m; value ranging: histogram equalization; mean value = 0.98). ).

Figure 10 .
Figure 10.Choropleth map of the mean G2met index values for Cu per basin in January (value ranging: histogram equalization; indicative example of Kouris basin, 1.156).

Figure 11 .
Figure 11.Choropleth map of the mean G2met index values for Cu per basin in May (value ranging: histogram equalization; indicative example of Kouris basin, 0.422).

Figure 12 .
Figure 12.Located in the south part of Cyprus, Kouris reservoir was used for a risk assessment verification study.

Figure 13 .
Figure 13.Monthly distribution of G2met index values for Ni, Cr, and Cu in Kouris reservoir.

Figure 14 .
Figure 14.Indicative overview graph of estimated risk values (logarithmic scale) of Ni in the Kouris basin and Ni concentrations in sediments (in mg•kg −1 ) in Kouris reservoir.

Table 1 .
Heavy metal surface mean concentrations in various geologic units of Cyprus (in mg•kg −1 ).

Table 2 .
Heavy metal reference concentrations in the various geologic units (in mg kg −1 ).

Table 3 .
Mean values of G2met layers per heavy metal and month (Hakanson values included for comparison).

Table 4 .
Mean values in mg kg −1 and standard deviation (SD) of heavy metals of sediment in Kouris reservoir (DL-Detection Level).