Regulation of Critical Capillary Number: High-Efficiency Displacement Mechanism of Deep Coalbed Methane Considering Wettability Heterogeneity of Proppants between Fractures

1  Introduction

Coal reservoirs, as the primary medium for the storage and migration of coalbed methane, the permeability of their internal fractures directly limits the flow capacity of gas-water two phases within the reservoir [1,2]. In their in situ state, coal reservoirs often face challenges such as limited effective flow channels and low permeability. Although technologies for increasing fracture density through artificial means, such as hydraulic fracturing, have matured [3,4], these fractures are frequently closed under reservoir pressure. This closure behavior results from the stress sensitivity inherent to coal reservoir fractures, which exhibit nonlinear closure behavior with increasing confining pressure; fracture compressibility is significantly correlated with effective stress [5,6]. This characteristic is particularly pronounced in the high-pressure environments of deep reservoirs [7,8]. To increase the degree of crack opening, quartz-based proppants are typically injected during the process of seam making. This results in fractures containing filling materials such as proppants and coal bridges. The effective flow channel of the fracture cannot be simplified as an open parallel plate, but tends to be more like a porous medium model [9,10]. Conventional fracture two-phase flow often utilizes unified wettability data to examine the characteristics of effective seepage paths and the residual phase distribution, ignoring the differing wettability among proppants, coal bridges, and other filling materials, leading to significant variations in effective seepage paths. These differences directly impact the judgement of real permeability of reservoirs [11]. Thus, investigation into the synergistic regulatory mechanisms of gas-water two-phase flow within mixed wettability fractures is warranted to uncover the actual seepage behavior of deep coal reservoirs, optimize the design of fracturing proppants, and enhance coalbed methane (CBM) recovery [12].

Gas-water two-phase flow directly determines reservoir relative permeability and CBM recovery efficiency, making it a core focus of seepage research. From the perspective of macroscopic experiments [13], Durucan et al. confirmed through multi-rank coal seepage tests that high irreducible water saturation severely inhibits gas flow, highlighting the restrictive effect of two-phase seepage mechanisms on production capacity [14]. Zhang et al. further revealed via NMR-displacement experiments that movable water saturation dominates relative permeability in low-medium rank coal, with bound water exhibiting stronger gas flow obstruction, resulting in narrow two-phase flow intervals [15]. In terms of microscopic simulations, Qian et al. established a fluid-solid coupled microscale two-phase seepage model by integrating Navier-Stokes equations, phase-field method and solid mechanics, which was verified to have a maximum relative error of only 6.19% against experimental results [16]. Wu et al. validated the feasibility of finite element simulations in pore-scale two-phase flow studies through glass etching experiments, confirming consistent laws of water displacement and saturation evolution [17]. Wettability is a key parameter regulating phase interface distribution and displacement efficiency. Chen et al. found that coal surface roughness enhances pore water retention, reducing bubble contact angle and prolonging adhesion via mixed wettability experiments [18]. Gerami et al. constructed coal-based microfluidic chips and revealed that mixed wettability dominates flow paths by regulating capillary pressure distribution [19]. However, existing experimental methods struggle to dynamically characterize cross-scale regulation of seepage paths by wettability heterogeneity, while numerical simulations mostly adopt uniform wettability parameters, failing to reflect the actual wettability differences between proppants and coal matrices. This research gap motivates the systematic investigation of synergistic regulatory mechanisms of gas-water flow in mixed wettability fractures.

In summary, the current research on gas-water flow in mixed wettability fractures is not yet comprehensive. It is difficult to directly observe the residual phase distribution and pressure field characteristics in gas-water seepage experiments. Numerical simulation methods are mostly based on uniform wettability parameters. Therefore, this study focuses on the Baicheng mining area (Xinjiang, China) as its research context. By constructing a geometric generalization model of mixed wettability fractures and combining the phase field method to track the dynamic evolution of the solid-liquid-gas three-phase interface, this work systematically reveals the regulation of proppant wettability, spatial distribution, and coal bridge contact angle on the seepage path. The research results aim to provide theoretical support for optimizing fracturing design and enhancing seepage efficiency in deep coalbed methane reservoirs, thereby promoting technological innovation in the development of clean energy from low-permeability coal reservoirs.

2  Material

2.1 Coal Sample Selection

The experimental specimens were obtained from the complete coal section in the middle of the A5 coal seam of the Lower Jurassic Taliqike Formation, Well BCS-15, Baicheng mining area, Xinjiang, at a sampling depth of 582.65 m, as illustrated in Fig. 1. Situated within the thrust-nappe belt of the southern foothills of the Tianshan Mountains, along the northern margin of the Tarim Basin, this region developed a complex tectonic framework dominated by high-angle thrust faults and associated secondary folds through multi-stage tectonic superposition during the Yanshanian to Himalayan orogenies. The interwoven network of natural fractures and bedding fractures forms a three-dimensional seepage system, establishing a typical geological prototype for the simulation of proppant accumulation law and fluid migration pathways [20]. Serving as the primary recoverable stratum in the Kubai coal field, the A5 coal seam exhibits medium-rank metamorphism characteristics with semi-bright coal predominance. The matrix properties (low porosity/low permeability) form a distinct coupling with well-developed fracture system, demonstrating high compatibility with fluid-proppant interaction mechanisms under mixed-wettability environments. The combined reservoir pressure gradient (9.50 kPa/m) and geothermal gradient (3.38°C/100 m) together constitute a stable atmospheric geothermal environment, providing a real physical boundary for numerical simulations of gas-water two-phase seepage. The samples were collected by strictly avoiding fault fracture zones and tectonically disturbed areas to ensure both the integrity and representativeness of the coal rock pore-fracture structures. Their geological parameters effectively support phase-field simulations in analyzing the mechanistic influence of quartz proppant wettability and flow velocity on gas-water migration within fractures, thereby providing a reliable experimental basis for elucidating the development law of coalbed methane under mixed-wettability conditions.

Figure 1: Sampling location of coal samples in Baicheng mining area.

2.2 In-Situ Wettability of Reservoir

To investigate the wettability evolution law of coal fracture media under in-situ temperature-pressure coupling, this study independently developed a wettability dynamic test system capable of replicating reservoir temperature and pressure environment. The device accurately simulates the conditions of in-situ pressure of 5.53 MPa and temperature of 307.84 K through pneumatic pressurization and PID temperature control module. It captures real time dynamic changes in droplet contact angles based on high-speed microscopic imaging technology, as illustrated in Fig. 2. The experimental results provide a quantitative characterization of the wettability of coal samples under various temperature and pressure combinations, offering essential data support for revealing the multiphase seepage mechanisms in reservoirs [21].

Figure 2: In-situ test system for wettability of deep reservoir. (a) Temperature-pressure coupling experimental device; (b) schematic diagram of contact angle dynamic test principle [28].

Dynamic reproduction of reservoir temperature and pressure. The coal samples were processed into cylindrical specimens with a diameter of 50 mm, and their surfaces were finely polished using sandpaper [22]. The treated coal samples were placed in a corrosion-resistant alloy high-pressure cavity (Tuochuang, China; V = 1 L, maximum pressure 40 MPa) [23]. A gas booster pump (STK, China; pressure range 0.7–40 MPa, pressure accuracy 0.5% FS) was employed to drive the multi-channel gas mixing system, which injected simulated formation gas (50.75 mol% N2, 6.47 mol% CO2, 42.72 mol% CH4) at a pressurization rate of 0.5 MPa/min until reaching 5.53 MPa. Subsequently, the closed-loop temperature control module (TES-300 A, temperature control ± 0.5 K) was activated, and the temperature was stepped up to 307.84 K, with temperature and pressure fluctuations maintained at less than 0.3% and 1.5%, respectively [24,25].

Precise transport of micro-nano droplets. Under constant temperature and pressure conditions, via a high-pressure nano-injection system (ISCO TC-100D, USA; flow rate range 0.1 μL/min-30 mL/min, accuracy 0.3% FS), deionized water (surface tension 72.1 mN/m, 25°C) was delivered to the coal samples surface through a 316L stainless steel microtube. A piezoelectric ceramic drive technology was employed to release droplets at a quasi-static flow rate of 5 nL/s, thereby avoiding the influence of inertial forces on the three-phase contact line.

Dynamic analysis of interface morphology. The dynamic spreading process of droplets on the coal surface was captured in real-time through the embedded sapphire window (transmittance > 92%) on the side wall of the cavity, utilizing a high-speed camera (M230 high-resolution digital camera, Revealer, Hefei) coupled with a stereomicroscope. After obtaining sequential images of the contact angle evolution over time, images corresponding to 8 to 15 s after the droplet’s release, which exhibited relatively stable morphology, were selected as the quantitative analysis objects.

High-pressure in situ sessile drop method was employed for contact angle quantification. The droplet images were processed using ImageJ to reconstruct three-dimensional morphology through edge detection and elliptical fitting. Based on the intersection points between the elliptical geometric model and the substrate reference line, the tangent values of the gas-liquid-solid three-phase contact angles were calculated by incorporating the Young-Laplace equation. During the experiment, a high-precision temperature and pressure control device was used to maintain the in-situ environmental parameters of the reservoir and to monitor temperature and pressure fluctuations in real time. This method captures the characteristics of the three-phase interface in situ and provides real wettability parameters, which are more accurate than atmospheric pressure tests for characterizing the coal-rock interface properties [26,27].

3  Mathematical Model Construction

3.1 Governing Equation

The phase field method treats interfacial tension and capillary forces using the continuous medium assumption, without the need for explicit tracking of the interface [29], making it suitable for capturing complex interfacial morphologies in mixed wettability fractures at the microscale. It accurately simulates capillary force gradients and the dynamic competitions of the gas-liquid-solid three-phase contact line, which are influenced by wettability differences. This provides a high-resolution numerical characterization of microscale multiphase seepage [30]. Under mixed wettability conditions, gas-water two-phase flow in coal microfractures adheres to the laws of mass conservation and momentum conservation.

The competition between capillary and viscous forces is inherently captured through the coupling of the phase-field equation (capillary/interfacial effects) and the Navier–Stokes equation (viscous effects), and is globally characterized by the capillary number Ca. We employ an improved Navier-Stokes equation to describe the flow of incompressible fluids [31,32]:

ρ(∂u∂t+u⋅∇u)=−∇p+∇⋅(η(∇u+(∇u)T))+σ∇ϕ+Fext(1)

∇⋅u=0(2)

where ρ is the fluid density, kg/m3; u is the velocity field, m/s; t is the time, s; p is the fluid pressure, Pa; η is fluid viscosity, Pa·s; σ is a parameter for controlling the surface tension of the interface; ϕ is the Phase field variables; Fext is the external volume force, N/m3.

The phase field variable ϕ [−1, 1] is introduced, and its governing equation is extended based on Cahn-Hilliard theory [33]. Eq. (4) introduces the auxiliary variable Ψ to transform the fourth-order Cahn–Hilliard equation into a system of two second-order equations (Eqs. (3) and (4)), which are solved simultaneously in the finite-element framework to track the phase interface.

∂ϕ∂t+u⋅∇ϕ=∇⋅(γλε2∇μ)(3)

Ψ=−∇⋅(ε2∇ϕ)+ϕ(ϕ2−1)+ε2cos⁡θw(x)(4)

where Ψ is the phase field auxiliary variable; γ is the mobility; μ is the chemical potential, defined as the variational derivative of the free energy functional, J/m3; λ is the mixed energy density; ε is the interface thickness parameter, m; θw(x) is a space-dependent contact angle distribution function.

Here, θw(x) is a spatially dependent contact angle distribution function that encodes the wettability heterogeneity: θw(x) = 27.3° for proppant regions and θw(x) = 75° for coal bridge regions, as determined from in-situ measurements. This spatial variation directly influences the capillary pressure field via the wall energy term in Eq. (8).

The mixing energy density k is related to the surface tension coefficient σ as: [34]

k=3ε⋅σ22(5)

The interface thickness parameter ε is generally defined as:

ε=hc2(6)

where hc is the size of the feature mesh in the region through which the interface passes.

γ=ε2(7)

In numerical implementation, the product γ·ε2 is treated as a numerical stability parameter, consistent with the phase-field literature, and does not affect the equilibrium contact angle or interfacial tension. The free energy density function quantitatively characterizes the mechanical state of the system by integrating the free energy of the body and the interface energy term [35]. The function is specifically expressed as:

f(ϕ)=ε4(ϕ2−1)2+σ2|∇ϕ|2(8)

3.2 Model Construction

3.2.1 Model Assumptions

Natural coal reservoir fracture networks exhibit a high degree of heterogeneity, including multi-scale pore throat structures, irregular filler distributions, and complex surface morphologies [36]. Direct modeling at full scale faces multiple challenges, including redundant geometric parameters and limited computing resources. Based on the percolation threshold screening and pore topological equivalence principle, the real pore throat structure is abstracted into a microscale periodic lattice model. Additionally, it is assumed that the contact interface between the proppant and the coal body is an ideal geometric shape, while the direct dynamic influence of surface roughness on wettability is disregarded. It is important to note that the wettability parameters of the proppant and coal bridge are measured under in-situ temperature and pressure conditions, which account for the effects of roughness. Consequently, this simplified assumption indirectly characterizes the influence of surface morphology through the contact angle boundary condition. This assumption achieves substantial computational complexity reduction while preserving the essential characteristics of the connectivity of the seepage channel, ensuring that the flow response law of the key pore throat structure can be accurately characterized.

The wettability of actual coal fracture fillers, including proppant, coal bridge, pulverized coal, etc., exhibits a multimodal distribution. However, the experiment needs to concentrate on the core regulation of the wettability difference between the proppant and the coal bridge in relation to two-phase flow. By assuming that the fracture filler consists solely of hydrophilic proppants and hydrophobic coal bridges, the filler is simplified into a binary proppant-coal bridge model. This simplification eliminates the interference of other secondary components and clarifies clear isolation of wettability gradient effects on capillary pressure distribution. Consequently, this establishes an idealized framework for elucidating the competition mechanism of “hydrophilic binding-hydrophobic breakthrough”.

Based on the dynamic characteristics of multiphase flow and the actual seepage conditions of deep coal reservoirs, the following flow condition constraints are proposed:

     i)  Fluid is an incompressible Newtonian fluid, and its flow process follows the Navier-Stokes equation;

    ii)  The evolution of the gas-water two-phase interface is described by the phase field method, and the interface thickness and surface tension adhere to the continuous medium assumption;

   iii)  The coal bridge is a rigid skeleton, and its pore structure does not undergo deformation or stress sensitivity during the flow process;

    iv)  The mixed wettability of the porous media wall is characterized by local contact angle distribution, while the dynamic wettability hysteresis is not considered.

This simplified binary model, while neglecting multimodal wettability distributions, clearly isolates the competing mechanisms of “hydrophilic adhesion vs. hydrophobic breakthrough.” It highlights the central regulatory role of wettability gradients in the spatial distribution of capillary pressure, providing an idealized analytical framework for understanding the evolution of air-water interfaces under heterogeneous wettability conditions.

3.2.2 Geometric Configuration

Based on the multi-physical field coupling numerical simulation method, a micro-scale digital model characterized by the porous structure of coal was constructed on the COMSOL Multiphysics platform. As illustrated in Fig. 3, the white area represents the solid matrix of coal rock, while the gray area is the microfracture space. The computational domain features a macroscopic anisotropic rectangular matrix measuring 6000 × 3080 μm, within which a periodic lattice system with a minimum characteristic spacing of lmin = 60 μm is established in the rectangular matrix [37]. A circular pore structure unit with a radius of r = 190 μm is positioned at the nodes of the system. Subsequently, the Non-uniform Boolean Subtraction was performed to exclude these geometric entities from the matrix domain, thereby constructing a rigid skeleton structure with a percolation threshold. The remaining continuous phase automatically forms the porous media seepage channels [38]. Verified through the spatial segmentation algorithm of Voronoi diagram theory, the equivalent porosity Φ of the obtained porous media model reaches 44.36%.

Figure 3: Geometric model and mesh generation of porous media.

3.2.3 Mesh Generation

Due to the significant scale variations in pore structure, the span between the maximum feature size lmax = 440√2 μm and the minimum feature size lmin = 60 μm is relatively large. Therefore, the COMSOL geometric adaptive algorithm is employed for meshing. In critical areas, such as throats that influence fluid flow, local encryption processing is applied to limit the mesh size to ≤5 μm, ensuring an accurate analysis of flow characteristics. The main pore area utilizes a relatively coarse grid to reduce the amount of calculation. The circular pore boundaries are discretized using curvature-adaptive triangular meshes to optimize the fit of complex geometric boundaries. After the above processing, the whole computational domain generates about 60,000 grid units.

To ensure the grid independence of the numerical results, this study verifies the convergence of key physical quantities based on multi-level grid density. Firstly, the local velocity distribution and the overall pressure drop in the throat area are selected as the sensitivity indexes, and the basic grid, refined grid and coarse grid are used for comparative calculation. The results show that when the throat grid size is refined from 5 to 3 μm, the maximum flow velocity deviation is 1.2%, and the global pressure drop difference is only 0.8%. When the mesh size is coarsened to 8 μm, the flow velocity and pressure drop errors increase to 4.5% and 3.7%, respectively, indicating that the scheme with the basic mesh size of 5 μm has met the convergence requirements. Further, the curvature adaptive algorithm is used to optimize the boundary mesh quality. The global unit distortion is less than 0.25, and the L2 norm error shows an exponential decay trend with mesh refinement, which verifies the balance rationality of the current mesh division between geometric resolution and computational efficiency, and can accurately characterize the multi-scale flow characteristics of the pore-throat system.

3.3 Boundary Condition

3.3.1 Model Boundary Setting

In constructing the boundary conditions of the model, different boundary types are configured to fit the actual physical scene. The left side of the model is set as a velocity inlet, featuring a 133 μm wide flow transition zone, where fluid is injected at a constant flow rate ranging from v = 1.27 × 10−5~7 × 10−1 m/s, covering the typical velocities of fluid in coal porous media. The pressure outlet boundary is adopted on the right side, maintaining a fixed static pressure of 0 Pa to simulate the pressure state under infinite boundary conditions. The non-slip symmetry condition is applied to the top and bottom boundaries of the model to ensure that the tangential component of the fluid velocity at the wall surface is zero, thereby preventing the non-physical phenomenon of fluid penetration. In the parameter settings of the phase field model, the interface mobility is set to 1 to regulate the diffusion and migration of the phase interface. The interface thickness is dynamically determined based on 50% of the maximum grid size in the computational domain, so as to optimize the computational efficiency while ensuring the computational accuracy. The solid wall boundary adopts the non-permeable condition to constrain the normal velocity component, which strictly limits the penetration behavior of the fluid at the wall surface.

To accurately represent the wettability characteristics of coal porous media surfaces, the porous media wall is designed as a mixed-wettability boundary, with the hydrophilic area representing proppant surfaces set to 25%. This ratio is designed to reflect the typical coverage scene of the local enrichment area of proppant in the fracture after hydraulic fracturing. To achieve spatial heterogeneity of wettability, we accomplish this by randomly generating seed points using a Poisson disk distribution within the coal matrix to simulate the positions of proppant particles, with an influence radius of r = 200 μm. By adjusting the surface density of the seed points, the total proportion of the wall area designated as hydrophilic in the model reaches the preset 25%, as shown in Fig. 4. This setting enables us to systematically investigate of how the differences in wettability and the spatial distribution of proppant within the coal bridges regulate two-phase flow, particularly under conditions where proppant occupies a significant proportion of the space. The wetting characteristics need to be defined on the wall surface (∂Ω), which can be achieved by introducing the surface energy density fw(ϕ) on the solid boundary and defining the boundary conditions [39]:

σn⋅∇ϕ+∂fw∂ϕ=0(9)

fw(ϕ)=−γwcos⁡θe(10)

where θe is the contact angle of the solid-liquid-gas three-phase equilibrium, which is used to define the wall wettability; γw is the wall energy coefficient, in the phase-field module as γw = σ·cos θe, where σ is the surface tension and θe is the prescribed contact angle for the wall material. n is the outward-pointing unit normal vector to the wall boundary.

Figure 4: Visualization characterization of coal fracture proppant injection and profile mixed wetting.

3.3.2 Mixed Wetting Implementation

To quantitatively characterize the regulation of wettability heterogeneity between the proppant and the coal matrix on the gas-water flow path, silica was selected as the representative hydrophilic proppant in this study. Based on in-situ temperature and pressure (5.53 MPa, 307.84 K), the surface contact angle is 27.3°, indicating strong hydrophilicity. The coal matrix (as a natural fracture component) exhibits hydrophobicity, with an in-situ contact angle of 75°. This constitutes a distinct “hydrophilic-hydrophobic” wettability gradient system. Notably, the above contact angles incorporate the comprehensive effects of surface roughness, accurately reflecting the reservoir’s real interface characteristics. The spatially heterogeneous distribution of proppants within fractures was achieved through rigorous numerical sampling: A Poisson disk distribution lattice, satisfying the minimum spacing constraint of 60 μm, was generated using the Bridson algorithm on the MATLAB platform. Circular solid elements with a radius of 190 μm in the porous medium were randomly identified. By controlling the surface density of the lattice, the hydrophilic area proportion in the model reached 27%, matching the typical scenario of local proppant enrichment after fracturing.

The circular units identified through sampling were assigned hydrophilic silica proppant properties (contact angle: 27.3°), whereas unsampled units maintained the hydrophobicity of the native coal matrix (contact angle: 75°). This setting strictly corresponds to the visual characterization in Fig. 4—the non-uniform distribution characteristics of white quartz proppant particles within the microscopic single fracture, depicted in the left diagram, are highly consistent with the Poisson disk distribution generated by the Bridson algorithm. In the right section, the randomly distributed white cylindrical units represent sampled proppants, while the black units correspond to unsampled coal matrix, clearly demonstrating the spatial gradient of wettability heterogeneity. This method not only quantitatively reproduces the mixed-wettability characteristics of fractures after proppant injection, but also more accurately captures the spatial gradients of capillary forces induced by wettability differences. It establishes a reliable computational foundation for simulating capillary-driven dynamic gas-water interface evolution processes, including liquid film formation, fragmentation, and gas-phase seepage path reconstruction.

4  Model Verification

To verify the reliability of the numerical model, this study systematically evaluated the phase-field method’s (PFM) capability in characterizing gas-liquid interface behavior under mixed wettability conditions through experimental-simulation comparisons [40]. Experiments were conducted using a self-developed in-situ temperature-pressure contact angle test device for deep reservoirs, replicating the actual temperature and pressure conditions of the reservoir (5.53 MPa, 307.84 K). This yielded contact angles of 27.3° for quartz proppants and 75° for coal bridges. The numerical model establishes boundary conditions based on the measured wettability parameters, which are modified by roughness. Although the model ignores the dynamic evolution of roughness, the measured contact angle inherently accounts for its comprehensive effect. As demonstrated in Fig. 5, the simulation exhibited a high degree of consistency with experimental observations in both interface dynamics and contact angle evolution [41], thereby verifying PFM’s reliability in characterizing wettability differences and gas-liquid displacement behavior [42].

Figure 5: Comparison of wettability experiment and simulation results under high pressure environment. (a) Quartz proppant (contact angle: 27.3°); (b) Coal bridge (contact angle: 75°).

5  Results and Analysis

5.1 Asynchronous Regulatory Mechanism of Capillary Force and Viscous Force Competition

The following results are obtained from post-processing the simulated fields. The gas volume fraction is calculated as (1 + ϕ)/2. Residual liquid clusters are identified via connected-component analysis of regions where ϕ > 0. Saturation is the volume fraction of liquid phase. The capillary number is defined as Ca = μv/σ, where v is the mean inlet velocity.

5.1.1 Low Ca Region: Capillary Force-Dominated Liquid Phase Morphology Pre-Reconstruction

Within the low capillary number region lgCa ≤ −4.32, the spatial distribution of residual liquid phase and gas phase flow paths exhibit significant progressive evolution with increasing Ca, as shown in Fig. 6. Capillary forces dominate this process, regulating the dynamics of the two-phase interface. When the Ca is extremely low (lgCa = −6.31), the gas phase is dispersed as isolated clusters within the fractured porous medium, and the residual liquid phase, under the local influence of capillary forces, can only sparsely adhere to the surface of hydrophilic proppants and cannot overcome the interfacial energy barrier formed by hydrophobic coal bridges, thus failing to form a continuous phase across pore throats; Fig. 7 shows that at this point, the residual liquid phase occupies 26 pores and 18 throats, with a strongly heterogeneous pressure field. The gas phase struggles to form effective connectivity pathways due to the dispersion of the liquid phase. As the Ca increases to lgCa = −5.56, capillary forces strengthen, driving droplets to initially coalesce on the hydrophilic support particles’ surfaces. Some throats are filled and locally connected, with the number of pores and throats increasing to 28 and 21, respectively. Continuing to increase to the critical Ca, at lgCa = −4.32, the residual liquid phase aggregates in the throats to form a continuous liquid film, and the gas phase escapes its isolated state and begins to connect. At this point, the number of pores remains at 28, while the number of throats increases to 23. The liquid phase coalesces to the saturation threshold, and the pressure field heterogeneity weakens, laying the foundation for the transition to a high-Ca region dominated by viscous forces.

Figure 6: Gas phase volume fraction and pressure distribution with low Ca. (a) Gas phase volume fraction and pressure distribution at lgCa = −6.31; (b) gas phase volume fraction and pressure distribution at lgCa = −5.56; (c) gas phase volume fraction and pressure distribution at lgCa = −4.32.

Figure 7: Residual liquid phase distribution in the throat under different low Ca.

It is important to note that the term “asynchrony” used here refers to a temporal lag and mechanistic separation between the pre-reconfiguration of liquid-phase morphology, governed by capillary forces in the low Ca region, and the transition to a flow state dominated by viscous forces in the high Ca region, occurring on either side of the critical Ca. In other words, morphological changes precede alterations in flow connectivity, reflecting the asynchronous response of competing forces under a wetting gradient. Pressure fields reflect capillary-dominated regimes: high heterogeneity at low Ca, transitioning toward more continuous gradients as liquid films coalesce.

The aforementioned phenomenon stems from the weak viscous force effect in low Ca regions, with capillary forces becoming the primary driving force for the evolution of liquid phase morphology. At extremely low Ca values where lgCa = −6.31, the strength of capillary forces is limited, capable only of driving droplets to undergo local coalescence within a single hydrophilic support agent region. They cannot overcome the energy barrier formed by the hydrophobic coal bridge interface, resulting in the liquid phase remaining highly dispersed with poor throat connectivity. As Ca values increase, enhanced capillary forces gradually gain the ability to overcome hydrophobic resistance, driving dispersed droplets to migrate across hydrophilic-hydrophobic boundaries and coalesce, ultimately forming a continuous liquid film network in the throat regions. This transition from a dispersed state to a continuous state is essentially a “morphological pre-reconstruction” driven by capillary forces. This reconfiguration process optimizes pore throat connectivity, expands effective flow channels, and lays a critical morphological foundation for efficient gas phase connectivity and flow state transition after the intervention of viscous forces in high Ca regions. At this stage, only liquid-phase morphological reconstruction occurs, with residual liquid transitioning from a dispersed state to a preliminary continuous liquid film, without triggering gas-phase flow state transition. This reflects the cross-stage asynchrony of “low Ca region morphological pre-reconstruction preceding high Ca region flow state transition”.

5.1.2 High Ca Region: Viscosity-Driven Flow Regime Transition and Morphological Secondary Reconstruction

In the low Ca region, where lgCa = −4.32, the continuous liquid film has completely covered the throat, with 23 throat channels occupied. This provides the “continuous phase foundation” for the viscosity-driven liquid film rupture and rapid release of throat channels in the high Ca region—compared to discrete liquid phases, the shear rupture efficiency of a continuous liquid film is higher, enabling the number of throat channels occupied by the liquid phase to directly decrease from 23 to 18 when the initial lgCa in the high Ca region is −3.64, thereby accelerating the release of gas phase channels. Upon entering the high Ca region (lgCa ≥ −3.64), the gas phase seepage pathways and residual liquid phase distribution undergo significant transformations, as shown in Fig. 8. When lgCa = −3.64, viscous forces begin to tear apart the continuous liquid film network, and the gas phase begins to form seepage pathways, but the throat channels still retain continuous liquid films. The pressure field exhibits a directed gradient with local distortions, with 24 pore channels and 18 throat channels at this point. As lgCa increases to −2.43, viscous forces completely dominate the liquid film breakage process, with liquid films shearing and breaking at the pore throat interface: larger-sized pore liquid films preferentially break into displaceable small droplets, reducing the number of pores occupied by the liquid phase from 24 to 21; The throat dimensions are small, with a minimum feature spacing of 60 μm. The broken droplets easily anchor in the corners of the throats covered by hydrophilic proppants, and the number of throats decreases only from 18 to 15. The gas phase connectivity increases, and the pressure field gradient becomes regular. Under the highest Ca condition (lgCa = −1.56), strongly anchored droplets within the pore channels continue to be displaced, reducing the number of pore channels to 13; When throat droplets reach the “minimum anchoring threshold” and the hydrophilic region can no longer be reduced, the number stabilizes at 15, as illustrated in Fig. 9. At this point, the gas phase dominates the flow path, with residual liquid phase droplets anchored in isolation. The pressure field is uniformly oriented, and flow resistance is minimized. The reduction in the number of liquid phase pores coupled with the increase in effective gas phase channels reflects improved gas phase flow efficiency. The pressure field evolves from highly heterogeneous (Fig. 6a) to more uniform (Fig. 8c) as viscous forces become dominant, indicating reduced capillary barriers.

Figure 8: Gas volume fraction and pressure distribution with high Ca. (a) Gas phase volume fraction and pressure distribution at lgCa = −3.64; (b) gas phase volume fraction and pressure distribution at lgCa = −2.43; (c) gas phase volume fraction and pressure distribution at lgCa = −1.56.

Figure 9: Residual liquid phase distribution in the throat under different Ca.

The aforementioned evolutionary process is primarily driven by viscous forces that overcome capillary forces. Enhanced viscous shear forces first initiate localized rupture of the continuous liquid film network near the critical point lgCa = −3.64, with the initial appearance of secondary morphological reconstruction, thereby creating initial channels for gas phase seepage. As the flow velocity continues to increase to lgCa = −2.43, viscous forces completely overcome capillary constraints, completely breaking the remaining liquid film into discrete droplets, completing the secondary morphological reconstruction, significantly reducing the blocking effect on the flow channel, and markedly enhancing gas phase transport efficiency. At the high Ca end with lgCa = −1.56, viscous forces dominate, continuously displacing and compressing the discrete droplets, confining them to hydrophilic regions, and establishing a gas phase flow network with no obstruction. At this stage, the viscous force first drives the secondary morphological reconstruction, that is, the liquid film breaking, and then drives the flow pattern change, that is, the gas phase connectivity, and triggers the ‘secondary morphological reconstruction’ of the liquid phase from the continuous film to the discrete droplets. Under the final maximum Ca condition, the discrete droplet morphology is completely synchronized with the flow state dominated by viscous force, which marks the end of “morphology-flow asynchrony”.

5.1.3 Macroscopic Phase Evolution of Residual Liquid Phase Clusters

An increase in the Ca drives significant morphological and distributional evolution in residual liquid phase clusters, revealing deep non-monotonic regulatory patterns. As shown in Fig. 10, as lgCa increases from −6.31 to −1.56, the number of clusters first decreases sharply to a minimum due to capillary force-dominated aggregation across hydrophilic regions, then increases significantly due to the liquid film rupture effect caused by enhanced viscous forces. The maximum characteristic radius of the clusters, which evolves synchronously, reaches its peak when capillary forces are at their maximum, as calculated by the Eqs. (11) and (12) to quantitatively calculate the equivalent radius of all residual liquid phase clusters, with the maximum value taken as the maximum characteristic radius. This value reaches a peak of 0.494 mm at lgCa = −5.56 and subsequently contracts to a minimum of 0.248 mm as viscous shear forces intensify. This inverse relationship between quantity and size reveals a dynamic phase transition mechanism. This mechanism arises from the competition between capillary and viscous forces: in the low Ca region, capillary forces maximize liquid phase continuity, optimize pore throat connectivity, and create pre-connected pathways for gas phase flow; whereas in the high Ca region, viscous forces tear open liquid films to release flow channels, yet cause discrete droplets to anchor at hydrophilic interfaces, forming new microscopic flow barriers.

ri=Aiπ(11)

rmax=max{r1,r2,⋯,rN}(12)

where A is a single cluster area, mm2; and r is the radius of the cluster, mm. N is the total number of residual liquid phase clusters under the current Ca.

Figure 10: The evolution of the number of residual liquid phase clusters and the maximum characteristic radius with the Ca.

When lgCa = −4.32, it represents the critical Ca, marking the transition from capillary force dominance to viscous force dominance: below the capillary force coalescence threshold, the liquid phase morphology undergoes pre-reconfiguration, forming a continuous liquid film covering the throat; while the viscous force fragmentation threshold triggers liquid film instability, driving the secondary dispersion of the residual phase. This non-monotonic phase evolution stems from the differential regulation of the two-phase interface by the spatial gradient of wetting heterogeneity—hydrophilic proppants attract liquid phase coalescence, while hydrophobic coal bridges promote shear fragmentation. Ultimately, the coupled effects of a sharp increase in cluster number and a sudden decrease in characteristic radius directly lead to fundamental changes in the occupancy range and distribution morphology of the liquid phase in the fracture space. This microscopic cluster “coalescence-breakup” phase evolution is precisely the core microscopic trigger for the non-monotonic response of macroscopic liquid phase saturation, and the two exhibit a clear intrinsic coupling relationship with changes in Ca.

5.2 Evolution Mechanism of Residual Liquid Phase Saturation

5.2.1 Non-Monotonic Saturation Response Dominated by Ca

As a macroscopic quantitative representation of the continuous or discrete microstructural morphology of residual liquid phase clusters, the evolution of liquid phase saturation with Ca is completely coupled with the cluster phase transition process from aggregation to fragmentation revealed in this study, and therefore exhibits significant non-monotonicity, as shown in Fig. 11. The critical Ca serves as the threshold for the transition of the dominant force mechanism, after which liquid phase saturation reaches a peak and then transitions to a linear decline. In the low Ca region (lgCa ≤ −4.32), capillary forces dominate liquid phase distribution, driving discrete water phases to aggregate, coalesce, and extend on hydrophilic support surfaces to form a more continuous liquid film network. This process aligns with the coalescence characteristics of reduced numbers and increased radii in the cluster phase, ultimately expanding the liquid phase’s occupancy in fracture spaces, with saturation rising continuously from 29.38% to a peak of 34.93%. While this process enhances liquid phase continuity, it also increases gas phase flow resistance due to the liquid film’s blocking effect on pore throats.

Figure 11: Non-monotonic evolution of liquid saturation with Ca.

When lgCa > −4.32, viscous forces surpass capillary forces to become dominant, with their shear action continuously tearing apart and fracturing the continuous liquid film network, synchronized with the fracturing characteristics of a sharp increase in the number of clusters and a sudden decrease in their radii; The fragmented microdroplets are displaced and anchored in the downstream hydrophilic region, causing the effective occupancy range of the liquid phase in the fissure space to significantly contract, resulting in a linear decrease in saturation to 20.60% and effectively releasing the gas phase flow channel. This two-stage response of first increasing then decreasing is the core manifestation of the dynamic competition between capillary forces and viscous forces in mixed-wetting fissures under a wetting gradient background. The critical Ca not only marks the transition point of the dominant force mechanism but also defines the transformation of liquid phase morphology from capillary coalescence to viscous fragmentation; The increase in saturation in the low Ca zone is driven by capillary-driven liquid phase reorganization and extension, while the linear decay in the high Ca zone quantitatively characterizes the separation efficiency of viscous forces. The final residual saturation is constrained by the capture effect of hydrophilic proppants on liquid droplets. This mechanism provides a theoretical basis for controlling Ca near the critical value through optimized injection strategies, efficiently fracturing confined liquid phases, and enhancing gas phase migration efficiency.

5.2.2 Configuration-Dependent Phase Transition Kinetics Controlled by Displacement Time

Under different Ca s, the dynamic process of gas phase displacement within mixed-wetted coal fractures exhibits significant differences. To uniformly quantify the displacement sequence and precisely compare the evolution characteristics of the residual liquid phase under different Cas, this study defines the breakthrough time tB as the moment when the vapor phase first reaches the model’s exit cross-section. The ratio of the displacement time at any given time to the breakthrough time (ti/tB) as the standardized intrusion time. Fig. 12 shows the variation in residual liquid phase saturation with this standardized invasion time during the gas phase displacement process, along with fitting results for different Cas. Based on the corresponding data, it can be observed that the evolution of residual liquid phase saturation with displacement time exhibits significant differences under different Ca, and is deeply coupled with the configuration phase transition process dominated by the competition between capillary forces and viscous forces. In the low Ca range where lgCa = −6.31, the residual liquid phase saturation decreases in stages as the standardized invasion time increases, dropping steadily from 100% to 29.38% in the early stage and then stabilizing with only a slight decrease to 27.42% in the later stage; At the critical Ca lgCa = −4.32, the saturation exhibits a “slow decrease–stable amplitude” transition characteristic, decreasing slowly from 100% to 34.93% in the early stage, then slowly decreasing to 30.12% in the later stage before stabilizing; Under high Ca conditions (lgCa = −1.56), the saturation evolution exhibits a ‘rapid decrease–early stabilization’ characteristic, with the saturation rapidly decreasing from 100% to 20.60% in the early stage and then stabilizing prematurely at 20.60%—i.e., the droplets have been displaced to the downstream hydrophilic proppant region and reached the minimum anchoring threshold that is the adsorption limit corresponding to the 25% hydrophilic area ratio—making further migration impossible.

Figure 12: Changes in residual liquid phase saturation under low, critical, and high Ca.

The core driving mechanism behind these differentiated evolution characteristics lies in the competitive relationship between capillary forces and viscous forces determined by Ca. Within the low Ca range, capillary forces dominate, only driving discrete droplets to locally coalesce along the hydrophilic support surface. The energy barrier formed by hydrophobic coal bridges restricts further migration and fragmentation of the liquid film, leading to incomplete pre-reconstruction of the liquid phase morphology. The regulatory effect of displacement time is limited to the initial stage of liquid phase contraction; The critical Ca corresponds to the transition point where capillary forces shift dominance to viscous forces. During the early displacement phase, capillary forces still dominate droplet coalescence to form a continuous liquid film covering the throat, resulting in a gradual decrease in saturation rate. In the later phase, viscous forces initiate local fragmentation of the liquid film but have not yet fully overcome the anchoring effect of capillary forces, leading to a dynamic equilibrium between liquid film fragmentation and droplet re-coalescence, with the configurational phase transition entering a transitional state; At high Ca, viscous forces absolutely dominate, with strong viscous shear forces rapidly tearing the liquid film and breaking it into small droplets that can be displaced. These droplets are displaced to the downstream hydrophilic support agent region and anchored, causing the liquid phase to rapidly reach the “minimum anchoring threshold,” achieving complete morphological secondary restructuring. Therefore, the saturation rate decreases rapidly and enters the stable stage prematurely. This mechanism reveals the intrinsic logic of how displacement time and Ca synergistically regulate liquid phase configuration phase transitions, providing a critical foundation for understanding the dynamic evolution of gas-water two-phase systems within mixed-wetting fractures.

6  Conclusions

This study developed a hybrid wettability fracture dynamics model by coupling the Navier-Stokes equations using the phase-field method, combined with in-situ temperature-pressure wettability test data from the Baicheng mining area, to systematically reveal the critical control mechanism of support agent wettability heterogeneity on gas-water two-phase flow paths in deep coal reservoirs. The main conclusions are as follows:

(1)   The critical value lgCa = −4.32 marks the transition point where the dominant force governing flow state changes. In the low Ca region, capillary forces drive droplet coalescence into a continuous liquid film, causing liquid phase saturation to peak at 34.93% and inhibiting gas phase connectivity. In the high Ca region, viscous forces dominate, causing liquid film fragmentation, with saturation linearly decreasing to 20.60%, significantly enhancing gas phase flow efficiency.

(2)   Hydrophilic proppants induce capillary forces to drive liquid phase coalescence, resulting in a sharp decrease in cluster numbers and an increase in maximum radius to 0.488 mm in the low Ca region; hydrophobic coal bridges synergize with viscous forces to trigger liquid film breakage, leading to a sharp increase in cluster numbers and a contraction in radius to 0.248 mm in the high Ca region. The wetting gradient dominates the asynchronous phase transition of the flow path.

(3)   The standardized invasion time ti/tB indicates: In the low Ca zone, displacement only regulates the initial contraction of the liquid phase, with saturation gradually decreasing in the early stage. In the later stage, due to capillary anchoring effects, further driving of liquid phase migration becomes difficult, ultimately stabilizing at 27.42%; In the high Ca zone, viscous forces dominate the rapid rupture of the liquid film, causing the liquid phase to quickly reach the minimum anchoring threshold. The saturation stabilizes prematurely at 20.60% and no longer changes with extended displacement time.

(4)   By actively controlling the gas injection rate to approach the critical Ca range, liquid-phase fracturing can be efficiently initiated. Proppant wettability design should adhere to the principle of “hydrophilic-hydrophobic dynamic equilibrium.” In fracturing design, surface modification can regulate the proportion of proppant hydrophilicity. When combined with injection rate optimization, this approach enables synergistic pre-reconstruction of the liquid-phase film and secondary fractures, thereby maximizing the development of gas-phase flow channels.

This study establishes a dynamic model for mixed-wettability fractures under coupled temperature-pressure conditions. Using random wettability point mapping, we reveal for the first time how the critical Ca regulates flow paths. The study elucidates how wettability heterogeneity drives the asynchronous phase transition of liquid “coalescence-breakage” through the competition between capillary and viscous forces, directly controlling the dynamic reconstruction of gas-phase flow channels. This research provides a theoretical benchmark and engineering optimization pathway for the efficient development of deep coalbed methane through critical Ca regulation.

Acknowledgement: Thanks to all colleagues in the research group for their helpful discussions and suggestions during the manuscript preparation.

Funding Statement: This study was funded by the Xinjiang “Two Zones” Science and Technology Development Program: “Occurrence Mechanism and Key Development Technologies of Deep Coalbed Methane in Xinjiang” (Grant No. 2023LQ01005).

Author Contributions: Zeliang Liang: Writing, Conceptualization-Original Draft and Resources. Jia Tan: Methodology and Supervision-Review & Editing. Jiachao She: Conceptualization and Resources. Fengnian Wang: Supervision and Editing. Yi Wang: Supervision and Resources. Donghao Li: Data Curation and Validation. Rugang Duan: Formal Analysis and Investigation. Haotian Chu: Visualization and Writing-Review & Editing. All authors reviewed and approved the final version of the manuscript.

Availability of Data and Materials: The data that support the findings of this study are available from the corresponding author upon reasonable request.

Ethics Approval: This study involves only secondary analysis of publicly available data without any personal identifiable information, and thus does not require ethical approval from an ethics committee.

Conflicts of Interest: The authors declare no conflicts of interest.

References

1. Zhou S, Liu D, Cai Y, Yao Y. Fractal characterization of pore-fracture in low-rank coals using a low-field NMR relaxation method. Fuel. 2016;181(4):218–26. doi:10.1016/j.fuel.2016.04.119. [Google Scholar] [CrossRef]

2. Ramandi HL, Mostaghimi P, Armstrong RT, Saadatfar M, Pinczewski WV. Porosity and permeability characterization of coal: a micro-computed tomography study. Int J Coal Geol. 2016;154:57–68. doi:10.1016/j.coal.2015.10.001. [Google Scholar] [CrossRef]

3. Heo W, Lee W, Lee DS. Hydraulic fracturing design for coalbed methane in Barito Basin, Indonesia. Geosyst Eng. 2015;18(3):151–62. doi:10.1080/12269328.2015.1023468. [Google Scholar] [CrossRef]

4. Lu Y, Meng Z, Su X, Yu Y. Experimental study of dynamic permeability changes in coals of various ranks during hydraulic fracturing. Nat Resour Res. 2022;31(6):3253–72. doi:10.1007/s11053-022-10115-x. [Google Scholar] [CrossRef]

5. Liu D, Zhao B, Cui J, Cai Y, Sun F, Wang B, et al. Pore-fissure compressibility and structural dynamic evolution of coal reservoir under confining pressure. Energy Fuels. 2024;38(20):19509–24. doi:10.1021/acs.energyfuels.4c03072. [Google Scholar] [CrossRef]

6. Liu J, Ding Y, Xue F, Wei J, Lin H, Dai H. Confinement pressure effect and influence mechanism of water injection-induced slip of shale fracture. Eng Geol. 2025;352(6318):108061. doi:10.1016/j.enggeo.2025.108061. [Google Scholar] [CrossRef]

7. Gao L, Zhao H, Liu C. Study on the influence of true triaxial principal stress loading and unloading on gas permeability characteristics of coal seams containing gangue. Results Eng. 2025;27:105667. doi:10.1016/j.rineng.2025.105667. [Google Scholar] [CrossRef]

8. Meng Z, Li G. Experimental research on the permeability of high-rank coal under a varying stress and its influencing factors. Eng Geol. 2013;162(10):108–17. doi:10.1016/j.enggeo.2013.04.013. [Google Scholar] [CrossRef]

9. Barbati AC, Desroches J, Robisson A, McKinley GH. Complex fluids and hydraulic fracturing. Annu Rev Chem Biomol Eng. 2016;7(1):415–53. doi:10.1146/annurev-chembioeng-080615-033630. [Google Scholar] [PubMed] [CrossRef]

10. Shu Y, Sang S, Zhou X, Zhao F. A coupled hydraulic-mechanical model with two-phase flow for fracturing development of undersaturated coalbed methane reservoirs considering permeability velocity-sensitive damage. Nat Resour Res. 2023;32(5):2053–76. doi:10.1007/s11053-023-10220-5. [Google Scholar] [CrossRef]

11. Huang F, Dong C, Shang X, You Z. Effects of proppant wettability and size on transport and retention of coal fines in saturated proppant packs: experimental and theoretical studies. Energy Fuels. 2021;35(15):11976–91. doi:10.1021/acs.energyfuels.1c01418.s001. [Google Scholar] [CrossRef]

12. Zhang JC, Chen SJ, Ji ST, Feng F, Aleksey AK, Zaki MM. Imaging dynamic water invasion behavior in microfractures based on microfluidics. J Cent South Univ. 2022;29(12):3986–4001. doi:10.1007/s11771-022-5202-7. [Google Scholar] [CrossRef]

13. Chen S, Zhang J, Mohammed MZ, Li F, Yan Z, Ding YS. Seepage characteristics of mixed-wettability porous media on the phase-field model. ACS Omega. 2022;7(34):30104–12. doi:10.1021/acsomega.2c03143. [Google Scholar] [PubMed] [CrossRef]

14. Durucan S, Ahsan M, Shi JQ, Syed A, Korre A. Two phase relative permeabilities for gas and water in selected European coals. Fuel. 2014;134(1):226–36. doi:10.1016/j.fuel.2014.05.040. [Google Scholar] [CrossRef]

15. Zhang Z, Yan D, Yang S, Zhuang X, Li G, Wang G, et al. Experimental studies on the movable-water saturations of different-scale pores and relative permeability of low-medium rank coals from the Southern Junggar Basin. J Nat Gas Sci Eng. 2020;83(1–4):103585. doi:10.1016/j.jngse.2020.103585. [Google Scholar] [CrossRef]

16. Qian C, Xie Y, Zhang X, Zhou R, Mou B. Study on numerical simulation of gas-water two-phase micro-seepage considering fluid-solid coupling in the cleats of coal rocks. Energies. 2024;17(4):928. doi:10.3390/en17040928. [Google Scholar] [CrossRef]

17. Wu F, Yao C, Cong L, Xi Y. Pore-scale gas-water flow in rock: visualization experiment and simulation. Open Geosci. 2020;12(1):532–46. doi:10.1515/geo-2020-0105. [Google Scholar] [CrossRef]

18. Chen Y, Xia W, Xie G. Contact angle and induction time of air bubble on flat coal surface of different roughness. Fuel. 2018;222:35–41. doi:10.1016/j.fuel.2018.02.140. [Google Scholar] [CrossRef]

19. Gerami A, Armstrong RT, Johnston B, Warkiani ME, Mosavat N, Mostaghimi P. Coal-on-a-chip: visualizing flow in coal fractures. Energy Fuels. 2017;31(10):10393–403. doi:10.1021/acs.energyfuels.7b01046. [Google Scholar] [CrossRef]

20. Guo Z, Dong C, Wang B, Li L, Pang Z, Chen X, et al. Geological characteristics of Jurassic coal rock gas and evaluation of favorable zones in the northern structural belt of the Kuqa Depression, Tarim Basin, China. J Nat Gas Geosci. 2025;10(4):219–38. doi:10.1016/j.jnggs.2025.07.002. [Google Scholar] [CrossRef]

21. Zheng S, Yao Y, Elsworth D, Liu D, Cai Y. Dynamic fluid interactions during CO2-ECBM and CO2 sequestration in coal seams. Part 2: CO2-H2O wettability. Fuel. 2020;279(7):118560. doi:10.1016/j.fuel.2020.118560. [Google Scholar] [CrossRef]

22. Siemons N, Bruining H, Castelijns H, Wolf KH. Pressure dependence of the contact angle in a CO2–H2O-coal system. J Colloid Interface Sci. 2006;297(2):755–61. doi:10.1016/j.jcis.2005.11.047. [Google Scholar] [PubMed] [CrossRef]

23. Iglauer S, Pentland CH, Busch A. CO2 wettability of seal and reservoir rocks and the implications for carbon geo-sequestration. Water Resour Res. 2015;51(1):729–74. doi:10.1002/2014WR015553. [Google Scholar] [CrossRef]

24. Zhang J, Lv D, Zhang JJ, Zhang X, Yu H, Li X, et al. Microfluidic diagnostics: evolution of gas critical paths based on dynamic alterations of cleats wettability. Chem Eng J. 2025;515:163662. doi:10.1016/j.cej.2025.163662. [Google Scholar] [CrossRef]

25. Wang K, Du F, Wang G. Investigation of gas pressure and temperature effects on the permeability and steady-state time of Chinese anthracite coal: an experimental study. J Nat Gas Sci Eng. 2017;40(6):179–88. doi:10.1016/j.jngse.2017.02.014. [Google Scholar] [CrossRef]

26. Meng Y, Li Z, Tang S, Chen W. Experimental study on gas-water seepage in high-rank coal and its research significance. Nat Resour Res. 2023;32(6):2867–81. doi:10.1007/s11053-023-10278-1. [Google Scholar] [CrossRef]

27. Zhang J, Lv D, Zhang JJ, Wang F, Yin D, Yu H. Dual-scale insights of two-phase flow in inter-cleats based on microfluidics: interface jumps and energy dissipation. Int J Min Sci Technol. 2025;35(3):451–65. doi:10.1016/j.ijmst.2025.01.010. [Google Scholar] [CrossRef]

28. Gang S, Kang I, Lee J, Jung J. Effect of surface-treated silica nanofluid on supercritical CO2 injection efficiency: application to deep saline aquifers. J CO2 Util. 2025;97(48):103103. doi:10.1016/j.jcou.2025.103103. [Google Scholar] [CrossRef]

29. Chen Y, Deng Z. Hydrodynamics of a droplet passing through a microfluidic T-junction. J Fluid Mech. 2017;819:401–34. doi:10.1017/jfm.2017.181. [Google Scholar] [CrossRef]

30. Santra S, Mandal S, Chakraborty S. Phase-field modeling of multicomponent and multiphase flows in microfluidic systems: a review. Int J Numer Meth Heat Fluid Flow. 2021;31(10):3089–131. doi:10.1108/hff-01-2020-0001. [Google Scholar] [CrossRef]

31. Guo Y, Zhang L, Zhu G, Yao J, Sun H, Song W, et al. A pore-scale investigation of residual oil distributions and enhanced oil recovery methods. Energies. 2019;12(19):3732. doi:10.3390/en12193732. [Google Scholar] [CrossRef]

32. Bai T, Chen Z, Aminossadati SM, Li L, Liu J, Lu H. Dimensional analysis and prediction of coal fines generation under two-phase flow conditions. Fuel. 2017;194(Part 2):460–79. doi:10.1016/j.fuel.2017.01.040. [Google Scholar] [CrossRef]

33. Zhou S, Zhuang X, Rabczuk T. A phase-field modeling approach of fracture propagation in poroelastic media. Eng Geol. 2018;240(1):189–203. doi:10.1016/j.enggeo.2018.04.008. [Google Scholar] [CrossRef]

34. Tardy E, Poirot-Crouvezier JP, Schott P, Morel C, Serre G, Bultel Y. Investigation of liquid water heterogeneities in large area proton exchange membrane fuel cells using a Darcy two-phase flow model in a multiphysics code. Int J Hydrogen Energy. 2022;47(91):38721–35. doi:10.1016/j.ijhydene.2022.09.039. [Google Scholar] [CrossRef]

35. Suzuki T, Ode M, Kim SG, Kim WT. Phase-field model of dendritic growth. J Cryst Growth. 2002;237:125–31. doi:10.1016/S0022-0248(01)01891-7. [Google Scholar] [CrossRef]

36. Hao J, Shu L, Huo Z, Fan Y, Wu S, Li Y, et al. Quantitative characterization of pore-fracture structure of medium and high-rank coal based on micro-CT technology. Int J Coal Prep Util. 2024;44(3):358–75. doi:10.1080/19392699.2023.2190102. [Google Scholar] [CrossRef]

37. Schmuck M, Pradas M, Pavliotis GA, Kalliadasis S. Derivation of effective macroscopic Stokes-Cahn–Hilliard equations for periodic immiscible flows in porous media. Nonlinearity. 2013;26(12):3259–77. doi:10.1088/0951-7715/26/12/3259. [Google Scholar] [CrossRef]

38. Liu S, Wang L, Jiang Y, Wang W, Yu M, Li H, et al. Random pore-network development methodology based on Voronoi and Delaunay tessellations for residual coal under axial stress. Fuel. 2023;353(4):129267. doi:10.1016/j.fuel.2023.129267. [Google Scholar] [CrossRef]

39. Chu H, Zhang J, Cao Y, Li X, Yu H, Li F, et al. Water invasion and residual gas distribution in partially filled fractures via phase-field method. Phys Fluids. 2024;36(11):112032. doi:10.1063/5.0236093. [Google Scholar] [CrossRef]

40. Zhang J, Lv D, Yin D, Zhang X, Li X, Fan K. Gas recovery and flowback in trans-coal-limestone fracture: an in-situ wettability microscale visualization insight. Gas Sci Eng. 2025;142(2):205707. doi:10.1016/j.jgsce.2025.205707. [Google Scholar] [CrossRef]

41. Chen S, Zhang J, Yin D, Li F, Lu J, Zhu P. Visualizing experimental investigation on gas-liquid replacements in a microcleat model using the reconstruction method. Deep Undergr Sci Eng. 2023;2(3):295–303. doi:10.1002/dug2.12044. [Google Scholar] [CrossRef]

42. Heider Y, Sun W. A phase field framework for capillary-induced fracture in unsaturated porous media: drying-induced vs. hydraulic cracking. Comput Meth Appl Mech Eng. 2020;359(6):112647. doi:10.1016/j.cma.2019.112647. [Google Scholar] [CrossRef]