Study on Fluorine Pollution in a Slag Yard

Abstract

This paper suggests a large-scale three-dimensional numerical simulation method to investigate the fluorine pollution near a slag yard. The large-scale three-dimensional numerical simulation method included an experimental investigation, laboratory studies of solute transport during absorption of water by soil, and large-scale three-dimensional numerical simulations of solute transport. The experimental results showed that the concentrations of fluorine from smelting slag and construction waste soil were well over the discharge limit of 0.1 kg/m³ recommended by Chinese guidelines. The key parameters of the materials used for large-scale three-dimensional numerical simulations were determined based on an experimental investigation, laboratory studies, and soil saturation of survey results and back analyses. A large-scale three-dimensional numerical simulation of solute transport was performed, and its results were compared to the experiment results. The simulation results showed that the clay near the slag had a high saturation of approximately 0.9, consistent with the survey results. Comparison of the results showed that the results of the numerical simulation of solute transport and the test results were nearly identical, and that the numerical simulation results could be used as the basis for groundwater environmental evaluation.

1. Introduction

Fluoride can be harmful to biology. Fluoride is toxic to vegetation and feeding animals [1]. Excess fluoride intake causes different types of fluorosis, primarily dental and skeletal fluorosis, depending on the level and period of exposure [2]. Fluorine concentration may reach contamination levels in the groundwater as a result of specific hydrogeochemical processes due to human activities. Elevated concentrations of fluorine in groundwater, exceeding the World Health Organization (WHO) limit of 0.0015 kg/m³ in drinking water, have been detected in many parts of the world, such as India [2], China [3], Tanzania [4], and others.

The finite element approach was first used for solute transport analysis by Voss [5]. To obtain an understanding of solute transport over longer time periods and under more complicated liner conditions, numerical modelling of solute transport in saturated–unsaturated fluid flow was conducted using the finite element method (FEM) [6,7,8,9]. However, research of FEM simulation of groundwater fluoride concentration was rare.

Most previous studies on fluorine pollution have focused on the experimental investigation of hydrogeochemical processes impacting fluoride in aquifers [1,10,11,12,13,14,15], such as analytical work of fluoride abundance based on experimental fluoride concentrations of samples [1,10,11,13], analytical work of mineralogical sources of groundwater fluoride in Archaen bedrock/regolith aquifers based on leaching experiments [12], controls of high fluoride groundwater and assessment of processes controlling the regional distribution of fluoride based on experiments [14,15]. Nevertheless, the identification of these processes remains challenging in areas with complex hydrogeochemical or geo-hydrogeological conditions. Research such as large-scale three-dimensional numerical simulations of groundwater fluoride concentration is rare. In this study, the numerical simulation method was used to predict the transport of fluoride pollutant of a hazardous slag yard in southwest China.

In this paper, a large-scale three-dimensional numerical simulation method to investigate the fluorine pollution near a slag yard was performed. The key parameters of the materials used for the large-scale three-dimensional numerical simulation were determined based on the field investigations, the laboratory studies, and back analyses. The results were compared to the experiment results. The research of fluorine migration in water near a slag yard was presented to provide a scientific theory basis for the management of the slag yard [16]. This method could be used as a reference for similar studies.

2. Materials and Methods

The study site is located at a hazardous slag yard in southwest China. Ground layers in the study area are mainly composed of clay and carbonate rocks. Slag has been stacked over a clay layer without disposal since 1958. The material in the slag yard was mainly composed of smelting slag and construction waste soil. The typical matters in the slag yard are shown in Figure 1.

The hazardous materials in the slag yard included smelting slag and construction waste soil. Three samples of leaching solution from these materials were collected, according to the China criterion for identification of extraction toxicity [17] and were tested.

To obtain the real leaching water of rainfall in the slag yard, some perforated PVC tubes were made and buried in the bottom of the slag yard. Two samples of the real leaching water, from sampling sites CK-1 and CK-2 in Figure 1, were collected and tested.

Water quality near the slag field was collected and analyzed to evaluate the slag yard effects on surrounding groundwater. The groundwater sampling sites were ZKJ-1, ZKJ-2, and ZKJ-3, shown in Figure 1.

Table 1. Concentration (C) of fluorine from the test results.

Sample	Plane Site	Depth	C/(kg/m3)
ZKJ-1	Distance away from the slag along y-direct-ion is 1457 m.	20 m	0.00005
ZKJ-2	Distance away from the slag along x-direct-ion is −105.67 m.	148 m	0.00006
ZKJ-3	Distance away from the slag along y-direct-ion is −286 m.	130 m	0.00007

also lists the groundwater sampling sites. The water samples were obtained from the wells.

3. Results of the Experimental Investigation

The three samples of leaching solution from the hazardous matters were tested, and the concentration of fluorine in each sample is shown in

Table 2. Concentration of fluorine.

Ground Layer	Test Results C (kg/m3)	Limit Value C (kg/m3)
Smelting slag	0.152	0.100
Construction waste soil	0.340	0.100
Clay in the bottom of the slag yard	0.0024	0.100

. The test results showed that the concentrations of fluorine from smelting slag and construction waste soil all exceeded the discharge limit of 0.1 kg/m³ of the China criterion [18]. The test results also showed that the clay in the bottom of the slag yard absorbed some fluorine.

The concentrations of fluorine in the leaching water samples from the slag yard after some days of rainfall are shown in

Table 3. Concentration of fluorine.

Sample	CK-2	CK-3
Measuring C (kg/m3)	0.00446	0.00447
Remarks	15 days after rainfall	30 days after rainfall

. The test results showed that the concentrations were stable approximately 15 days after rainfall. The results also showed that the concentrations of fluorine were well over the World Health Organization (WHO) limit of 0.0015 kg/m³ in drinking water and the China health limit of 0.001 kg/m³ in drinking water [19].

Table 1 shows the test results of the groundwater samples. The concentrations of fluorine did not reach the China health limit of 0.001 kg/m³ in drinking water.

The test results showed that the concentrations of fluorine from smelting slag and construction waste soil exceeded the discharge limit of 0.1 kg/m³ of the China criterion [18] and the concentrations of fluorine of the real leaching water of rainfall in the slag yard. Thus, fluorine was a contaminant, and it was necessary to study the scope of fluorine pollution of the slag yard.

4. Laboratory Studies of Solute Transport During Absorption of Water by Soil

The experimental studies of hydrodynamic dispersion of low concentration solutions during absorption into horizontal columns of soil were made for the clay under the slag yard. The main objective of these studies was to obtain the molecular diffusion coefficient of the clay for numerical simulation.

4.1. Theory for the Laboratory Studies

The conventional approach to hydrodynamic dispersion during one-dimensional water movement in soils involves the equation [20]:

$$\frac{\partial \left(\theta C\right)}{\partial t}=\frac{\partial }{\partial x}\left[\theta D_{sh}\left(\theta \right)\frac{\partial C}{\partial x}\right]-\frac{\partial \left(qC\right)}{\partial x}$$

(1)

in which C is the solute concentration in the water, θ is the volume fraction of fluid within the soil, q is the volume flux density of the solution, x is the space coordinate, t is time, and D_sh is the longitudinal dispersion coefficient.

The experiments were designed to realize the initial and boundary conditions:

$$\theta ={\theta }_i,C=C_i,x>0,t=0$$

(2a)

$$\theta ={\theta }_0,C=C_0,x=0,t>0$$

(2b)

In horizontal non-hysteretic systems, the flow of water may be described by the equations:

$$q=-D\left(\theta \right)\frac{\partial \theta }{\partial x}$$

(3)

$$\frac{\partial \theta }{\partial t}=\frac{\partial q}{\partial x}$$

(4)

where D(θ) is the soil–water diffusivity.

Setting $\lambda =x{t}^{-1/2}$, the following equations are derived [20]:

$$D\left(\theta \right)=-\frac{1}{2}\frac{d\lambda }{d\theta }{\displaystyle {\int }_{{\theta }_i}^{\theta }\lambda }d\theta$$

(5)

$$g\left(\theta \right)=\theta \lambda +2D\left(\theta \right)\frac{d\theta }{d\lambda }$$

(6)

$$D_{sh}(\theta )=-\frac{1}{2\theta }\frac{d\lambda }{d{C}^{\prime }}{\int }_0^{C^{\prime }}g(\theta )d{C}^{\prime }$$

(7)

$$C^{\prime }=\frac{C-C_i}{C_0-C_i}$$

(8)

It can be assumed that D_sh(θ) is a function principally of θ and that any dependence on q could be neglected. Thus, the following equations are used:

$$\mathsf{\theta }{D}_{sh}(\theta )=\mathsf{\theta }{D}_m{\mathsf{\theta }}^{\frac{7}{3}}/{\mathsf{\theta }}_r$$

(9)

$$D_m={\mathsf{\theta }}_{\mathrm{r}}{D}_{sh}(\theta )/{\mathsf{\theta }}^{\frac{7}{3}}$$

(10)

in which D_m is the molecular diffusion coefficient and θ_r is the residual water volume fraction of soil.

4.2. Tests and Results

Two clay samples designed to have the same dry bulk density as a clay sample in site were used for experiments. The test conditions are shown in

Table 4. Test conditions.

Soil Sample Number	1	2
Test time (min)	720.45	670.73
Total wetting front distance (cm)	24	24
Dry bulk density (g/cm3)	1.22	1.22

. The main experiment results are shown in Figure 2, Figure 3 and Figure 4 and

Table 5. Test results.

Number	Θ (cm3/cm3)	D (cm2/min)	Dsh (cm2/min)	Dm (m2/s)
1	0.5172	0.4747	0.7013	1.36 × 10−7
2	0.5388	1.1872	1.6810	2.96 × 10−7
3	0.5491	1.8559	5.2774	8.89 × 10−7
4	0.5572	2.1421	10.6023	1.73 × 10−6
5	0.5651	3.1772	18.4831	2.91 × 10−6
6	0.5692	4.6549	28.7940	4.46 × 10−6
7	0.5726	5.3263	36.1603	5.53 × 10−6
8	0.5757	5.7641	44.3241	6.69 × 10−6
9	0.5786	6.3421	60.4441	9.02 × 10−6
10	0.5812	7.5682	80.4299	1.19 × 10−5
11	0.5833	8.4745	89.8375	1.31 × 10−5

.

5. Numerical Modelling Investigation

The fate and transport of pollution were simulated by COMSOL Multiphysics, which is a groundwater flow and solute transport visual software simulating the movement of solutes. Based on the program COMSOL Multiphysics, a model was built for solute transport in porous media, and the following equations were used in this model [21].

5.1. Fluid Flow

Richards’ equation governs the saturated–unsaturated flow of water in the soil. The soil air is open to the atmosphere, so it can be assumed that pressure change in air does not affect the flow. Here, Richards’ equation is used for single-phase flow. Richards’ equation can be expressed as:

$$(C_h+S_e\times S)\frac{\partial H_p}{\partial t}+\nabla ·(K\nabla (H_p+E))=Q_s$$

(11)

where $C_h$ denotes the specific moisture capacity (m⁻¹); S_e is the effective saturation of the soil (dimensionless); S is a storage coefficient (m⁻¹); H_p represents the dependent variable pressure head (m); K equals the hydraulic conductivity function (m s⁻¹); E is the coordinate (for example, x, y, or z) that represents vertical elevation (m); and Q_s is a fluid source defined by volumetric flow rate per unit volume of soil (s⁻¹). In this equation, S = (θ_s − θ_r)/(1 m) where θ_s denotes the volume fraction of fluid under saturation.

For unsaturated flow, the fluid moves through the soil but may or may not completely fill the pores in the soil. The coefficients C_h, S_e, and K vary with the pressure head, and H_p, together with θ, makes Richards’ equation nonlinear. The specific moisture capacity, C_h, relates the variations in soil moisture to pressure head as C_h = ∂θ/∂H_p. Because C_h, the first term in the time coefficient, goes to zero at saturation, time change in storage depends on the compression of the aquifer and water under saturated conditions. The saturated storage varies with the effective saturation, as represented by the second term in the time-coefficient. Furthermore, K is a function that defines how readily the porous media transmits fluid. The relative permeability of the soil, kr, increases with increasing fluid content given by K = K_sk_r, where Ks (m/s) is the constant hydraulic conductivity at saturation.

Groundwater flow and solute transport can be linked by fluid velocities. To follow the form of the transport equation, the fluid velocities need to be deduced from Darcy’s law:

$$u=-k_r{K}_S\nabla (H_p+D)$$

(12)

5.2. Solute Transport

The equation that governs advection and dispersion is:

$$\frac{\partial }{\partial {\mathrm{x}}_{\mathrm{i}}}\left[\theta {\mathrm{D}}_{\mathrm{Lij}}\frac{\partial \mathrm{c}}{\partial {\mathrm{x}}_{\mathrm{j}}}\right]-\frac{\partial }{\partial {\mathrm{x}}_{\mathrm{i}}}\left({\mathrm{cu}}_{\mathrm{i}}\right)+\theta {\displaystyle \sum _{\mathrm{k}=1}^{\mathrm{n}}{\mathrm{R}}_{\mathrm{k}}}=\frac{\partial \left(\theta \mathrm{c}\right)}{\partial \mathrm{t}}$$

(13)

where D_L is the hydrodynamic dispersion tensor and is referred to as D_sh(θ) in the preceding paragraph; R_k is the added solute per unit volume of water per unit time; and x_i is the Cartesian coordinate. So, θD_L in the three-dimensional case can be expressed as [20]:

$$\{\begin{array}{l}\theta {\mathrm{D}}_{\mathrm{Lxx}}={\alpha }_1\frac{u_{}^2}{\left|u\right|}+{\alpha }_2\frac{v_{}^2}{\left|u\right|}+{\alpha }_3\frac{w_{}^2}{\left|u\right|}+\theta D_m{\theta }^{7/3}/{\theta }_s{}^2\\ \theta {\mathrm{D}}_{\mathrm{Lyy}}={\alpha }_2\frac{u_{}^2}{\left|u\right|}+{\alpha }_1\frac{v_{}^2}{\left|u\right|}+{\alpha }_3\frac{w_{}^2}{\left|u\right|}+\theta D_m{\theta }^{7/3}/{\theta }_s{}^2\\ \theta {\mathrm{D}}_{\mathrm{Lzz}}={\alpha }_1\frac{w_{}^2}{\left|u\right|}+{\alpha }_3\frac{v_{}^2}{\left|u\right|}+{\alpha }_3\frac{u_{}^2}{\left|u\right|}+\theta D_m{\theta }^{7/3}/{\theta }_s{}^2\\ \theta {\mathrm{D}}_{\mathrm{Lxy}}=\theta {\mathrm{D}}_{\mathrm{Lyx}}=({\alpha }_1-{\alpha }_2)\frac{uv}{\left|u\right|}\\ \theta {\mathrm{D}}_{\mathrm{Lxz}}=\theta {\mathrm{D}}_{\mathrm{Lzx}}=({\alpha }_1-{\alpha }_3)\frac{uw}{\left|u\right|}\\ \theta {\mathrm{D}}_{\mathrm{Lzy}}=\theta {\mathrm{D}}_{\mathrm{Lyz}}=({\alpha }_2-{\alpha }_3)\frac{vw}{\left|u\right|}\end{array}$$

(14)

5.3. Geological Model for the Hydrogeology Study in the Study Area

The ground layers were divided into six materials mainly based on their permeability difference, as follows [22]: material 1 was mainly composed of slag and construction waste soil; material 2 was clay; material 3 was unsaturated limestone located a certain distance away from groundwater; material 4 was saturated limestone within the scope of groundwater; material 5 was base rock with very low permeability; material 6 was sandy gravel. Material 3 and material 4 offered different permeability coefficients, which are shown in

Table 6. Model parameters for flow simulation.

Material	mat.1	mat.2	mat.3	mat.4	mat.5	mat.6
kv/m/s	0.0001		0.00001	0.001	10−6	0.001
kh/m/s	0.0001		0.00001	0.001	10−6	0.001
Pore volume fraction	0.32	0.571	0.25	0.25	0.20	0.35
Residual water volume fraction	0.06	0.158	0.01	0.01	0.01	0.01

.

Two simulation models based on exploratory boring are shown in Figure 5 and Figure 6 [22]. The 2D simulation model in Figure 5 was 10.462 km in length and 1.058 km in height and jointed by a lake. The 3D simulation model in Figure 6 was modelled by 8 profiles and was 1.840 km in length along the y direction, 0.89 km in width along the x direction, and 0.439 km in height. The layout of the ground layers was shown in Figure 5a,b and Figure 6a. Materials 1–4 of the 2D-model and the 3D-model were the same, but the 3D-model did not include material 5 because the 3D-model was shallower than the 2D-model.The 2D model was divided into 42,236 elements with 21,648 nodes, while the 3D model was divided into 67,968 elements with 14,620 nodes.

5.4. Boundary and Model Parameters for Flow Simulation

The seepage field must be simulated as accurately as possible to obtain reasonable simulation results of contaminant migration. Here, the seepage field was simulated as saturated–unsaturated stable seepage. For the model in Figure 5, the boundary conditions were set as follows [22]. The upper surface accepted rainfall infiltration and belonged to the flux boundary. The flux at the slag surface was 3.2 × 10⁻⁸ m/s, which was equal to the local average annual rainfall of 1012 mm. The flux at the clay surface was 1.6 × 10⁻⁹ m/s, which was 5% of the rainfall. The flux of the slag surfaces was 3.2 × 10⁻⁸ m/s because of the following reasons. First, the saturated permeability coefficients of the slag (0.0001 m/s and 0.0001 m/s) were very high so the rain water on the slag surface easily flowed into the slag and could be stored in the lower slag. Second, the clay surface under the slag was sunken so it could keep the water from the slag surface and some water from the clay surface. Finally, the water stored in the lower slag was covered by thick slag, so the evaporation of the rain water to atmosphere was very little. The flux of the clay surfaces was 1.6 × 10⁻⁹ m/s because of the following reasons. The saturated permeability coefficients of the clay were very low, so most of the rain water on the clay surface was easy to run off. Some of the water evaporated off, so the water penetration through clay surface would be very low. The two border surfaces vertical to the x-axis were set as boundaries with constant hydraulic head: one of the constant hydraulic heads was the hydraulic head in the exploratory hole; the other was the hydraulic head of the lake at the boundary. The bottom surface was set as a boundary with a constant flux of zero.

The seepage field of the 3D model in Figure 6 was also simulated as saturated–unsaturated stable seepage, and the following boundary conditions of seepage simulation were used [22]. The upper surface accepted rainfall infiltration, and the border belonged to flux boundary. The flux of the slag surfaces was 3.2 × 10⁻⁸ m/s while that of the clay surfaces was 1.6 × 10⁻⁹ m/s. The two border surfaces vertical to the x-axis were set as symmetric boundaries because they belonged to watersheds. The two border surfaces vertical to the y-axis were set as boundaries with constant hydraulic head: one was the hydraulic head in the exploratory hole on the boundary; the other was −61 m based on the simulation results of 2D FEM model. The bottom surface was set as a boundary with a constant flux of zero.

Soil–water characteristic curves for material 1 and material 2 were obtained from a volumetric pressure plate test, and the pore volume fractions were estimated from the corresponding soil–water characteristic curves. The residual water volume fractions for material 1 and material 2 were obtained from the drying method. The saturated permeability coefficients for material 1 and material 4 were estimated from a water injection test. Other parameters were determined according to the experience values in geotechnical engineering. All the saturated permeability coefficients of the materials, except those of material 2, are shown in Table 6 for simulation. The permeability coefficients of material 2 were important for determining the seepage field. Therefore, they were determined by inverse analysis based on the test results. Materials 1 and 2, whose soil–water characteristic curves were obtained from the test results, were set to be unsaturated materials. The soil–water characteristic curves of material 2 are shown in Figure 7, and the experimental curve of sample 82 was used for simulation. The test results indicated that the saturated permeability coefficient of material 2 in the vertical direction ranged from 5 × 10⁻⁹ m/s to 1 × 10⁻⁷ m/s. Simulation analysis showed that the simulation results of the seepage field could fit the actual seepage field better if the saturated permeability coefficient in the vertical direction was 8 × 10⁻⁸ m/s and that in the horizontal direction was 3 × 10⁻⁸ m/s.

5.5. Model Definition for Dissolved Solute Transport

The solute component identification test showed that the concentration value of fluorine in the slag solute was approximately 0.0045 kg/m³, while that of the solute on the clay surface was approximately 0.00005 kg/m³.

The following conditions were set for solute migration simulation. The initial value of slag solute was set to 0.0045 kg/m³, while that of other solutes were set to be 0.00005 kg/m³. Boundary conditions were as follows. Flux boundary was set to the upper surface of the model. The flux of the slag yard surfaces was (3.2 × 10⁻⁸) × 0.0045 kg m⁻² s⁻¹, while that of the clay surfaces was (3.2 × 10⁻⁸) × 0.00005 kg m⁻² s⁻¹. The concentration of fluorine in the water of the border surface with the higher constant hydraulic heads was set to 0.00005 kg/m³. The advective flux model was set to the border surface with the lower constant hydraulic heads. The flux of the bottom surface was zero. The fluxes of the two border surfaces vertical to the x-axis were set to be zero for the 3D FEM model.

The parameters of the materials used in the dispersion calculation are shown in

Table 7. Model parameters for pollutant migration.

Material	mat.1	mat.2	mat.3	mat.4	mat.5	mat.6
dm, m2 s−1	10−9	4.85 × 10−6	10−9	10−9	10−9	10−9
${\alpha }_1$, ${\alpha }_2$, m	1	0.01	2	2	2	3
${\alpha }_3$, m	2	0.02	2	2	2	3

. The molecular diffusion coefficient of material 2 was determined as D_m = 4.85 × 10⁻⁶ m² s⁻¹ based on the unsaturated dispersion test results shown in Table 5 and the clay saturation of survey results shown in

Table 8. The clay saturation of survey results.

Soil sample No	Depth/m	Saturation/%
1754	7.50–7.70	100
1755	17.90–18.10	99
1756	14.30–14.50	98
1757	12.90–13.10	100
1758	12.00–12.20	95
1759	8.60–8.80	86
1760	18.20–18.50	97
1761	18.20–18.40	99
1762	15.70–15.90	95
1763	20.70–20.90	95
1764	21.60–21.80	98
1765	25.40–25.60	93

. Other parameters, shown in Table 7, were determined according to the experience values in geotechnical engineering.

5.6. Results

Survey results showed that clay had a high saturation of approximately 0.95, as shown in Table 8. The simulation results of the 2D model showed that the saturation of clay near the slag in Figure 8 ranged from 0.89 to 0.95. The simulation results of the 3D model showed that the saturation of clay near the slag in Figure 8 ranged from 0.89 to 0.91. The saturation of clay near the slag in Figure 8 was approximately 0.9. The simulation results were consistent with survey results. The results of the calculations showed that the relative error between the simulation results and the survey results was less than approximately 10%. According to the Technical Requirements for Groundwater Numerical Simulation [23], a relative error below 10% is acceptable. Therefore, the formation model established for this study was feasible for the simulation.

Groundwater quality near the slag field was analyzed to evaluate the slag field effects on the surrounding groundwater.

Table 9. Comparison between test results and simulation results.

Sample	Test Result C/(kg/m3)	Simulation Result C/(kg/m3)
ZKJ-1	0.00005	0.000053
ZKJ-2	0.00006	0.000051
ZKJ-3	0.00007	0.000054

shows the groundwater sample test results and corresponding numerical simulation results. It can be seen that the numerical simulation results and test results are nearly identical. Hence, the numerical simulation results could be used as the basis of groundwater environment evaluation.

Figure 9 provides some solute migration simulation results of the 2D model. Over time, the solution with the fluorine concentration in the slag field migrated outward slowly. The dispersion velocity near the groundwater was so large that the fluorine concentration in the groundwater did not increase with time. The final fluorine concentration in the groundwater, 100 years after disposal started, did not exceed 0.001 kg m⁻³, which is the legislative upper limit of the second-grade water of China. The concentration in the groundwater did not exceed the legislative upper limit of the second-grade water of China. The maximal pollution water range, from the slag yard to the south, 52 years after disposal started, was approximately 123 m, and was approximately 172 m 100 years after disposal started.

Figure 10, Figure 11 and Figure 12 provide some solute migration simulation results of the 3D model. Over time, the solution with the fluorine concentration in the slag field migrated outward slowly. The dispersion velocity near groundwater was so large that the fluorine concentration in the groundwater did not increase with time. The final fluorine concentration in the groundwater 100 years after disposal started did not exceed 0.001 kg m⁻³. The concentration in the groundwater did not exceed the legislative upper limit of the second-grade water of China. The maximal pollution water range, from the slag yard to the south, 52 years after disposal started, was approximately 185 m, and the range 100 years after disposal started was approximately 330 m. The maximal pollution water range, from the slag yard to the west, 52 years after disposal started, was approximately 75 m, and the range 100 years after disposal started was approximately 135 m. The maximal pollution water range, from the slag yard to the east, 52 years after disposal started, was approximately 60 m and the range 100 years after disposal started was approximately 80 m.

5.7. Discussion

Comparing to the study in the literature, a large-scale three-dimensional numerical simulation of groundwater fluoride concentration was performed, and the maximal fluoride pollution water range was attempted to be determined based on a large-scale three-dimensional numerical simulation of groundwater fluoride concentration for the first time.

5.8. Limitations

The limitation of this study was that the sorption of fluoride in clay and the mineralogical sources of groundwater fluoride from the matters out of the slag were not considered. The fluorine concentrations would be less if the sorption of fluoride in clay and the mineralogical sources of groundwater fluoride from the matters out of the slag were considered [16].

6. Conclusions

In this paper, a large-scale three-dimensional numerical simulation method to investigate the fluorine pollution near a slag yard was implemented. The large-scale three-dimensional numerical simulation method included experimental investigation, laboratory studies of solute transport during absorption of water by soil, and a large-scale three-dimensional numerical simulation of solute transport. Experimental investigations were performed, and the test results showed that the concentrations of fluorine from smelting slag and construction waste soil exceeded the discharge limit of 0.1 kg/m³ of the China criterion. The key parameters of the materials used for numerical simulation were determined based on the experimental investigation, the laboratory studies, the soil saturation of survey results, and back analyses. A large-scale two-dimensional numerical simulation and a large-scale three-dimensional numerical simulation of solute transport were implemented, and their results were compared to the experimental results. The following conclusions were made.

(1) The simulation results showed that the clay near the slag had high saturation of approximately 0.9, which was consistent with survey results.

(2) Comparison of the results showed that the numerical simulation results of solute transport and test results were basically identical. Thus, the numerical simulation results could be used as the basis of groundwater environment evaluation.

(3) Over 100 years of simulation, the pollution plume from the slag field migrated outward slowly. The dispersion velocity near the groundwater was so large that the fluorine concentration in the groundwater did not increase with time. The final fluorine concentration in the groundwater 100 years after disposal started did not exceed 0.001 kg m⁻³, the legislative upper limit of the second-grade water of China. The maximal pollution water range from the slag yard to the south, 52 years after disposal started, was approximately 185 m, and the range 100 years after disposal started was approximately 330 m. The maximal pollution water range from the slag yard to the west, 52 years after disposal started, was approximately 75 m. The range 100 years after disposal started was approximately 135 m. The maximal pollution water range from the slag yard to the east, 52 years after disposal started, was approximately 60 m, and the range 100 years after disposal started was approximately 80 m.

To compare with the study in the literature, a large-scale three-dimensional numerical simulation of groundwater fluoride concentration was performed, and the maximal fluoride pollution water range was attempted to be determined based on a large-scale three-dimensional numerical simulation of groundwater fluoride concentration for the first time.