Seepage Characteristics of Mixed-Wettability Porous Media on the Phase-Field Model.

Clarifying the microscale gas-water flow behaviors in a mixed wettability reservoir is of great importance for underground engineering. A numerical model of mixed wettability based on circular particles was constructed using the MATLAB stochastic distribution program, and the gas-water flow was simulated based on the phase-field method. The Navier-Stokes equations were solved by the finite element method. The work analyzed the effects of the content of heterogeneous wetting particles, wettability, and inversed wettability of the matrix on the flow path and pressure distribution of the mixed wettability model. Besides, the two-phase flow behaviors were evaluated in microscale mixed-wettability porous media. The simulation results revealed that (i) the residual saturation of the gas phase showed a positive correlation with the hydrophobic particle content, and closed gases only existed in isolated pore channels with small content. Isolated closed gases gradually connected as the content increased. (ii) The residual gas content in the corner and tail end increased as the hydrophobicity of particles increased in hydrophilic matrices. Hydrophobic matrices showed a negative correlation, with the greatest pressure drop due to capillary resistance and step changes in the neutral-hydrophobic transition zone. (iii) Water-phase breakthrough time and gas-phase residual saturation showed a negative correlation change. The more space occupied by the gas phase, the faster the water-phase breakthrough. Moreover, the saturation no longer changes after the breakthrough. The work provides a guideline for determining the dominant flow path of phase displacements and the distribution of residual phases.

Introduction

The study on gas-water flow in the reservoir is significant for unconventional gas extraction, geological storage of carbon dioxide, and geothermal energy utilization.[1−5] There are pore and fracture structures with different development degrees in the reservoir.[6,7] A high-permeability pressure environment forms under engineering disturbances such as mining or fracturing. Unsaturated flow such as the water-driven gas and gas-driven water occurs in the reservoir (see Figure a), and the flow process is influenced by the pore distribution of media,[8] roughness,[9] and wettability.[10] Wettability can characterize the solid–gas–liquid interaction and determine the gas–water flow. However, geological bodies are generally composed of multiple components, and different components exhibit different hydrophilic and hydrophobic properties.[11] Quartz, feldspar, and mica in rock reservoirs generally exhibit hydrophilic properties, while talc and metal sulfides exhibit hydrophobic properties.[12] The wettability of the coal reservoir matrix is related to the degree of metamorphism,[13] and its inclusions of clay exhibit strong hydrophilic properties, while hydrocarbon organic matter is hydrophobic.[14] The different affinities of different wettability media for water directly affect the flow orientations of water in the media, which in turn affects strength degradation of rocks, mineral dissolution, and unconventional gas recovery. Therefore, the two-phase flow behaviors of gas–water in mixed-wettability media should be studied for the stability analysis of large underground projects and evaluation of oil-and-gas resource extraction.

Figure 1

Water-driven gas in porous media. (a) Reservoir exploitation and (b) pore throat structure.

Most of the existing studies have used the characteristic wettability components to represent the overall wettability characteristics, and the direct application of the results to reservoirs with complex mineral compositions causes large errors.[15] Cassie and Baxter[16] proposed a model for calculating the average contact angle of rocks by the volume fractions and contact angles of rock compositions. Besides, the relationship between surface wettability and single-component mixed-wettability media is characterized by the equation cos θ = ∑(f·cos θ). The model has been widely used in cases with a high degree of homogeneity of components and the low influence of heterogeneous wetting particles.[17−19] For mixed-wettability reservoirs, the diversity of rock components leads to a complex surface characterized by mixed wettability.[20] Mixed-wettability particles control the fluid trapping behavior under multiphase fluid coexistence, displacement relationships, residual phase distribution, and hydration reactions of rocks in reservoirs.[21] However, it is difficult for the quantitative formulation of different wettability components. The migration process of the air–water interface at the microscale cannot be fully understood and described due to the limited experimental conditions of indoor seepage.

The study of multiphase flow in the pore scale of the inhomogeneous wet reservoir[20,22] has become possible with advanced technical means. Gerami et al.[23] silanized and plasma-treated microfluidic chips etched with coal microfractures to obtain multiple network-channel structures containing hydrophilic, hydrophobic, and mixed wettability. Bright-zone coal (gas-wet) and dark-zone coal (water-wet) are compared by the water injection test. Harrison et al.[24] made a microscopic model of a porous medium containing hydromagnesite. Strong water absorption in the hydrophilic region of the model at the pore scale l increases the precipitation of hydration products in this region, with the change in wettability of the ″dry″ and ″wet″ regions. Then, the solid–gas–water interface gradually migrates with the hydration reaction, which in turn changes the pore structure of the model. Xie et al.[25] conducted a wettability test of shales containing nonclay minerals using the contact-angle measurement system based on the seated-drop method and the suspended-drop method. The effects of brine types and mineralization on mineral wettability are analyzed to reveal the effects of different wettability of different minerals within rocks on the shale reservoir.

More and more researchers have investigated it by employing numerical simulations due to a large number of uncontrollable and expensive factors in indoor tests.[26] Mahmud et al.[27] constructed a mixed wettability model in Eclipse by varying the wettability of the oil phase in contact with the solid surface in a hydrophilic matrix. JBN simulations are used to investigate the effects of recovery, saturation, and pressures on the relative permeability of the model. The nonwetting phase is enclosed in the small pores of the mixed hydrophilic medium at high capillary numbers, and the relative permeability values are lower than those of the homogeneous model.

The work focused on the interaction of the solid–gas–liquid interface in mixed-wettability reservoirs. Gas–water displacements in different mixed-wettability reservoirs were performed to explore the factors influencing the dominant channel of the water phase, breakthrough time, local confined gas distribution, and residual gas saturation under the pore scale of the reservoir. It provides strong theoretical and technical support for the in-depth study of the pore-scale gas–water distribution, long-term water–rock interaction, and gas recovery.

Governing Equations

Phase-Field Method

The existing means do not allow direct detection of flow and pressure changes in the channel and do not visualize the restricted behavior of the two phases in the microchannel; simulations were performed in COMSOL using the phase-field approach. The model is based on the Navier–Stokes equations to simulate the mass and momentum transfer of the fluid, and surface tension is added as a mass force to the Navier–Stokes equations in the two-phase flow:[28]where u is the velocity of the fluid, m/s; t is the time, s; ρ is the density of the fluid, kg/m3; μ is the dynamic viscosity of the fluid, mPa·s; Fst is the surface tension at the gas/water interface, N/m3; p denotes the pressure, Pa; I is the unit matrix.

The water displacement process requires tracking of the displacement interface, and the Cahn–Hilliard equation[29] is coupled with the Navier–Stokes equation using the phase-field approach to describe the diffusion interface of an immiscible two-phase fluid; the expression of the Cahn–Hilliard equation iswhere φ is the factor phase-field variable; ψ is the phase-field assistant variable; λ is the mixing energy density; ε is the interface migration thickness; γ is related to ε through γ = χε2; χ is the mobility parameter; fext is the external free energy density.

The relationship between mixing energy density, interface thickness, and surface tension can be expressed as[30]

The volume fraction of the two-phase fluid is expressed aswhere Vf1 is the repulsing phase; Vf2 is the repulsed phase; and −1 ≤ φ ≤ 1.

The volume fraction-weighted average method was used to calculate the gas–liquid phase mixing properties, and the fluid density and viscosity were expressed as[31]

Coupling wall wettability to the fluid equation:where n is the interface normal vector.

In the phase-field method, the surface tension can be calculated with the help of the diffusion interface representation according to the following expressions:[32]where G is the chemical potential at the interface.

Bringing the surface tension into the N–S equation yields:

Capillary Pressure

Figure b shows the conceptual model of water-driven gas in porous media, assuming that the matrix part is round particles and the interparticle part can be divided into pores and throats, and the throat is the key factor affecting the percolation performance of the pore medium. It is assumed that the pore throat is saturated by the gas phase at the initial moment, and the radius of the throat suddenly decreases during the water-driven process, increasing the local capillary force. In the hydrophilic environment, capillary force is expressed as a dynamic force, and in a hydrophobic environment, capillary force is expressed as resistance. The capillary force of a single capillary can be expressed as[33]where Pc is the capillary force; Pnw is the gas-phase pressure, Pw is the water-phase pressure, and Rh is the characteristic scale of the pore throat.

Validation of the Simulation Method

To verify the accuracy of the model, we first simulated the three-phase solid–liquid contact angle using the phase-field method and compared it with the spread of droplets on the solid surface when measuring the contact angle of rock (coal) by the seat-drop method. The experimental parameters and conditions are shown in Table . The droplet method is the most common and direct method of measuring the contact angle, using optical means to obtain the solid–liquid–gas three-phase contact line and then processing the image to obtain the contact angle value. This method can be used to measure flat surfaces, curved surfaces, and even solids with nonuniform surfaces. A comparison of the numerical simulation results with the experimental results of the coal–water static contact angle measured by Zhu et al.[34] at atmospheric pressure is shown in Figure and with the rock–water contact angle measured by Maribel et al.[35] in Figure , which shows a good agreement between the arch height and the three-phase contact line.

Table 1

Experimental Parameters and Conditions

fluid parameters	numerical value
water density, ρw (kg/m3)	1000
water viscosity, μw (Pa·s)	0.001
methane density, ρm (kg/m3)	0.648
methane viscosity, μm (Pa·s)	1.107 × 10–5
gas density, ρg (kg/m3)	1.169
gas viscosity, μg (Pa·s)	1.845 × 10–5
temperature, T (K)	298.15
pressure, P (Pa)	1.01 × 105

Figure 2

Coal–water static contact angle: (a) θ = 48.5°; (b) θ = 71°, the experimental photograph courtesy of Zhu et al.[34] Copyright 2020.

Figure 3

Rock–water static contact angle: (a) θ = 32°; (b) θ = 93°; (c) θ = 117°; (d) θ = 133°, the experimental photograph courtesy of Maribel et al.[35] Copyright 2018.

Implementation of Numerical Models

Model Assumption

This work aimed to study the effect of mixed wetting components on two-phase flow, so the microstructure of the reservoir is simplified to a homogeneous porous medium consisting of round particles and pores (see Figure a), and the following assumptions are made on the model:

Figure 4

(a and b) Simplification of the microstructure and geometric-model meshing.

The porous medium structure is a uniformly distributed circular one without considering the nonhomogeneity of the geometry.

No dissolution occurs between the liquid and gas phases during the two-phase displacement.

Both phases in the porous medium flow at a low velocity and fall within the range of laminar flow.

The gas phase under study satisfies the ideal gas equation of state.

Heterogeneous particles are randomly distributed in porous media.

Establishment of the Model

Figure b shows the geometric model of the porous media, where the white part indicates the solid part of the porous media and the gray part indicates the pore space. A rectangular area with a side of length 6790 × 4840 μm is used as the matrix. Circles and semicircles with radius R = 190 μm are arranged within the rectangular area at intervals lmin = 60 μm, and the calculated porosity is φ = 0.452. The left end is the injection end of the water phase, and an injection buffer with a width of 200 μm is constructed. Besides, the inlet boundary is set to injection v = 1 × 10–3 m/s at a constant flow rate. The right end was set as the outflow end of the boundary at a constant pressure, and the static pressure of the outlet was set to 0 Pa in turn. The upper and lower boundaries are symmetrically no-slip.[36] The mobility adjustment parameter of the phase-field interface is set to 1, and the interface thickness control parameter is set to half of the maximum number of grids in the interface flow area. A wet boundary is used on the surface of the circular particles, and the fluid-flow velocity is normal to the wall component of 0. The pore space varies widely in different regions of the model (lmax = 440√2 μm and lmin = 60 μm), so the embedded meshing method in COMSOL was used to refine the meshes within different pores by adjusting the number of meshes with narrow regional resolutions. The meshing in the throat region should be dense enough to ensure computational accuracy, and the meshes in the pore region should be coarsened appropriately to reduce the computational time. A triangular mesh was chosen to fit the boundary of the circular pore space. The dichotomous method was used to debug the density of grids. The calculation time of the model was shortened as much as possible by setting different grid divisions for the exact and the rough calculation areas. The whole model was divided into about 150,000 grid cells on the premise of ensuring calculation accuracy and convergence.

Random Distribution of Particles

Different reservoirs contained different amounts of heterogeneous wetting particles, and the contents of heterogeneous wetting particles varied greatly in different regions of the same reservoir. A total of 180 circular media in 12 rows and 15 columns were constructed in the model to characterize the feature in rock reservoirs. The RAND function of MATLAB was used to randomly select the circular media in porous media, which characterized the effects of inhomogeneous wetting particles on flow in porous media. Also, four combinations of inhomogeneous wetting particles were constructed with the contents of 18.9, 33.9, 43.3, and 52.8% (see Figure for the distribution of heterogeneous wetting particles).

Figure 5

Random distribution of heterogeneous wetting points in porous media; the content is case (a) 18.9%, case (b) 33.9%, case (c) 43.3%, and case (d) 52.8%.

Red was the repelled phase (water) and blue was the repelled phase (gas) in the simulation. The flow characteristics and residual-phase distribution were studied by changing the volume fractions of the two phases. A stable flow channel was formed after the breakthrough of the water-phase exit. The volume fraction was binarized using ImageJ software, and the threshold was adjusted to enhance the contrast of the gas–water region. Besides, the Analyze Particles algorithm was used to obtain the areas of the two phases after stabilization as well as saturation. The residual saturation of the gas phase in each pore column was normalized, and the area in the last pore column was removed in the calculation to reduce the capillary end effect on the results. Total residual saturation was plotted against breakthrough time as a change function of heterogeneous particles.

Results and Discussion

Homogeneous Scenarios

The water phase was injected into the pore and throat space from the left inlet of the homogeneous wetting model at a constant flow rate (see Table for the simulation conditions and Figure for the results). The displacing phases of the homogeneous model were uniformly distributed in the porous medium during the changing process in flow over time. Neither an obvious finger-in phenomenon nor large area around the flow phenomenon occurred. The whole process of water drive showed a uniformly layered displacement, and the flow pattern did not change during this process.

Table 2

Mixed Wetting Models with Different Wetting Particle Contents

cases	heterogeneous wetted particle content	matrix contact angle	heterogeneous contact angle
case 1c	18.9%	70°	120°
case 2c	33.9%	70°	120°
case 3c	43.3%	70°	120°
case 4c	52.8%	70°	120°

Figure 6

Volume fraction and pressure distribution in homogeneous wetting models.

The pressure decreased from the left inlet to the right outlet, and the velocity was larger at the narrow throat. t = 2.88 s when the water phase breaks through at the outlet. Only residual gases existed near the outlet, and the residual saturation of the internal gas phase was 0. Since little difference was marked between the hydrophilic and hydrophobic scenarios in terms of single-water saturation flow, both of them were categorized as the same type. The homogeneous model that acted as a control in the work was therefore not described here in detail.

Heterogeneous Particle Contents

Four different contents of heterogeneous wetting particles were randomly distributed in a hydrophilic matrix (70°) at a contact angle of 120°. Table and Figures , 8 show the model conditions and the simulation results, respectively. The interfaces advanced along the porous media wall with time. The water phase first occupies the smaller throat space of the matrix under capillary force and then gradually filled the larger pore space. The filling order was reversed at the heterogeneous particles due to the influence of capillary force between hydrophobic heterogeneous particles, and flow showed hysteresis. The model pressure dropped and residual saturation of the gas phase increased with the increased particle content. The hydrophobic particles in the scenario were isolated from each other and difficult to connect in the first half of the zone at a low content. Besides, residual gas saturation was small, and the hydrophilic particles in the matrix maintained good connectivity with each other. The aqueous phase advanced in the form of layered repetition and was delayed by capillary resistances when flowing through the heterogeneous particles. The heterogeneous wet particles had less influence on the overall repetition, so the model cases 1c and 2c had the lowest residual gas contents of 17.9 and 19%, respectively, distributed around a few isolated pore spaces and outlets.

Table 3

Heterogeneous Particle Models with Different Contact Angles

cases	particle content	matrix contact angle	particle contact angle
case 2a	33.9%	70°	60°
case 2b	33.9%	70°	90°
case 2c	33.9%	70°	120°
case 2d	33.9%	70°	150°

Figure 7

Variation of volume fraction and pressure distribution with a heterogeneous particle content.

Figure 8

Gas saturation and breakthrough time in different heterogeneous particle content scenarios: (a) local gas-phase residual saturation; (b) gas saturation and breakthrough time.

The hydrophobic-particle spacing became smaller in the high content scenario, and the residual gases were constrained by capillary force to accumulate at the hydrophobic points. Flow channels were compressed by local capillary force, and isolated residual gases were easily connected to the regions; therefore, residual gas saturation in the first half of the region was generally higher than in the low content scenario. The formation of dominant flow channels became more apparent as more hydrophobic particles were present. As the content of heterogeneous wet particles continued to increase (see case 4c), the aqueous phase dominant channels were compressed by capillary resistances, and residual gas saturation increased at a slower rate. The change in the breakthrough time of the aqueous phase at the outlet was opposite to the change in residual saturation. The high capillary resistance zone formed by hydrophobic particles caused the aqueous phase to form a dominant flow channel in the porous medium and breakthrough prematurely due to larger residual gas saturation in the latter half of the zone. Hydrophobic particles were more likely to form a finger at a high content. The greater the content of hydrophobic particles, the shorter the breakthrough time and the larger the amount of residual gas formed in the model.

Heterogeneous Particle Wettability

Heterogeneous particles with a content of 33.9% were randomly distributed in a hydrophilic (70°) matrix to characterize the effects of heterogeneous wetting particles on flow in mixed-wet reservoirs. Wettability was set to 60, 90, 120, and 150°, respectively (see Table for the model conditions and Figures and 11 for the simulation results). When the contact angles of the heterogeneous particles and matrix were similar (see case 2a), the water drive behaved as full displacement, and the residual saturation of the gas phase was 0. When the contact angle of particles increased to neutrality (see case 2b), the capillary force was 0. Closed gases gradually appeared inside the porous medium but mostly in isolated pore units.

Table 4

Hydrophobic Matrix Model

cases	heterogeneous wetted particle content	matrix contact angle	contact angle
*case 2a	33.9%	110°	60°
*case 2b	33.9%	110°	90°
*case 2c	33.9%	110°	120°
*case 2d	33.9%	110°	150°

Figure 9

Evolution of volume fraction and pressure distribution with inhomogeneous wetting.

Figure 11

Gas saturation and breakthrough time in water-wet matrix scenarios: (a) local gas-phase residual saturation confined gas; (b) gas saturation and breakthrough time.

Residual gas saturation increased from 0 to 12.8% when the contact angle was hydrophobic (see case 2c). The water-phase breakthrough of the pore containing isolated closed gases increased with resistant capillary force. The upper boundary with a large number of gas phases failed to repel, and residual gas saturation was 17.5%. As the contact angle continued to increase (see case 2d), the residual saturation of the gas phase reached a maximum of 18.8%. Heterogeneous particles exhibited super-hydrophobicity, and the water phase needed to break through a greater pressure when flowing through the throat. The water phase was forced to bypass the strong resistance zone to flow from the matrix, which formed a continuous gas-phase closure zone around these particles (see Figure a). The first half of the case had the most homogeneous flow, and the finger-in phenomenon was not obvious.

Figure 10

Microscale gas–water distribution display: (a) confined gas; (b) liquid film on the surface of water-wet particles.

Case 2d showed a sharp decrease in the pressure of the water phase after the breakthrough because the aqueous phase breaks through to form a conduction channel along which the subsequent injected phase flowed with minimal pressure and without overcoming additional capillary pressures. The residual saturation of the gas phase in the first half of all cases was less than 30%; the residual saturation varied with the increased hydrophobicity of heterogeneous particles in the second half. The greater the hydrophobicity, the greater the residual saturation. The whole model exhibited the maximum pressure drop. Based on the simulation results of Figures b and 13b, using the nonwetting phase to displace the wetting phase contributed to a more residual displaced phase, which was unfavorable for gas escapes. This was attributed to the affinity of the wetting phase for coal, which made the wetting phase form thin films adhered to porous-medium surfaces (see Figure b).

Figure 13

Gas saturation and breakthrough time in gas-wet matrix scenarios: (a) local gas-phase residual saturation; (b) gas saturation and breakthrough time.

Gas-Wet Matrix Scenarios

Chemical reagents were often injected into unconventional oil and gas exploration in geological reservoirs to improve wettability. The wettability of the matrix was changed from hydrophilic to hydrophobic in this section. Table and Figures , 13 present the simulation conditions and results, respectively. The capillary force was expressed as resistance in the hydrophobic matrix. The inlet to outlet exhibited a greater pressure drop than in the hydrophilic model, with more residual gases. The regional saturation curves overlapped more when heterogeneous particles were hydrophilic and neutral (see *cases 2a and 2b), and the saturation curves overlapped more when wettability was hydrophobic (see *cases 2c and 2d). The saturation curves also had a high degree of overlap when wettability expressed hydrophobic (see *cases 2c and 2d).

Figure 12

Evolution of volume and pressure distribution of the wettability inversion for the matrix.

Figures a and 13a show that residual saturation in the first half of the gas-wet matrix scenario was greater, and more residual gases were connected into connected regions in the first half of the scenario. The water drive was blocked and stalled at the hydrophobic throat. When the leading edge of the water phase came into contact with hydrophilic particles (see *case 2a), the local flow rate increased, and bypassing occurred to form locally closed gases, which exhibited maximum residual saturation (42.4%). When the particle contact angle was neutral (see *case 2b), the capillary force was 0 and the model contained a large amount of circled closed gases. However, residual gas saturation was slightly reduced (41.6%) compared to *case 2a, with the obvious finger-in phenomenon.

As the particle contact angle continued to increase, both the matrix and the particles exhibited hydrophobicity (see case 2c). Only a small amount of closed gas remained in the pore channel in the front part of the model under capillary resistances, and more closed gases remained in the blind end of the model wall after the water-phase breakthrough. Meanwhile, the trapped gas could be finally expelled or permanently retained depending on the size of the microstructures and the pressure exerted by the subsequent injection of fluids. Deadend microfractures tended to hold residual water permanently. When the particle’s contact angle increased to 150° (see case 2d), the capillary resistance at particles was 2.5 times that of the matrix, and the inlet pressure reached its maximum. Residual gases in the matrix were enclosed in a smaller space by capillary force. The water phase exited the breakthrough time. The residual saturation of the gas phase showed an obvious negative correlation and increased with the increased particle contact angle. The increase rate was consistent with decreased residual saturation.

Conclusions

The effects of inhomogeneous wetting particle content, wettability, and matrix wettability on displacements were investigated by simulating two-phase displacements in mixed wetting porous media. Besides, the work analyzed the relationships between the mixed wetting characteristics of the reservoir and the dominant channel of the water phase, breakthrough time, local confined gas distribution, and residual saturation of the gas phase.

Different mixed wetting models were constructed by randomly aligning the wetting points with the MATLAB random distribution function. The water drive behaved as a stratified flow in the homogeneous wetting model, and both pro- and hydrophobic media could reach the full replacement. The local change in particle wettability after reaching saturation did not affect flow.

When the content of hydrophobic wetting particles in the hydrophilic matrix was small, the effect on the residual phase was small, and residual gases were distributed in the corner area and the closed pore space.

Heterogeneous particles’ contact angle and matrix contact angle were similar to the full displacement in the hydrophilic matrix. The pore containing closed gas increased with the increased particle contact angle under capillary resistance. Residual gases were connected at the trailing end, which compressed the water phase to flow to the outlet in narrow channels.

The leading edge of the water phase flowing through the hydrophilic particles increased local flow velocity around the flow to produce closed gases, and the formation of hydrophilic channels in the region was occupied by hydrophilic particles. As the contact angle of heterogeneous particles increased, the residual saturation of the gas phase nonlinearly decreased, and the local pressure difference increased under capillary resistances. The model showed a higher pressure drop.

The changing trend of the breakthrough time of the water phase exited and the change of gas-phase residual saturation showed good regularity (except for the full displacement state). More gas-phase residuals meant a shorter breakthrough time of the water phase. The breakthrough time was always between 1.5 and 3 s.

The mixed wetting characteristics of the reservoir caused a nonhomogeneous distribution of the gas and water phases. The nonhomogeneous flow of the water phase in the thermal reservoir caused nonhomogeneous heat transfer at the solid–liquid interface, leading to local nonthermal equilibrium. We will consider the effects of changes in temperature difference-induced surface tensions (Marangoni effect) on nonhomogeneous flow in microscopic pores based on this model in the next work.