Identification of Environmental Damage Process of a Chromium-Contaminated Site in China

Abstract

Identifying the source and impact pathways of soil heavy-metal pollution is critical for its assessment and remediation. Numerical simulation has been widely used to simulate soil heavy-metal pollution processes and predict risks. However, traditional numerical simulation software requires a large number of parameters, which are difficult to obtain in site-scale studies. This study proposes a rapid method for identifying soil heavy-metal pollution processes using the TOUGH2/EOS7 software. It has automatic calibration and uncertainty analysis capabilities, which can effectively reduce the demand for parameters. This study established a method, including model selection, simulation, validation, and error analysis, to verify the effectiveness of the proposed method. This study identified the most realistic scenario for chromium pollution and simulated its release over 20 years, and the results met accuracy requirements with a best-case fit of 0.9998. The results showed that the method can quickly identify the source and impact pathways of soil heavy-metal pollution, providing strong evidence for environmental damage assessment.

1. Introduction

In the field of environmental damage assessment, the results are often used as judicial expertise in court, implying the necessity of a complete chain of evidence for the damage event [1]. Therefore, in the case of environmental damage, it is crucial to determine the cause–effect relationship by identifying the source and process of the damage [2,3,4,5]. The identification of the damage process requires a clear understanding of the complete pathway between the source and the receptor, encompassing both spatial continuity and temporal accuracy [6,7,8]. Historically, the lack of rigor in identifying the damage process has frequently been the weakest and most susceptible aspect of environmental assessments. This deficiency leads to an incomplete chain of evidence for cause–effect relationships and, thus, significantly reduces the credibility and acceptance rate of the overall environmental damage assessment [9,10].

To reconstruct the process of soil heavy-metal damage, a comprehensive understanding of the sources, behavior, bioavailability, migration, and transformation mechanisms of heavy metals in the soil is essential along with an assessment of their potential risks to the ecosystem. Research methods can generally be divided into laboratory simulation studies, field investigations and monitoring, and model simulation analysis [11].

Although laboratory simulation and field investigation methods can help to understand the migration and transformation process of pollutants and the identification and analysis of pollution sources [12,13,14,15,16], they cannot reproduce the pollution process of heavy metals in soil. At present, numerical simulation is mainly used to predict and evaluate the migration and transformation process of pollutants in groundwater [17,18,19,20]. However, there are relatively few studies on the migration of pollutants in soil porous media. There is currently a lack of rapid methods to reproduce the migration process of pollutants in soil porous media, which is an important evaluation basis for environmental damage assessment.

Due to the spatial heterogeneity of soil and the medium’s non-conformity with fundamental hydrodynamic assumptions, various tracing techniques and migration simulation methods often fail to be effectively utilized, adding to the difficulty of damage process identification. The HYDRUS V5.04, developed by the U.S. Department of Agriculture’s Soil Laboratory, is widely used for vadose zone migration simulation. Its derivatives, HYDRUS-1D, HYDRUS-2D, and HYDRUS-3D, can simulate the movement of water, heat, and solutes in the unsaturated zone under different dimensions and scenarios [21,22]. Recently, some scholars have also conducted research in the analysis of pollutant migration pathways in the vadose zone using HYDRUS, such as Anlauf et al. using HYDRUS-1D to simulate the migration process of pesticides in soil [23]; Shi and Ren using HYDRUS-1D to simulate the migration risk of manganese [21]; and Zeng and Guo developing a new model to simulate the transport and risk of perfluoro and polyfluoroalkyl substances in the vadose zone [24]. A few studies have involved using migration simulation methods for source tracing, such as combining vadose zone numerical simulation to identify potential source locations based on the spatial variability of VOCs [25] and Dafny using migration simulation of trichloroethylene in the vadose zone to compare incidental and one-time pollution sources [26]. However, these software require a large number of measured parameters [27], which are often difficult to obtain in site-scale studies. Especially in the practice of environmental damage assessment in China, due to legal constraints on detention periods and investigation custody after arrest, it is often necessary to complete the entire environmental damage assessment in a short time. Therefore, it is not feasible to conduct extensive experiments to obtain the parameters required for HYDRUS simulation. To address this issue, we need to adopt more suitable simulation software to depict the damage process and provide some data support for the damage assessment.

TOUGH2 V3.0 is a general-purpose numerical simulation software developed by Lawrence Berkeley National Laboratory in the United States. It is used to simulate the transport of multiphase, multicomponent fluids and heat in porous media [28]. TOUGH2 features a modular design and wide applications. It can be used to simulate various complex problems in the fields of geothermal energy development, groundwater remediation, CO₂ geological storage, natural gas hydrate exploitation, and nuclear waste disposal [29]. However, there are no cases of using TOUGH2 software for simulation and analysis in the identification of environmental damage processes. In this study, the TOUGH2 simulation software will be used to address challenges such as lack of basic data, scarcity of parameters, and the many uncertainties in the pollution process. Compared to the HYDRUS series of software [30], TOUGH2 allows for appropriate parametrization of underground structures and their properties. It can simulate multiple parameter sets, and through automatic model calibration, local and global sensitivity analysis, data analysis, and uncertainty analysis, it effectively resolves issues related to the lack of simulation parameters and the highly nonlinear problem-solving process [31]. Additionally, through model verification and the design of multiple scenarios, the most realistic damage processes are identified in this study. The research method provides scientific references and data support for a comprehensive assessment method in contaminated soil.

2. Methodology

2.1. Study Area

The study area, situated in a city in Southwest China, is primarily used for industrial purposes with a history of leather and phosphate fertilizer production, covering an area of approximately 24,000 m². Characterized by a generally flat terrain, the study area is adjacent to the River X, with the groundwater mainly consisting of loose rock pore water (Figure 1). The depth of the groundwater varies from 3 to 7 m, flowing perpendicular to the river with a hydraulic gradient of approximately 1.1 to 5‰. The aquifer is relatively abundant in water with a single well yield of approximately 300 to 365 m³/day, primarily recharged by atmospheric precipitation.

2.2. Sampling and Analysis

2.2.1. Sampling

Drawing upon the identification of production functional areas, facilities, and key zones within the study area and adhering to the relevant monitoring technical guidelines [32,33,34], this study established 44 sampling points in the production functional areas. Since one of the monitoring points was located in a gravel area and no soil was collected, there were a total of 43 sampling points. Taking into account the historical production and production facilities distribution in the study area, the soil investigation points were densely arranged in areas where pollution might occur. Therefore, the sampling points were placed at a density of one per 400 m² in key zones and one per 1600 m² in general zones. This was performed to ensure the representativeness of the points, reflect the overall soil quality in the area, and assist in determining the pollution boundaries of the key investigation zones.

According to the stratigraphic distribution and on-site rapid monitoring results in the study area, the soil sampling depths were set at 0.5 m, 2 m, 4 m, 6 m, and 8 m. Due to the proximity to the riverbed, some boreholes did not reach 8 m before encountering a cobble layer. Additionally, the soil layers below 8 m were mainly gravel-bearing aquifers, which were not the primary focus of this study. Therefore, the main research depth of the study area was determined to be 8 m.

2.2.2. Spatial Distribution Analysis

Employing classical statistical methods and IBM SPSS 22.0 software, the distribution of chromium content in the soil of the study area was meticulously analyzed across each layer (

Table 1. Description and statistics of Cr content in soil of the study area.

Pollutant	Depth (m)	Content Range (mg/kg)	Mean (mg/kg)	Median (mg/kg)	Standard Deviation (mg/kg)	Variance (mg/kg)	Coefficient of Variation (%)
Cr	0	43.00–2230.00	227.71	82.00	453.93	206,056.69	199.35
	2	42.00–206.00	89.74	78.00	39.58	1566.57	44.11
	4	31.00–190.00	75.95	72.00	29.07	844.82	38.27
	6	27.00–341.00	69.36	63.00	43.37	1881.00	62.53
	8	26.00–501.00	75.97	59.00	78.70	6193.36	103.59

). As shown in Table 1, the content of Cr in the soil of the study area tended to be higher in the surface layer (0~0.5 m) and gradually decreased with depth, with a slight increase at a depth of 6 m. This study chose the IDW (inverse distance weighting) method to predict the spatial distribution of pollutants in the soil of the study area. The spatial distribution of Cr across 120 soil samples was mapped using the spatial analysis module within Arc GIS 10.3 (Figure 2).

Figure 2 shows that, in the study area, the overall trend of chromium, a pollutant, exhibits a gradual decrease in concentration with increasing depth and a corresponding reduction in the area of high-concentration zones. Given that the study area predominantly consists of unconfined aquifers, it remains uncertain whether the pollution process has affected the groundwater at present. In order to more clearly restore the vertical pollution process of Cr in the past 30 years, this study will select the severely polluted drilling site S43 as a typical pollution point for analysis of its vertical profile.

3. Numerical Model

To establish a regional unsaturated zone groundwater flow model, the TOUGH2 V3.0/EOS7 module was used for the numerical simulation. This module, a part of the TOUGH2 software, can simulate multiphase fluid flow and phase transition processes, effectively reconstructing the migration paths and processes of pollutants in the soil and exploring the simulation of heavy metal damage in site soil.

3.1. Conceptual Model

Based on the vertical stratification observed in the soil of the study area, the primary soil composition consists of silt and silty clay with relatively high permeability. To better illustrate the groundwater level changes and the soil moisture content in the unsaturated zone, a columnar model was established at the S43 site for the study. Vertically, the unsaturated zone was refined into 16 layers, with an average thickness of 0.5 m per layer, reaching a depth of 8 m. Horizontally, the column was divided into a single grid of 20 m × 20 m. Based on soil properties, it was divided into three layers with thicknesses of 2 m, 4 m, and 2 m, respectively (Figure 3). The top 143 boundary of the model is set as an atmospheric pressure boundary, defined as a constant temperature and 144 pressure (1 atm) boundary. The bottom surface of the model is conceptualized as an impermeable boundary. Vertically, the model is divided into 16 layers. In addition, it is assumed that the initial condition of the model is a steady flow field, and the temperature in the study area remains constant. The simulation period covers 1992–2021, encompassing a total of 30 years.

3.2. Numerical Model

The core control equations of the TOUGH2/EOS7 module include fluid flow, heat transfer, mass transfer, and phase change processes. They typically describe the transport of multiphase fluids and heat in porous media, with the relevant control equations as follows (

Table 2. Governing equations for TOUGH2/EOS7.

Description	Governing Equations
Fluid Flow Equation	$q=-k\frac{\nabla p}{\mu }$
Heat Transfer Equation	$\rho C_p\frac{\partial T}{\partial t}+\rho C_p\left(\nu \cdot \nabla T\right)=\nabla \cdot \left(k_T\nabla T\right)+Q$
Mass Transfer Equation	$\frac{\partial \left(\varphi \rho C\right)}{\partial t}+\nabla \cdot \left(\rho C_v\right)=\nabla \cdot \left(D\nabla C\right)+S_C$
Phase Change Equation	f (P,T,composition) = 0

q is the fluid velocity [m/s], k is the permeability coefficient [m²], p is the pressure [Pa], $\mu$ is the fluid viscosity [Pa·s], T is the temperature [K], $\rho$ is the density [kg/m³], c_p is the specific heat capacity [J/(kg·K)], v is the velocity [m/s], k_T is the thermal conductivity [W/(m·K)], Q is the heat source term [W/m³], C is the concentration of substances [mol/m³], $\varphi$ is the porosity (dimensionless), D is the diffusion coefficient [m²/s], S_C is the source term [mol/(m³·s)], P is the pressure [Pa], and T is the temperature [K]; composition represents the composition of components, typically in terms of mole fraction or mass fraction (dimensionless).

).

3.3. Parameter Settings

To ascertain the permeability of typical soil layers within the investigation site and to assess the pollution prevention performance of the vadose zone, soil samples were taken from typical boreholes on the site during the investigation period and sent to the local geotechnical material testing center for permeability coefficient testing. The test results show that both the vertical and horizontal permeability coefficients in the study area are below 10⁻⁵ cm/s, indicating medium to strong pollution prevention performance. The soil layers have a relatively good blocking effect on the migration of pollutants, with specific parameters shown in

Table 3. Characteristic parameters of typical drilling soil layers of the study area.

Soil Sample Number	Depth (m)	Color	Geotechnical Classification	Permeability Coefficient
				Vertical (cm/s)	Horizontal (cm/s)
S43-0.5	0.5	Dark Brownish	Silt	3.43 × 10−6	2.80 × 10−6
S43-2	2.0	Brownish Red	Silty Clay	1.52 × 10−6	1.24 × 10−6
S43-4	4.0	Brown	Silty Clay	1.47 × 10−6	1.33 × 10−6
S43-6	6.0	Brownish Yellow (Brown)	Silty Clay	1.68 × 10−6	8.41 × 10−7
S43-8	8.0	Brownish Yellow	Silt	5.75 × 10−5	4.47 × 10−5

.

3.4. Sources

The groundwater replenishment in the model mainly considers precipitation infiltration, and the discharge of groundwater includes lateral runoff outside the area. This study obtained annual rainfall data from 1992 to 2021 as shown below (Figure 4). The data were obtained at monthly scale using remote-sensing methods combined with a literature search [35]. It shows that the rainfall in the study area is significant, with the maximum annual rainfall reaching 1224 mm/a and the minimum being 686 mm/a. In this study, based on empirical values, the infiltration recharge coefficient for the unsaturated zone in the model is set to 0.1.

4. Results and Discussion

4.1. Model Identification

4.1.1. Analysis and Validation of Liquid Saturation Changes in the Unsaturated Zone

The identification of the unsaturated zone liquid saturation variation and verification was carried out by simulating liquid saturation in the unsaturated zone. Figure 5 shows the soil vertical saturation changes in the years 1992, 1997, 2002, 2007, 2012, and 2021 after the model ran for 1, 5, 10, 15, 20, and 30 years, respectively. The calculation shows that, when the saturation is 1, the water level depth obtained from the model simulation is consistent with that calculated from the saturation.

Furthermore, the vertical liquid saturation of the boreholes conforms to the variation characteristics of the soil moisture characteristic curve (Figure 5). Overall, owing to the ample rainfall prevalent in the study area, the soil surface layer of the boreholes reaches a saturated state in approximately 15 years, and an obvious inflection point appears at a depth of 6~8 m. This is due to the presence of a perennial groundwater fluctuation zone at this depth, suggesting that the concentration increase is related to the accumulation of pollutants caused by water level fluctuations.

4.1.2. Verification of Vertical Migration Process of Pollutants

The vertical migration process of pollutants was verified based on the actual measured values of chromium concentrations at different depths in a single borehole in 2021. For the damage process simulation, this study chose an initial concentration injection of 5 × 10⁻⁸ kg/s. The simulation results when compared with the actual concentration measurements (Figure 6) show that, although the fit of pollutant concentrations at depths of 6~8 m is poor, the simulation concentrations elsewhere are very close to the measured concentrations. The possible reason is that this study only considered the migration path of pollutants in the unsaturated zone, ignoring other potential factors.

The analysis of the fitting results suggests that the regional unsaturated zone solute migration model established in this study exhibits high accuracy and reliability, accurately reflecting the pollution characteristics of the study area. This indicates that the model used in this study is suitable for the migration of pollutants in the unsaturated zone and can provide valuable references for related research in this field.

4.2. Scenario Design

According to the distribution and verification of pollutant concentrations, when chromium is injected at an initial concentration of 5 × 10⁻⁸ kg/s for 20 years, the simulation results match well with the actual measured values in the soil of the study area, suggesting that this injection scenario is close to the actual contamination process. Scenario analysis and simulation should be conducted around this scenario to obtain the damage process closest to the actual situation.

Taking into account the production history of the surveyed area, we note that it transitioned from being a phosphate fertilizer plant to being a leather factory in 1980, its initial production scale remained relatively small. Between 1980 and 1998, there were several changes in management, maintaining a low output for a long time. In 1999, the leather factory expanded its production capacity and started large-scale leather production but did not install corresponding environmental protection facilities. No professional anti-seepage work was carried out in the areas prone to pollution, such as the production and raw material areas, which could lead to soil contamination in case of pollutant leakage.

To identify the damage process, scenario design is carried out. By simulating the migration of chromium in the vadose zone under different scenarios, the model’s simulated values are fitted with the actual site measurements. The scenario with the best fit is identified as the one closest to the actual damage process. To obtain a damage scenario closer to reality, before carrying out formal scenario design and simulation, simulations of different orders of injection concentrations and pollution durations around the 5 × 10⁻⁸ kg/s injection for 20 years are conducted. The calculations show that injecting chromium in the range of 5 × 10⁻⁷ kg/s to 5 × 10⁻⁹ kg/s for 5 to 30 years results in a relatively good fit between the model simulations and actual measurements, while scenarios outside this range show a poorer fit, indicating that the actual damage scenario should be within the aforementioned range.

Therefore, combining the actual production history of the study area, a cross-scenario design is carried out with five different durations and three different injection concentrations (

Table 4. Simulation schemes for 15 environmental damage scenarios.

Injection Scheme	5 Years	10 Years	15 Years	20 Years	30 Years
Injection Concentration (kg/s)	5 × 10−7	5 × 10−7	5 × 10−7	5 × 10−7	5 × 10−7
	5 × 10−8	5 × 10−8	5 × 10−8	5 × 10−8	5 × 10−8
	5 × 10−9	5 × 10−9	5 × 10−9	5 × 10−9	5 × 10−9

) to simulate the chromium damage process in the soil of the study area under different scenarios.

4.2.1. Damage Process Identification

Through model calculations, the simulated concentrations of pollutants in boreholes for the above 15 pollution scenario schemes are compared and analyzed with the actual measured concentrations (as shown in Figure 7). Comparing the shapes of the simulation results with the actual measurements, when the pollutant release concentration reaches 5 × 10⁻⁷ kg/s, the simulated values for all five scenarios are higher than the actual values, especially near the surface where the pollutant concentration error is quite large, with a maximum error of up to 3300 mg/kg. This indicates that the pollutant release concentration in this simulated scenario does not match the actual pollution situation.

When the pollutant release concentration is 5 × 10⁻⁸ kg/s, the simulated values of the five different scenarios are close to the actual values, and the trend of pollutant concentration change with depth matches the trend of the actual measurements. In the 0–2 m underground interval, there is a greater range of pollutant variation due to the predominance of silty and silty clay soils at the surface, resulting in a lower vertical permeability rate of pollutants and, thus, slower infiltration into the ground. In the 2–8 m underground range, the trend of pollutant concentration change with depth is relatively stable. However, at approximately 6 m underground, an actual peak in pollutant concentration is observed, while the simulated peak appears at 7–8 m. This is because the numerical simulation only considered the movement of pollutants in the unsaturated zone and did not account for changes in pollutant concentration caused by groundwater fluctuations. In wet years in the study area, the groundwater fluctuation zone can reach 6–8 m underground, and the rise in the groundwater level promotes the appearance of pollutant peaks at shallower depths than the simulated values. Looking at the five different injection time scenarios, the longer the injection time, the more stable the trend of the simulated values. When continuously injected for approximately 20 years, the curve shape of pollutant concentration change with depth matches the actual measured values most closely.

When the pollutant release concentration is 5 × 10⁻⁹ kg/s, the simulated values for all five scenarios are less than the actual values. The error in pollutant concentration near the surface is quite significant, reaching approximately 200 mg/kg, and as the injection time increases, the impact on the simulated values is relatively small but still significantly different from the actual values, indicating that the pollution in this simulation scenario cannot replicate the actual situation.

In summary, the comparative analysis demonstrates that, when the pollutant release is at a concentration of 5 × 10⁻⁸ kg/s over a period of approximately 20 years, the simulation results match well with the actual results in terms of pollutant concentration magnitude, and the curve shape of pollutant concentration change with depth also matches the actual measured values.

4.2.2. Simulation Error Analysis

To more accurately perform error analysis, this study chooses four indicators to measure the accuracy of the simulation results: the coefficient of determination (R²), mean absolute error (MAE), root mean squared error (RMSE), and mean absolute percent error (MAPE). The larger the R² and the smaller the other three indicators indicate higher calculation accuracy. The formulas for calculating these evaluation indicators are as follows:

$${\mathrm{R}}^2=1 -\frac{\sum _i^n{(\widehat{y_i}-y_i)}^2}{\sum _i^n{\left(y_i -\stackrel{\mathrm{-}}{y}\right)}^2}$$

(1)

$$\mathrm{M}\mathrm{A}\mathrm{E}=\frac{1}{n}{\sum }_{i=1}^n\left|y_i-\widehat{y_i}\right|$$

(2)

$$\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}=\sqrt{\frac{1}{n}{\sum }_{i=1}^n{(y_i-\widehat{y_i})}^2}$$

(3)

$$\mathrm{M}\mathrm{A}\mathrm{P}\mathrm{E}=\frac{1}{n}{\sum }_{i=1}^n\left|\frac{y_i-\widehat{y_i}}{y_i}\right| \times 100\%$$

(4)

Here, i represents different depths, n represents the number of sample points, y_i represents the actual measured value, y_i represents the simulated value, and y represents the average value of the sample points.

Based on the formulas of the four evaluation indexes, the accuracy calculation results between the actual measured values and the simulated values of each scenario are shown in

Table 5. Calculation results of evaluation indexes of measured values and simulated values.

Injection Scheme	5 × 10−7 kg/s
	5 Years	10 Years	15 Years	20 Years	30 Years
R2	0.74469	0.64669	0.58536	0.39776	−0.37785
MAE	323.00	378.00	409.00	493.00	763.00
RMSE	684.53	805.26	872.36	1051.34	1590.23
MAPE (%)	33.04	37.25	40.38	51.11	115.45
Injection Scheme	5 × 10−8 kg/s
	5 Years	10 Years	15 Years	20 Years	30 Years
R2	0.98562	0.99965	0.99976	0.99982	0.99730
MAE	100.00	17.80	14.40	11.00	44.00
RMSE	162.48	25.43	21.18	18.34	70.34
MAPE (%)	153.44	5.58	4.73	4.58	13.97
Injection Scheme	5 × 10−9 kg/s
	5 Years	10 Years	15 Years	20 Years	30 Years
R2	0.99454	0.99476	0.99498	0.99518	0.99538
MAE	54.60	54.00	53.40	52.80	52.20
RMSE	100.06	98.03	96.02	94.02	92.05
MAPE (%)	8.30	8.38	8.45	8.52	8.59

and Figure 8.

5. Conclusions

This study presents a method for identifying the timing and process of damage in the unsaturated zone using numerical simulation with the TOUGH2 model. This approach addresses the lack of site parameters and the challenges of highly nonlinear solutions in the vadose zone. By combining parameter analysis, we established a vadose zone water flow model and solute migration model. Through multi-scenario simulation, this study identifies the most realistic damage process. The simulation of the chromium migration process in the vadose zone of the typical pollution point in the study area suggests that a release concentration of 5 × 10⁸ kg/s for 20 years closely approximates the actual heavy metal damage process in the soil. The simulation results meet the model’s accuracy requirements, with the best scenario achieving a fit of 0.9998. Its mean absolute error, root mean squared error, and mean absolute percent error are the smallest among the scenarios, and the simulated values correspond with the actual pollution intensity in the soil of the study area. This reliability of the simulation results effectively resolves the difficulty of identifying the process of soil heavy-metal damage. However, in this study, due to the small simulation area, the impact of groundwater flow velocity was ignored. Additionally, the model also ignored the effects of microbial degradation and changes in temperature and humidity.

Existing studies and current technical guidelines for soil environmental damage assessment in China lack specific methods for identifying soil damage processes. Through method identification, application, and verification, this study constructs a method for identifying the process of soil heavy-metal damage, effectively resolving the challenges associated with process identification. In future research, plans are in place to conduct simulation studies at a regional scale, considering the establishment of complex simulations of aquifer heterogeneity across the area. Additionally, if conditions permit, it is hoped to include various influencing factors, such as microbial degradation, climate change, and groundwater level fluctuations, in the model to achieve more realistic simulation results and improve accuracy. Identifying the damage process is a crucial step in determining causality, and this innovative research provides scientific reference and data support for assessing the causal relationship of soil heavy-metal pollution.