Effect of Wettability Heterogeneity on Water-Gas Two-Phase Displacement Behavior in a Complex Pore Structure by Phase-Field Model

Abstract

Understanding the immiscible displacement mechanism in porous media is vital to enhancing the hydrocarbon resources in the oil and gas reservoir. Improving resource recovery requires quantitatively characterizing the effect of wettability heterogeneity on the immiscible displacement behaviors at the pore scale, which can be used to predict the displacement distribution of multiphase fluids and evaluate the optimal wettability strategy in porous media. The heterogeneity of fluid wettability in a natural rock makes it extremely hard to directly observe the fluid displacement behaviors in the reservoir rocks and quantify the sensitivity of preferential displacement path and displacement efficiency to wettability distribution. In this study, the phase-field model coupling wettability heterogeneity was established. The gas-water two-phase displacement process was simulated under various wettability distributions and injecting flux rates in a complex pore structure. The effect of wettability heterogeneity on immiscible displacement behavior was analyzed. The results indicated that wettability heterogeneity significantly affects the fluid displacement path and invasion patterns, while the injecting flux rate negatively influences the capillary–viscous crossover flow regime. The continuous wetting patches enhanced the preferential flow and hindered displacement, whereas the dalmatian wetting patches promoted a higher displacement efficiency. The results of the fractal dimensions and specific surface area also quantitatively show the effects of wettability distribution and heterogeneity on the complexity of the two-phase fluid distribution. The research provides the theoretical foundation and analysis approach for designing an optimal wettability strategy for injecting fluid into unconventional oil and gas reservoirs.

1. Introduction

Accurate understanding and description of multiphase fluid flow in reservoir rocks are vital for exploiting unconventional oil and gas and improving recovery efficiency in energy engineering. Despite significant advancements in the multiphase flow mechanism in porous media [1,2,3,4,5,6], our understanding of the multiphase flow behavior at pore scale under versatile conditions is still in its infancy. Pore-scale analysis and characterization of the multiphase fluid displacement behaviors will contribute to predicting the preferential flow of oil and gas resources and evaluating the recovery efficiency of the subsurface hydrocarbon resources. However, because unconventional oil and gas reservoirs often have a burial depth of more than 3000 m, the complex rock properties and geological conditions make it extremely difficult to visually observe and quantitatively analyze the dynamic characteristics and migration paths of multiphase fluids in porous media [7]. Moreover, the development of fluid displacement paths in porous media is susceptible to flow conditions, and several physical and chemical properties affect the dynamic behaviors of immiscible fluids. Recent studies mainly focused on the properties containing the capillary number, viscosity ratio, pore geometry, and wettability heterogeneity [8,9,10,11,12,13,14]. Specifically, fluid wettability on the surface of porous media shows a substantial impact on fluid saturation during the immiscible phase displacement process and, therefore, significantly affects the fluid displacement efficiency in the hydrocarbon recovery of the oil and gas reservoirs [1,15,16].

Wettability in porous media is generally classified as homogeneous or heterogeneous [17]. Owing to the complex mineral compositions of natural rock reservoirs, the rock surface wettability, which changes significantly with the mineral properties of matrices or the reservoir environment [18,19,20], is heterogeneous at the core and pore scales and has a broad spectrum ranging from water-wet to mixed-wet and oil-wet [21]. The pore structures of shales have dual surface wettability with gas-wet organic pores and water-wet inorganic pores [22,23]. Radke et al. [17] stated that two types of heterogeneous wettability, mixed wettability and dalmatian wettability, are generally recognized. Mixed wettability refers to two continuous and distinct types of wet surfaces that coexist and span a porous medium. Dalmatian, also speckled or spotted, wettability refers to one continuous wet surface enclosing regions of another discontinuous wet surface or vice versa [24]. Heterogeneous wettability can result in different heterogeneity scales within a reservoir [25]; thus, the interior pore structure shows versatile physical and chemical properties in the pore walls and the presence of residual fluid bubbles and clusters, known as a mixed-wet system. In such systems, the stability and movement of residual fluids significantly affect the contacting properties of invading fluid and lead to an unknown development of preferential flow path occurring at the pore scale. The complex wettability heterogeneity of rock reservoirs makes it more challenging to understand the multiphase displacement behaviors in such a mixed-wet pore system, which further requires quantitative investigations.

Previous studies have found that the wettability characteristics of solid surfaces and the spreading characteristics of fluid systems play vital roles [26,27,28,29] and have identified their importance for pore-scale fluid behaviors in geological engineering [21,30,31,32,33]. The wettability influence on multiphase flow behavior, capillary pressure, and relative permeability were also explored numerically and experimentally in earlier stages using the two-dimensional (2D) optical method or three-dimensional (3D) imaging method. Zhao et al. [30] performed microfluidic experiments to study the immiscible displacement process in the 2D micromodels with vertical posts and revealed a critical affinity of water on the solid surface for a more effective displacement of the non-wetting phase. Avendano et al. [34] investigated the wettability effect on immiscible displacement in glass microfluidic models under oil-wet and water-wet conditions. They found water-wet pore system exhibits more and larger trapped oil clusters, leading to a lower displacement efficiency compared to the results in the oil-wet device. These studies analyzed immiscible displacement mechanisms in porous media with different wettability conditions. At the same time, most of them focused on wettability systems in homogeneous pore structures [35,36] or small-scale heterogeneous wettability distribution in individual pores [37].

Averaging contact angles described the effects of spatial wetting distribution over several pores [23,38], which cannot sufficiently characterize the non-spreading systems [39] or flow heterogeneity [40]. The wettability effect on the immiscible displacement process is more difficult to understand when the porous medium exhibits complicated wettability distribution, namely in the mixed-wet system. Chang et al. [41] investigated and discussed the influence of mixed-wettability on brine drainage for a mixed-wet micromodel. Murison et al. [38] reconstructed a mixed-wet system by randomly packing functionalized glass beads with well-defined spatial wettability correlations and visually performed two-phase invasion experiments in such porous media and concluded that the wettability heterogeneity controls the flow dissipation. More recently, Rücker et al. [42] also experimentally indicated that the flow regime in a mixed-wet rock-oil-brine system showed significantly more ganglion dynamics than in water-wet rock.

Moreover, the mixed-wettability distribution was probed experimentally in the reservoirs, and the researchers also attempted to explore the effect of wettability distribution on capillary pressure, relative permeability, and fluid displacement efficiency in porous media. Foroughi et al. [43,44] proposed an algorithm to estimate the randomly assigned distribution of fluid wettability in the reservoir rocks by matching the macroscopic properties experimentally obtained and discussed the effect of spatial correlation in wettability on the predicted capillary pressure and relative permeability. Garfi et al. [45,46] used X-ray imaging to characterize the wetting state of fluids within the pores of reservoir rocks and visualized the spatial distribution of fluid coverage on the mineral surface. Their results revealed variant wettability properties for the mineral compositions, leading to mixed wettability in the rock. These studies mainly focused on investigating and characterizing the spatial wettability distribution in the natural rock pore structures and demonstrated the vital role of wettability heterogeneity in the multiphase flow. However, the sensitivity of immiscible displacement behaviors to wettability heterogeneity in porous media has not been well understood, especially in a heterogeneous pore structure. Little attention has been given to quantifying the effect of wettability heterogeneity and types on fluid flow in heterogeneous structures at the pore scale, resulting in a lack of accurate knowledge of the multiphase flow mechanism in rocky pore structures.

As an alternative effective research tool, many numerical simulations have also been carried out to investigate the multiphase flow behavior in the mixed-wet rocky structures. The pore network modeling approach [47,48] was developed to investigate the fluid flow in heterogeneous pore structures by assigning different contact angles to different pores in the network. Aghaei and Piri [47] presented a new dynamic pore network model capable of up-scaling two-phase flow processes from pore to core and carried out the core-flooding experiment validations. Song et al. [23] reconstructed a 3D shale pore network model to study nanoscale confined gas and water transport behavior in dual wettability nanoporous shale, while the model was restricted in the piston displacement during injected water flow in a process according to water advancing contact angle measurement results in the literature. The results also concluded that accurate modeling of two-phase flow behavior should consider the significantly different surface wettability in the dual rocky pore structures. Cha et al. [49] performed the direct pore-scale modeling to investigate the effects of mixed wettability and stressed the great impact of wettability heterogeneity on the fluid distribution. The wettability heterogeneity in rock pore structures has been experimentally and numerically characterized and studied. However, a quantitative analysis is still inadequate concerning the effect of wettability heterogeneity on the multiphase flow pattern and paths in the complex pore structure. Furthermore, the influence of dalmatian wettability on the two-phase flow behavior in porous media has rarely been discussed.

The phase-field method is popular in modeling multiphase flow in porous media. To simulate the interface dynamics of the immiscible fluids, the phase-field method utilizes order parameters such as the concentration or volume fraction of the fluids to represent different fluid phases. Thus, the sharp interface of immiscible fluids is treated as a transitional region with finite thickness. The order parameter is typically evolved by the Cahn–Hilliard (C–H) equation, the Allen–Cahn equation, or other types of dynamic equations to capture the interfacial motion of the immiscible fluids. Various numerical methods have been developed for modeling the displacement process of immiscible fluids based on the phase field method, such as the phase-field-lattice-Boltzmann method [50,51] and the phase-field-Navier-Stokes method [52,53]. In the former method, the C–H equation controlling the fluid interface is transformed into the corresponding LBE (Lattice Boltzmann Equation), and the combination of LBE equations controlling the interfacial dynamics and interior flow of immiscible fluids are solved simultaneously. Moreover, in the latter method, the combined equations of the phase-field model and the Navier–Stokes equations are solved directly.

In the present study, the phase-field model coupling wettability heterogeneity was established based on the pore networks of natural reservoir rocks. A series of simulations with monotonically increasing gas-wet fraction and two wettability heterogeneity types were performed to investigate the effect of wettability heterogeneity on the dynamic flow behavior of immiscible two-phase fluids in the complex pore structure. The sensitivity of flow patterns and paths to the wettability heterogeneity were quantitatively analyzed. Our results intuitively revealed the significant effect of wettability heterogeneity and types on the flow paths and displacement patterns during the two-phase displacement process. The present study provides a research basis for quantitatively analyzing the multiphase flow behavior in natural rock structures.

2. Numerical Models

2.1. Establishment of Heterogeneous Wetting Pore Models

The porous network used in the present study was composed of nine identical substructures in a 3 × 3 array, which could exclude the effects of the pore structures, as indicated by the red dashed lines in Figure 1a. For the same pore structure, different wettability distributions and heterogeneity were set to examine the effects of wettability distribution and heterogeneity on two-phase flow patterns. Similar to our previous relevant studies [53,54], the substructure was also based on studies by Sirivithayapakorn and Keller [55] and Auset and Keller [56]. The distance map in Figure 1b shows the pore radius, defined as the maximum value of the local distance of each void point away from the solid wall. In ImageJ, it was computed using the maximum ball algorithm and is equal to the inscribed radius of the pore space. The pore radii varied from 800 μm to 3000 μm. As shown in Figure 2a, the pore radius has a wide distribution range, and the corresponding pore volumes at various pore radii were not concentrated at a specific value. Therefore, the pores with various pore radii had approximate volume fractions, and their structural importance was equivalent, showing the heterogeneity of pore structures. The geometric heterogeneity of pore structure was also computed via the self-programming code in Matlab, see Figure 1c. Here, the heterogeneity of pore distribution is quantified using the relative standard deviation of each block’s porosity [57]. The heterogeneity stabilized at around 0.65 with the increasing number of domain divisions, which quantitatively proves the heterogeneity of this pore structure.

From the porosity tests, the porous structure had a representative element volume (REV), and its porosity stabilized at 0.69 (see Figure 2b); therefore, it was applicable for integral rock cores to analyze the two-phase flow behavior and mechanism in such porous structures. It is essential to understand the influence mechanism of wettability and structure heterogeneity on the two-phase flow behavior in such heterogeneous pore structures.

2.2. Case Studies and Boundary Conditions

Twelve cases with different wettability distributions were designed to simulate the gas-water flow process in porous structures by assigning different gas/water contact angles based on the desired wettability pattern, as shown in Figure 3. We used the ratio of the gas-wet area $r_{g-w}$ to quantify the wettability distribution:

$$r_{g-w}=\frac{A_{g-w}}{A}$$

(1)

where $A_{g-w}$ is the gas-wet area and A is the total wet area. The final distribution of gas and water in porous structures is established by simulating a two-phase displacement, corresponding to drainage in the gas-wet region and imbibition in the water-wet region. Cases A1, A4 and B1, B4 ($r_{g-w}$ = 0, 1) were used to study the effect of homogenous wettability on the two-phase flow. Cases A2, A3 and B2, B3 ($r_{g-w}$ = 1/3, 2/3) were used to study the effect of continuous wettability heterogeneity on the two-phase flow. Cases A5, A6 and B5, B6 were obtained by adjusting the gas-wet wall positions of A2, A3 and B2, B3 under the same wettability ratio, $r_{g-w}$ = 1/3, 2/3, and were used to examine the effect of alternating wettability heterogeneity on the immiscible two-phase flow paths.

Owing to the higher viscosity and density ratio of gas and water, the phase-field method using the Cahn–Hilliard equation [52,58] was used to solve the gas-water fluid systems, formulated as follows:

$$\frac{\partial \varnothing }{\partial t}+u·\nabla \varnothing =\nabla ·\frac{\gamma \lambda }{{\epsilon }^2}\nabla \Psi$$

(2)

and

$$\Psi =-\nabla ·{\epsilon }^2\nabla \varnothing +\left({\varnothing }^2-1\right)\varnothing$$

(3)

where $\varnothing$ is the phase-field parameter, $u$ is the fluid velocity (m/s), $\gamma$ is the mobility (m³·s/kg), $\lambda$ is the mixing energy density (N), and $\epsilon$ is the interface thickness parameter (m). The variable $\Psi$ is referred to as the phase-field help variable. The compressive Navier–Stokes equations [59] were used to describe the transport of mass and momentum for fluids:

$$\rho \frac{\partial u}{\partial t}+\rho \left(u·\nabla \right)u=-\nabla p+\nabla ·\left(\mu \left(\nabla u+{\left(\nabla u\right)}^T\right)-\frac{2}{3}\mu \left(\nabla ·u\right)I\right)+F_{st}$$

(4)

$$\frac{\partial \rho }{\partial t}+\nabla ·\left(\rho u\right)=0$$

(5)

where, $\rho$ denotes the density (kg/m³), $\mu$ denotes the dynamic viscosity (Pa $·$ s), $u$ represents the velocity (m/s), $p$ denotes the pressure (Pa), and $F_{st}$ is the surface tension force acting at the two-phase interface [52]. Herein, we performed the simulation of static contact angles for immiscible fluids on a solid surface under different wettability conditions. Figure 4 shows various interfacial configurations of immiscible fluids with contact angles from 15° to 165°. The used phase model has been verified, and more detailed theories and governing equations can be found in our previous studies [53,54].

During a simulation of the gas-water flow process in porous media, the fluids inside the pores did not penetrate the solid parts (black area in Figure 1a). The displacing phase fluid was gas, and the displaced phase was water. The flow rate q_v = 3 mL/min and 10 mL/min were set as a constant at the inlets, and the gas volume fraction was one at the inlets in the simulations.

Table 1. Parameters used for the case studies.

No.	qv (mL/min)	${\mathit{r}}_{\mathit{g}-\mathit{w}}$	Wettability	μw (×10−3 Pa·s)	μgas (×10−5 Pa·s)	ρw (kg/m3)	ρa (kg/m3)	θ (ο)	σ (×10−2 N/m)
A1	3	0	Mixed	1	1.81	1000	1.28	65.5/120	7.21
B1	10
A2	3	1/3
B2	10
A3	3	2/3
B3	10
A4	3	1
B4	10
A5	3	1/3	Dalmatian	1	1.81	1000	1.28	65.5/120	7.21
B5	10
A6	3	2/3
B6	10

summarizes the basic data and fluid properties used in the simulations. The capillary number Ca (= (μ × v)/σ) and mobility ratio M (=μ_gas/μ_w) were used to determine the various flow conditions, where μ and v are the viscosity and Darcy velocity of the fluid, respectively, and σ (=0.0721 N/m) and θ are the interfacial tension and contact angle of the two-phase fluids, respectively. The mobility ratio is 0.0181 and the capillary numbers are 6.28 × 10⁻⁵ and 1.26 × 10⁻⁴ for cases A1–A6 and B1–B6, respectively. According to the phase diagram (logCa−logM) for the interfacial stability areas indicated by Zhang et al. [60] and Lenormand et al. [61], these scenarios are expected to fall within the crossover range from a capillary to a viscous fingering regime. Figure 5 depicts the computational domain and the initial and boundary conditions used in the simulation. The model used tetrahedral and adaptive meshes to optimize the two-phase interface to ensure the accuracy of the calculations. The detailed results are discussed in Section 3.

3. Results and Discussion

Figure 6, Figure 7 and Figure 8 provide snapshots of the two-phase fluid distributions for all cases at the quasi-steady state (when the gas saturation remains unchanged mainly). Figure 9 gives the variation of gas saturation as a function of time for all cases. In all cases, the highest displacement efficiency was obtained in the uniformly water-wet medium A1, coinciding with homogenous pore models from Laroche et al. [28]. Figure 6 gives the phase distributions in porous structures with homogenous wettability. More branches, especially the displacement path along the vertical direction, were formed in the entire gas-wet structures (A4 and B4) during the gas-water displacement process. Although A1 and B1 reached a quasi-static state earlier than A4 and B4, respectively, A1 and B1 presented a uniform horizontal displacement and had no obvious preferential paths during the initial displacement period. Models A4 and B4 had three finger paths in the initial horizontal displacement, followed by vertical displacement to invade the porous regions in the side boundary. The immiscible interface in models A1 and B1 shows higher gas saturation than A4 and B4, respectively. The interface evolution indicates that the gas phase prefers to uniformly displace the water phase in entire water-wet structures, compared to the cases with a complete gas-wet system. Gas entered the broader channels and exhibited pronounced preferential flow during the immiscible two-phase displacement process in the complete gas-wet structures A4 and B4 (see Figure 6b). At this stage, the flow paths were mainly determined by the pore patterns, which have been discussed and verified in a previous study [53]. Within the range of flow regimes used in this study, for the homogenous wetting structures, two kinds of wettability types led to totally different flow paths in the same porous structures, showing that the wettability significantly influences the preferential paths and displacement patterns during the immiscible two-phase flow process.

For mixed wettability structures (A2/B2 and A3/B3), gas first passed through the gas-wet channels and caused developed finger patterns (Figure 7). In these models, the gas-wet walls play a guiding role in the immiscible interface evolution. In addition, gas flowed through the lower channels rather than the top channels, despite having the same pore structure and wetting distribution. It can be attributed to the pore morphology, indicating that the development of a preferential flow path is simultaneously affected by the pore geometry heterogeneity.

For the dalmatian wetting structures A5/A6 and B5/B6, complicated preferential paths were formed during the immiscible fluid displacement process, as shown in Figure 8. A6/B6 reached a quasi-static state earlier than A5/B5 and obtained lower gas saturation. In the assembled heterogeneous pore structures, a more significant fraction of the gas-wet zone facilitated the formation of preferential fluid flow and hindered a compact displacement pattern. The fluids’ distribution in the pore models contains both viscous fingering and capillary fingering simultaneously. Compared to the invasion patterns in models A5/B5, the vertical fingers are suppressed, and the immiscible interface tends to forward along the fluid injecting direction in models A6/B6.

In this capillary–viscous crossover flow regime, the spatial wettability correlation also showed a necessary effect on the displacement paths and fluid patterns of the two-phase system, comparing the interface dynamics between the pore models A2/B2 and A5/B5, as well as A3/B3 and A6/B6. Figure 9 shows that the gas saturation and breakthrough time gradually decreased with an increased wettability ratio $r_{g-w}$. Contrast groups A2 and A5 (B2 and B5), and A3 and A6 (B3 and B6) in Figure 7 and Figure 8, respectively, had the same wettability ratio $r_{g-w}$. Murison et al. [38] performed immiscible displacement experiments in a bead pack consisting of monodisperse wetting and non-wetting spheres. Their experimental results indicate that the spatial wettability correlation impacts the irreducible water saturation, and the pore structure with a weak correlation of the same-type wettability shows a higher saturation of the nonaqueous phase. We also obtained a similar trend for the effect of spatial wettability correlation on the gas saturation; namely, the ultimate gas saturation in models A5 and A6 (the dalmatian wetting structures) is higher than the value in models A2 and A3 (the mixed wetting structures), respectively. For the same wettability ratio, the dalmatian wetting patches more easily enhanced displacement than the mixed wetting patches. Continuous displacement phase wet patches guided the formation of preferential flow paths and hindered the gas-phase invasion into neighboring pore space, whereas alternating wet patches improved flow sweep efficiency in complex porous structures. Similarly, Cha et al. [49] proposed a heterogeneous index (HI) to quantify the wettability heterogeneity and demonstrated that the case with lower HI (i.e., less heterogeneous wettability distribution) leads to a higher saturation of injecting fluid at breakthrough time. Although the displacement conditions of immiscible fluids are different between these studies and our work, a common understanding is that the spatial correlation of the wettability heterogeneity state affects the fluid distribution, and the pore structure with less wettability heterogeneity is beneficial to a higher displacement efficiency.

Moreover, this study also explored the wettability influence on immiscible displacement behaviors under different injecting flux rates. Comparison in the two-phase distributions between A-models and B-models indicates similar displacement patterns and gas saturation. In this capillary–viscous crossover flow regime, the injecting flux rate has a weaker impact on the immiscible displacement process than the influence resulting from the wettability heterogeneity.

Figure 10 shows the temporary variation of the gas front position from A1–A6 in the horizontal direction as a function of time. In six cases, A1 did not exhibit obvious preferential paths with the entire water-wet walls and had the highest saturation variation ratio with the x-position and the most uniform displacement pattern. Therefore, the capillary flow promoted better displacement. However, A4 exhibited apparent preferential flows from entire gas walls and obtained a poorer displacement effect. For the structures with the same wettability ratio but two different types of heterogeneity, A2 and A3 exhibited more obvious preferential paths and obtained poorer displacement than A5 and A6, respectively. The mixed wetting patches enhanced the preferential flow in the displacing phase wetting channels and hindered the displacement, yet the continuously displaced phase wet patches promoted a uniform flow pattern and improved the displacement efficiency. In actual engineering, it can be attempted to suppress such preferential paths to improve oil and gas production.

To quantify the effect of wettability heterogeneity on flow patterns, we used fractal dimension D and specific surface area A_ss to characterize the shape of the flooded channels. The fractal dimension D of all cases at the quasi-steady state was computed using the fractal box counter method in ImageJ software. Referring to previous studies [62,63], the specific surface area A_ss is defined as the area per unit length (m⁻¹) and is calculated using the following formula:

$$A_{ss}=\frac{P}{A}$$

(6)

where, A and P are the area and edge length of the injected fluids, respectively.

Figure 11 shows that the fractal dimension of the flooded channels decreased with the wettability ratio $r_{g-w}$ despite the minor differences; therefore, the wettability heterogeneity affected the complexity of the two-phase fluid distribution. In essence, the fractal dimension reflects the complexity of the spatial distribution of the invading fluid, which is determined by both wettability heterogeneity and geometry heterogeneity. In this study, the wettability heterogeneity for the gas-wetting ratio induced essential changes in the complexity of immiscible fluids distribution at the quasi-static state. However, the coupling effect between wettability heterogeneity and geometry heterogeneity is not fully known yet and needs further investigation. Comparing cases A2 and A5 (B2 and B5), and A3 and A6 (B3 and B6), the types of wettability heterogeneity (mixed wettability or alternating wettability) had little effect on the fractal dimensions of the flooded channels. Similarly, in Figure 12, the specific surface area A_ss increased with the wettability ratio $r_{g-w}$ and then decreased, and the types of wettability heterogeneity slightly affected the specific surface area A_ss. This weak impact was also demonstrated in Foroughi et al.’s work [43] (2021), which indicated that spatial correlation in wettability is not significant in the rock samples.

To analyze the effect of wettability heterogeneity further in the multiphase flow, the relative permeability curves under different wettability conditions are shown in Figure 13. Most points of the water relative permeability are close to zero in cases A2 and A3 with continuum wetting-state distribution. Moreover, cases A5 and A6 with the same wetting surface fraction show higher fluid permeability, which is attributed to the trapping of the water phase as the large cluster, as shown in Figure 7. The 3D pore characteristics of pore structure are restricted in our models, such as pore connectivity and pore number. Therefore, the relative permeability curve only qualitatively demonstrated the impact of spatial wettability heterogeneity, and further research is required to quantify the effect of wettability heterogeneity concerning the spatial distribution and surface fraction for different contact angles.

4. Conclusions

We reconstructed a pore structure model with heterogeneous pore morphology and heterogeneous wettability distribution in the present study. The phase-field method coupling wettability heterogeneity was established to simulate the gas-water two-phase displacement under different injecting flux rates in the reconstructed pore structures. The effect of wettability heterogeneity and types on the flow characterization and paths during the immiscible displacement process were quantitatively discussed. The results indicated the significant role of wettability heterogeneity in the dynamic behaviors of two-phase fluids in a complex pore structure, particularly in the capillary–viscous crossover flow regime. The key findings are summarized as follows:

(I)

The surface fraction of gas-wetting and alternating wettability heterogeneity all present an essential effect on the fluid displacement path and invasion patterns, while the injecting flux rate has less influence in the capillary–viscous crossover flow regime.

(II)

Among various surface fractions of gas-wetting, a uniformly water-wet medium prefers to produce a higher gas-phase displacement efficiency.

(III)

For heterogeneous porous structures with the uniform gas-wetting distribution, the dalmatian wetting structure has a higher displacement efficiency than the mixed wetting structure during the immiscible two-phase displacement process.

(IV)

The complexity of the two-phase fluid distribution, characterized by the fractal dimension and specific surface area, is sensitive to the wettability heterogeneity, while the spatial wettability correlation has a weaker impact.

This work established a pore structure model with various wettability distributions and quantitatively illustrated the effect of wettability heterogeneity on the two-phase displacement behaviors in a complex pore structure. The present study leads to more understanding of the two-phase displacement behaviors in porous media with mixed wettability. We expect to further characterize the complex 3D wettability distribution in natural rock and reveal its quantitative effect on the oil and gas recovery efficiency.