Probabilistic ecological risk assessment of heavy metals in western Laizhou Bay, Shandong Province, China

Considering the serious land-based pollution and the weak water exchange ability of western Laizhou Bay, it is essential to conduct an ecological risk assessment of the pollutants in this area. In this study, the ecological risk caused by heavy metals deposited in the surface sediments and those resuspended in the seawater of western Laizhou Bay was evaluated using probabilistic approaches. First, the concentrations of seven heavy metals, namely As, Cd, Cr, Cu, Hg, Pb, and Zn, in the surface sediments and seawater of western Laizhou Bay were detected during the spring and autumn of 2016. The concentrations of As, Cd, Cr, Cu, and Pb were found to be at levels comparable to those in the other global coastal systems, while those of Hg and Zn were lower than those in other coastal areas. Next, an ecological risk assessment of heavy metals in the surface sediments was performed using a typical potential ecological risk index and refined by using a Monte Carlo simulation. The results suggested low risk for the heavy metals detected in the sediments of western Laizhou Bay, with the exception of Hg in September 2016, which showed a probability (0.03%) of moderate risk. Meanwhile, the aquatic ecological risk assessment of the heavy metals was performed by applying a combination of hazard quotient (HQ) and joint probability curve. While the ecological risk of Cd, Hg, and Pb was found to be acceptable, the HQs for Cr, Cu, and Zn were greater than 1, and the overall risk probability of their adverse effects was higher than 0.05, suggesting certain ecological risk. Specifically, in the case of As, the overall risk probability was lower than 0.05, suggesting that its ecological risk was acceptable, although its HQ was greater than 1. Thus, by applying the probabilistic approaches, the ecological risk of the heavy metals in western Laizhou Bay was better characterized in this study, avoiding both overestimation and underestimation of ecological risk.


Introduction
Due to their poor biodegradability, easy bioaccumulation, and high toxicity, heavy metals discharged into the sea from different sources may pose serious threats to marine organisms. For

Sample collection and analysis
The field survey in this study was conducted based on the Marine Environmental Impact Assessment Project of Guangli Port Logistics Park, which was approved by Dongying Marine and Fisheries Bureau. Twenty sampling stations were set up in western Laizhou Bay (1185 5 0 03.08@-119˚25 0 28.78@ E, 37˚16 0 50.22@-37˚36 0 11.78@ N) in May (spring) and September (autumn) 2016 (Fig 1 and Table 1).
Surface sediment samples (0-5 cm sediment layer) were obtained by a grab sampler and collected using glass jars for Hg analysis and polythene bags for As, Cd, Cr, Cu, Pb, and Zn detection. Surface seawater samples (at a depth of 0.5 m) were collected using glass bottles (for Hg detection) and plastic bottles (for analysis of the other heavy metals). Three seawater samples and three sediment samples were collected from each station. The sealed sediment and water samples were sent to our laboratory for treatment and analysis, in accordance with the Specifications for Marine Monitoring (GB 17387. 4-2007 andGB 17387.5-2007). The analytic techniques and the detection limits are shown in Table 2.

Toxicity data collection
Chronic toxicity data for the seven heavy metals with respect to their impact towards marine species were collected from the ECOTOX database (http://cfpub.epa.gov/ecotox/) and screened according to the criteria of reliability, relevance, and adequacy [27]. No observed effect concentration (NOEC) was adopted as the primary endpoint representing chronic toxicity, with maximum acceptable toxicant concentration (MATC) and lowest observed effect concentration/level (LOEC/LOEL) serving as a supplement. Data were adopted only from exposure experiment with adequate duration. To be specific, the exposure period should be  Table 3). The data on the toxicity values for each pollutant with respect to all six major functional groups of the marine ecosystem were involved, meeting the requirement by the US EPA of at least eight families in three classes of tested organisms.

Ecological risk assessment approach
Potential ecological risk index. According to Hakanson [28], the potential ecological risk of a given substance in the sediments was calculated as follows: where E i r is the potential ecological risk factor of substance ''i", T i r is the toxic response factor of substance ''i" (which is 10 for As, 30 for Cd, 2 for Cr, 5 for Cu and Pb, 40 for Hg, and 1 for Zn [28]), C i f is the contamination factor of substance ''i", C i 0 is the measured concentrations in the sediments of substance ''i", and C i r is the background reference level for substance ''i". Grade I of the Marine Sediment Quality (GB 18668-2002) was adopted as C i r in this study. The following grades were used for the E i r value: (I) low risk: E i r < 40; (II) moderate risk: 40 � E i r < 80; (III) considerable risk: 80 � E i r < 160; (IV) high risk: 160 � E i r < 320; (V) very high risk: E i r � 320. RI represents the ecological risk for the sediment. It was calculated as the sum of E i r and categorized into the following four classes: (I) low risk: RI < 150; (II) moderate risk: 150 � RI < 300; (III) considerable risk: 300 � RI < 600; and (IV) high risk: RI � 600: Monte Carlo simulation. The probability distribution of E i r and RI values was obtained by using Monte Carlo simulation [29]. The measured environmental concentrations of each metal in the sediments were used as a data set comprised of random variables that conform to a certain probability distribution. The commonly used cumulative probability distribution functions, which mainly include Weibull, log-normal, log-logistic, and Burr III methods, were all applied for the fitting of the data set. The most suitable model was selected based on the Kolmogorov-Smirnov test (S1 Table): the closer the P value is to 1, the better is the fitting effect. When the Burr III distribution was used, its limit distribution appeared, and the fitting effect was poor. In addition, Burr III distribution was often recommended for the species sensitivity distribution (SSD) model fitting [30], but not for that of environmental monitoring data. Therefore, Burr III distribution is not included in S1 Table. Therefore, the log-logistic distribution was found to be the most suitable; the parameters of this distribution are presented in S2 Table. Monte Carlo simulations were performed for 100,000 times using the MATLAB 2017b software. Hazard quotient. Hazard quotient is the quotient of the environmental exposure concentration (EEC) and predicted no effect concentration (PNEC). While HQ > 1 indicates potential ecological risk, HQ < 1 suggests that the ecological risk is at an acceptable level [26,31,32]. In this study, the geometric mean of the heavy metal concentrations detected in the seawater were used as the EEC to calculate the HQ in general. The PNEC value was calculated as HC 5 (hazardous concentration affecting 5% of species) divided by the safety factor (SF = 5) [33,34]. The value of HC 5 was derived from the SSD [26,31]. The above-mentioned cumulative probability distribution functions were applied to derive SSD. The most suitable model was selected based on the Anderson-Darling test (S3 Table): the closer the P value is to 1 and the smaller the Akaike information criterion (AIC), the better is the fitting effect [31,35,36]. Overall, the log-logistic distribution was found to be the most suitable, and its parameters are presented in S4 Table. Joint probability curve. The same methods and criteria used for selecting the probability distribution model applied to measured concentrations of heavy metals in the surface sediments were also applied to that of the seawater. The results were similar as well, i.e., the fitting results of the log-logistic distribution model were better than those of the log-normal and Weibull methods (S5 and S6 Tables), which were adopted in the JPC construction. The JPC was generated using the cumulative probability of the toxicity data from the SSD as an independent variable and the reverse cumulative probability of the exposure data (or exceedance probability, EXP) as the dependent variable to describe the probability of a certain proportion of species expected to be adversely affected [17]. The distance between the generated curve and the axes positively indicated the risk level, and the area under the curve showed the overall risk probability (ORP) of the adverse effects: where x is the proportion of adversely affected species, and EXP(x) is the exceedance probability of the exposure data associated with 100x% of the adversely affected species.

Measured concentrations of heavy metals in the surface sediments of western Laizhou Bay
The concentrations (mg/kg) of heavy metals in the surface sediments of western Laizhou Bay were in the range of 9.20-12.70 (average 11.01) for As, 0. , and were comparable to (for As, Cd, Cr, Cu, and Pb) or lower than (for Hg and Zn) those of other coastal systems around the world. The results of the matched-pair t-test (pair of May-Sep) showed significant differences for all the seven elements. Thus, the two data sets were separately analyzed in the following ecological risk assessment (S7 Table).

Ecological risk of heavy metals in the surface sediments of western Laizhou Bay
By applying a Monte Carlo simulation, the ecological risks of each heavy metal and the mixture in the surface sediments of western Laizhou Bay were expressed as a probability distribution of E i r and RI values instead of single-point estimates. Fig 2 shows the cumulative probability curves of E i r for each heavy metal. Apparently, the E i r curves of Hg, Cd, and As are towards the right compared to those of Cr, Cu, Pb, and Zn; however, all curves are to the left of the straight line E i r ¼ 40, indicating low risk. The Monte Carlo simulation demonstrated that only Hg in September 2016 showed a probability (0.03%) of moderate risk ( Table 6). The sources of Hg in this area were mainly land-based human activities, including factory discharge and combustion of fossil fuels, and river transportation is possibly the main means by which Hg enters Laizhou Bay [12,14].
The RI value was 28.50 and 41.04 for May and September 2016, respectively, suggesting low risk for the sediments of western Laizhou Bay. The Monte Carlo simulation also showed that the combined ecological risk caused by these seven metals is 100% low (Fig 3). Generally, the average, conservative, or maximum values were adopted in typical risk assessment indices, to estimate the average, conservative, or the worst case of ecological risk [45,46]; by doing this, the risks might be either under-or overestimated. For example, in this study, although the average estimation of E i r for Hg in September 2016 was 20.24, which is only 0.51 times of 40 (the upper limit of the low risk), a probability (0.03%) of moderate risk still existed (Table 6); similarly, although high risk was identified in the Xiangjiang River and Dianchi Lake according to the average RI values, the Monte Carlo simulation indicated that the probabilities of considerable risk level reached only as high as 43.3% in the Xiangjiang River and 47.1% in the Dianchi Lake [29]. The combined approach using PERI and Monte Carlo simulation joint approach may therefore avoid either under-or overestimation of the ecological risk and provide a more objective scientific evidence for the environmental management of  the polluted aquatic bodies. Furthermore, when the potential ecological risks of heavy metals in soil in a study area in the urban-rural transition zone of the Wuhan City, China was characterized, the spatial distribution of heavy metal was simulated using sequential Gaussian simulation (SGS), which is also a Monte Carlo method, and then the simulated realizations was fed into the Hakanson PERI computation equation to obtain the response maps of the PERI for each metal [47]. Combining SGS or other geostatistical stochastic simulation methods and the Hakanson PERI may be a way for assessing the spatial distribution and uncertainty of the potential ecological risk of heavy metals in sediments.

Measured concentrations of heavy metals in the surface seawater of western Laizhou Bay
The heavy metal levels in the surface seawater of western Laizhou Bay were 3.01-3.87 μg/L for As, 0.  Table), the concentrations of Cd and Cu in September were significantly higher than those in May (P < 0.05). Hence, the ecological risks were separately assessed for the two time points.

Ecological risk of heavy metals in the surface seawater of western Laizhou Bay
The ecological risk of heavy metals to the marine ecosystems should include the effect of their deposition in the sediments and resuspension into the surface seawater. However, compared to the various index approaches available for the ecological risk assessment of heavy metals in the sediments, methods for those in the seawater are very limited. Typically, only Single Factor and Nemerow index methods are available. In recent years, new approaches were developed based on the SSD theory [48,49], and HQ, which is the simplest among them, was adopted in this study. The results showed that HQ was less than or equal to 1 for Cd, Hg, and Pb, indicating that their potential ecological risk was acceptable, whereas HQ was greater than 1 for As, Cr, Cu, and Zn, suggesting unacceptable ecological risk (Table 8).
Although the HQ method has its advantages such as simplicity and low data requirements [25,50], it is a screening-level method applied to focus on the most dominant pollutants. HQ > 1 does not necessarily demonstrate the real risk of As, Cr, Cu, and Zn. Refinement with higher level methods should follow to alleviate uncertainty to an acceptable degree. Therefore, in this study, the JPC method was employed to refine the aquatic ecological risk assessment. For As, Cr, Cu, and Zn, the distances between their JPC curves and the axes were further than those for the other three elements, indicating their higher risk level (Fig 4). Generally, the ORP of the adverse effects is considered to be acceptable when it is not higher than 0.05 [51]. However, the ORP of Cr, Cu, and Zn ranged from 0.086-0.087, 0.092-0.096, and 0.070-0.073, respectively, suggesting certain ecological risk. This was consistent with the results of HQ. In the case of As, although its HQ was also greater than 1, its ORP was < 0.05, suggesting an overestimation of its ecological risk in the water column by HQ ( Table 9). The sources of Cu and Zn were mainly from natural contribution, including coastal erosion, weathering products carried by the surrounding short rivers, and loess matters carried by the Yellow River, while those of Cr were likely from anthropogenic discharges such as metal processing, fuel burning, and domestic sewage, in addition to natural inputs [13][14][15].
According to the Marine Functional Zoning of Shandong Province (2011-2020), Grade I (for site 1, 2, 3, and 6 in the protection zone) or Grade II (for most sites except for those in the protection zone) of the Sea Water Quality Standard (GB 3097-1997) should be adopted in the study area. Chemical monitoring showed that concentrations of Cr and Cu met Grade I, while the Zn levels met Grade II. However, according to the results of HQ and JPC, Cr, Cu, and Zn posed potential ecological risk to the aquatic environment of western Laizhou Bay, while Cr and Cu were considered as the main aquatic pollutants. The HC 5 value of Zn was higher than Grade I and lower than Grade II, while that of Cr and Cu was lower than Grade I. The criterion continuous concentrations (CCCs) suggested by US EPA [52] were 50 μg/L for Cr, 3.1 μg/L for Cu, and 81 μg/L for Zn, respectively, and the HC 5 values of Cr, Cu, and Zn were 0.02, 0.28, and 0.36 times each CCC, respectively, showing certain differences. Mu et al. [53] reported a HC 5 value of 1.8 μg/L for Cd using chronic toxicity data from local marine organisms in Bohai Bay, which was slightly less than the value of 2.57 μg/L obtained in this study. Thus, the thresholds for environmental protection derived from the chronic toxicity data in this study are observed to be relatively low, and the risk assessment using these values can better characterize the ecological risk posed by the heavy metals in western Laizhou Bay.  On the one hand, uncertainty is inevitable in ecological risk assessment even when high tier approaches are conducted. This could be due to the following reasons. First, physical and chemical parameters (such as hardness, pH, and suspended solid) may affect the distribution and bioavailability of heavy metals in the water environment, and ultimately affect their toxicities to aquatic organisms [54]. Second, the toxicity data used for the construction of SSDs are from marine organisms all over the world, instead of a pool of the local species in the study area. For example, 38 d LOEC of H 2 CrO 4 to the sensitive species Palaemon elegans (0.003 μg/ L) and 14 d LOEC of H 2 CrO 4 to the sensitive species Palaemonetes varians (0.001 μg/L) were adopted when the HC 5 of Cr was calculated [55]. Thus, a high HQ was derived for Cr, and the fitting effect of the SSD curve (especially the tail) was also influenced. On the other hand, the spatial distribution of pollution risks was not mapped in this study. Since sampling in marine environments is difficult and expansive, 20 stations were set in this study, which were comparable to those in other studies focusing on heavy metals or other pollutants in western or southwestern Laizhou Bay [12,13,16,56]. However, the sample data is limited when the variogram estimation is conducted. If more sample data is available especially when a wider sea area is studied, the SGS can be used to conduct a stochastic spatial simulation and map the pollution risks of heavy metals [25,47,57,58].

Conclusion
The presence of the heavy metals namely As, Cd, Cr, Cu, Hg, Pb, and Zn was detected in the surface sediments and seawater of western Laizhou Bay during the spring and autumn of 2016; and their concentrations were found to be comparable to or lower than those in other coastal areas around the world. The typical potential ecological risk index and Monte Carlo simulation suggested low risk for the sediments of western Laizhou Bay, with the exception of Hg during September 2016, which showed the probability of a moderate risk. The HQ and JPC indicated certain ecological risk for Cr, Cu, and Zn, and acceptable risk for Cd, Hg, Pb, and As in the surface seawater. The ecological risk of heavy metals in western Laizhou Bay was better characterized in this study by applying the probabilistic approaches.