Abstract
Adsorption reaction in unsaturated porous media is of great importance for soil and groundwater remediation. In this study, the influence of the Peclet number (Pe) and water saturation on adsorption behavior at liquid-liquid interfaces was quantitatively investigated. The pore-scale reactive transport in unsaturated porous media was directly simulated. The Navier-Stokes equations, the surface transfer and adsorption reaction equations, and the advection-diffusion equation (ADE) were coupled to obtain the flow and concentration fields. The results showed that water saturation had a significant influence on the complexity of the flow field. A nonmonotonic relationship was found between water saturation and the uniformity of the flow field. Peclet number had little influence on the maximum adsorption. On the other hand, the adsorption time showed a nearly linear relationship with the Peclet number and increased with increasing Peclet number. Additionally, a nonlinear relationship was found between water saturation and the maximum adsorption. As water saturation increased, the maximum adsorption tended to increase to a peak and then decrease. The peak of the maximum adsorption occurred at , , and the shortest adsorption time was observed at . However, the difference in adsorption times for saturations of 0.458 and 0.698 was not significant and was only about 20 PV difference.
1. Introduction
Soil is an important resource for human survival and a site for various reactions in groundwater. Surface soils are mostly unsaturated porous media subject to natural influences. In recent years, it has become evident that uncontrolled use and immature environmental regulations have led to a variety of ecological and environmental problems. Pollutants such as organic matter, heavy metals, pesticides, and fertilizers led to long-term soil pollution. In addition, soil pollution caused significant pollution of groundwater due to adsorption reactions at the soil surface and hydrodynamic influences on the groundwater system [1–5]. Therefore, understanding and applying reactive solute transport in unsaturated porous media and their internal adsorption reactions are crucial for soil and groundwater remediation.
Solute transport in porous media has been widely studied [6–9]. As in previous studies, solute transport research focused on the macroscopic scale. Transport processes in porous media were studied by constructing unit bodies and using averaging parameters to make macroscopic generalizations about their relevant properties. Wu et al. [10] determined the topology of spatial structure between particles of agar bead-filled beds under water-saturated conditions using magnetic resonance imaging (MRI). Wang et al. [11] studied the influence of soil anisotropy on solute transport. Pore-scale studies have made significant progress in recent years due to advances in computer technology and the refinement of numerical simulation techniques. Current methods for determining microscopic pore structure in porous media include physical imaging techniques and numerical methods. Zhou and Xiao [12] proposed a method to reconstruct 3D multiphase rock models using a combination of computerized tomography (CT) and multiple features to study simulated fracture formation and obtained better results compared to experimental data. Scanning electron microscope (SEM) is also widely used; Liu et al. [13] and Zhao et al. [14] produced gray-scale images to refine the microstructure and quantify the characteristic parameters of the pore structure. Luo et al. [15] combined process-driven three-dimensional (3D) digital rock modeling with fractal analysis to create a digital rock model of simulated sandstone. Zhao et al. [16] determined the porosity of residual granite beds under external loading by CT scanning and derived the microstructural evolution process. However, due to the economic and technical limitations of the physical imaging method, the numerical simulation method has been used more frequently.
Although considerable efforts have been made in the field of solute transport in recent years, many studies have addressed only partial aspects that have been covered, especially in reactive solute transport problems [17–21]. Unsaturated brucite [Mg(OH)2] bearing column experiments have been used to evaluate the influence of the change in reactive mineral surface area, water content, precipitation of hydrous Mg carbonate, and gas distribution on the reaction of brucite [Mg(OH)2] with CO2 gas [22, 23]. The results confirmed that conventional geometric surface updating models do not adequately represent the observed reaction process during the carbonation of brucite. Raoof et al. [24] presented a pore network modeling tool capable of simulating fluid flow and reactive and adsorptive transport of multiple components under saturated and differentially saturated conditions. Jiménez-Martínez et al. [25] pointed out mechanisms by which the mixing of an invading fluid with the resident fluid is considerably enhanced as saturation decreases, and the effective reactivity is much greater than under saturated conditions. Mixing-driven reactions in porous media are ubiquitous, but predicting where and how fast the reactions occur has been a tremendous challenge due to the complex and nonuniform nature of porous media flows [26]. Dou et al. [27] developed a pore-level numerical simulation model to investigate the saturation dependence of mass transfer in unsaturated porous residual media and the influence of saturation topography on mass transfer. They discovered a nonmonotonic relationship between dispersivity and residual saturation. However, the reactive properties at the liquid-liquid interface are neglected. Especially for adsorption, Yue et al. [28] studied the adsorption and desorption of the herbicide atrazine on soils, which is dominated by surface adsorption at low concentrations and influenced by adsorption sites on the adsorption surface. According to Lyu et al. [29], the spherical structure and ultrahigh surface ratio of the extremely efficient catalyst were studied to make it more efficient and stable. Thus, interfacial adsorption is the key to soil remediation, and the research combining unsaturated porous media with interfacial adsorption will receive more attention.
In this study, the complex pore structure was reconstructed to simulate the reactive solute transport in unsaturated porous media considering adsorption reactions at the interfaces at the pore scale. An immiscible two-phase transport (ITPT) model was used to construct the geometric model of the unsaturated porous media with different water saturations. The Navier-Stokes equation, the surface transfer and adsorption reaction equations, and the advection-diffusion equation (ADE) were coupled and applied to the geometric model using COMSOL Multiphysics. We investigated the influence of water saturation on the flow field and quantitatively evaluated the uniformity of the flow velocity distribution. The influence of the Peclet number and water saturation on adsorption behavior at liquid-liquid interfaces was analyzed.
2. Materials and Methods
2.1. Unsaturated Porous Media Generation
A series of regular circles with different radii were simulated as solid particles in the porous media. The microscopic pore structure of saturated porous media particles was generated using the method described by Dou et al., i.e., a single porous medium with target porosity is generated in a specific region [30]. It should be noted that the size of the solid particles must be appropriate. To this end, the radius of the particle size should not only be within the range where the mean particle size and the standard deviation of the particle size distribution were usually controlled but also followed a truncated log-normal distribution. In this study, the radius ranged from 0.3 mm to 0.9 mm and was randomly distributed in a size range of mm. The coefficient of variation (COV) was equal to 0.25, and the porosity of the generated porous media was set to 0.42. Due to the random distribution of particles and heterogeneous radius sizes, the internal spatial structure of the porous media was complex. Therefore, the porous media were considered heterogeneous in this study.
Our numerical simulations were based on four porous medium topologies with different water saturations (, , , and ). The topology was determined by an ITPT model consisting of immiscible two-phase flows at the pore scale of Dou et al. [27]. Transport at the immiscible boundary was tracked using the phase field method (PFM). Trichloroethylene (TCE) was assumed to be the nonwetting phase fluid, and water was injected as the wetting phase fluid. The porous medium was first saturated with TCE and then injected with water at three rates to displace the TCE. In the end, we obtained models with different water saturations. As shown in Figure 1, the solid grains were shown in white. The liquid in the nonwetting phase was shown in gray, the blue areas showed the flow paths, and the interface between gray and blue is the liquid-liquid interface, where the adsorption and desorption reaction occurred.
(a)
(b)
(c)
(d)
2.2. Fluid Flow and Adsorptive Transport Modeling
To solve the mass transfer problem in unsaturated porous media, it is necessary to couple the groundwater flow field with the mass transfer concentration field. The mathematical fields of the reactive adsorption process must also be coupled with the two aforementioned fields to account for the adsorption reaction of solutes with the internal residual non-wetting phase.
To simulate water flow in unsaturated porous media, we used the Navier-Stokes equation as the control equation for water flow. where is the density of water, is the acceleration of gravity, and represents the water head. is the dynamic viscosity coefficient of water, and is the velocity.
The adsorption/desorption reaction at the interface of unsaturated porous media can be described as where the forward reaction rate is , the reverse reaction rate is , is the molecular material in the solute, and is the molecular material adsorbed on the surface.
In this study, adsorption in a single molecular layer was used for the reaction interface, and molecular forces between adsorbed molecules were neglected to account for adsorption at the interface during solute transport. According to the Langmuir hypothesis of adsorption in a single molecular layer, the adsorption capacity at all adsorbable interfaces was the same because the interface was homogeneous throughout. Therefore, the adsorption and desorption processes were in a dynamic equilibrium state during solute transport.
Assuming that is the percentage of the surface covered by adsorbed molecules, the percentage of the free zone is 1-. According to Langmuir’s theory of monolayer adsorption, the desorption rate of molecules, , was proportional to and equal to , where was constant at a constant temperature. The adsorption rate of molecules was proportional to the percentage of free zone 1- and the rate at which the molecule collided with the surface. The collision rate was in turn proportional to the partial pressure of the substance . The adsorption rate (mol/(m2·s)) was therefore expressed as , which was equal to (1-). The was constant at a constant temperature.
To establish the equation by the concentrations, the following substitutions are made. where we set the in-domain concentration of to and the surface concentration to . is the total concentration of the surface that can be covered by the adsorbed molecules, is the gas constant, and is the temperature.
Eq. (3) and Eq. (4) were used to derive the adsorption and desorption rates. where the rate constant for the desorption reaction and the rate constant for the adsorption reaction .
This leads to the equations for surface transfer and adsorption.
The nonuniform velocity field of groundwater flow in unsaturated porous media and the surface reaction rate can be derived from Eq. (5) and Eq. (6), which are inserted into ADE to obtain the following equation. where the surface diffusion coefficient and the reactant diffusivity .
2.3. Numerical Modeling Strategy and Setup
COMSOL Multiphysics finite element software (COMSOL Inc., Burlington, MA, USA) was used to simulate water flow, solute transport, and adsorption at the interface using. The model was set up with a constant pressure outlet boundary on the right and a constant pressure inlet boundary on the left. The particle boundary and the upper and lower boundaries of the geometric model were set as slip-free boundaries. The slip-free boundary condition is a central part of the boundary conditions used in fluid dynamics. It assumes that the fluid in direct contact with a solid surface adheres to that surface and has the same velocity as the surface. In other words, the velocity of the fluid directly at the boundary is zero. The inclusion of the slip-free boundary condition greatly simplifies the mathematical representation of the interactions between different materials and their corresponding surfaces. Different pressure differences were set for the different models to characterize the distribution of water pressure, and the water flowed under the initial pressure difference.
As a fundamental requirement for solute transport, it was important to determine the principles controlling the fluctuation of the flow field in porous media. The Navier-Stokes equation was coupled with the advection-diffusion equation to regulate the flow field. Water flow is controlled by applying the pressure difference (5 Pa) between the inlet and outlet. To simulate adsorption behavior at liquid-liquid interfaces in unsaturated porous media in the pore structure, the solute with the specific concentration () was injected along the inlet boundary. At the same time, the corresponding reaction rate constants were determined to represent the adsorption reaction at the interface. The specific parameters are given in Table 1.
3. Results and Discussion
3.1. Impact of Water Saturation on the Flow Fields
Figure 2 shows the distribution of the flow field at different water saturations. It can be seen that the distribution of the flow field was complex, and there was a clear phenomenon of preferential flow. This preferential flow generally occurred between two groups of solid particles that formed the dominant channel and were almost parallel to each other. Compared with the flow velocity in the surrounding area, the dominant channel had a higher velocity, and the flow direction was approximately parallel to the water flow. At the same water saturation, the distribution of dominant channels was complex and random due to the influence of pore structure and the distribution of water. Multiple dominant channels may be connected at one location to form a larger dominant channel or separated at another location to form multiple smaller dominant channels. This random separation creates the heterogeneity of the flow field in porous media. In Figure 2, it can be seen that water saturation has a significant influence on the occurrence of connectivity and separation of the dominant flow. The strongest connectivity and separation occurred in the least saturated porous media.
To investigate the phenomenon of preferential flow, the location of the peak velocity of the flow field was selected for magnification. In Figures 2(a)–2(d), it can be seen that there was a clear phenomenon in the velocity distribution of the flow field in the pore structure: the peak velocity was located on the line connecting the centers of the circles of two groups of nearly parallel solid particles, and the velocity decreased in all directions with the peak velocity as the center. In Figure 2(e), the water flowed in the dominant channel. As the distance between the two groups of solid particles decreased, the flow velocity peaked. However, due to the random distribution of the pores, the particles differed greatly in size, and the dominant channel was divided into two parts, forming two tributaries that continued the flow.
In this study, the profiles of flow velocity were determined to quantitatively evaluate the uniformity of flow velocity. The coefficient of variation [31] and uniformity index [32] were used for the comprehensive analysis.
The coefficient of variation (CV) is dimensionless. It reflects the degree of dispersion of the flow velocity compared to the mean flow velocity and exhibits a negative relationship with the flow field uniformity. where is the standard deviation. [m∙s-1] is the velocity value at , and [m∙s-1] is the mean velocity. represents the total number of sample points obtained.
As the crucial evaluation parameter of the flow field uniformity, the uniformity index () shows the flow velocity distribution across the entire section. The range of values is 0 to 1.
According to Table 2, the CV were all close to 1.000 and did not significantly change with water saturation. Meanwhile, the values at four water saturation were all more than 0.5000. Poor uniformity of the flow field was discovered after a thorough investigation of the two indicators. We plotted the relationship between uniformity and water saturation, as seen in Figure 3, in order to clearly analyze the relationship. Since the CV was inversely related to uniformity, we used the inverse of CV () and the value to effectively indicate the magnitude of the flow field uniformity. The and value showed the same trend, increased, and then decreased as water saturation decreased. In other words, the flow uniformity showed a nonmonotonic relationship as water saturation decreased.
3.2. Impact of the Peclet Number on Adsorption Behavior at Liquid-Liquid Interfaces
Four different Peclet numbers (, , , and ) were set at the same water saturation to investigate the influence of different Peclet numbers on reactive transport in unsaturated porous media.
It is evident from Figure 4 that the Peclet number has a significant influence on reactive transport in unsaturated porous media. One or more preferential fronts with high transport rates occurred in the porous media. These fronts were distributed at higher concentrations when the solutes were still at low concentrations in other areas. During solute transport, the solutes first occupied the larger pores and then diffused toward the smaller pores, with several larger pores together forming the preferred channel for solute transport. Large pore channels had a better connection and higher flow velocity than small pore channels, which allowed large pores to load up with solutes at higher head pressures. However, the velocity in small pores was significantly lower than in large pores that form the preferred channels because small pore channels were mainly filled by diffusion due to weak advection. The direction of movement of the fronts formed by solute transport was affected by the distribution of the larger pore spaces, and the fronts moved along the channels formed mainly by the larger pore spaces.
Considering four groups of models of unsaturated porous media with different Peclet numbers at the same water saturation, we find the following: at 0.5 PV, all the solute was injected, and the concentration near the left inlet was in the higher range, but some of the ranges lagged behind the concentration changes in the surrounding area. Then, the injection of water is continued. At 2 PV, the solutes had essentially flowed out of the right outlet, and the entire region was in a state of zero concentration. However, some areas still had a high residual solute concentration and lagged behind the change in the concentration center. Because of this phenomenon, these areas were referred to as “stationary zones.” Stationary zones were usually made up of smaller pores or dead-end pores with a slow flow velocity and weak advection that relied on diffusion filling of solute molecules, so these areas lagged behind changes in the surrounding area as if they were stopped. In contrast, there were also “mobile zones” (see Figure 5).
In the following simulation period, the solutes in the stationary zone were largely dispersed by the continuous input of water. However, due to adsorption, a large number of solutes were adsorbed on the surface of the nonwetting liquid phase in the initial phase. The desorption reaction at the interface gradually became visible due to the repulsion of water, as shown by the results in Figures 4(i)–4(p). It can be observed that the desorption is strongly affected by the Peclet number. For the models with unsaturated porous media geometry and smaller Peclet number, significantly more solute was desorbed from the surface of the nonwetting liquid, while less was desorbed for the larger Peclet number.
In this study, the outlet boundary was used as the location for calculating the concentration of the breakthrough curves. The simulation time was set to 5 PV in the COMSOL Multiphysics software (0-0.5 PV for solute injection and 0.5-5 PV for water injection), and this was used as the horizontal coordinate. The expression for the dimensionless concentration was as follows. where the solute concentration is [mol∙L-1] and the initial concentration of solute is [mol∙L-1]. The length of the outlet boundary is . The average velocity of water flow is .
The breakthrough curves for different Peclet numbers at four different water saturations are shown in Figure 6. All BTCs at the same water saturation showed typical characteristics of anomalous transport behavior (e.g., long tails of BTCs).
(a)
(b)
(c)
(d)
As can be seen in Figure 6, breakthrough curves follow a similar trend at different Peclet numbers. As the Peclet number increased, the arrival time was earlier. This is because advection was the primary mechanism of reactive solute transport. And at larger Peclet numbers, the flow was faster. It should be noted that the breakthrough curves exhibited left-right asymmetry and a significant tailing phenomenon. For the tailing phenomenon, the saturated and unsaturated curves were diametrically opposed. At unsaturated curves, the trailing phenomenon became more pronounced as the Peclet number decreased. In this study, solute chemistry and adsorption in saturated conditions were ignored, so the above phenomena were caused by the microscopic morphology of the pore structure and the inhomogeneous particles within the porous media. In unsaturated conditions, it was the result of a combination of adsorption reactions at the surface of the nonwetting liquid phase and the influencing variables in saturated media.
In order to quantitatively describe the influence of Peclet number on adsorption behavior at liquid-liquid interfaces, 500 PV of contaminate was injected to reach adsorption saturation. The influence of adsorption within 10 PV was exhibited to clearly highlight the trend, as shown in Figure 7. The adsorbed amount at different Peclet numbers varied regularly with time and gradually increased to equilibrium as time increased. At early times, the adsorbed amount increased in a very similar way for all cases, and the adsorption rate decreased with increasing Peclet number. As time increased, the adsorption rate decreased, and the adsorption reached saturation. Table 3 shows the maximum adsorption and adsorption time at different Peclet numbers and water saturations. Further, we also plotted the maximum adsorption at different Peclet numbers as shown in Figure 8. It indicated that the Peclet number had little influence on the maximum adsorption. The adsorption time at saturation, on the other hand, increased with increasing Peclet number and showed a nearly linear relationship with the Peclet number. This is due to the fact that at higher velocities, the ability of solutes to penetrate the interface and experience adsorption is less than the transport capacity in the direction of water flow. Therefore, the adsorption rate is lower at higher flow velocities, and thus, more adsorption time is required.
3.3. Impact of Water Saturation on Adsorption Behavior at Liquid-Liquid Interfaces
The relevant parameters were determined as described above. Numerical simulations with constant PV were performed in four groups of geometric models with different water saturations (, , , and ).
At the same Peclet number, as shown in Figure 9, it became clear that the same phenomenon of preferential flow occurred in reactive transport in the four different water saturation groups. In the less saturated geometric model, more fluids were in the nonwetting phase, resulting in a greater number of fronts with significant separation, compared to the more saturated geometric model in which there were fewer fronts without significant separation. Similarly, there were stationary and mobile zones as the water spread, resulting in strong stagnation.
Figure 10 clearly depicts the influence of the nonwetting phase fluid (TCE) and of the heterogeneity; it adds to the pore structure on reactive transport. The solute shows larger spreading, with earlier arrival times and a more pronounced tailing as water saturation decreases. This is due to the fact that with decreasing water saturation, more fronts existed in the unsaturated porous media, and reactive solute transport along the dominant channel became more obvious. Although the internal heterogeneity resulted in strong lateral dispersion, the change in transport channels dominated in comparison, resulting in a decrease in transport time.
(a)
(b)
(c)
(d)
Similarly, we investigated the direct influence of saturation on adsorption behavior at liquid-liquid interfaces. The adsorbed amount at different water saturations is shown in Figure 11. In the early stage, the adsorbed amount increased in a very similar way, and the adsorption rate increased and then decreased with increasing water saturation. Figure 12(a) also shows a nonlinear relationship between the maximum adsorption and water saturation. The maximum adsorption tended to increase and then decrease as saturation increased. We explained this trend by finding a link with the area in porous media available for adsorption. It can be seen that the adsorption area tends to increase and then decrease with saturation as shown in Figure 12(b). This was consistent with the change in the maximum adsorption with saturation. The maximum adsorption at exceeded that of (the largest adsorption area), and this was thought to be due to the combined influence of the stationary zone and the preferential front. According to Table 3, the shortest adsorption times occurred in porous media with . Due to the microscopic morphology of the pore structure and the inhomogeneous particles within the porous media, and were substantially longer than those for . However, the difference in adsorption times for saturations of 0.458 and 0.698 was not significant and was only about 20 PV difference.
(a)
(b)
4. Conclusions
(1)The complexity of the flow field was influenced by water saturation, and the uniformity of the flow showed a nonmonotonic relationship with water saturation. Quantitative evaluation of the CV and values showed that the uniformity of the flow field first increased and then decreased with decreasing water saturation(2)The Peclet number had little influence on the maximum adsorption. On the other hand, the adsorption time showed a nearly linear relationship with the Peclet number and increased with increasing Peclet number(3)A nonlinear relationship was found between water saturation and the maximum adsorption. With increasing water saturation, the maximum adsorption tended to increase to a peak and then decrease. The peak of the maximum adsorption occurred at , , and the shortest adsorption time was observed at . However, the difference in adsorption times for saturations of 0.458 and 0.698 was not significant and was only about 20 PV difference
Data Availability
The underlying data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.
Acknowledgments
The study is financially supported by the National Key Research and Development Program of China (No. 2019YFC1804300), the National Natural Science Foundation of China (Grant No. 42272278), and the China Scholarship Council (CSC202106710083).