Migration characteristics of oil pollutants into soil pores influenced by surface wettability

ABSTRACT Pore-scale modeling plays an indispensable role in unveiling the migration mechanism of oil pollutants into the water-saturated soil. Navier Stokes equations and conservation of mass, coupled with the phase field method, are employed to investigate the effect of soil wettability on oil pollutants migration into the soil porous media. The results show that the oil-water interface in the migration process does not advance evenly in three wetting states. When the model system is in oil-wetting state, oil pollutants obviously migrate at a higher speed. Through monitoring the evolution of oil front location and analyzing the velocity distribution vertical to the flow direction at the central axis of the model, the channel corresponding to the peak velocity is the dominant flow channel for the oil pollutants. This study will have the potential ability in obtaining the pore-scale migration regulation of oil pollutions as well as more accurately predicting the pollution range.


Introduction
The leaked oil from the buried pipelines can adversely affect soil ecosystems, groundwater quality and surroundings near the leakage sites, leading to serious environment pollution and health problems of local mankind and other living beings and this problem has become conspicuous day by day [1][2][3]. The exhaustive literature study reveals that most oil leakage around the world occurs onshore due to pipeline distribution features. The leaked oil may act as contaminant or may carry other pollutants to enter into the subsurface soil system, displacing the original fluid in soil pores and eventually forming the immiscible phase fluid flow. In some areas with high rainfall or abundant groundwater, the soil water content may be extremely high or reach the saturation state [4]. Therefore, considering the serious consequences of oil leakage, a comprehensive and systematic research of the leaked oil migration characteristics in water saturated soil is required to weigh up the disadvantages to environment.
Researchers have conducted a large number of studies on the migration characteristic of the leaked oil based on the experimental and simulation methods. Pu et al [5] experimentally simulated the porous soil with glass spheres to conduct the resistance coefficient experiment and demonstrated the migration ability of fluid in soil porous medium. Liu and Zhang [6] established a sandbox model and obtained the migration law and distribution characteristic of the non-aqueous phase liquid in porous media. But due to long experimental period, high cost and scale effects, these disadvantages make experiment method rather limited. Accordingly, simulations have been widely applied and made great progress in describing oil migration regularity in soil. During the long-distance migration of oil leaked from the pipeline, the migration laws under different oil densities and leakage hole sizes were determined by Zhu et al. [7] through computational fluid dynamics method. He et al. [8] considered that oil went through such a process: pipe flow, leakage hole outflow and seepage diffusion. Thus a coupled model including the pipe flow and seepage-diffusion processes was established to study the effects of leakage flow rate, pipe diameter and pipe depth on the leakage range and volume. Although these macroscopic investigations have revealed a series of information such as the leakage volume, pollution range and global motion trend of the leaked oil, this approach cannot accurately obtain the micro-level behavior, mechanism and dynamics of pollutants migration process which are only available from the microcosmic modeling. Lattice Boltzmann Method (LBM) was employed to simulate the mesoscopic scale migration of the leaked oil in the reconstructed soil porous structure based on Quartet Structure Generation Set algorithm (QSGS). The influence of different porosities on the velocity field, pressure field and streamline distribution of oil pollutants in the computational domain was found [5]. The leaked oil migrated in the subsurface soil, forming the immiscible CONTACT Yu Pu puyunepu@163.com; Di Wang wangdinepu@163.com This article has been republished with minor changes. These changes do not impact the academic content of the article. two-phase flow in pores. A simplified model with linear porosity was used to simulate soil porous media. Based on the oil trajectory in soil porous media, oil leakage process and the existence state of leaked oil was analyzed [9]. Jiang et al. [10,11] used volume of fluid (VOF) method [12] to study the migration characteristics of leaked oil in turbulent state. These studies find the filling state of oil in soil model and the pollution spread range, but do not reflect the influence of pore structure characteristics of the real soil itself on oil pollutants migration. Considering complex pore geometries and topological changes in soil [13], oil through the porethroat will form the sophisticated moving interfaces and the interface evolution in space and time [14]. Phase field method [15,16] possesses higher potential and becomes a very popular approach to investigate the two phase motion at pore scale [17,18]. Through combining Cahn-Hilliard phase field method, which is a relatively simple fourth-order partial differential equation (PDE), with Navier Stokes equations, the water-oil two phases flow in a complexly heterogeneous porous media was modeled to finally find the flow regimes controlled by viscous and capillary forces [19,20]. Thus phase field method can realize to obtain the oil pollutants migration mechanism at pore scale. Because the leaked oil flows in soil pores for a long time, soil properties will inevitably affect oil migration process. Wettability, as an important hydraulic parameter of soil, is closely related to the formation of preferential flow, the flow state and occurrence state of fluid at pore scale [21,22]. However, very few literatures specifically focus on the effect of soil wettability on leaked oil migration. Therefore, it is necessary to further investigate the effect of the soil wettability on the oil pollutants migration behavior, retention and distribution at pore scale.
In this study, pore-scale model of oil pollutants migration is established and validated by comparing with the previous experimental results. The next step consists then of investigating the effect of three typical wetting states on the oil pollutants migration behavior.
Finally, the observed pore-scale migration characteristics and underlying mechanisms of some typical phenomena are discussed, which is conductive to accurately predict the pollution level and subsequently control the oil pollutants in soil.

Physical model
The leaked oil flows along the leakage hole under the pressure gradient and subsequently enters into the soil. In the research, the focus is put on the oil pollutants migration in soil porous media. Figure 1 shows a schematic of the realistic porous domain in the present study. The dimension of the pore-scale physical model is 350 μm × 250 μm. Meanwhile, the physical properties of the soil porous media are 17.5 D for permeability and 38.5% for porosity. Considering the abundant lakes and swamps around the oil pipeline in study area, the soil is rich in water content, so air can be neglected to simplify the model [23]. Accordingly, the oil pollutants migration can be assumed to be an oil-water two-phase flow. The viscosity and density of oil are 8.9 mPa·s and 867 kg/m 3 , respectively. The interfacial tension of the two phases is 2.8 × 10 -2 N/ m. Before running each simulation, water is initially filled in the soil model, while oil is uniformly distributed in the inlet buffer zone on the left side of the model. Then the oil continuously migrates into the soil pore under the inlet and outlet pressure of 2.6242 × 10 5 Pa and 2.1242 × 10 5 Pa, respectively. The wettability condition at the soil particle walls is set to the wetted walls with a specified contact angle according to the following research requirement.

Mathematical model
Due to the constant inflow of oil, the interface location of two phases flow will continuously change during the pollutants migration process. Therefore, a moving interface should be taken into consideration to conserve mass. The model of the immiscible two-phase flow during the pollutants migration process is established by the Cahn-Hilliard phase-field method coupled with Navier-Stokes and continuity equations. Assuming that these two kinds of fluids are incompressible, the correlative equations in the micro channel are given as follows [14]: where p is pressure, I is the identity matrix, u denotes the fluid velocity, g is the gravitational acceleration, μ equals the dynamic viscosity. ρ denotes the density.
The surface tension F st in the momentum equation plays an important role in the phase field interface tracking simulation, which can be obtained as follows: The chemical potential G is approximated as follows: In the interfacial dynamics simulation, the migration interface of finite thickness is defined as the region where the dimensionless phase field variable (Φ) can change rapidly from −1 to 1 across the interface. In this definition, Φ = ±1 represents the bulks of the two phases. Correspondingly, −1 < Φ < 1 represents the interface. The Cahn-Hilliard equation can be decomposed into two second-order partial differential equations: where γ is the mobility, and λ is the mixing energy density. ε is the interface thickness parameter, which is typically defined as a half of the largest value of the mesh size when the interface possibly passes through the region. ψ denotes the phase field auxiliary variable. Based on the mixing energy density and the interface thickness, as ε→0, the surface tension coefficient can be defined by the following equation [14]: All the fluid physical properties between the two phases can be calculated by means of the relative volume fractions and dimensionless phase field variable. The viscosity and density are then as follow: where subscripts 1 and 2 denote two different phases, respectively. Symmetry boundary conditions are established on the lateral sides of the model. As previously discussed, the inlet and outlet boundary conditions are the real pressure conditions from the three dimensional simulation. On the surfaces of the grains, the wetted walls with a specified contact angle are employed as the boundary conditions as follows: in which n is the unit normal to the wall and θ is the surface contact angle.

Model validation
In order to validate the accuracy of this model, the phase field interface tracking method needs to be confirmed further. In view of the fact that there are few pore scale oil pollutants migration study, the new developed model is used to simulate the process of oil displacement and compared with the previous experimental results by Sun et al. [24]. For this purpose, a CAD micro-fluidic chip model identical to that employed in Sun's experiment is constructed. In this way, our numerical simulations for water-oil displacement processes in micro-fluidic chip model are qualitatively compared with Sun's available experimental observations. As shown in Figure 2, our numerical approach almost synchronously reenacts the flow trends observed during the experiments for the water/oil displacement processes in micro-fluidic chip model. It can be concluded that the excellent agreement of our model with the previous experimental results are proved and the numerical approach can be used to investigate oil leakage migration characteristics in soil.

Results and discussions
In this section, several simulations are performed to investigate the oil pollutants migration behavior and pollution level in the real soil porous media with three wetting states. Three different contact angles frequently found in the soil porous structure are considered such as 45°, 90° and 135°, which correspond to water wet, intermediate wet and oil wet [25].

Migration process evolution of oil pollutants
The evolution of oil migration process and oil-water interface position at different contact angles are shown in Figure 3. At initial phase, the flow rate and flow trend of two-phase interface are nearly same. But due to the complex pore throat structure of real soil, the oil-water interface does not advance evenly. By 7 × 10 -4 s, the front interface position in oil wet state is farthest from the entrance. After that, it is the first to break through the model outlet with the shortest time of 1.227 × 10 -3 s. It can be concluded when the system is in oil-wetting state, the oil obviously migrates at a higher speed, resulting in the extension of the oil pollutants range. From the analyses of forces acting on the pore-throat, the direction of capillary force in oil wet is identical to that of flow velocity, behaving as the driving force. Correspondingly the direction of capillary force in water wet is opposite to that of the seepage velocity, acting as the resistance force, which will hinder the migration. In the oil wet state, the generated shear force is insufficient to peel off the water droplet, finally retaining in the back slope of the soil particle and forming the some cluster stagnant water. Figure 3 illustrates that the channel to first breakthrough the outlet of the model is same in three wetting states. This indicates there may exist a predominant migration pathway for the oil pollutants. Thus, the oil front positions as function of time are monitored and plotted at contact angle of 45°, 90° and 135°, as shown in Figure 4. In the early stage of migration, the slope of oil front position curves along the X-direction is significantly higher than that in the middle and late stage of migration for three wetting conditions, which explains that the oil pollutants migration ability is abate in the late stage. At the  beginning, the oil front positions have not significant difference. After 3 × 10 -4 s, the values of distance from the entrance end of the model along the X direction at the contact angle of ~ 135° are always higher than that of the other two, which indicates that the migration rate of oil pollutants along the flow direction is faster in oil wetting state. In addition, through observing these three breakthrough timelines (the lines that are perpendicular to the flow direction), it can also be found that the oil front in oil-wetting state first break through, followed by intermediate wetting, and finally water wetting. This is consistent with the results obtained from curve slope analysis. However, under three contact angles, the y-coordinate values mainly fluctuate between 148 and 177 (as shown by the blue dotted line in Figure 4), which further explains that the corresponding channel in this range is the main migration channel. The velocity at the central axis of the model (i.e. x = 180 μm) under three contact angles is shown in Figure 5. The velocity profiles are found to be five-peak shaped, which correspond to five oil migration channels on the central axis. Velocity curve presents ascension, descension and interruption. This means that the velocity in the center of each channel is higher, while the flow velocity near the particle surface is lower due to the shear effect near the channel surface. In addition, the interruption of the velocity also indicates that the pores in real soil are discontinuous. It is also found that the velocity in the pore throat continuously decreases with time, which indicates that the oil migration level is weakened in the late stage. In terms of the flow velocity, there exists a higher value of oil flow velocity in the channel between 158.0623 μm and 179.4124 μm, which means that this channel is the dominant flow pathway for the leaked oil. With the increase of contact angle, the velocity peak distribution sequence in every channel is consistent, which shows that wettability does not change the dominant channel. In case of the occurrence of preferential flow, the leaked oil flows preferentially along this channel, which will eventually lead to a sharp increase of pollution distance just in the locality of the channel. In the actual process, this phenomenon often occurs in high permeability areas composed of connected macro pores.

Behavior of stripping and coalescence of water droplet
The dynamic process of stripping and coalescence of water droplet in the oil migration process under the contact angle of ~45° is shown in Figure 6. Note that the process has been observed by analyzing a series of dynamic micro-images. Although the leaked oil initially passes through the pore throat, a few tiny film-shaped water droplets still adhere to the surface of soil particle. The main reason for this phenomenon is that the shearing force is now not enough to separate the two small water droplets from the wall of soil grains. Due to the two kinds of oil behaviors including 'pull and push', the two water droplets become deformed, gradually approach. For the droplets with different sizes, the forces with smaller order of magnitude (such as additional mass force, pressure gradient force and the like) are ignored, and hence the drag force is the main force. The drag force of large droplets is obviously greater than that of small droplets, so the big droplet with a higher speed gradually approaches the adjacent small droplet until they contact and coalesce to form a bigger water droplet, eventually breaking away from the soil particle surface.
This phenomenon of the entrained droplets chasing and merging with the advancing interface is perfectly presented in Figure 7. It can be seen that the velocity of the entrained water droplet is slightly higher than that of the advancing interface in this pore throat, which eventually makes it catch up with the oil-water advancing interface. The friction loss of the flowing fluid in porous media mainly originates from two parts. One is the friction resistance at solid-liquid interface between the flowing fluid and the pore wall. The other is the interface friction resistance between two phase fluids in the pore throats. For the entrained droplet, it is only affected by the interface friction resistance between itself and the surrounding carrier oil. As for the surrounding carrier oil, it is not only influence by the  interface friction resistance between the surrounding carrier oil and the entrained droplets, but also by the friction resistance at the wall of the pore throats. It is obvious that the entrained droplets suffer less resistance and therefore have a higher velocity. Under the effect of the velocity difference, the distance between the entrained water droplet and the advancing interface decreases with time until they contact. In an instant of contact and merge, the liquid bridge connection is formed. Under the effect of the interfacial tension, the liquid bridge rapidly expands and its volume increases. Due to the restriction of the pore wall and the interfacial action of fluid, the entrained water droplet penetrates into the water phase in front of the advance interface. After reaching the steady state, the meniscus advancing interface continues to move forward along the pore throat.

Analysis of typically unswept area
The lower left corner of the model is unswept by the leaked oil under three wettability modes, and the location and shape of the unswept area are almost similar. Further, the reasons for the formation of the unswept soil pores are also discussed in detail, which are shown in Figure 8. Figure 8a clearly displays the inverted T-shape pore-throat is unswept by the leaked oil. When the leaked oil arrives point M and does not point N, the pressure difference between these two points continues to increase. However the maximum pressure difference between the two points does not invariably overcome the threshold capillary force at the narrowest pore throat, which leads that the leaked oil cannot pass through the pore throat. In oil arriving in point N, the pressure difference between two points descends rapidly, in which the stagnant water is still detained. Figure 8 b and c indicates that the velocity in this area is very low which may explain the stagnant water is still entrapped in the pores. In order to further analyze the fluid flow direction in the pore throat, the streamlines at different times are displayed in Figure 8c. It can be found that there are vortices in this area, which makes it difficult for fluid to flow in and out. Thus the leaked oil cannot migrate this area. This reveals the reason for the emergence and existence of stagnant water at the lower left of the model.

Conclusions
In this study, the effect of the soil wettability on the leakage oil migration at pore scale is examined through a series of numerical simulations. The observed porescale behaviors provide a detail understanding of the oil pollutants migration mechanisms and interfacial phenomena, which will contribute to precisely predict the environmental fate of the leaked oil pollutants. Several important conclusions are drawn as follows: (1) Due to the complex pore throat structure of real soil, the oil-water interface does not advance evenly under three wet states. In oil-wetting state, oil obviously migrates at a higher speed and becomes first to break through the model outlet with the shortest time. In the practical process, this eventually results in the extension of the oil pollutants range. In case of leakage from the buried pipeline, a large amount of oil enters into the soil, which easily causes a wide range of pollution. Pore scale simulation of oil pollutants migration process is helpful to clarify the flow characteristics of pollutants, so as to predict the flow range of pollutants with high accuracy. On this basis, emergency response shall be made to reduce the negative influences to the lowest level as far as possible.

Disclosure statement
No potential conflict of interest was reported by the author(s).