Integrated Ecological Risk Assessment of the Agricultural Area under a High Anthropopressure Based on Chemical, Ecotoxicological and Ecological Indicators

Abstract

Agricultural land is often located close to highly urbanised/industrialised areas and is subject to continuous anthropogenic pressure associated with the emission of many pollutants, ultimately deposited in the soil. Most studies on ecological risk assessment have only analysed the total contaminants’ concentration, which does not reflect their bioavailability or toxicity and often leads to an overestimation of risk. Therefore, in our study, we used an interdisciplinary approach, whereby the final conclusions about the risk in a given area are based on the integration of detailed data from chemical, ecotoxicological and ecological analysis. The research was carried out on agricultural land exposed to high levels of anthropopression for more than 100 years. Chemical measurements comprised both the total and bioavailable PAH content. A battery of bio-assays describing effects on soil retention and habitat function was used for ecotoxicity testing, and ecological indicators included enzymatic activity, respiration, microbial biomass, carbon mineralisation and nitrification. The integrated IntRisk index ranged from 0.19 to 0.94, and this was mainly due to high values of the chemical risk index, while the ecotoxicological and ecological results indicated no or low risk. The majority of the area (almost 90%) had acceptable risk levels, no/low risk (IntRisk < 0.5) at 57% of the sites and medium risk at 28% of the area. Very high unacceptable risk (IntRisk 0.77–0.94) was only at three sampling sites. The integration of data from a set of 15 indicators allowed us to derive quantitative risk indexes and delineate the limited area which needs additional action.

1. Introduction

Soils are fundamental and non-renewable natural resources providing goods and services essential for ecosystems and human life. The condition of soils is significantly affected by ongoing anthropogenic pressures associated with agricultural (use of organic and mineral fertilisers, pesticides) and non-agricultural activities (industrial emissions), which result in the release of a range of contaminants (e.g., polycyclic aromatic hydrocarbons, PAHs) that ultimately reach the soil. This may affect soil biodiversity, thereby reducing resistance to stress, influencing soil retention and habitat functions and, finally, deteriorating soil quality, fertility and productivity [1]. The Intergovernmental Technical Panel on Soil has identified soil contamination as one of the most worrying threats to soil health and the provision of key soil ecosystem functions and services [2]. Therefore, the protection of agricultural soils is one of the most important issues in current European policy, as reflected by many soil-related initiatives and strategies, such as the European Green Deal, the Biodiversity Strategy and the new European Union Soil Strategy. Given the objectives of the Green Deal, it is crucial to identify contaminated sites, and the final decision to remediate them should be based on site-specific risk assessment procedures to determine when risk is unacceptable.

A lot of studies on ecological risk assessment are focused on the prediction of potential risk based on the chemical analysis of the total pollutants concentrations and chemical risk indexes derived from such data [3,4,5,6,7]. Various indexes were used, e.g., the hazard quotient and hazard index [6,7], contamination factor and pollution load index [7], toxic units for individual PAHs and PAHs mixtures [6,7] and toxic equivalent concentration [5,6,7]. In each of these cases, the total contaminant content was compared with the generic reference values from the applicable soil quality standards. The final result of the assessment often depended on the type of chemical index [6,7], the model used for the calculation [6] and the assessment criterion adopted [5,6]. This traditional approach is not sufficient for fully assessing the response of soil organisms, does not provide information on the bioavailability of pollutants and their actual toxicity to soil biota and often leads to an overestimation or underestimation of risk [8,9]. Therefore, chemical analyses should be complemented by the results of ecotoxicological and ecological studies, preferably using organisms representing different groups of biological organisation.

Triad is one of the weight-of-evidence methods in which conclusions are based on the quantitative analysis of three independent assessment tools (lines of evidence, LoEs), i.e., chemical measurements of pollutant concentrations, results of ecotoxicological testing and ecological indicators [10]. Combining the three LoEs assessments is necessary because no single line of evidence provides comprehensive information, whereas the combined efforts and assessment of results from three independent disciplines provide a pragmatic reduction in the conceptual uncertainty associated with risk assessment. So far, the Triad method has mainly been used to assess the risk of sites contaminated by metals [11,12,13,14,15,16,17,18,19,20]. There are almost no reports on other groups of contaminants; only a few studies on the use of this method for the risk assessment of soils contaminated with organic pollutants can be found in the literature [21,22]. This method usually consists of several tiers and is recommended for site-specific risk assessment [10,23]. In practice, Tier 1 comprises screening studies to determine whether further stages of ecological risk assessment are required at a site. The screening phase is often based on a generic assessment, i.e., a simple comparison of environmental concentrations of pollutants with available threshold values [11,21,22], as well as on simple biotests, which, at this stage, should be quick, sensitive, easy to perform and inexpensive [12,22,23]. Higher tiers of this assessment (refined and detailed assessment), levels 2 and 3, are carried out to verify and reduce uncertainty and should take into account a wider set of more relevant and sophisticated indicators.

This research is the next step forward in the ecological risk assessment study focusing on the ecotoxicological and ecological effects in the agricultural area long-time-contaminated with organic pollutants (PAHs). Our previous studies characterised the contamination status of the area [22,24], a screening risk assessment based on chemical indexes [6] and simple bioassays implementation [22]. The aim of the current study was the deployment of a detailed integrated risk assessment procedure including a multidisciplinary approach based on chemical, ecotoxicological and ecological lines of evidence. We hypothesised that organic pollutants present in soils at high concentrations may deteriorate their retention and habitat functions and that the inclusion of parameters characterising the availability of pollutants in the assessment increases the reliability of the procedure. To the best of our knowledge, this study presents the first application of a higher-tier Triad approach to the ecological risk assessment of agricultural land with a long history of PAH contamination.

2. Materials and Methods

2.1. Study Area and Sampling Campaigns

The study was carried out in the area located in the south-western part of the Silesian region in Czerwionka, Rybnik district, Poland. The territory is very diverse; it includes agricultural and urban parts as well as industrial and post-industrial sites (Figure 1).

The high anthropopressure in this region results from diverse diffuse sources: the long-time coke production, former coal mining activity, fossil fuel combustion, waste recovery and road transport. Moreover, high values were found for anthropogenic indices such as total dust emissions (6827 kg year⁻¹ km⁻²) and dust emissions from industrial sources (1084 kg year⁻¹ km⁻²). The detailed characteristics of the region can be found elsewhere [22,24]. The climate in the region is influenced by three colliding subtropical, Arctic and continental air masses. The mean annual precipitation and temperature during the sampling period were 719 mm and 8.9 °C, respectively. The wind distribution was dominated by south-westerly winds with an average speed of 3.8 m s⁻¹. Soil samples were collected from agricultural land after the growing season during two sampling campaigns; initially, soil samples were collected to characterise the contamination status of the area (115 km²; one sampling site per ca. 4.8 km²) and to conduct a screening risk assessment [6,22]. Tier 1 assessment [22] showed that the acceptable risk was exceeded, or the uncertainty was too high in 46% of investigated soils, and it was concluded that, for 11 sampling sites, the risk analysis needed to be refined with a broader battery of tests. Therefore, the next phase—the second campaign (current research)—focused on a region of possible high risk limited to the northern and central part of the Czerwionka region (38 km²). Finally, 31 soil samples (1 sample per 1.2 km²) from the 0–0.3 m layer were collected after the growing season; each sample (approximately 5 kg) consisted of 6 subsamples taken from a 1 m² area. After initial manual mixing and homogenisation in the field, the samples were divided in the laboratory. Part of the sample was then air-dried and prepared for chemical analyses. Directly after collection, fresh soil samples were sieved through a 2 mm sieve, and ecotoxicological and ecological analyses were successively performed. If analyses could not be carried out immediately, the samples were stored under ISO 10381-6 standard conditions (temperature 4 ± 2 °C, in the dark with access to air) for a maximum of 7 days.

2.2. Soil Characteristic and Chemical Analysis

Soils were characterised in terms of texture (laser diffraction method), total organic carbon content (TOC, sulfochromic oxidation method), total carbon and nitrogen content (TC and TN, dry combustion) and pH in 1 mol L⁻¹ KCl (potentiometric method) (

Table 1. Summary statistics for the selected physical-chemical properties and PAH concentration (n = 31).

Indicator	Mean	Median	SD	Min	Max	CV%
Sand (%)	72.2	70.7	7.4	58.1	90.6	10
Silt (%)	25.6	26.8	6.8	9.2	39.0	27
Clay (%)	2.1	2.06	0.7	0.2	3.5	33
TOC (g kg−1)	17.4	11.5	20.0	5.2	106.9	115
TC (g kg−1)	24.9	14.6	32.7	6.3	155.0	131
pHKCl	5.4	5.4	0.9	3.9	7.2	16
TN (g kg−1)	1.0	0.7	0.9	0.2	4.1	88
16PAH bioav. (µg kg−1)	97.1	33.6	239.6	22.6	1183.8	247
16PAH (µg kg−1) a	12,620	1412	56,461	311	316,085	447

^a Adapted from Smreczak [26] and Klimkowicz-Pawlas and Debaene [6].

). The total content of 16 PAHs (16PAH) was determined by gas chromatography–mass spectrometry (GC-MS) using the Agilent GC-MS system (Agilent Technologies, Santa Clara, CA, USA) after the extraction of soil samples with dichloromethane in an ASE200 Accelerated Solvent Extractor (Dionex Co., Sunnyvale, CA, USA) [22,25]. Moreover, the bioavailable fraction of the 16 PAH compounds (16PAH bioav) was measured by GC-MS after extraction with Tenax-TA as a solid-phase sorbent according to Smreczak [26].

Soil physical-chemical properties (including quality assurance (QA)/quality control (QC) parameters) were described in detail by Klimkowicz-Pawlas et al. [22,24] and Ukalska-Jaruga et al. [25]. The study area was dominated by loamy sand soils (sand 72%, silt 26%, clay 2%, on average; CV = 10, 27 and 33%, respectively) with acidic reaction (mean pH_KCl value equal to 5.4). TOC and TN contents were in the ranges of 5.2–106.9 g kg⁻¹ (CV = 115%) and 0.2–4.1 g kg⁻¹ (CV = 88%), respectively. Considerable differences were observed for the total (CV = 447%) and bioavailable PAH concentration (CV = 247%). 16PAHs averaged 12,620 µg kg⁻¹ and reached levels as high as 316,085 µg kg⁻¹, but the content of bioavailable PAHs was significantly lower (Table 1).

The reference site was selected taking into account the following criteria: no contamination, pedological properties (organic carbon and nitrogen content, pH and texture) similar to those of the tested soils, the same land use and similar biological properties [15,27]. The reference site had the following characteristics: loamy sand texture (70% sand, 28% silt, 2% clay), TOC content of 12.9 g kg⁻¹, pH of 5.4 and TN contents of 0.6 g kg⁻¹ and 408 µg Σ16 PAHs kg⁻¹. A more detailed analysis of the difficulties involved in selecting a reference site is given elsewhere [22,28,29].

2.3. Soil Ecotoxicity Testing with Bioassays Battery (Ecotoxicological Line of Evidence)

Ecological risk was assessed towards organisms from different taxonomic groups: microorganisms, invertebrates and plants (Tables S1 and S3). Soil elutriates were prepared according to the ISO 21268-1 [30] leaching procedure: soils were mixed with 0.001 mol L⁻¹ CaCl₂ in a ratio of 1:2 (w/v), shacked through 24 h (20 ± 2 °C, 125 rpm) and filtered with the 0.45 µm syringe filters. We applied two kinds of tests: a water-phase to assess the soil retention function and a solid-phase to assess the soil habitat function.

2.3.1. Water-Phase Ecotoxicity Testing

An acute toxicity assay with the luminescent bacteria Aliivibrio fischeri was performed according to the Microtox 81.9% Screening test procedure [31]. Duplicated samples of soil leachates were incubated with bacteria at 15 °C for 15 min, and then the light output of the samples was recorded with a Microtox Model 500 analyser and analysed using the MicrotoxOmni software. The results were expressed as a percentage of the toxic effect (PE).

Rapidtoxkit is an indicator test for detecting the presence of toxic substances using crustacean larvae Thamnocephalus platyurus [22,32,33], which take up their food by filtration from the surrounding solution. The rate and level of filtration decrease or, in extreme cases, stop completely after the exposure of T. platyurus to toxic substances. Prior to the test, the cysts were rehydrated, transferred to a hatching container and incubated for 34 h at 25 ± 2 °C under continuous illumination (3000–4000 lux). T. platyurus larvae were then exposed to soil elutriates (or a culture medium as a control) for 1 h (25 °C, in the dark), followed by feeding with artificial food (suspension of red latex microspheres) for 30 min. The presence or absence of coloured microspheres in the digestive tract of the crustacean larvae was observed under a stereomicroscope, and the results were expressed as the mean percentage inhibition of particle uptake.

The Phytotestkit was conducted in flat test plates using two plant species: white mustard (Sinapis alba L.) and garden cress (Lepidium sativum L.), for which germination inhibition and root length inhibition were determined. The test was performed according to the procedure recommended by the manufacturer [34]. The lower part of the plate (21 × 15.5 × 0.8 cm) was filled with a foam pad and covered with a white thick cellulose filter, on which 25 mL of soil elutriate (in 0.001 mol L⁻¹ CaCl₂) was slowly spread in order to hydrate the filter paper completely. Then, the lower part of the plate was covered with a black paper filter, and 10 seeds of the test plants were placed on it. The test plates were incubated vertically at 25 ± 2°C in darkness for 3 days. A photo of the plate was taken using a digital camera, and the root length was measured in the ‘Image Tools’ image analysis program ver. 3.0. The test control was a 0.001 mol L⁻¹ CaCl₂ solution used to prepare the soil extracts. The test was carried out in three replicates for each of the soil/plant combinations tested.

2.3.2. Solid-Phase Ecotoxicity Testing

A direct contact toxicity assay with the ostracods Heterocypris incongruens was performed according to ISO 14371 [35]. The assay was carried out in the six-well polystyrene multi-well plates. Ten neonates of the H. incongruens hatched from cysts were added to each well of the microplate containing the tested soil sample (2 g), standard nutrient medium (2 mL) and algal suspension (1 mL) as food. After 6 days of exposure to contaminated soil samples (25 ± 2 °C, in the darkness), the percentage mortality and the growth of the H. incongruens were determined and compared to those of the control (non-toxic sediment).

Two plant growth tests were carried out: seed germination and an early growth assay with white mustard [36] and a screening test for the emergence of lettuce seedlings [37]. Approximately 90 mL of the soil sample was placed on the bottom part of the test plate, hydrated with distilled water to 100% of the water holding capacity and covered with a paper filter. Ten seeds of the same plant species (Sinapis alba L. or Lactuca sativa L.) were sown in one line at the same distance from one another; the test was carried out in triplicate. After closing with a transparent cover, the test plates were placed vertically in a holder and incubated at 25 °C for 3 days (S. alba) or 5 days (L. sativa). After incubation, a photo was taken with a digital camera, and the length of the roots was measured with the use of the image analysis program (ImageTool, ver. 3.0). The test control was OECD soil, the standard recommended for use in the plant growth test [36].

2.4. Soil Microbial and Biochemical Parameters (Ecological Line of Evidence)

Soil enzymatic activity and respiration, microbial biomass, nitrification and carbon mineralisation were considered as indicators within the ecological line of evidence (Tables S1 and S3).

Dehydrogenases activity (DH) was assessed after 24 h of incubation of the soil samples (37 ± 2 °C) with 3% triphenyltetrazolium chloride (TTC) as an electron acceptor [38]. The resulting triphenylformazan (TPF) was extracted from soils with ethanol and determined colorimetrically (wavelength of 485 nm) using a Lambda 45 UV-VIS spectrophotometer (Perkin Elmer, Waltham, MA, USA).

The basal (BR) and substrate-induced (SIR) soil respiration were analysed according to the ISO 16072 and ISO 14240-1 methods, respectively. Field-moist soil subsamples were incubated for 24 h (BR) or 6 hours (SIR) at 20 ± 2 °C in darkness. Glucose at 10 g kg⁻¹ was applied as a substrate in SIR measurements. CO₂ released during incubation (in both BR and SIR analyses) was entrapped in 0.05 mol L⁻¹ NaOH and precipitated as barium carbonate by adding 0.5 mol L⁻¹ BaCl₂ solution. The non-consumed NaOH was titrated with HCl in the presence of phenolphthalein as an indicator. BR and SIR were expressed as a µg of released CO₂ per gram of dry weight of soil per hour. The microbial biomass (MB) was based on the substrate-induced respiration ISO 14240-1 method and was expressed in µg of microbial C per g of dry weight of soil.

The carbon mineralisation test (CMIN) was carried out according to the OECD 217 [39] standard. Soil samples (100 g) of 60% water holding capacity were incubated for 28 days (20 ± 2 °C, darkness). The soil moisture was checked daily by the weighting of each experimental pot, and the water content was replenished as needed. C mineralisation was assessed at the beginning of the experiment (time 0 days) and after 28 days of soil incubation. At each time interval, 20 g of the soil subsamples was amended with 200 mg of glucose powder, and CMIN was calculated by CO₂ evolution after the next 24 h, as described in the ISO 16072 method.

The nitrification potential (NIT) was determined according to the method described by Maliszewska-Kordybach et al. [40]. Briefly, soil samples were incubated (24 h, 125 rpm, 20 ± 2 °C) with a mineral medium containing ammonium sulphate as the substrate, potassium phosphate buffer and sodium chlorate (to prevent the further oxidation of NO₂⁻ to NO₃⁻). The sample slurry was then mixed with 4 mol L⁻¹ KCl and filtered, and the colour reagent (containing sulphanilamide and N-(1-naphthyl)ethylenediamine dihydrochloride) was added. The amount of NO₂⁻ formed by the nitrification process was determined after 1 h using a Lambda 45 UV-VIS spectrophotometer (Perkin Elmer, Waltham, MA, USA) at 543 nm and expressed as μg NO₂⁻ per gram of soil dry weight per 24 h. All microbiological and biochemical measurements were made in triplicate, and each series of analyses included control samples without soil material. The accuracy of the biological methods, expressed in terms of the relative standard deviation, ranged from 5 to 25% depending on the parameter analysed.

2.5. Statistical Analysis and Integrated Risk Calculation

Statistical parameters—mean, median, standard deviation (SD), range (Min and Max) and coefficient of variation (CV)—were analysed using the Statgraphics Centurion software (ver. XVIII, Statpoint Technologies). Differences between the response of test organisms and ecological parameters in contaminated and reference soil were assessed via a one-way ANOVA variance analysis, followed by the Duncan’s test (p ≤ 0.05).

The method described by Jensen et al. [23] and Niemeyer et al. [14] was used to calculate risk, including scaling the results for each singular parameter, integrating the scaled results for each assessment line (chemical—IR_Chem, ecological—IR_Ecol and ecotoxicological—IR_Ecotox) and calculating an overall integrated risk index (IntRisk) (Table S2). All parameters were normalised with respect to the reference soil and expressed on a scale from zero (no effect) to one (maximal ecosystem effect). In the chemical line of evidence, the toxic pressure (TP′) for each PAH (expressed as a potentially affected fraction of species) was calculated according to Equation (1) [23].

$$T{P}^{\prime }=\frac{1}{1+exp\left\{\frac{logMPC-logPAH}{0.4}\right\}}$$

(1)

The PAH content was related to the maximum permissible concentration (MPC) given by Verbruggen [41].

Next, the scaled TP′ were integrated into one TP index according to the formula [6,23,28]:

TP_16PAH = 1 − (1 − TP′₁) (1 − TP′₂) … (1 − TP′_n)

(2)

The chemical, ecotoxicological and ecological lines of evidence included two, seven and six indicators (mentioned in Section 2.2, Section 2.3, Section 2.4 and Tables S1 and S3), respectively, and IR_Chem, IR_Ecotox and IR_Ecol were calculated as follows [14,19,23,29]:

$${\mathrm{I}\mathrm{R}}_{\mathrm{C}\mathrm{h}\mathrm{e}\mathrm{m}} \mathrm{o}\mathrm{r} {\mathrm{I}\mathrm{R}}_{\mathrm{E}\mathrm{c}\mathrm{o}\mathrm{t}\mathrm{o}\mathrm{x}} \mathrm{o}\mathrm{r} {\mathrm{I}\mathrm{R}}_{\mathrm{E}\mathrm{c}\mathrm{o}\mathrm{l}}=1-\left({10}^{averageof log(1-A)}\right)$$

(3)

where A is the scaled value of each chemical, ecotoxicological or ecological parameter.

In the final step, the integrated risk index was determined according to Equation (4) [23,42], and uncertainty was assessed based on SD calculation for three separate IR indexes.

$$IntRisk=1-\left({10}^{\left[\log \left(1-I{R\_}_{Chem}\right)+\log \left(1-I{R\_}_{Ecotox}\right)+\log \left(1-I{R}_{\_Ecol}\right)\right]/n}\right)$$

(4)

3. Results and Discussion

3.1. Chemical Line of Evidence

The site-specific toxic pressure of the 16PAH mixture, reflecting the potentially affected fraction of species, was calculated based on a response addition model [3,23] using the total measured concentration of PAHs (Table 1) and MPC values reported by Verbruggen [41]. Maximum permissible concentrations were derived from the no-observed-effect concentration (NOEC) or 10% effect concentration (EC₁₀) value of the most sensitive terrestrial species and should be protective for all species [3,41]. For 30% of the sampling sites [6], TP_16PAH values ranged from 0.01 to 0.20 and were below the lowest threshold (0.25) according to the classification given by Jensen et al. [23] and Dagnino et al. [28], indicating no risk. Furthermore, the Σ16PAH content in these soils ranged from 311 to 591 µg kg⁻¹ and was at the level recorded in the reference soil; such contaminant content <600 µg kg⁻¹ corresponds to the class of uncontaminated soils according to the classification of Maliszewska-Kordybach [43]. The permissible limits (100 or 200 µg kg⁻¹) for 10 compounds from Polish soil guidelines [44] were also not exceeded, and these nine sampling points were therefore excluded from further stages of the assessment.

The TP_16PAH values in the remaining area were rather high, ranging from 0.36 to 1.00 (Table 6); 63% of the soils can be classified as medium and high risk, with TP values of 0.56–0.74 and 0.86–1.00, respectively. The high toxic pressure is due to the high total concentration of PAHs (956–316,085 µg kg⁻¹ (Table S3). As highlighted in previous studies by Semenzin et al. [45], Hong et al. [18] and Grassi et al. [46], such a conservative approach considering only the total content of contaminants may lead to a biased risk assessment and should only be used at screening stages. According to Ortega-Calvo et al. [8] and Peijnenburg [9], determining the bioavailability of contaminants is crucial for a proper and realistic assessment of contaminated sites. In our study, the bioavailable PAHs ranged from 22.6 to 1184 µg kg⁻¹ (Table 1), representing only 0.4 to 7.5% of the total 16PAHs. It was in line with the previous data of Smreczak [26] and Ukalska-Jaruga and Smreczak [47], who reported that the potentially bioavailable PAHs in historically contaminated soils did not exceed 8%. Therefore, both total and bioavailable 16PAH concentrations were used in our study for risk calculation.

3.2. Ecotoxicity Analysis

The battery of bioassays with various organisms was utilised to assess the toxicity of soils and soil elutriates (

Table 2. Response of organisms (PE%) for water-phase and solid soil samples—ecotoxicological line of evidence (n = 31).

Organism	Endpoint	Mean	Median	SD	Min	Max	CV%
Water-phase (retention function)
A.fischeri	LUI_Scr	13.9	13.7	15.5	–10.9	45.7	111
T. platyurus	IF	25.2	33.6	21.9	–10.0	56.7	87
L. sativum	IRG_LEP	7.9	6.8	12.7	–9.2	39.8	161
S. alba	IRG_SIN	–4.6	–8.6	22.1	–48.1	37.8	482
Solid-phase (habitat function)
H. incongruens	GI	30.0	30.2	17.0	2.5	63.8	57
S. alba	IRG_SINS	15.6	11.1	18.5	–12.6	60.6	118
L. sativa	IRG_LAC	21.5	19.1	30.9	–24.9	72.9	144

CV%—variation coefficient; LUI_Scr—luminescence inhibition of A. fischeri; IF—inhibition of food ingestion by T. platyurus; IRG_SIN and IRG_LEP—inhibition of root growth of S. alba and L. sativum, respectively (in soil leachates); GI—growth inhibition of H. incongruens; IRG_SINS and IRG_LAC—inhibition of root growth of S. alba and L. sativa, respectively (soil samples).

and Table S1). Liquid-phase tests with soil elutriates were applied to obtain information about the fraction of contaminants potentially reaching the groundwater by the water path (retention function of soils), whereas solid-phase tests were used to assess the soil habitat function [4,14]. The response of the test organisms varied significantly according to the type of samples tested and the type of test species (Table 2, Figure 2).

The soil and soil leachates toxicity was evaluated using the toxicity assessment method developed by Persoone et al. [48], where PE < 20% means no toxic effect, PE from 20 to 50% means a low toxic sample, PE in the range of 50–100% means a toxic sample and PE = 100% means a very toxic sample. In general, a weaker response of the organisms was observed for the water-phase tests. The sensitivity of the organisms can be ordered as follows: crustacean T. platyurus > bacteria A. fischeri > plants L. sativum = S. alba (Table 2, Figure 2).

Food uptake by T. platyurus was inhibited by an average of 25.2%; for most samples, it was within the range of up to 50%, indicating slight toxicity [48]. The A. fischeri luminescence inhibition was at a level of 14%; for a majority of samples, no toxic effect was observed, and only 24% of samples exhibited low toxicity (Figure 1), which was correlated with soil pH_KCl (r = −0.40, p < 0.05). Both plants from the Phytotestkit (L. sativum and S. alba) test showed a similar response, averaging 7.9% (IRG_LEP) and –4.6% (IRG_SIN), respectively, with the stimulation of root growth more frequently observed in mustard. For extracts from several samples (17, 18 and 35), the inhibition of root growth was approximately 30–40%.

Significantly higher effects were observed in the solid-phase tests (Table 2, Figure 2). The most sensitive was H. incongruens; no crustacean mortality was observed, while the growth inhibition ranged from 2.5 to 64% and averaged 30%. Only seven soil samples were non-toxic; the remaining samples were low-toxic (61% of sampling points) or toxic (9% of soils), corresponding to toxicity classes II and III, according to Persoone et al. [48]. The inhibition of plant root elongation (IRG_SINS and IRG_LAC) was 15.6% and 21.5% for white mustard and lettuce, respectively (Table 2). For both plants, 50–60% of the samples were non-toxic (Figure 2), with root growth stimulation effects often observed. Lettuce was the more sensitive plant; based on the IRG_LAC endpoint, 23% of the samples can be classified as slightly toxic, and 26% of the samples can be classified as toxic.

The no toxicity or low toxicity for test organisms observed in our study (Table 2, Figure 2), particularly in the aqueous phase tests, indicates the limited availability and leaching of PAHs from the soil system and thus a low risk of their transfer to the water environment. Soil contamination in the Czerwionka region is mainly the result of long-term deposition from several sources, i.e., coking plant, hard coal mining activity, storage of mine wastes and bituminous substances production [22,24,25]. Several authors [25,26,47,49,50,51,52] indicated that, in historically contaminated soils, PAHs compounds are subjected to ageing processes leading to significant PAH sequestration or binding, resulting in the reduction in their mobility, extractability and availability. Soil organic matter plays a major role in these processes [25,47,50,51]. Although most of the soils from the region were characterised by a low TOC content (Table 1), according to a previous study by Ukalska-Jaruga et al. [25], these soils showed a significant contribution of humin (up to 123 g kg⁻¹) and black carbon (up to 45 g kg⁻¹) fractions, further limiting the desorption and bioavailability/bioaccessibility of PAHs.

3.3. Ecological Indicators

Soil contamination can affect both the abundance and activity of soil microorganisms, which play a key role in many soil biochemical processes (e.g., decomposition of organic compounds, organic matter synthesis, nutrient cycling, soil detoxification, etc.) essential for the soil ecological functions [40,53,54]. Furthermore, due to their direct contact with the soil, soil microorganisms can respond quickly to changes caused by natural and anthropogenic factors and can therefore be a direct measure of soil function [55,56]. Our study included parameters reflecting the soil’s functional diversity: enzymatic activity (DH), soil respiration (BR and SIR), microbial biomass (MB), nitrification (NIT) and carbon mineralisation (CMIN) (

Table 3. Statistical evaluation of biological indicators used to calculate risk within the ecological line of evidence (n = 31).

Indicator	Unit	Mean	Median	SD	Min	Max	CV%
DH	µg TPF g−1 d.w. 24 h−1	33.0	32.6	13.7	13.1	86.9	41
NIT	µg NO2− g−1 d.w. 24 h−1	3.0	2.4	2.4	0.8	13.7	80
BR	µg CO2 g−1 d.w. h−1	2.9	2.5	1.7	0.3	7.8	58
SIR	µg CO2 g−1 d.w. h−1	10.6	8.4	6.1	2.6	26.3	57
MB	µg Cmic g−1 d.w.	21.9	17.5	12.4	5.7	54.0	56
CMIN	µg CO2 g−1 d.w. h−1	15.9	14.0	9.3	1.9	38.2	58

and Table S3), indicators commonly used in both soil quality and contaminant exposure assessments [14,22,53,57,58].

Overall, it can be concluded that all soils were characterised by relatively low biological activity, which, for most indicators, did not vary significantly between sampling sites (CV 41–58%). Slightly greater variability (CV = 80%) was found for the nitrification potential, which reflects the activity of the ammonia-oxidising bacteria responsible for the first stage of the nitrification process and is considered one of the most sensitive microbiological parameters responding strongly to the natural factors [59] and the presence of pollutants such as PAHs [40] or metals [14]. NIT values ranged from 0.8 to 13.7 µg NO₂⁻ g⁻¹ d.w. 24h⁻¹, with a median value of 2.4 µg NO₂⁻ g⁻¹ d.w. 24 h⁻¹, and were strongly related to both soil properties (r for TOC = 0.72, p < 0.000; r for pH_KCl = 0.67, p < 0.000) and the pollutants content (r for 16PAHs = 0.82, p < 0.000). For the other parameters (BR, SIR, MB), these relationships were much weaker, and the correlation coefficients with TOC, soil pH and PAH content were 0.36–0.54, 0.39–0.61 and 0.35–0.45 (p < 0.05), respectively. The dehydrogenase activity, intracellular enzymes related to the oxidative phosphorylation process, was 33 µg TPF g⁻¹ d.w. 24 h⁻¹, on average (Table 3). According to Chakravarty et al. [7], low enzymatic activity (e.g., dehydrogenases, urease, phosphatase) in contaminated soils could be attributed to the poor microbial growth and deficiency of available substrates. Soil microbial respiration and soil microbial biomass are indicators closely related to organic matter degradation and synthesis; soil respiration reflected soil carbon availability to microorganisms, and MB is an important reservoir of nutrients in ecosystems [53,60]. Sites with a high microbial biomass can store and recycle more nutrients [53]; both soil respiration and biomass, in our study, were rather low at the levels of 2.5 µg CO₂ g⁻¹ d.w. h⁻¹ (for BR) and 17.5 µg Cmic g⁻¹ d.w. (MB), indicating that these functions are impaired. It was also confirmed by the low carbon mineralisation (CMIN) rate (Table 3). Under stress conditions (e.g., presence of pollutants), microorganisms usually exhibit a lower metabolic efficiency and need more energy for survival, resulting in less TOC incorporation into the microbial biomass [53].

3.4. Integrated Environmental Risk Assessment

The results from all analyses (chemical, ecotoxicological and ecological) were used to determine the environmental risks associated with PAHs in soils. Initially, the results were normalised to the reference soil, and risk indices were calculated for the individual indicators measured (

Table 4. Scaled and corrected risk values (RI) calculated for single ecotoxicological test endpoints and combined risk for retention (IR_reten) and habitat function (IR_habit), respectively.

	Retention Function					Habitat Function
Sample	LUI_Scr	IF	IRG_SIN	IRG_LEP	IR_Reten	GI	IRG_SINS	IRG_LAC	IR_Habit
3	0.00	0.00	0.11	0.09	0.05	0.00	0.00	0.48	0.20
4	0.00	0.00	0.06	0.16	0.06	0.26	0.00	0.00	0.09
5	0.00	0.00	0.00	0.07	0.02	0.34	0.21	0.60	0.32
6	0.20	0.00	0.00	0.01	0.06	0.59	0.00	0.30	0.67
7	0.10	0.00	0.00	0.00	0.03	0.03	0.04	0.59	0.28
8	0.24	0.00	0.00	0.07	0.08	0.18	0.22	0.50	0.25
16	0.07	0.00	0.21	0.16	0.11	0.26	0.02	0.73	0.34
17	0.01	0.00	0.29	0.20	0.14	0.21	0.11	0.57	0.26
18	0.07	0.00	0.16	0.29	0.14	0.28	0.04	0.67	0.55
19	0.27	0.00	0.00	0.00	0.08	0.17	0.00	0.00	0.05
20	0.07	0.00	0.00	0.00	0.02	0.37	0.05	0.35	0.34
21	0.18	0.00	0.00	0.00	0.05	0.36	0.18	0.00	0.15
21a	0.04	0.00	0.00	0.07	0.03	0.17	0.17	0.14	0.38
22a	0.00	0.00	0.00	0.02	0.01	0.00	0.00	0.02	0.01
26	0.39	0.00	0.00	0.00	0.12	0.37	0.14	0.00	0.38
27	0.00	0.00	0.15	0.26	0.11	0.00	0.34	0.00	0.10
29	0.00	0.00	0.00	0.00	0.00	0.58	0.38	0.39	0.37
30	0.00	0.00	0.00	0.04	0.01	0.24	0.28	0.14	0.70
32	0.00	0.00	0.06	0.19	0.07	0.35	0.43	0.00	0.22
33	0.00	0.00	0.00	0.02	0.00	0.00	0.37	0.57	0.28
35	0.00	0.00	0.38	0.40	0.22	0.06	0.23	0.00	0.08

Explanation of abbreviations under Table 2.

,

Table 5. Scaled and corrected risk values (RI) calculated for single biological activity indicators according to the ecological line of evidence.

Sample	DH	NIT	BR	SIR	MB	CMIN
3	0.20	0.62	0.13	0.50	0.49	0.01
4	0.38	0.45	0.06	0.10	0.10	0.61
5	0.23	0.58	0.00	0.54	0.53	0.00
6	0.10	0.78	0.00	0.31	0.30	0.67
7	0.46	0.68	0.41	0.26	0.25	0.44
8	0.12	0.44	0.09	0.33	0.32	0.00
16	0.00	0.00	0.00	0.00	0.00	0.00
17	0.01	0.02	0.36	0.00	0.00	0.04
18	0.10	0.52	0.84	0.71	0.70	0.59
19	0.00	0.44	0.00	0.00	0.00	0.00
20	0.09	0.52	0.00	0.07	0.07	0.13
21	0.00	0.23	0.42	0.00	0.00	0.00
21a	0.00	0.55	0.40	0.36	0.36	0.18
22a	0.44	0.00	0.00	0.00	0.00	0.00
26	0.00	0.61	0.00	0.00	0.00	0.47
27	0.06	0.00	0.00	0.00	0.00	0.43
29	0.00	0.00	0.00	0.00	0.00	0.05
30	0.21	0.72	0.00	0.13	0.13	0.74
32	0.33	0.00	0.00	0.19	0.19	0.28
33	0.37	0.00	0.00	0.00	0.00	0.30
35	0.26	0.15	0.00	0.00	0.00	0.00

and Table S2); the results were then integrated across the different lines of evidence (IR_Chem, IR_Ecotox, IR_Ecol in

Table 6. Toxic pressure for total (TP_16PAH) and bioavailable PAHs (TP_{16PAH bioav.}), risk indexes (IR) for individual LoEs and final integrated risk index (IntRisk).

Sample	TP16PAH	TP16PAH bioav.	IR_Chem	IR_Ecotox	IR_Ecol	IntRisk	SD	Evaluation *
3	0.39	0.01	0.23	0.12	0.36	0.24	0.12	no risk
4	0.59	0.01	0.36	0.07	0.32	0.26	0.16	low risk
5	0.67	0.00	0.43	0.18	0.36	0.33	0.13	low risk
6	0.63	0.01	0.40	0.44	0.44	0.43	0.03	low risk
7	0.36	0.00	0.20	0.14	0.44	0.27	0.16	low risk
8	0.98	0.00	0.85	0.17	0.23	0.55	0.38	moderate risk
16	0.99	0.00	0.92	0.23	0.00	0.61	0.48	moderate risk **
17	0.71	0.00	0.46	0.20	0.08	0.27	0.19	low risk
18	0.89	0.01	0.67	0.38	0.63	0.58	0.16	moderate risk
19	0.74	0.01	0.49	0.06	0.09	0.24	0.24	no risk
20	0.97	0.00	0.83	0.19	0.17	0.52	0.38	moderate risk
21	0.92	0.00	0.72	0.10	0.12	0.40	0.35	low risk
21a	0.99	0.00	0.92	0.22	0.33	0.65	0.37	moderate risk
22a	0.65	0.00	0.41	0.01	0.09	0.19	0.21	no risk
26	1.00	0.08	0.98	0.26	0.23	0.77	0.42	high risk **
27	0.56	n.d.	0.56	0.10	0.10	0.29	0.27	low risk
29	0.69	0.07	0.46	0.21	0.01	0.25	0.23	no risk
30	0.86	0.10	0.64	0.45	0.40	0.51	0.13	moderate risk
32	0.59	n.d.	0.59	0.14	0.18	0.34	0.25	low risk
33	1.00	0.53	1.00	0.15	0.13	0.94	0.50	high risk **
35	1.00	0.08	1.00	0.15	0.07	0.94	0.51	high risk **

n.d.—not determined; SD—standard deviation; * evaluation of risk level based on Jensen et al. [23] and Dagnino et al. [28]; ** SD value > 0.40 indicates high uncertainty and no acceptable risk.

and Table S2), with the combined risk for retention (IR_reten) and habitat function (IR_habit) additionally calculated for the ecotoxicological assessment line (Table 4). Finally, after considering all available data, an integrated risk index (IntRisk) was calculated for each of the study sites (Equation (4), Table 6 and Table S2).

The combined IR_Chem index included the total and bioavailable fraction of the 16PAH content and was slightly lower compared to TP based on total PAHs only (Table 6). High chemical risk values (0.67–1.00) were found in nine sampling sites (8, 16, 18, 20, 21, 21a, 26, 33 and 35) (Table 6). On the contrary, the risk index (TP_{16PAH bioav.}) calculated for the bioavailable fraction was at a level below 0.1 (no-risk), and for one sampling site (sample 33), it was equal to 0.53. This sample had an extremely high 16PAH content (316,085 µg kg⁻¹, Table S3) and, thus, high 16PAH bioav. (1184 µg kg⁻¹); however, 16PAH bioav. accounted for only 0.4% of the total PAH content. Although the consideration of bioavailability in environmental risk analyses is recommended by other authors [8,9,10,46,50], so far, only a few papers describing the Triad method have considered availability analysis, and this has been for metals only [12,14,20,61].

As described in Section 3.2, organisms reacted differently to the presence of contaminants, resulting in different risk index values in relation to ecotoxicological parameters (Table 4). The risk values for the soil retention function (IR_reten) did not exceed 0.22, thus indicating no risk. When considering single parameters, low risk was found in only two (RI for LUI_Scr and IRG_SIN) or three cases (RI for IRG_LEP). Interestingly, for T. palyurus, characterised by the highest sensitivity (Figure 2), the risk index was 0.00, undoubtedly related to the higher toxicity found in the reference soil; similar observations were also reported by Gutierrez et al. [15] and Niemeyer et al. [14], indicating that the selection of a reference site is a critical point in the risk assessment process.

The RI values for tests describing the effect of pollutants on the soil’s habitat function were significantly different; IR_habit was three times higher and resulted mainly from high indices for H. incongruens and L. sativa (Table 4), for which the highest toxicity was observed. The risk index value finally calculated for all seven ecotoxicological endpoints (IR_Ecotox) was in the range of 0.01–0.45 (Table 6) and indicated no risk in 80% of the sampling sites; only for four sampling points (6, 18, 26 and 30) could the risk be described as low.

Upon analysing the risk indices for the ecological line of evidence (Table 6), no risk was found for 50% of the samples (IR_Ecol from 0.00 to 0.23). The remaining area showed low risk (0.32–0.44), and only soil 18 showed medium risk (IR_Ecol = 0.63), closely linked to the lower activity of NIT, BR, SIR, MB and CMIN in this case (Table 5), due to the acidic soil reaction (pH_KCl = 4.1).

The integrated IntRisk (calculated from the three independent LoEs using Equation (4)) ranged from 0.19 to 0.94 (Table 6). The risk level was determined according to the classification given by Jensen et al. [23] and Dagnino et al. [28], according to the following: IntRisk ≤ 0.25 no risk; 0.25 < IntRisk ≤ 0.50 low risk; 0.50 < IntRisk ≤ 0.75 moderate risk and 0.75 < IntRisk ≤ 1.00 high risk. Furthermore, in order to determine the uncertainty of the results obtained, the SD between LoEs was calculated, assuming that IntRisk values with a standard deviation greater than 0.4 [23] indicate high uncertainty in the risk assessment and unacceptable risk.

The lowest IntRisk below 0.25 (no risk) was found in sampling points 3, 19, 22a and 29, located at a distance of about 5–8 km from the coking plant (the main pollution source in the area), with 16PAHs from 989 to 1620 µg kg⁻¹. A low risk level (IntRisk 0.26–0.43) was associated with eight locations (4, 5, 6, 7, 17, 21, 27 and 32) unevenly distributed over the study area with a similar 16PAH concentration (not exceeding 1500 µg kg⁻¹), and due to the low discrepancy in the results between the three lines (SD values < 0.4), the risk should be considered acceptable. For the other six soils, the IntRisk ranged from 0.51 to 0.65 (moderate risk) and was mainly linked with high IR_Chem values; based on adopted criteria [23,28], risk in five sampling sites (8, 18, 20, 21a and 30) can be assessed as acceptable. Only for one site (sample 16 located near the coal recovery plant) were the differences between three LoEs high, indicating that additional evaluation in this area is needed. The IntRisk had very high and unacceptable values (0.77 and 0.94, respectively) in only three sampling points (26, 33 and 35), representing less than 10% of the study area. These sites were in close vicinity to the coke plant (site 33 approx. 400 m and site 35 about 100 m from the source) and to the asphalt processing plant (site 26). In addition, extremely high concentrations of 16 PAHs were recorded at these points (Table S3): 7604 µg kg⁻¹ (site 26), 316,085 µg kg⁻¹ (site 33) and 21,768 µg kg⁻¹ (site 35). However, it should be emphasised that the ecological and ecotoxicological risk indices in this sampling point were low and did not exceed 0.25. In the previous study of Klimkowicz-Pawlas et al. [22], describing the Tier 1 assessment, almost 50% of the area was classified under medium and high risk. Moreover, high differences between chemical, ecological and ecotoxicological RI indexes were observed.

In our detailed risk assessment, a total of 15 different indicators giving information on the PAH content, bioavailability, toxicity and soil biota activity were included. This comprehensive approach is very valuable, as previous studies mostly reported a screening risk assessment [12,17,20,22] and analysed only chemical and ecotoxicological data [4,17,20,61], and only a few papers [14,19] described higher tiers of the Triad procedure. Most often, these studies focused on post-industrial sites contaminated by metals [11,12,13,17,18,19], and only a few papers focused on sites contaminated by PAHs [21,22]. A very important issue in the Triad methodology is the weighting of results from different disciplines. Equal weighting is recommended by default [10,46], and this approach was used in our study. However, when we assume that chemical analyses (despite taking into account the bioavailability of PAHs) are less important than the results of the ecotoxicological and ecological tests and we adopt the different weights respectively, e.g., 1, 1.5 and 2 suggested by Dagnino et al. [28] or Terekhova et al. [13], the final result of the risk assessment may be quite different, especially in cases where the contribution of IR_Chem to the total IntRisk is very high. This also applies to our study, where the unacceptable risks found at sampling points 26, 33 and 35 are mainly due to extremely high chemical indices (Table 6). Recalculating the data with different weights gave a lower IntRisk index (0.41, 0.32 and 0.31 for sites 26, 33 and 35 respectively), within the range of low risk values. Nevertheless, due to the high uncertainty of the results (SD > 0.4), the risk is still unacceptable.

4. Conclusions

The detailed assessment using the Triad approach was conducted for the agricultural area contaminated with organic pollutants, with PAHs as an example. A comprehensive set of 15 indicators was analysed and subsequently scaled and integrated to obtain risk indicators for the individual assessment routes and the final integrated IntRisk. The values for IR_Chem, IR_Ecotox and IR_Ecol were 0.20–1.00, 0.00–0.45 and 0.00–0.63, respectively, while IntRisk ranged from 0.19 to 0.94 and was mainly dependent on high chemical risk index values. Almost 90% of the study area had no or low risk (IntRisk was < 0.5). Very high unacceptable risk (IntRisk 0.77–0.94) was found at only three sampling sites located in close proximity (100–400 m) to a coke plant or asphalt processing plant. Our study revealed that the integration of data on the organism’s response and ecological indicators into the risk assessment process provides a realistic assessment and allows for deriving quantitative risk indexes and for delineating the area of the high risk which needs additional action, e.g., further monitoring, changes in land use or remediation.

Several important points should be highlighted, which are crucial for the final outcome of an environmental risk assessment. Our research confirms that the risk assessment cannot be based solely on the total content of contaminants; the determination of the available fraction of compounds should also be included in the chemical line of evidence. In addition, a more realistic assessment is provided by the consideration of toxicity tests and ecological indicators. However, it has to be kept in mind that no species is sensitive to all contaminants, so bioassays representing different taxonomic groups of organisms should be included in the studies. A certain challenge is the selection of a suitable reference site, which should be uncontaminated and have pedological characteristics similar to those of contaminated sites. Sometimes, however, despite meeting these criteria, higher toxicity or lower biological activity may be observed for the reference site, which may be related to the presence of other contaminants not considered in the study. Another issue is the use of the appropriate weighting of results obtained using equal weights (assuming an equal influence of different LoEs on the final result) or different weights (assigning higher weights to indicators that better reflect the state of the soil ecosystem). In each case, the final value of the risk index may be different.

However, it is important to emphasise that the calculation of integrated risk indexes provides a more comprehensive and realistic assessment, enables the prioritisation of contaminated areas and allows for a decision on the actions needed at a site. Such a methodology can be used in the future to assess other pollutants that are a threat to agricultural lands, including emerging contaminants such as pharmaceuticals or microplastics.