Rizhao city (119°0′43″E-119°39′22″E and 35°4′55″N-35°36′25″N) is located in the southeastern part of Shandong province, East China (Figure 1), and extends for about 1800 km². The study area consists of an alluvial plain in the eastern part and mountainous area in the western part. The elevation presents a decreasing trend from southwest to northeast with a range of 1.3-656.9 m. Skeleton and brown soils are located in the western mountainous area, which are commonly originated from granite and metamorphic rocks. Alluviation and marine deposition are the main parent materials as the origin of moisture soil and cinnamon soil. The area has a temperate continental climate with a mean annual temperature of 13℃. The average annual precipitation in the area is 870 mm, and the annual average relative humidity is 83%. With the rapid urbanization and industrialization, currently the urban area has increased to 116 km² from 13 km² in the early 1980. The industry concentrates on the coastal area in the eastern part, while the agriculture is mainly located in the western part of the study area.

The study area was divided into 2 km×2 km cells, and then the 445 sites are selected at the centre of each cell (without sites in the water body). At each sampling site, four to six subsamples of topsoil (0-20 cm) were taken and mixed thoroughly to obtain a composite sample. In the case of soil being unavailable at the centre of the cell (i.e. urban area, road, etc.), an alternative location is selected as close as possible to the centre of the cell to find the natural soils. Locations of sampling sites are shown in Figure 1.

Samples were air-dried at room temperature (25℃), and were conducted to remove stones, coarse materials, and other debris. The samples with 10g weight were ground to a fine 0.15 mm powder. The samples for analyzing Cd, Cr, Cu, Mn, Ni, Pb, and Zn were digested using H₂SO₄-HNO₃-HF, the levels of Cr, Cu, Mn, Ni, Pb, and Zn concentrations in the samples were determined by flame atomic absorption spectrophotometer, while Cd was determined by graphite furnace atomic absorption spectrophotometer. The samples for Hg determination were digested by H₂SO₄-HNO₃-KMnO₄, and were determined with an atomic fluorescence spectrometer. The recovery of chemical measure was 10±100%, indicating that the accuracy was satisfactory for this study.

Descriptive statistical indicators including range, mean, median, standard deviation (SD), variation of coefficient (CV), skewness and kurtosis were computed to examine the data structures. Correlation analysis, principal component analysis (PCA), and cluster analysis (CA) were used to evaluate the data to acquire the sources of heavy metals, and were performed using SPSS 16.0 software.

Indicator Kriging (IK) is a non-parametric geostatistical method for estimating the probability of exceeding a specific threshold value at a given location. IK includes Univariate Indicator Kriging and Multiple-variable Indicator Kriging.

In Univariate Indicator Kriging, the original value I(u; z_k) was transformed into indicator variables (Lark and Ferguson, 2004; Lin et al., 2010), which is written as follows:

The following calculation step is concomitant with Ordinary Kriging.

A linear combination of indicator variables can be used to estimate the probability in unsampling location, according to The Indicator Kriging system can be solved using where I^*(u₀; z_k) represents the values of the indicator variable at the measured locations; η is the Lagrange multiplier; γ_j is the variogram of the indicator variable at the respective lag distance, λ_i is a weighting factor in estimation.

Multiple-variable Indicator Kriging can combine indicator variables resulted from respective cut-off together into one integrative variable. Based on the results of various variables calculated in Eq. 1, the multiple-variables integration is obtained by Eq. 4, and the following calculation steps are the same to Univariate Indicator Kriging.

where the weights of various heavy metals W_k are determined by the contribution of respective toxic response factors (Hakanson 1980) according to Eq. 5. The toxic response factors of Hg, Cd, As, Pb, Co, Cu, Ni, Cr, Zn, and Mn are 40, 30, 10, 5, 5, 5, 5, 2, 1, and 1, respectively, which can be used to assess the comprehensive risk from the combination of elements.

where W_k represents the weight of each heavy metal, R_k is toxic response factors; m is the number of elements in the comprehensive evaluation on potential risk.

In this study, the potential hazardous risk derived from Indicator Kriging was conducted using Geostatistical analyst in ArcGIS software (Esri inc., USA).

The descriptive statistics of 445 heavy metals samples are summarized in Table 1. The average contents of As, Cd, Co, Cr, Cu, Hg, Mn, Ni, Pb and Zn were 5.04, 0.20, 10.87, 54.09, 17.57, 0.043, 597.08, 23.52, 27.78 and 63.10 mg/kg, respectively. The mean concentrations of As, Co, Cr and Cu were lower than the background values of eastern Shandong province (Dai et al., 2011), but their maximum concentrations were much higher than background values. The mean concentrations of Pb, Zn, Hg, Cd, Ni and Mn exceeded background values; in particular, the mean contents of Hg and Cd almost reached twice time their responding background values. The coefficient of variation (CV) is an independent measure of relative dispersion, and it is a useful tool to demonstrate the degree of discrete distribution of different elements. Commonly, the lithogenic element has relatively low CV, while the CV of the elements affected by anthropogenic sources tends to be quite high. Based on CV classification, As, Co, Mn, Pb and Zn (26.5%, 32.2%, 21.7%, 21.8% and 31.8%) had moderate variability (15%<CVs<0.36%) , while Cd, Cr, Cu, Hg and Ni (196.9%, 50.6%, 50.6%, 385.0% and 54.6%) had high variability (CVs>36%). It could be suggested that the CVs of Hg and Cd were higher than other elements, and the two elements may be dominated by anthropogenic sources. The skewness values of heavy metal contents showed larger variability of heavy metal contents, with positively skewed frequency distribution. This is common for heavy metals as a result of its common presence of a point source contamination. The skewness of various elements followed a decreasing trend of Cd > Hg > Ni > Cu > Cr > Pb > Zn > Co > As> Mn, the skewness of Cd and Hg were quiet higher than others, indicating that the distribution of Cd and Hg in the study area was highly asymmetric.

Table 2 illustrates the heavy metals levels in Rizhao compared with other regions in the world. The mean content of As in the study area was the lowest among various regions mentioned above. The average content of Cd was lower than Alicante, Ebro, Zagreb, Galicia, Kavadarci, Ireland, and Yangzhong, and higher than other regions. For Co, the mean content was lower than other regions expect for Alicante and Ireland. The mean content of Cr was close to Luhe, Zagreb and Galicia, and lower than Ju, Zhengding and Yangzhong. The level of Cu was relatively low, merely higher than Huizhou, Ebro, Ireland and Duero. The mean value of Hg contents sampled in Rizhao was close to Ebro and Juxian county, and lower than other regions. Mn showed the average content close to Juxian county and Zagreb. The mean content of Ni, similar to Galicia, was lower than Juxian county, Piemonte, Luhe, Zagreb, Zhengding, Kavadarci, and Yangzhong. The mean value of Pb contents, except lower than Juxian county, Huizhou, and Yangzhong, was higher than other regions. The average content of Zn in the study area was close to Piemonte and Ireland, but lower than Juxian county, Luhe, Zagreb, Galicia, Shunyi, Zhengding and Yangzhong. Overall, in the presented study, As, Cu and Hg exhibited relatively lower level compared with other regions, Cd, Co, Ni, Mn, and Zn showed moderate position in all the regions, while Cr and Pb have an upper position among various regions.

3.3.1 Correlation analysis results

Correlation analysis can provide the interesting information on the relationship among various elements, and commonly high correlation between elements could indicate the common source. Pearson’s correlation coefficients of heavy metals in soils are summarized in Table 3. Correlation coefficients of Co-Cr, Mn-Co, Ni-Co, Mn-Cr, Ni-Cr and Ni-Mn were 0.739, 0.719, 0.700, 0.403, 0.947 and 0.367, respectively, and have passed the 0.01 level significance tests. The strongly positive correlation among Co, Cr, Mn and Ni commonly indicates their natural sources. As, Cd, Cu, Hg, Pb and Zn had little correlation with other elements (r<0.25), and these relationships among these heavy metals should be further assessed by multivariate analysis. The correlation analysis results suggested that there existed some original relationship between heavy metals, and indicated different possible sources of elements.

3.3.2 PCA analysis results

Principal component analysis is a powerful tool to identify effective sources of heavy metals. The results of PCA in the topsoil are shown in Table 4 and Figure 2. Four principal components (PC1, PC2, PC3 and PC4) were obtained with an Eigen value higher than 1, together explaining 73.7% variance of the total data. The 36.1% variance was explained by PC1, showing elevated loadings of Co, Cr, Mn, Ni, Zn, and moderate loading of Cu and Pb. PC2 amounted on 15.4% variance of total data, and showed significant loadings of Cd (0.881) and Pb (0.726). PC3 with the 12.1% variance was dominated by the loading As and Cu. PC4 with a variance loading of 10.1% was merely dominated by Hg.

3.3.3 Cluster analysis results

Although CA is not notably different from PCA, this method is an alternative to confirm the results. The cluster analysis results for the heavy metal contents in the studied soils are illustrated in Figure 3. It can be seen that cluster analysis coincided with the PCA. The total data can be divided in to four clusters. Cluster 1 included Cr, Ni, Co, Mn and Zn, which was consistent with the PC1. Cluster 2 contained As and Cu, which were grouped in the PC 3. Pb, Cd and Hg were contained in the cluster 3, which could reflect the PC2 and PC4. It can be demonstrated that the relationship of heavy metals was consistent with PCA.

3.3.4 Sources identification of heavy metals in soils

Based on correlation analysis, PCA and CA results, four elements groups could be identified. As, Cr, Co, Cu, Ni, Mn, partially Pb and partially Zn associated with PC1 were mainly controlled by natural sources. Cd and partially Pb originated from human activities. Hg, as an isolated element, was dominated by human inputs.

The heavy metals in Group 1, Cr, Ni, Co and Mn, with strongly positive correlations in the correlation analysis, had highly positive loading in PC1 and were classified together in CA. The mean contents of Co, Cr, Mn, Ni and Zn were close to or lower than the background value of eastern Shandong province. Sajn et al. (2011), Facchinelli et al. (2001), and Boruvka et al. (2005) suggested a lithogenic origin of Co, Cr, Mn, Ni and Zn in soils, and basic and ultrabasic rocks tend to present high contents than other parent materials. Therefore, it seemed reasonable to infer that PC1 represented the natural origin.

Pb and Cd are commonly influenced by human activity. Pb mainly resulted from the vehicle exhaust and coal combustion. Cd came from the application of fertilizers and pesticides; meanwhile waste water discharge from the electroplating and metallurgy also contribute to the Cd contents (Mico et al., 2006; Rodríguez Martín et al., 2006). Furthermore, the mean contents of Cd and Pb in the study area were much larger than the background values. Therefore, Cd and Pb were seemed to result from the anthropogenic sources.

As and Cu could also be considered lithogenic elements. The sources of As and Cu depend on the human activities in different areas (Sun et al., 2013). Some previous works found that As and Cu were affected by human inputs, such as coal combustion and the application of fertilizers. In the present study, the mean content of As and Cu were lower than background values, and they were negatively correlated with human elements. Therefore, based on the multivariate analysis of topsoil samples, we concluded that the natural sources dominated the As and Cu contents.

In this study, Hg is an isolated anthropogenic element in multivariate analysis, which also has been demonstrated in the works of Cai et al. (2012) and Davis et al. (2009). If an element has equivalent loads in two PCs, this element has the common properties of both PCs (Wang and Lu, 2011). Zn and Pb had the high load in both PC1 and PC2, and it could be concluded that the parent materials and human activities determined the contents of these elements. Cu with the loads in the two natural PCs could be considered as the natural sources.

Geostatistics were commonly employed to examine the spatial structure and variability of heavy metals. In the present study, Ordinary Kriging was applied to map the spatial patterns of four PCs derived from principal component analysis. The potential hazardous risk of elements is defined here as probability of exceeding the respective background value in eastern Shandong province (Dai et al., 2011), which could be delineated by Indicator Kriging.

3.4.1 Variograms fitting

The variogram could depict the variance of the sample values at various distances of separation, and provides the input parameters for the spatial interpolation of Kriging. The types of variograms mainly include Spherical, Exponential and Gaussian models, with parameters including nugget (Co), sill (Co+C), range values, coefficients of determination (R²). Nugget represents the experimental error and field variation within the minimum sampling spacing. The ratio of nugget to sill can be used to reflect the extent of spatial autocorrelations of the factors. If the ratio is less than 25%, the variable has strong spatial dependence; between 25% and 75%, the variable has moderate spatial dependence; and greater than 75%, the variable shows only weak spatial dependence. Strong spatial dependence of soil properties can be contributed to intrinsic factors (soil formation factors), and weak spatial dependence can be contributed to extrinsic factors. Coefficients of determination (R²) reflect the precision that elements are fitted by variogram.

The attributes of the variogram fitting are summarized in Table 5. The variograms of PC1, PC4, Cd, Co, Cr, Cu, Hg, Mn, Ni, Co-Cr-Mn-Ni-Zn, Cd-Pb, As-Cu and all elements integration variables were adjusted best to an exponential model and their ranges were between 1660m to 45900m. PC2, PC3, As and Pb variables could be fitted with spherical model, and Gaussian model was used to fit Zn variables. The values of R² were greater than 0.553, showing that the variogram models gave good descriptions of spatial structure. The nugget/sill ratios of PC1, PC3, As, Co, Cr, Cu, Mn, Ni, Co-Cr-Mn-Ni-Zn and As-Cu were 0.192, 0.234, 0.197, 0.145, 0.191, 0.123, 0.125, 0.124, 0.212 and 0.125 respectively, indicating high spatial dependence due to the effects of natural factors such as parent materials and topography. The nugget/sill ratios of PC2, PC4, Cd, Hg, Pb, Zn, Cd-Pb and all elements integration ranged from 0.409 to 0.818, suggesting low spatial dependence related to anthropogenic factors.

3.4.2 Spatial patterns of PCs

Figure 4 indicates the spatial patterns of various principal components. The spatial pattern of PC1 was characterized by high contents in western mountainous region and low content along the coastal areas, which was similar to the pattern of parent material. The soils from granite and metamorphic rock have higher geological background than alluviation and marine deposit along the coastal area. Therefore, the distribution of PC1 was controlled by the parent materials, which is consistent with the results of variogram analysis.

The high values of PC2 were located in western agricultural area and eastern urban area (Figure 4b). The developed economy and intensive industry resulted in the high values in eastern urban area. There is serious soil intrusion in the western part; however the soil fertility in the western part is higher than that in the eastern part, which could be attributed to the excessive application of chemical fertilizer in the western part. The application of chemical fertilizer in the western part is 2.39 times higher than the mean level in the whole study area (RMBS, 2012), which contributed to the high PC2 scores in the western part.

The higher region of PC3 is located in the coastal area (Figure 4c). The As content in urban area was higher than others, but it always was lower than the background value. The high geological background in alluviation and marine deposit was the dominating factor on the content of As, which is consistent with the tendency of the whole Shandong province (Dai et al., 2011). The spatial pattern of PC4 followed the decreasing trend from the eastern to the western part, and the high scores were consistent with the location of urban area, indicating its human source (Figure 4d).

3.4.3 Delineating high risk areas by Indicator Kriging

Figures 5 and 6 illustrate the spatial distribution of hazardous risk of various elements and elements association in the topsoils. The kriged maps showed that the single element, elements association and their corresponding PCs had some similar features, indicating that it is reasonable to apply Indicator Kriging to delineate hazardous risk of heavy metals.

The spatial distributions of Co, Cr, Mn, Ni, Zn and Co-Cr-Mn-Ni-Zn were characterized by an overall tendency downward from western to eastern part, and the areas with high risk probabilities (Prob>0.8) were mainly located in the western part originated from granite and metamorphic rocks (Figures 5c, 5d, 5g, 5h 5j and 6a). For Cd-Pb, the high risk probability accounting for 5.5% of the total area was associated with the eastern urban part and the western study area (Figure 6b and Table 6). The map of As-Cu had the highest proportion of non-risk level and showed no area of high risk level (Figure 6c and Table 6). The probability map of Hg showed that high risk probability of Hg was evidenced in the eastern urban area, and covered 18503 hm² of the total area (Figure 5f and Table 6).

For the comprehensive risk of all elements integration, no area was found to be in the highest risk level, and low and moderate risk level covered more than 90% of the study area (Table 6). The considerable risk was mainly located in western and eastern parts, which could be attributed to the high geological background in the western part and strong human activity in the eastern part.