Heavy Metals in Bottom Sediments of Reservoirs in the Lowland Area of Western Poland: Concentrations, Distribution, Sources and Ecological Risk

Abstract

The paper presents the results of a study of heavy metals (HMs) concentrations in six retention reservoirs located in the lowland area of western Poland. The objectives of this study were to analyze the Cd, Cr, Cu, Ni, Pb and Zn concentrations, assess contamination and ecological risk, analyze the spatial variability of HM concentrations and identify potential sources and factors determining the concentration and spatial distribution. The bottom sediment pollution by HMs was assessed on the basis of the index of geo-accumulation (I_geo), enrichment factor (EF), pollution load index (PLI) and metal pollution index (MPI). To assess the ecological risk associated with multiple HMs, the mean probable effect concentration (PEC) quotient (Q_m-PEC) and the toxic risk index (TRI) were used. In order to determine the similarities and differences between sampling sites in regard to the HM concentration, cluster analysis (CA) was applied. Principal component analysis (PCA) was performed to assess the impact of grain size, total organic matter (TOM) content and sampling site location on HM spatial distribution. Additionally, PCA was used to assess the impact of catchment, reservoir characteristics and hydrological conditions. The values of I_geo, EF, MPI and PLI show that Cd, Cr, Cu, Ni and Pb mainly originate from geogenic sources. In contrast, Zn concentrations come from point sources related to agriculture. The mean PEC quotient (Q_m-PEC) and TRI value show that the greatest ecological risk occurred at the inlet to the reservoir and near the dam. The analysis showed that the HMs concentration depends on silt and sand content. However, the Pb, Cu, Cd and Zn concentrations are associated with TOM as well. The relationship between individual HMs and silt was stronger than with TOM. The PCA results indicate that HMs with the exception of Zn originate from geogenic sources—weathering of rock material. However, the Ni concentration may additionally depend on road traffic. The results show that a reservoir with more frequent water exchange has higher HMs concentrations, whereas the Zn concentration in bottom sediments is associated with agricultural point sources.

1. Introduction

Reservoirs play important roles in water supply, irrigation, power generation and flood control. However, dam construction changes the hydrological regime and sediment transport [1,2,3]. The suspended sediments are stored in the reservoir and with them pollutants including heavy metals (HMs) [4]. Bottom sediment plays an important role in monitoring the aquatic environment [5], especially evaluating the contamination levels [6,7], and ecological risk assessment [8].

HMs have natural, i.e., geogenic (weathering and erosion), and anthropogenic origin (urban domestic and industrial waste and the use of chemical fertilizers, roadside soil) [5,7,9,10,11]. They are classified as serious environmental pollutants due to persistence, tendency to bioaccumulate and high toxicity [12].

The spatial distribution in HM concentrations among reservoirs is associated with different anthropogenic activities [13,14]. The distribution depends on hydrodynamic conditions, type of sediment and metal sources [15]. Moreover, shape and morphology of the reservoir [16], reservoir operation [17] and the biochemical processes [18] modify the HM deposition. The highest concentrations have been detected generally in the vicinity of polluted water inflow [10]. River-reservoir interaction might modify the impact of HM point sources [19]. Intensified accumulation of sediments and pollutants has been observed especially in the vicinity of dams [17,20]. The concentrations of HMs in bottom sediments varied in different sampling locations and layers, but the concentrations in the surface layer are higher than those in the deeper layers [21]. Generally, the concentration of HMs was higher in reservoir sediments than in river sediments [16,22]. However, Frémion et al. [23] indicated that the HM concentrations are lower in the reservoir than upstream and downstream, which is caused by the high fresh organic matter inputs, diluting the contamination. Seasonal variations of metal concentrations depend on hydrological conditions, the impact of non-point sources and lithogenic sources [24].

The distribution of HMs in the sediments was mainly connected with grain size and organic matter content [25,26]. Fine-grain sediments contain higher HM concentrations [27], which is related to higher magnetic susceptibility [8]. Farhat and Aly [28] suggested that the organic matter was more important in controlling HM distribution. Organic matter was a key mediator for ecological risk [29]. Metals released into reservoirs are generally bound to the sediment bottom [30,31]. However, HMs may be released into the water column via sediment resuspension, due to changing chemical and hydrological conditions [23,32,33,34,35] and accumulate in plants and animals [2,26,36,37]. Higher concentrations of HMs can cause toxicity risk to the biota and human health [38,39,40,41]. Birch and Apostolatos [42] showed that anthropogenic metals have higher mobility and bioavailability than metals from geogenic origin.

Flow regulation, local geomorphological characteristics and reservoir operations determine the redistribution of HMs [39,43]. In dam reservoirs the water-level fluctuation zone is important for the accumulation and redistribution of HMs [29,44]. Moreover, dredging operations, the emptying of reservoirs and flood events might lead to release of HMs to the water environment [45,46,47,48].

The transport and enrichment of HMs in reservoir bottom sediments have been intensively studied, including concentration, spatial distribution, source identification and pollution assessment [13,49,50,51], temporal variability [52], bioavailability [53], ecological risk and organism toxicity [4,54,55]. Moreover, some studies relate to HM concentrations in reservoir tributaries [16,56]. Many methods have been used for the assessment of bottom sediment HM pollution and to determine the potential risk of heavy metal contamination [29,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71].

However, for the assessment of spatial distribution, identification of pollution sources and factors affecting their content in bottom sediments, multivariate statistical have been used [7,26,71,72,73,74,75,76].

Despite the numerous research results available, there are some deficiencies that still need attention. The studies usually focus on a certain reservoir, and less information is available about the HMs’ variation in a group of reservoirs [14]. There is little information available about heavy metal pollution in reservoirs and their changes after prolonged exploitation [77]. Moreover, Wu et al. [78] suggested that future climate change will aggravate the ecological risk of HMs in the water environment due to the release of HMs from sediments to the water environment.

Poland is characterized by high seasonal and spatial variability of water resources. Therefore, in order to increase the efficiency of water management and ensure flood protection, retention reservoirs have been constructed for over 50 years. The largest number of retention reservoirs were constructed in lowland areas of Poland. Their location in the agricultural landscape means that they are exposed to a high supply of nitrogen and phosphorus [79,80,81]. In the reservoirs, seasonal algal blooms and an overgrowing process are observed [82,83,84]. Studies conducted by Baran et al. [85], Wiatkowski [86] and Kasperek and Wiatkowski [87] show that the heavy metal concentrations in reservoir bottom sediments vary across Poland. Concentrations of heavy metals in bottom sediments depend on the catchment land use and the occurrence of pollution point sources [25,88,89]. Therefore, each reservoir should be considered individually in order to expose the factors affecting the bottom sediment pollution. In this study, six reservoirs located in the lowland area of western Poland were selected for analysis. Their selection was influenced by the time of their construction, different morphometric parameters, catchment land use and hydrological conditions. The objectives of this study are to: (1) analyze the concentrations of HMs in bottom sediments, (2) show the spatial variability of heavy metal concentrations, (3) analyze the bottom sediment contamination, (4) assess ecological risk, (5) identify potential sources and factors determining the content and spatial distribution of HMs in the reservoirs.

2. Materials and Methods

2.1. Study Area and Sampling Sites

The research was carried out on the example of six reservoirs located in lowland areas in western Poland (Figure 1).

The reservoirs are located in the areas most vulnerable to drought in Poland [90]. The main functions of the reservoirs are to provide water for irrigation purposes. The reservoirs were built between 1967 and 2013 (

Table 1. Basic parameters of the reservoirs and their catchments.

Parameter	Jeżewo	Jutrosin	Pakosław	Rydzyna	Środa	Września
Reservoir characteristics
Start of functioning—AG (year)	2003	2011	2006	2013	1971	1967
X coordinate of centroid	51°57′15.523”	51°39′51.876”	51°35′25.904”	51°48′15.419”	52°14′33.014”	52°20′27.316”
Y coordinate of centroid	17°14′9.603”	17°10′29.229”	17°3′35.704”	16°39′10.306”	17°17′38.652”	17°32′28.053”
Area—AR (ha)	73	91	54	41	39	33
Mean depth—MD (m) 1	2.9	2.7	1.9	2.3	2.3	0.9
Volume—VO (million m3)	2.1	2.4	1.0	1.0	0.9	0.3
Shoreline development—SD (-)	2.8	1.6	0.8	1.6	2.7	2.9
Hydrological conditions
Retention time—RT (d)	77	122	122	109	23	4
Catchment characteristics
Catchment area—CA (km2)	104.1	517.6	773.5	24.3	148.5	272.5
Mean elevation—ME (m a.s.l.)	126.3	128.7	123.9	104.5	104.5	114.8
Mean slope—MS (°)	0.4	0.4	0.5	0.4	0.5	0.4
Urban/Industry—UI (%)	4	6	5	7	2	2
Agriculture—A (%)	79	84	78	75	81	68
Forest—F (%)	17	10	17	18	17	30
Stream density—DD (km·km−2)	1.9	2.1	2.3	1.2	1.6	1.7
Road density—RD (km·km−2)	1.0	1.0	0.9	1.1	1.0	1.0
Road and river crossings—RRC (no.· km−2)	0.8	0.7	0.6	0.7	0.5	0.6

¹ Mean depth was calculated as ratio of volume and area of the reservoirs.

). Their surface area ranges from 33 to 91 ha and their volume from 0.3 to 2.4 million m³. They are shallow reservoirs with mean depths ranging from 1.3 to 2.1 m. The greatest depths from 3.5 to 5.5 m are found near the dam. The reservoirs Jeżewo, Jutrosin, Rydzyna, Środa and Września have an elongated shape and their shoreline development ranges from 1.6 to 2.9. However, the Pakosław reservoir has a rounded shape and the shoreline development is 0.8. In the catchment land use structure, agriculture area dominates, its percentage ranging from 68 to 81%. The proportion of urban and industrial areas is small and ranges from 2 to 7% in the Środa, Rydzyna and Jutrosin catchments, respectively. The stream and road density in the catchment range from 1.2 to 2.3 km·km⁻² and from 0.9 to 1.1 km·km⁻² respectively. The number of river and road crossings where HM inflow may occur ranges from 0.5 to 0.8 no.·km⁻².

2.2. Sample Collection and Preparation

The sampling and field surveys took place during September 2017. The total 28 surface sediment samples were collected from six reservoirs (Figure 1). The sampling locations were recorded using a hand-held Global Positioning System (GPS) device (GARMIN OREGON 300). Due to the elongated shape of the Jeżewo, Jutrosin, Radzyny, Środa and Września reservoirs, five samples were taken along the potential flow path. In the case of the Pakosław reservoir, with a rounded shape, three samples were taken. The top 0–5 cm depth of sediment samples was manually collected by means of Nurek and Czapla devices and carried within zip-mouthed PVC packages. In the laboratory, the bottom sediment samples were subjected to grain size and total organic matter (TOM) analysis. The analysis was carried out in order to examine the relationship between bottom sediment characteristics and total HM concentration. In the first step the bottom sediment samples were air dried at room temperature and sieved through a 2 mm nylon sieve to remove coarse debris. Next the samples were homogenized and passed through 2 mm, 1 mm, 0.5 mm, 0.25 mm, 0.10 mm and 0.05 mm sieves. The analysis was done by means of an analytical sieve shaker (Model, AS 200). The grain size fraction less than 0.05 mm was determined according to the Casagrande aerometric method modified by Prószyński [91]. The sand (0.05–2.0 mm), silt (0.002–0.05 mm), and clay (<0.002 mm) fractions were determined. The total organic matter (TOM) was measured by the Tiurin method with wet oxidation, followed by ferrous ammonium sulfate titration [91]. Finally, the samples were prepared for chemical analysis. The bottom sediment samples were extracted with hydrochloric acid (Merck, Darmstadt, Germany) at a 1:4 ratio at 95 °C ± 5 °C in a Mars 5 Xpress microwave digestion system (CEM, Matthews, NC, USA).

2.3. Chemical Analysis

The concentrations of Cd, Cr, Cu, Pb, Ni, and Zn in the bottom sediment samples with a fraction of <2 mm were determined by inductively coupled plasma mass spectrometry (ICP-QQQ), model 8800 Triple Quad (Agilent Technologies Tokyo, Japan). The analytical procedure is described in previous papers [16,73,75]. The reagents used were analytically pure, and the water was purified to the resistivity of 18.2 MΩ·cm (at 25 °C) in a Direct-Q 3 Ultrapure Water Systems apparatus (Millipore Molsheim, France). During the determinations by the ICP-QQQ techniques, standard solutions produced by VHG Labs, Inc. (Manchester, UK) were used. We measured the Certified Reference Material CRM no. LGC 6187 used for River sediments (Manchester, UK), and high compliance with reference values was found.

2.4. Quantification of Metal Pollution

The assessment of bottom sediment pollution by HMs was performed on the basis of four indices: the index of geo-accumulation (I_geo), the enrichment factor (EF) index, the pollution load index (PLI), and the metal pollution index (MPI). I_geo and EF are single indices, while PLI and MPI are integrated indices.

2.4.1. Index of Geo-Accumulation

The index of geo-accumulation I_geo was determined by the equation proposed by Müller [92] as follows (1):

$${\mathrm{I}}_{\mathrm{geo}}{=\log }_2\left(\frac{{\mathrm{C}}_{\mathrm{n}}}{1.5{\mathrm{B}}_{\mathrm{n}}}\right)$$

(1)

I_geo is a single factor pollution index where C_n and B_n are the concentration of the selected HMs in bottom sediments and the geochemical background value, respectively (mg.kg⁻¹). The factor 1.5 is used for the possible variability of the background data due to lithological conditions. Values of I_geo were used to classify sediment samples into 7 purity classes: class 0 (I_geo ≤ 0)—uncontaminated, class 1 (0 < I_geo ≤ 1)—uncontaminated to moderately contaminated, class 2 (1 < I_geo ≤ 2)—moderately contaminated, class 3 (2 < I_geo ≤ 3)—moderately to heavily contaminated, class 4 (3 < I_geo ≤ 4)—highly contaminated, class 5 (4 < I_geo ≤ 5)—heavily to extremely contaminated, class 6 (5 < I_geo)—extremely contaminated [64].

2.4.2. Enrichment Factor

The enrichment factor (EF) is an efficient index for quantification of the human impact on the concentration of a given metal [93]. The EF is based on a ratio of measured metal concentration to the metal geochemical background value. To account for natural heavy metal concentrations, EF is normalized to sediment Al or Fe content. In this study the normalization with Fe was selected, as proposed by Deely and Fergusson [94]. The EF was calculated as follows [95] (2):

$${\mathrm{EF}}_{\mathrm{n}}=\frac{\frac{{\mathrm{C}}_{\mathrm{n}}}{{\mathrm{C}}_{\mathrm{Fe}}}}{\frac{{\mathrm{B}}_{\mathrm{n}}}{{\mathrm{B}}_{\mathrm{Fe}}}}$$

(2)

where (C_n/C_Fe) is the ratio of the concentration of heavy metal (C_n) to iron (C_Fe) in the sediment sample and (B_n/B_Fe) is the same ratio for the geochemical background [96,97]. Values of EF were used to assess the pollution of bottom sediment samples into the following classes: 0 (EF ≤ 1) no enrichment; 1 (1 < EF ≤ 3) is minor enrichment; 2 (3 < EF ≤ 5) is moderate enrichment; 3 (5 < EF ≤ 10) is moderately severe enrichment; 4 (10 < EF ≤ 25) is severe enrichment; 5 (25 < EF ≤ 50) is very severe enrichment; and 6 (EF > 50) is extremely severe enrichment [98,99].

2.4.3. Pollution Load Index

The pollution load index (PLI) is an experimental formula developed by Tomlinson et al. [100] (3):

$$\mathrm{PLI}={\left(\frac{{\mathrm{C}}_{\mathrm{n}1}}{{\mathrm{B}}_{\mathrm{n}1}}·\frac{{\mathrm{C}}_{\mathrm{n}2}}{{\mathrm{B}}_{\mathrm{n}2}}·\dots ·\frac{{\mathrm{C}}_{\mathrm{ni}}}{{\mathrm{B}}_{\mathrm{ni}}}\right)}^{\frac{1}{\mathrm{i}}},$$

(3)

where C_n/B_n is the ratio of selected elements to their background values and i is the number of elements. The empirical index provides simple comparisons of site average heavy metal pollution. PLI value = 0 denotes perfection, PLI < 1 denotes no pollution, and PLI > 1 indicates pollution [101].

2.4.4. Metal Pollution Index

The metal pollution index (MPI) was also used to compare the total metal content between different sampling sites within reservoirs. The MPI was calculated as follows (4):

$$\mathrm{MPI}={\left({\mathrm{C}}_{\mathrm{n}1}·{\mathrm{C}}_{\mathrm{n}2}·\dots ·{\mathrm{C}}_{\mathrm{ni}}\right)}^{\frac{1}{\mathrm{i}}}$$

(4)

where C_n are concentrations of HMs in the sample and i is the number of analyzed metals [96,102]. The background heavy metal concentrations in the study area were provided by Bojakowska and Sokołowska [103]. The geochemical background values for bottom sediments were as follows 1.0 for Cd, 10.0 for Cu, 10.0 for Cr, 10.0 for Ni, 25.0 for Pb, 100.0 for Zn and 10,000.0 mg^·kg⁻¹ for Fe.

2.5. Ecological Risk Assessment

A consensus-based sediment quality guideline (SQG) was introduced by Macdonald et al. [104] to predict the toxicity of a sediment samples. For each chemical contaminant there are two consensus-based values—threshold effect concentration (TEC) and probable effect concentration (PEC). The concentration of chemical contaminant below the TEC means that undesirable effects are not expected to occur, while above the PEC adverse effects are expected to occur more often than not [39,105,106]. To assess the effects of multiple HMs, the mean PEC quotient (Q_m-PEC) was proposed, which can be calculated using Equation (5):

$${\mathrm{Q}}_{\mathrm{m}-\mathrm{PEC}}=\overline{{\mathsf{\Sigma }}_{\mathrm{i}=1}^{\mathrm{n}}\frac{{\mathrm{C}}_{\mathrm{n}}}{{\mathrm{PEC}}_{\mathrm{n}}}},$$

(5)

where n is the number of HMs, Cn is the measured concentration of a heavy metal and PEC_n is the corresponding PEC. The PEC benchmark values for Cd, Cr, Cu, Ni, Pb and Zn are 4.98, 111, 149, 48.6, 128, and 459 mg·kg⁻¹, respectively. The categories for Q_m-PEC are classified as: not toxic (Q_m-PEC < 0.5) and toxic (Q_m-PEC > 0.5).

The toxic risk index (TRI) was developed by Zhang et al. [107] to assess the toxic risk for certain HMs in sediment samples. The TRI is based on threshold (TEL) and probable (PEL) level effects [108]; it can be calculated using Equation (6):

$${\mathrm{TRI}}_{\mathrm{i}}=\sqrt{\frac{{\left(\frac{{\mathrm{C}}_{\mathrm{i}}}{\mathrm{TEL}}\right)}^2+(\frac{{\mathrm{C}}_{\mathrm{i}}}{\mathrm{PEL}}{)}^2}{2}}$$

(6)

To assess integrated risk of multiple HMs, the TRI can be calculated using Equation (7):

$$\mathrm{TRI}={\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{n}}}{\mathrm{TRI}}_{\mathrm{i}}$$

(7)

where n is the number of HMs, C_i is the content of each heavy metal in the sediment sample, and TRI_i is the toxic risk index of each heavy metal. Five categories of the TRI are classified as: no toxic risk (TRI ≤ 5), low (5 < TRI ≤ 10), moderate (10 < TRI ≤ 15), considerable (15 < TRI ≤ 20), and very high (TRI > 20) [43].

2.6. Distribution and Source Identification

In order to determine the similarities and differences between sampling sites in regard to the HM concentration, cluster analysis (CA) was applied. The CA was carried out using the Ward method with square Euclidean distance as a measure of similarity between clusters. Groups and subgroups were distinguished using the cut-off criteria of 66% and 25% respectively [109]. For the purpose of comparison of HM concentrations between sample sites, the non-parametric Kruskal-Wallis test (K-W, p ≤ 0.05) followed by Dunn’s test as a post-hoc procedure was used. The K-W test was used for verification of the hypothesis of significance of differences between concentrations of HMs in distinguished groups.

Subsequently, to choose the appropriate ordination methods for future analysis, detrended component analysis (DCA) was conducted. The obtained gradients of the first axis of the DCA were shorter than 3.0 SD, which indicates that the HMs were in linear distribution. Glińska-Lewczuk et al. [110] recommended in such cases principal component analysis (PCA). PCA was employed to identify a potential source of HMs and factors that have the greatest impact on their concentration and spatial distribution. The PCA analysis was done in two variants. First the impact of grain size, total organic matter content and sampling site location on HMs’ spatial distribution was evaluated. The percentage content of sand (SA), silt (SI), clay (CL) and total organic matter (TOM) was calculated for each sample. In addition, the distance between sampling site location with respect to inlet (ID) and outlet (OD) of the reservoir was calculated. Second, the impact of other factors was analyzed. The following factors were distinguished: catchment characteristics (catchment area—CA, mean slope—MS, mean elevation—ME, urban and industrial—UI, agriculture—A, forest—F, stream density—DD, road density—RD, number of road and river crossings—RRC), reservoir characteristics (age—AG, area—AR, mean depth—MD, volume—VO, shoreline development—SD), and hydrological conditions (retention time—RT). The above factors were referred to the mean HM concentration in the samples located in the reservoir (J, Ju, P, R, S, W). Before starting the PCA analysis the factors were tested in the context of the occurrence of outliers. The analysis was performed using the two-sided Grubbs test at the significance level of 0.05. Subsequently the HM concentration and factors that may have an impact on HM content and distribution were log-transformed to obtain a normal distribution. In order to avoid misclassification due to differences in data dimensionality, multivariate statistical techniques were applied on standardized data through z-scale transformation [73,74]. The significant principal components were selected based on a Kaiser criterion with eigenvalues higher than 1. The correlation between principal components and analyzed data was classified according to the Liu et al. [111] criterion where values >0.75, 0.75–0.50, and 0.50–0.30 show a strong, moderate and weak relationship, respectively. CA and PCA were conducted using Statistica 13.1 and Canoco 5.0 respectively.

3. Results

3.1. Heavy Metal Content

The characteristic concentrations of HMs (Cd, Cr, Cu, Ni, Pb and Zn) in bottom sediments from 28 sites are shown in

Table 2. HM concentrations in bottom sediments of reservoirs (mg·kg⁻¹) ¹.

Reservoir	Cd	Cr	Cu	Ni	Pb	Zn
Jeżewo	0.18-0.72 0.4	2.79-11.0 6.5	4.80-15.7 10.1	2.68-9.66 5.9	9.29-26.1 17.6	143.6-1639.6 903.7
Jutrosin	0.08-0.29 0.2	0.60-5.01 3.1	0.45-7.48 3.6	0.38-6.61 3.7	3.12-13.6 6.2	6.08-60.1 23.1
Pakosław	0.01-0.14 0.1	0.66-4.66 2.0	0.33-4.67 2.0	0.44-4.80 2.0	0.91-5.28 2.6	8.23-510.7 221.7
Rydzyna	0.06-0.09 0.1	1.79-5.48 3.5	2.22-6.27 4.2	2.11-6.42 4.4	3.05-9.23 5.9	42.2-1355.9 436.5
Środa	0.05-0.30 0.2	1.73-9.25 4.8	1.36-8.73 4.6	1.01-5.99 3.5	2.81-13.0 7.4	50.8-1131.7 357.5
Września	0.19-0.71 0.4	3.91-9.32 6.1	3.74-18.7 9.4	2.47-9.82 5.5	7.04-32.8 15.2	27.8-1990.4 678.4

¹ Upper values: minimum-maximum. Lower values: mean.

. Generally, the mean concentrations of HMs in reservoirs follow a descending order of Zn > Pb > Cu > Cr > Ni > Cd. The Zn average concentration in reservoirs’ bottom sediments was 452.2 mg·kg⁻¹. The lowest Zn concentrations occurred in the Jutrosin and the highest in the Jeżewo reservoir, 23.1 and 903.7 mg·kg⁻¹ respectively. It was observed that in the Jeżewo, Rydzyna, Środa and Września reservoirs, Zn concentrations higher than 1000 mg·kg⁻¹ were observed at the inlet to the reservoir or near the dam. Pb concentrations at the highest level were in the Jeżewo reservoir (average 17.6 mg·kg⁻¹), while the lowest values were observed in Pakosław (average 2.6 mg·kg⁻¹). A similar tendency was observed in the case of Cu, Cr and Ni concentrations. The average concentrations of Cu, Cr and Ni in the Pakosław reservoir were at a similar level and was on average 2.0 mg·kg⁻¹. In the Jeżewo reservoir, the average concentrations of these elements in bottom sediments were 3 to 5 times higher. The Cd concentrations in the reservoirs were the lowest and ranged from 0.01 to 0.72 mg·kg⁻¹. In the Września reservoir, the average concentration of Cd was about seven times higher than in the Pakosław reservoir. It was observed that Zn concentrations are characterized by the highest variability between the sampling sites and reservoirs; however, Cr concentrations have the lowest variability. The non-parametric Kruskal-Wallis test showed significant differences in Cd, Cu, Zn and Pb concentrations between bottom sediments at individual reservoirs. Differences were not significant only for Cr and Ni concentrations.

3.2. Bottom Sediment Pollution Assessment

The I_geo was applied to assess the HM contamination level in bottom sediments. The I_geo values for Cd, Cr, Pb and Ni were below zero and could be classified as unpolluted to moderately polluted. The I_geo values for Cu were also lower than 0 except for two sampling sites, J1 and W3. At these points, I_geo values were slightly higher than 0, which indicates a small local contamination of bottom sediments with Cu. The largest contamination of Zn was noted in the reservoirs. I_geo values for Zn lower than 0 occurred only for the Jutrosin reservoir. In 23 measurement sampling sites located in other reservoirs, I_geo values higher than 0 were recorded fourteen times. The sediment bottom at sample sites R2, R3 and S2 demonstrates weak contamination. The values of I_geo were slightly higher than 0. The sediments in sampling sites P5, R5 and S4 demonstrate moderate contamination with the value of I_geo higher than 1. The I_geo for sampling sites J1, J2, J3, J5, R1, S5, W1 and W5 was between 2 and 4, which demonstrated moderate to severe pollution of bottom sediments (Figure 2a).

The enrichment factor (EF) was applied to assess the sources of HMs. The results of the EF calculations for Cd, Cr, Cu, Pb and Ni showed that the values generally did not exceed a value of 1 (Figure 2b). EF values ranging from 1 to 3 for Pb occurred once, for Cd and Cr twice, for Ni four times and for Cu nine times. In most cases, the values were less than 1.5, which may suggest that the HMs may come from natural weathering processes. In the case of Cr, at points S2 and W5 and for Cu at points S5 and W5, the EF values were higher than 1.5. However, an EF value higher than 1.5 indicates that HMs were delivered by other sources, such as point and non-point pollution [112,113]. The EF value of Zn was highest, and it ranged from 0.22 to 43.52. In fourteen sampling sites the EF values were higher than 3. Only in the Jutrosin reservoir were the EF values lower than 1. EF values for Zn higher than 10 were observed in the Jeżewo, Pakosław, Rydzyna, Środa and Września reservoirs, generally at sampling sites located at the inlet to the reservoir and at the dam. Two global indices of PLI and MPI were used to assess the bottom sediments’ contamination with HMs. PLI values ranged from 0.06 to 1.54. PLI values higher than 1 were noted four times at points J1, J5, W1 and W3 (Figure 3). These values may indicate the inflow of HMs from anthropogenic sources. The MPI index values ranged from 0.72 to 17.97. The highest MPI values occurred in the Jeżewo and Września reservoirs, on average 11.19 and 8.75 respectively. The lowest MPI values occurred in the Jutrosin and Pakosław reservoirs, 2.95 and 2.45 respectively (Figure 3). The highest values of the PLI and MPI indices were recorded at sampling sites located at the inlet to the reservoir and near the dam.

3.3. Ecological Risk Assessment

The toxic risk index was used to assess the toxic risk for HMs in bottom sediment samples. The values of the TRI index ranged from 0.17 to 13.56. In general, the TRI values were lower than 5 in 20 out of 28 sampling sites, which indicates a lack of ecological risk associated with HM concentrations in bottom sediments. Values lower than 5 occurred at all sampling sites located in Jutrosin reservoirs. TRI values ranging from 5 to 10 occurred six times (J1, J2, J3, R1, S5 and W5), which indicated low ecological risk. Only in sampling sites J5 and W1 were the TRI values higher than 10, which indicated moderate ecological risk. The highest contribution in the TRI index were Zn concentrations (Figure 4a), on average 59%. The average TRI value for the Jutrosin reservoir was the lowest at 0.64, while more than ten times higher values were found in the Jeżewo reservoir. The reservoirs according to TRI values can be arranged in ascending order Jutrosin < Pakosław < Środa < Rydzyna < Września < Jeżewo.

For the toxicity assessment of combined HMs, a consensus-based sediment quality guideline (SQG) was used. The Q_m-PEC values ranged from 0.07 to 4.88. The average Q_m-PEC values of the reservoirs ranked as follows: Jeżewo > Września > Rydzyna > Środa > Pakosław > Jutrosin. Q_m-PEC values higher than 0.5 were noted up to seventeen times, which may indicate that HMs can be toxic to certain sediment-dwelling organisms [104]. Values higher than 0.5 occurred at all sampling sites located in the Jeżewo reservoir, four times in the Rydzyna reservoir and three times in the Środa and Września reservoirs (Figure 4b). The contribution of Zn was the highest in the Q_m-PEC index of the analyzed sampling sites. However, at the sampling sites located in the Jutrosin reservoir and points W2, W3 and W4 (Września Reservoir), R4 (Rydzyna Reservoir) and S3 (Środa Reservoir), the percentages of HMs were at a similar level. It was observed that in the reservoirs, the greatest ecological risk associated with the pollution of HM deposits occurred at the inlet to the reservoir and near the dam.

3.4. Distribution and Source Identification

The CA analysis allowed the sampling sites to be divided into two groups, A and B, depending on the HM concentrations (Figure 5a). A total of 19 sampling sites were classified in group A and 9 in group B. The concentrations of all HMs in group A were lower than those observed in group B. In group B there are four measurement points located in the Jeżewo reservoir, three in the Września reservoir and two in the Środa reservoir. Group B was characterized by low variability of concentrations of individual HMs. Higher variability occurred in group A. Therefore, this group was divided into two subgroups: A1 and A2. Four sampling sites—S1, P1, Ju1 and P3—were classified in subgroup A1, in which the lowest concentrations of Cr, Ni, Cu and Pb were observed. However, in the case of Cd and Zn concentrations the differences between subgroups A1 and A2 were not statistically significant. Another division of the sampling sites into groups was obtained depending on the sediments’ bottom texture and TOM (Figure 5b). Samples were divided into two groups, I and II, which included 3 and 25 sampling sites, respectively. Sampling sites in group I were characterized by the highest content of TOM and silt fraction, on average 8.5% and 49.7%. At these sampling sites, the lowest sand fraction was observed, on average 47.0%. Within group II, 3 subgroups—IIa, IIb and IIc—were separated, into which 9, 3 and 13 sampling sites were classified respectively. In subgroup IIa there were sampling sites with the highest sand content (mean 92.9%) and the lowest content of silt and clay (average 6.7% and 0.4%). Also, the content of TOM was at the lowest level: 0.4%. Subgroup IIb was characterized by the highest content of clay (average 7.7%) and high content of sand (84.7%). In subgroup IIc the clay content was at a low level (on average 1.2%) and silt content at a high level (on average 21.8%). Dunn’s test showed that in groups I and II there are significant differences between sand, silt, clay and TOM content. Differences in all analyzed parameters were also present between subgroups IIa and IIc. However, between subgroups IIa and IIb there were only differences in the content of sand and clay and between subgroups IIb and IIc in the content of silt and clay. It was observed that the sample sites with the highest content of HMs (group B—Figure 5a) belonged to group I and subgroup IIc (Figure 5b). In these groups, high content of silt and TOM in group I was recorded.

3.5. Relationships between HM Concentrations and Sediment Texture and Organic Matter Content

In the first stage, the aim of the PCA analysis was to identify factors that determine the concentration of HMs in the bottom sediments. The PCA analysis allowed two significant factors to be distinguished, PC1 and PC2, whose eigenvalues were higher than 1. The PC1 explains as much as 75.9% of the variance. The PC1 was strongly negatively correlated with Cd, Cr, Cu, Ni and Pb concentrations (Figure 6).

PC2, explaining about 11.7% of the variance, was strongly negatively correlated with Zn concentrations. The projection of sampling sites perpendicular to individual HM vectors or their extension beyond the origin of the coordinate system allows estimation of their content in bottom sediments. For example, the highest Zn values occur at sampling sites W1, J5, J1, S5, W5 and R1 and the lowest at sampling sites Ju3, Ju2, Ju1 and P1. The PC1 factor was moderately positively correlated with sand content and negatively with silt and TOM content. The angle created by the vectors of individual HMs informs about the relationship between them. The smaller the angle is, the stronger is the relationship; however, the angle value of 90° means no correlation. The analysis showed that the concentrations of HMs in bottom sediments are positively correlated with the silt content and negatively with the sand content (except Zn), while no correlations were observed with clay. In addition, the studies showed a positive correlation between Pb, Cu, Cd and Zn concentrations and TOM content. However, the relationship between individual HMs and silt was stronger than with TOM. A positive correlation was observed between the sampling site distance from the inlet to the reservoir (ID) and the silt and TOM content and a negative correlation with the sand content. This indicates grain size segregation in reservoirs. Thicker, sandy fractions are deposited at the inlet of the reservoirs, while the silt and TOM fractions are deposited near the dam.

The second objective of the PCA was the identification of external factors that may affect the concentrations of HMs in reservoirs. The PCA analysis made it possible to distinguish two factors, PC1 and PC2, which explain 79.9 and 14.5% of variance respectively. The average concentrations of Cd, Cr, Cu, Ni and Pb were strongly negatively correlated with PC1, whereas Zn concentrations were strongly correlated with PC2 (Figure 7). Among the other factors that may have an impact on the content of HMs in bottom sediments, PC1 was strongly positively correlated with mean catchment slope and strongly negatively correlated with reservoir shoreline development. In addition, a positive correlation was observed with the catchment area and the water retention time in the reservoir. Only the share of agricultural areas was negatively correlated with PC2. The results indicate that HMs with the exception of Zn have the same source. Low concentrations in bottom sediments indicate that they originate from geogenic sources—weathering of rock material. The positive correlation with reservoir shoreline development indicates a possible impact of adjacent areas. The content of Ni in bottom sediments may additionally result from road traffic, which confirms the correlation with number of road and river crossings. The negative correlation of Cd, Cr, Cu, Pb and Ni with retention time indicates that in the reservoir with more frequent water exchange, HM concentrations are higher, which is related to the suspended sediment deposition, whereas the Zn in bottom sediments is correlated with agricultural land use. High concentrations of Zn in bottom sediments may result from the existence of anthropogenic sources and a potential supply with domestic wastewater. The data provided by the Central Statistical Office in Poland indicate that the number of people using the water supply network in the reservoir catchment varied from 93 to 98%. On the other hand, the number of people using the sewerage network varied from 36 to 65% for the Rydzyna and Września catchment respectively. The ratio of the length of the water supply to the length of the sewage system varied from 0.21 to 0.35. The above data indicate possible supply of Zn with domestic wastewater.

4. Discussion

The HM content in bottom sediments of retention reservoirs and relations between them result mainly from the impact of anthropogenic sources [14] as well as sediment texture and organic matter content [26]. Meteorological, hydrological and geological conditions of the catchment affect the content of sand, silt, clay and TOM in water. These factors influence the transport of HMs in the river-reservoir system. In the world, bottom sediments of retention reservoirs are characterized by relatively high variability of HM concentration. The analyzed reservoirs are located in agricultural catchment; therefore, HM concentrations in bottom sediments are at a lower level. Only Zn concentrations are at a higher level with respect to those noted by other authors in Polish reservoirs [25]. Probably Zn content may result from potential supply with domestic wastewater. The results presented in this paper confirmed the results obtained by Haziak et al. [114] which indicated that Zn values exceeding the geochemical background were associated with anthropogenic activities (run-off from farmland and domestic wastewater).

High values of the I_geo, EF, PLI and MPI indices at individual sampling sites of the reservoirs may suggest their origin of HMs from anthropogenic sources. High values at the inlet to the reservoir result from their inflow from point sources. Bing et al. [43] suggested that local anthropogenic activities increase contamination of specific HMs in the bottom sediments, especially of upper regions. Saleem et al. [55] observed relatively high content of HMs at sites which were adjacent to urban and semi-urban areas. The HMs may also come from the discharges of untreated urban/industrial wastes, agricultural runoffs and automobile emissions. High HM concentration near the dam is a consequence of changes in hydrological conditions in reservoirs during floods and high flow periods. HMs absorbed on silt and clay and TOM are transported and deposited near the dam. Also, Bing et al. [43] and Frémion et al. [17] reported relatively high values of HM concentrations close to the dam. Sojka et al. [25] noted that near the outlet pipe from the bottom sediments there can be discharged silt, clay and TOM, which are responsible for the transport of HM downstream from the reservoir.

The results have shown that concentrations of Cd, Cr, Cu, Ni, and Pb were strongly positively correlated with the silt and TOM contents in bottom sediments. However, the relationship with silt was at a higher level. Palma et al. [52] obtained slightly different results which showed stronger relationships between HMs and silt and clay content. However, the relationship between HMs and silt was stronger. Farhat and Aly [28] suggest that organic matter was more critical than grain size in controlling HM distribution in sediment. Different results were obtained by Frémion et al. [23], who found no relationships between trace elements and organic fraction.

The content in the bottom sediments of HMs involves the possibility of ecological risks. It was observed that the ecological risk in the reservoirs occurred at the inlet and near the dam, where the highest Q_m-PEC and TRI values were obtained. Slightly different results were obtained by Bing et al. [43], who found that the potential ecological risk level increased towards the dam. Results obtained by Zhao et al. [50] suggested that the concentration and potential ecological risk of HMs in sediments near the dam were higher compared to the upstream sampling sites. The probable reason for the very high HM content at the reservoir’s inlet was the date of sampling during the very low flow period. The obtained values of Q_m-PEC and TRI were the most strongly correlated with TOM. Our research confirmed the results obtained by Lin et al. [29] that the organic matter-bound fraction of heavy metal was a key mediator for ecological risk. Saleem et al. [55] demonstrated that organic matter may retain heavy elements in sediments and play an important role in bottom sediment quality assessment.

The present results showed that there is a strong correlation between Cd, Cr, Cu, Ni, and Pb. Suresh et al. [8] and Wang et al. [7] indicate that the HMs have common geochemical behaviors and originated from similar pollution sources when the correlation coefficient between them is higher. The absence of a correlation among the HMs suggests that the contents of these metals are not controlled by a single factor. Low concentrations of Cd, Cr, Cu, Ni, and Pb indicate their origin mainly from natural sources. In addition, during the research, the impact of the areas immediately adjacent to the reservoirs is marked; moreover, in the case of Ni the number of road and river crossings (RRC) plays an important role. Studies have shown that reservoirs with a shorter water retention time are more likely to supply HMs. In the case of Zn, it was shown that it may come from anthropogenic sources related to wastewater management in the reservoirs’ catchment.

5. Conclusions

The conducted analyses led us to draw the following conclusions:

• The content of HMs in bottom sediments is an individual feature of reservoirs.
• The PCA analysis and the values of I_geo, EF, MPI and PLI indices show that Cd, Cr, Cu, Ni and Pb in bottom sediments originate from geogenic sources – weathering of rock material. In contrast, Zn comes from anthropogenic sources.
• High variability of Zn concentrations between sampling sites and reservoirs confirms their origin from anthropogenic sources related to wastewater management.
• The PCA analysis indicates that the areas adjacent to reservoirs may have an impact on HM distribution, while the Ni concentration may additionally be affected by road traffic.
• In reservoirs with higher frequency of water exchange, higher HM concentrations were observed.
• The highest concentrations of HMs are observed at the inlet to the reservoir and near the dam, which causes the greatest ecological risk.
• The CA and PCA analysis show that the concentrations and spatial distribution of HMs depend on silt content. The Pb, Cu, Cd and Zn concentrations in reservoirs’ bottom sediments are associated with TOM content. However, the relation between individual HMs and silt was stronger than with TOM.
• Segregation of sediments along the reservoirs and the TOM content affect the concentrations and spatial distribution of HMs.