Abstract

By embedding a nonlinear nonequilibrium adsorption model considering the adsorption history (i.e., Bai model), a governing equation for the three-dimensional model is extended, which discussed the influence of varying seepage velocity and injection concentration. Compared with the previous linear adsorption model, the concentration peak value of the nonlinear nonequilibrium adsorption history model is higher than that of the linear model, and the time to reach the peak concentration is slightly earlier. In the case of horizontal seepage water, a pollution point source tends to migrate along the direction of water flow and has less ability to diffuse in the vertical direction. Compared with the adsorption history model, the linear model has a stronger blocking ability for pollutant migration, and the longer the time, the greater the gap between the two models. The longer the decay period, the wider the spread of contaminants, and the longer it takes for them to migrate out of the model completely. The more significant the head difference, the larger the diffusion area of pollutants in the main seepage direction, but it has no promotion effect on the lateral diffusion. The pollutant concentration is higher than that of the point source case, and the diffusion range in each section is wider. The closer it is to the center of the pollution source, the weaker the dispersion effect on the diffusion of pollutants in the main seepage direction, and the pollutants will spread in the countercurrent direction, causing pollution upstream.

1. Introduction

The migration of pollutants in the soil is affected by many conditions, such as seepage velocity, dispersion conditions, soil adsorption capacity, chemical factors, and biological factors [13]. The adsorption mechanism of soil is complicated for its complex nature [4, 5], and it also has specific characteristics for desorption. Previous studies have shown that the adsorption capacity of porous media does not change after the initial adsorption [68], but in practice, the residue of pollutants in porous media will weaken its adsorption capacity [911]. Diverse soil conditions will result in complex adsorption-desorption processes, which in turn affect the migration of pollutants. To better investigate the migration law of pollutants in soil, a more realistic adsorption-desorption model is needed to support the theoretical research [12, 13], which is the basis for the development of pollutant purification technology [8, 14].

Kim et al. [15] studied the influence of hydraulic conditions on the adsorption properties of porous media through experiments. According to different chemical conditions, it was found that the magnitude of particle separation increased with the improvement of particle size. Hosseini et al. [16] developed as an economical and alternative adsorbent to treat dye wastewater, in which the magnetic chitosan nanocomposites possessed of excellent adsorptive property for the recycle of wastewater contaminated by dye. García-García et al. [17] studied the effect of temperature on particle aggregation dynamics under different pH values and ionic strengths. They found that the particle coagulation rate decreases with the temperature increase when the environment is alkaline. Lake and Rowe [18] found through laboratory tests that the diffusion coefficient is related to the concentration of the pollution source. Cui et al. [19] explored the effects of particle sizes, particle concentrations, penetration velocities, and penetration directions on the dispersion and deposition mechanisms of suspended particles through cylindrical penetration tests. Porubcan and Xu [20] discussed the migration characteristics of colloidal particles with different particle sizes in the porous media of quartz sand mixed with different particle sizes through experiments. They found that when the proportion of fine sand in uneven quartz sand increased, the adsorption of colloidal particles increased. Mahmoodi et al. [21] investigated the hydration inhibitive ability of surfactants by various experiments, including zeta potential, XRD, FTIR, particle size, and SEM measures for a dye removal by the surfactant-modified montmorillonite.

Cui et al. [22] found that different metal ions have different migration characteristics through laboratory research. The difference is mainly reflected in the adsorption properties, resulting in apparent differences in the migration process of different metal ions. Cui et al. [23] considered that the attenuation of soil adsorption performance is affected by its adsorption historical and temperature effect and conducted a sensitivity analysis on relevant parameters [24]. Some researchers [9, 25, 26], based on the hydrodynamic dispersion coefficient and flow velocity under unstable conditions, established their pollutant transport equations in semi-infinite space concerning space-time changes and obtained the analytical solution under the action of the Gaussian pulse pollution source. Cui et al. [22] proposed a short-time pulse pollution source by studying pollution sources, established a pollutant migration model, and calculated the migration and deposition characteristics of pollutants in porous media. Leij et al. [27] studied the migration law of pollutants during equilibrium adsorption, established a semi-infinite equation, and obtained the analytical solution of a unidirectional stable seepage field through Laplace and Fourier transform. Dahaghi et al. [28] deduced a new mathematical model of pollutant migration, which considers the complex migration process of particles. Rowe and Bookeret [29] proposed a two-dimensional finite layer calculation method considering vertical and horizontal directions and obtained the migration law of layered soil under the condition of linear equilibrium adsorption.

Miracapillo and Ferroni [30] considered the nonuniform groundwater flow to simulate the migration process of nuclides in the geological repository. They discussed the influence of distribution coefficient and hydraulic gradient on pollutant migration, respectively. Trofa et al. [31] used the arbitrary Lagrange–Eulerian method for numerical analysis to solve the governing equation of pollutant transport and numerically simulated the two-dimensional migration of suspended particles at an infinite Reynolds number. Villone et al. [32] used the arbitrary Lagrange–Eulerian method to numerically analyze the model from a three-dimensional perspective and obtained the corresponding numerical analysis results. Mironenko and Pachepsky [33] considered the exchange process and biochemical reaction process of solutes in immobile and movable water and established a pollutant migration model under the condition of a one-dimensional semi-infinite soil column. Tiraferri et al. [1] first determined the functional relationship between particle filtration coefficient and pore water ions through the deposition test of particles in siliceous sand and then established the corresponding model.

To explore the transport mechanisms of pollutants under regular boundary conditions in three-dimensional spatial, based on the in-depth analysis of typical adsorption models, a nonlinear nonequilibrium adsorption model (i.e., Bai model) [2] considering the adsorption history is embedded to carry out the theoretical calculation of pollutant migration. Classified by pollution source diffusion, the three-dimensional point source model and three-dimensional volume source model were discussed, respectively. This research analyzed the pollution source changes and parameter changes and obtained pollutant migration characteristics under the corresponding model.

2. Theoretical Model

2.1. Governing Equations

The governing equation has been derived as [2, 25]where C is the concentration of pollutants in porous media (ML−3), θ is the porosity of porous media, t is the time (T), D is the hydrodynamic dispersion coefficient (T2T−1), u is the average velocity of cross section (LT−1), is the volumetric dry density of porous media (ML−3), and Cd is the concentration of adsorption (MM−1).

Equation (1) refers to the amount of solute adsorbed by the porous medium per unit mass. In the governing equation, the expression of the dispersion term is . The mechanical dispersion coefficient tensor and the molecular diffusion coefficient tensor are combined into the hydrodynamic dispersion coefficient tensor, i.e.,where is the hydrodynamic dispersion coefficient tensor (L2T−1); is the mechanical dispersion coefficient tensor (L2T−1); , represents the molecular diffusion coefficient tensor (L2T−1), is the molecular diffusion coefficient scalar, is the curvature tensor of porous media.

For three-dimensional cases, assume that the seepage direction of the fluid is horizontal, that is,  = 0, D can be obtained , , so the hydrodynamic dispersion coefficient matrix is

The expression of hydrodynamic dispersion coefficient tensor can be obtained when the main direction of seepage velocity in porous media [3436] is consistent with the selected coordinates, and molecular diffusion is ignoredwhere is the longitudinal dispersion of porous media (L), and is the transverse dispersion of porous media (L),

It is assumed that the velocity in the plane direction is mainly caused by pressure difference, and the velocity in the vertical direction is mainly caused by gravity [25]. The adsorption term of the governing equation represents the adsorption/desorption of pollutants in porous media [37], which has been discussed in previous literature [38, 39]. Actually, several important isothermal adsorption models such as Langmuir, Freundlich, and Dubinin–Radushkevich were reported previously, which was already applied to dye adsorption/desorption processes in environmental, biological, and chemical industries [40].

Combined with the famous adsorption history model proposed by Bai et al. [2, 3] (i.e., Bai model), we can getwhere is the rate of penetration (LT−1), R is the retardation factor, kd is the adsorption equilibrium coefficient (L3M−1), kr is the desorption equilibrium coefficient (L3M−1), kNF is the equilibrium coefficient (L3M−1), is the attenuation coefficient of initial adsorption (L3M−1), is the attenuation coefficient of desorption process (L3M−1), Cl is the characteristic concentration (MM−1), and Cp is the maximum concentration before desorption (MM−1).

When  =  = 0, it is a linear nonequilibrium adsorption process (Figure 1(a)). On the contrary, the adsorption history effect can be considered as a nonlinear nonequilibrium model.

In the past, the transport of pollution is generally characterized by the first-order attachment kinetics [25]. Here, a nonlinear attachment-detachment model with hysteresis is used. The attachment/desorption rate in equation (1) can be written as [2, 3]where λ is the reaction rate (T−1), and S is the equilibrium concentration (MM−1).

The parameter λ denotes the mass conversion between moving and attached pollution under water flow and actually indicates the hysteretic effect of pollution attached onto the matrix. The determination method of some parameters was discussed in some previous references [26, 41], which is a very promising method for measuring migration parameters.

2.2. Model Characteristics and Numerical Calculation

Figure 1 shows that when the concentration in the solution increases (e.g., from point O to point B and further to point A), the adsorption rate decreases. There is a maximum adsorption concentration (Cd, max). When the concentration decreases from a certain point (e.g., point B), a desorption reaction occurs, and its magnitude is related to the maximum concentration of CP. When adsorption occurs again due to the increment of solution concentration, the readsorption process is reversed along the desorption process (i.e., CB), and the readsorption capacity is weakened compared with the initial adsorption process (i.e., OB), indicating that the model is related to the adsorption history of porous media.

It is assumed that the concentration at a certain point in the pollution space changes with time, as shown in Figure 1(b). Then, in the stage OA, the pollutant concentration increases, and the porous medium undergoes an initial adsorption reaction. During stage AB, the concentration of pollutants decreases, and the porous medium undergoes a desorption reaction. In the later stage (e.g., stage BC), when the concentration increases again, the adsorption reaction occurs again, and the change path is opposite to that of the desorption reaction in stage AB. During the numerical calculation, it is necessary to record the concentration peak (i.e., deposition history) during the concentration change process. Of course, when the concentration is higher than point A in the CD stage, it continues along the initial adsorption process [4143].

The coupling calculation is carried out using the physical field of “dilute matter transfer in porous media” and the mathematical field of “domain ordinary differential and differential algebraic equations” in COMSOL Multiphysics 5.4. Special solution settings complete the numerical model calculation. The time step algorithm uses the backward difference formula (BDF). According to existed literature [2, 3, 6, 24], the selected calculation parameters can be seen in Table 1.

3. Adsorption History Effect of Pollutant Migration under a Point Source

3.1. Model Establishment

Taking a point source case as an example, the influence of adsorption history on the pollutant migration process is analyzed. It is specified that the positive direction of water flow is the positive direction of x, and it is assumed that the model is a part of the semi-infinite soil layer. Therefore, the six surface boundaries of the cuboid model are set as the outflow boundary. The pollutants can flow freely from the boundary during the migration process. The pollution source will be set in the cuboid formation and vary according to schemes. To improve the calculation accuracy, the relative tolerance is set to 0.0001, the duration is from 0 to 200 d, and the time step is 1 d.

For convenience, the point source of the pollution source is set at the center of the plane, that is, the coordinates (25 m, 25 m). The concentration can be set to exponential attenuation, i.e.,where C is the concentration (in mol/l), and T represents the attenuation period of pollutants.

The initial value of concentration is set as 10 mol/l, and the attenuation period is assumed to be 100 d. The migration of point source pollutants in cuboid formation can be obtained by calculating the model (Figure 2). Figure 2 indicates that the diffusion trajectory of the point source is approximately a cone, which gradually expands outward from the pollution source. Here, we will explore the characteristics of the adsorption history model from the two perspectives of line and surface.

3.2. Case of a Line Segment

Take a line segment to analyze the concentration change, and the coordinates at both ends of the line segment are (0, 25, 15), (100, 25, 15). The change of concentration with time at this line segment is taken and compared with the previous linear model (Figure 3). From the concentration change of the nonlinear adsorption history model, the pollutant distribution range is 30 m–40 m at 20 d, and the concentration value is small. With the increase of time, the diffusion range of pollutants gradually expands, and the maximum concentration first increases and then decreases. Thus, in the early stage of pollutant (point source attenuation) migration, the peak concentration is higher, but the pollutant distribution range is smaller. The peak concentration will decrease when the pollutant distribution range increases. The concentration distribution of the nonlinear model is similar to that of the adsorption history model. However, the pollutants in the adsorption history model are more significant, and the gap becomes larger and more prominent with time. This shows that the adsorption capacity of the adsorption history model is weaker than that of the nonlinear nonequilibrium model.

Similarly, take another line segment in the model, and the coordinates at both ends of the line segment are (0, 25, 35), (100, 25, 35). The concentration value at this curve has decreased, which proves that the pollutant concentration at the bottom of the model is higher than that at the top under the action of gravity. From the perspective of concentration changes, the peak concentration changes less, and the concentration changes are more smoothly with time.

3.3. Changes of Pollutants at the Section

The z = 10, 30 m plane is intercepted to study the changes of concentration distribution with time. Compare the adsorption history model with the linear model, as shown in Figure 4. For the adsorption history model, at 30 d, the pollutant migrates to the cross section at z= 10 m. At 60 d, the pollutant propagation area has reached 1316.7 m2, the propagation range in the x direction is 29.5 m–73.4 m, and the maximum concentration in the center is 0.01 mol/l. With increased time, the transmission area of pollutants gradually increases. Due to the attenuation of the pollution source, the propagation range moves with the seepage direction. In the case of continuous seepage, pollutants will eventually flow out of the model. Compared with the adsorption history model, the linear model has a more vital ability to block the migration of pollutants. For example, at 35 d, the pollutants diffuse to the cross section, and at 60 d, the pollutant propagation area reaches 1173.4 m2, the propagation range in the x direction is 28.5 m–68.3 m, and the maximum concentration in the center is 0.01 mol/l. Compared with the adsorption history model, the pollutant propagation range in the linear model is more diminutive. For section x = 30 m, its fundamental law is similar to that of section x = 10 m.

3.4. Influence of Head Difference Change

The head difference mainly controls the seepage velocity for the porous geomaterials such as loose soil layer, soft rock, and fractured rock mass. Hence, the calculation scheme considered the effect of different seepage velocities (Figure 5). For example, comparing the distribution of pollutants in 60 days vertically, when the head difference is 10 m, the distribution range of pollutants is 1620.2 m2. When the head difference is 25 m, the pollutant distribution range is 2759.8 m2, and the concentration value basically remains unchanged. The head difference (seepage velocity) significantly impacts the migration of pollutants. The larger the head difference, the larger the pollutants migration and diffusion area. From the comparison result, the head difference only affects the pollutant migration in the seepage direction but has no significant impact on the horizontal diffusion. By horizontal comparison, the concentration value decreases with the increase of pollutant transmission area. It should be noted that the physical mechanisms can be attributed to the head difference and gravity.

4. Volume Pollution Source

The established geometric model is 800 m long, 400 m wide, and 400 m high, which is set as isotropic saturated porous media. The volume source pollutants are simplified as a sphere with a diameter of 20 m, and its central coordinates are 200 m, 200 m, 200 m, located in the middle of the stratigraphic model near the left.

Assuming that there is a pollutant leakage in all directions of the volume source, the boundary condition setting and meshing are consistent with the point source model, and each side of the cuboid is set as the outflow boundary condition. The specific geometric modeling is shown in Figure 6. The spherical boundary is set as the concentration boundary condition. To reflect the concentration attenuation of pollutants, the concentration value is set as

The concentration attenuation period is set as 3000 d. The diffusion trajectory of pollutants obtained by calculation is shown in Figure 6. The characteristic sections z = 100, 300 m are intercepted, and the variation of the concentration of each section with time is shown in Figure 7. The area of pollutants is greatly expanded compared with the point source model. The diffusion rate of pollutants is also significantly increased. For example, at z = 200 m, the pollutant area is 97729 m2 at 300 d, the maximum concentration value is 7.41 mol/l, and the reverse migration distance of pollutants along the x-axis is 53 m. At 800 d, the pollutant area is 234920 m2, the maximum concentration is 4.49 mol/l, and the reverse migration distance of pollutants along the x-axis is 70 m. This section is in the center of the volume source, with a larger diffusion area and higher concentration of pollutants. It is evident that although the diffusion area of this section is larger, the propagation distance becomes shorter. The closer to the center of the pollution source, the weaker the influence of dispersion on the diffusion of pollutants in the main seepage direction. Under the volume source model, pollutants can propagate in the reverse seepage direction, causing pollution upstream. Therefore, in the actual process of pollutant prevention and control, we should also pay attention to the pollution of the groundwater upstream.

The characteristic sections y = 100, 300 m are intercepted, and the variation of the concentration of each section with time is shown in Figure 8. The maximum concentration value at this section shows a trend of increasing first and then decreasing. The maximum concentration value will also decrease with the increase of pollutant area in the volume source model. For example, at y = 100 m, the pollutant area is 50605 m2 at 300 d, and the maximum concentration value is 0.2 5 mol/l. At 600 d, the pollutant area is 161300 m2, and the maximum concentration is 0.31 mol/l. At 800 d, the pollutant area is 211270 m2, and the maximum concentration is 0.29 mol/l.

Comparing the two sections, it can be found that the changes of pollutants in the symmetrical sections on both sides of the center of the pollution source are the same. The diffusion of pollutants in the direction perpendicular to the main seepage is mainly affected by dispersion and is symmetrically distributed. The effect of countercurrent propagation of pollutants can be seen from the migration of pollutants in the section at y = 300 m.

The characteristic sections x = 300 m, 400 m are intercepted, and the variation of the concentration of each section with time is shown in Figure 9. In the direction perpendicular to the main seepage flow, the diffusion of pollutants is mainly affected by dispersion. The range and concentration of pollutants increase first and then decrease. Dispersion leads to the increase of pollutant area. However, because the seepage effect is stronger than the dispersion effect, the pollutant range and concentration value on the x section will decrease with the attenuation of the concentration of the pollution source. For example, at x = 400 m, the pollutant area is 7653 m2 at 150 d, and the maximum concentration value is 0.05 mol/l. At 300 d, the pollutant area is 72047 m2, and the maximum concentration is 1.18 mol/l. At 800 d, the pollutant area is 77171 m2, and the maximum concentration is 1.28 mol/l. By comparison, the longer the section is from the pollution source, the later the pollutants appear, the faster the diffusion speed of pollutants, and the lower the concentration of pollutants. In the volume source model, the propagation of pollutants is still conical. The closer the pollution source, the smaller the propagation range, and the higher the concentration value. The farther away from the pollution source, the larger the propagation range, and the smaller the concentration value.

The three-dimensional nonlinear nonequilibrium adsorption model that takes into account the adsorption history can be used to reveal the diffusion law of pollutants under regular three-dimensional boundary conditions. That is, this model considers the probable increase or decrease of seepage velocity as well as pollutant concentration in practical engineering. The transport laws of pollution sources in three-dimensional spatial are also of interest in real pollution problems such as changes in the diffusion rate and concentration field of pollutants. The problem of pollutant diffusion under complex boundary conditions (e.g., the variation of seepage direction and gravity effect) should be discussed in the future.

5. Conclusions

A governing equation is obtained by embedding a nonlinear nonequilibrium adsorption model (i.e., Bai model) considering the adsorption history in the three-dimensional situation, which can be easily applied to solve practical engineering problems. The concentration peak of the nonlinear nonequilibrium adsorption history model is higher than that of the linear model, and the peak concentration appears slightly in advance. The pollution source tends to migrate horizontally along the water flow direction, and the diffusion capacity in the vertical direction is small. Due to the influence of gravity, the concentration changes at the symmetrical position of the vertical distance are different, and the bottom concentration is higher than the upper concentration.

In the early stage of pollutant migration, the peak concentration is high, but the distribution range of pollutants is small. When the distribution range of pollutants increases, the peak concentration will decrease. Compared with the adsorption history model, the linear model has a more vital ability to block the migration of pollutants. The longer the time, the more significant the gap between the two models. Close to the center of the pollution source, the gap between the two models increases. By changing the attenuation period of the pollution source, the longer the attenuation period is, the wider the diffusion range is, but the longer the time it takes for the pollutants to migrate out of the model entirely, and different sections have different regular distributions.

For the volume source model, the pollutant concentration is higher than that of the point source model, and the diffusion range of each section is more comprehensive. The closer to the center of the pollution source, the weaker the influence of dispersion on the diffusion of pollutants in the main seepage direction, and the pollutants will spread in the countercurrent direction, causing pollution upstream. The diffusion of pollutants in the direction perpendicular to the main seepage is mainly affected by dispersion and is symmetrically distributed.

Notations

C:Pollutant concentration in porous media
Cd:Adsorption concentration
Cl:Characteristic concentration
Cp:Maximum concentration before desorption
D:Hydrodynamic dispersion coefficient
kd:Adsorption equilibrium coefficient
kr:Desorption equilibrium coefficient
kNF:Equilibrium coefficient
L:Length of soil column
r:Soil column radius
R:Retardation factor
t:Time
T:Pollutant attenuation period
u:Average velocity of cross section
:Penetration rate
β1:Attenuation coefficient of initial adsorption
β2:Attenuation coefficient of desorption process
∇:Vector gradient operator
ρb:Volumetric dry density of porous media
(Dd)ij:Mechanical dispersion coefficient tensor
Dτij:Molecular diffusion coefficient tensor
τij:Curvature tensor
αL:Longitudinal dispersion of porous media
αT:Transverse dispersion of porous media
αx:Tortuosity degree
θ:Porous media porosity.

Data Availability

The data supporting the current study are available from the corresponding author upon request.

Conflicts of Interest

WC, SZ, and WL are employed by the Beijing Guodaotong Highway Design & Research institute Co., Ltd. GZ, PZ, and SY are employed by the Beijing Uni.-Construction Group Co., Ltd. The authors declare that they have no conflicts of interest.

Authors’ Contributions

The main contribution of ZM in this paper is methodology, and the other authors (WC, GZ, PZ, SZ, WL, and SY) contributed to investigation and analysis.

Acknowledgments

This research was funded by the Beijing Natural Science Foundation (8222023).