A Pore-Scale Simulation of the Effect of Heterogeneity on Underground Hydrogen Storage

Abstract

Using underground hydrogen storage technology has been recognized as an effective way to store hydrogen on a large scale, yet the physical mechanisms of hydrogen flow in porous media remain complex and challenging. Studying the heterogeneity of pore structures is crucial to enhance the efficiency of hydrogen storage. In order to better understand the pore-scale behavior of hydrogen in underground heterogeneous porous structures, this paper investigates the effects of wettability, pore–throat ratio, and pore structure heterogeneity on the behavior of the two-phase H₂–brine flow using pore-scale simulations. The results show that the complex interactions between wettability, heterogeneity, and pore geometry play a crucial role in controlling the repulsion pattern. The flow of H₂ is more obstructed in the region of the low pore–throat ratio, and the obstructive effect is more obvious when adjacent to the region of the high pore–throat ratio than that when adjacent to the region of the medium pore–throat ratio. In high-pore–throat ratio structures, the interfacial velocity changes abruptly as it passes through a wide pore and adjacent narrower throat. Interfacial velocities at the local pore scale may increase by several orders of magnitude, leading to non-negligible viscous flow effects. It is observed that an increase in the pore–throat ratio from 6.35 (low pore–throat ratio) to 12.12 (medium pore–throat ratio) promotes H₂ flow, while an increase from 12.12 (medium pore–throat ratio) to 23.67 (high pore–throat ratio) negatively affects H₂ flow. Insights are provided for understanding the role of the heterogeneity of pore structures in H₂–brine two-phase flow during underground hydrogen storage.

1. Introduction

Hydrogen is gradually becoming one of the most important energy storage carriers for replacing fossil fuels as the energy system transitions to cleaner energy sources [1,2]. In the process of generating electricity from renewable energy sources (e.g., wind and solar), there is an imbalance between supply and demand. The surplus power produced especially during the peak energy period needs to be stored in an effective way for energy storage, and hydrogen prepared by electrolyzing water is becoming an ideal medium for energy storage. Utilizing hydrogen to store surplus renewable energy and releasing the energy again during peak demand periods not only balance supply and demand but also provide a stable energy supply for future clean energy systems [3,4,5,6]. Using underground hydrogen storage (UHS) technology is an effective way to store hydrogen on a large scale [7,8]. As a widely distributed underground porous medium, saline aquifers have enormous storage potential [9,10,11]. Their geological closure and natural pressure provide the necessary conditions for long-term hydrogen storage and reduce the risk of hydrogen leakage. However, the migration and storage of hydrogen in underground porous media still face many technical challenges; the high diffusivity, low density, and low viscosity properties of hydrogen make it especially difficult to effectively manage its flow behavior in the reservoir [12,13]. Therefore, further research on the seepage pattern of hydrogen in underground porous media and the related physical properties is the key to promoting the mature application of hydrogen underground storage technology.

When hydrogen is stored in porous media such as saline aquifers, the flow and storage characteristics of hydrogen are affected by a combination of multiple factors. These factors include the pore structure of the underground reservoir, the pressure and temperature conditions of the reservoir, the wettability of the rock, and the interaction between hydrogen and brine [14,15,16,17]. The pore structure of underground reservoirs is a heterogeneous system consisting of pore throats of different scales and shapes. This complex structure determines the flow paths and flow efficiency of hydrogen, and yet the mechanism governing the formation of preferential flow paths as well as the factors influencing them are still not fully understood [18,19,20]. To further advance the understanding and optimization of underground hydrogen storage (UHS), it is critical to explore the pore-scale mechanisms that control the behavior of hydrogen in complex porous media such as saline aquifers. Current studies have shown that parameters such as the pore–throat ratio, rock wettability, and interfacial tension significantly affect hydrogen seepage pathways, saturation distribution, and breakthrough behavior in reservoirs [10,21,22]. Further simulation and experimental studies at the microscale are essential, focusing on how these factors interact to form preferential flow paths and influence hydrogen retention in heterogeneous reservoirs. This will help to develop more efficient storage strategies and improve the overall performance of underground hydrogen storage systems. Current research on the effect of heterogeneity on underground hydrogen storage has mainly focused on the reservoir scale and core scale, with studies focusing on the characterization of reservoir heterogeneity through permeability and wettability [18,23,24,25]. However, research on how pore-scale heterogeneity, especially pore structure heterogeneity, affects the distribution of hydrogen in reservoirs is relatively lacking.

Research on the heterogeneity of pore structures is crucial to enhance hydrogen storage efficiency, which affects the flow path of hydrogen in reservoirs and the final storage efficiency [26,27]. The study of the microscopic pore scale is one of the important means for understanding the flow behavior of hydrogen in reservoirs. Pore-scale modeling provides detailed insights into how factors such as pore heterogeneity and pore throat structure affect hydrogen storage efficiency and flow dynamics. Researchers can observe the flow behavior of hydrogen in porous media in a more detailed way at the microscopic scale. Lysyy et al. [28] used an irregular pore model of natural sandstone to observe that H₂ saturation after primary drainage increased with capillary number, and the residual trapped H₂ showed a good ability to reconnect near the large pores of a wide pore throat. The effect of capillary number on the preferential flow path of hydrogen and H₂ saturation after drainage was also investigated [27]. Jangda et al. [29] performed H₂ flow get visualization experiments on porous rocks with different pores and throats and showed that capillary pressure effects complicate fluid replacement when capillary forces dominate, especially near heterogeneous boundaries, and significantly affect H₂ saturation. Bagheri et al. [30] investigated the effect of different flow rates, operating pressures, and several cycles on underground hydrogen storage using pore-scale simulations. In another paper, the effect of pore heterogeneity (e.g., curvature, coordination number, and pore geometry) on the flow characteristics was investigated, and it was observed that an increase in heterogeneity decreased the viscous fingering during the drainage period [31]. Based on the current study, the importance of heterogeneity in pore structures for hydrogen storage is evident. These studies provide key insights into the microscale for understanding and optimizing hydrogen storage, especially for reducing residual hydrogen capture and improving hydrogen recovery. Although existing experiments and simulations have revealed some phenomena, a limited number of experimental and simulation studies have been conducted on the pore-scale behavior of hydrogen in underground porous structures, and the physical mechanisms of hydrogen flow in complex porous media remain complex and challenging.

In this study, three basic pore structures with different pore–throat ratios were designed, and four heterogeneous pore models with different pore throat arrangement characteristics were constructed by combining these basic pore structures. Based on the OpenFOAM (version 9) platform, the two-phase flow process of H₂–brine in the four heterogeneous pore media was simulated by the Volume of Fluid (VOF) method. We investigated the differences in the H₂–brine flow distribution under different heterogeneous structures and wettability conditions. The effects of different pore–throat ratios and pore arrangement structures on H₂ flow were studied, especially the differences in preferential H₂ flow paths. This study provides insights into understanding the role of pore structure heterogeneity in H₂–brine two-phase flow during underground hydrogen storage.

2. Methodology

2.1. Geometry Models

In this study, the geometric models were designed to simulate pore media with different pore–throat ratios (PTRs). We drew three sets of basic geometrical models based on AutoCAD^® software (version 2020). Figure 1a–c represent the structure of high, medium, and low pore–throat ratios, respectively. To improve the accuracy of the flow path selection, two additional homogeneous inlet and outlet flow channels were attached to the left and right sides of the main geometric model.

Table 1. Physical parameters of basic pore structure.

		Low PTR	Median PTR	High PTR
Porosity		0.27	0.20	0.18
Pore–throat ratio		6.35	12.12	23.67
Specific surface area (103 m2/m3)		5.84	3.20	1.64
Size (μm)	Pore	317.54	606.22	1183.57
	Throat	50
	(Length, Width, Depth)	(4550, 4600, 50)

presents the detailed parameters of the three sets of models, which specifically include porosity, pore–throat ratio, Kozeny–Carman permeability, specific surface area, and pore throat size. To further investigate the effect of heterogeneity on the pore flow characteristics, we arranged these geometric models in three layers in the longitudinal direction and constructed four different nonhomogeneous structures by changing the order of the layers, as shown in Figure 2. The four heterogeneous models all contain three basic models of the pore–throat ratio. The four models are called the HML, HLM, MHL, and LMH models in order from top to bottom. These structures were designed to simulate the differences in the flow behavior of the two phases of H₂ and brine inside the pore medium under different heterogeneity structures. The simulations in this paper use a simplified two-dimensional model of the pore structure, specifically setting the depth of the model to 50 µm, which is about 1/100 of the length and width. Due to the relatively small depth, the flow effect in the depth direction can be neglected, so this two-dimensional approximation is reasonable.

After completing the design of the heterogeneity models, these models were imported into the computational environment of OpenFOAM using the STL file format [32]. Meshing was performed using OpenFOAM’s blockMesh and snappyHexMesh commands [33,34]. A background hexahedral mesh was generated by blockMesh, and these background meshes were further refined using the snappyHexMesh command. Especially in the pore and pore throat regions, a tetrahedral mesh was used to delineate the detail regions with higher precision to accurately capture the complex structures between the pores and the pore throats. Also, snappyHexMesh was responsible for excluding meshes outside the computational domain. The quality of the generated mesh was verified by the checkMesh tool to ensure that it satisfied the convergence and accuracy requirements of the numerical simulation. The final geometric model and mesh features ensured the accuracy of the physical design, as well as the stability and accuracy of the subsequent numerical simulations. The mesh and geometric model features are shown in

Table 2. Pore structure grid characterization parameters.

Parameter	Value
Size (μm) (L × W × H)	4550 × 4600 × 50
Original voxels	360 × 360 × 4
Number of cells	125,016
Number of points	180,120
Porosity	0.25

.

2.2. Governing Equations

The two-phase immiscible, compressible flow model is used to characterize underground hydrogen storage in saline aquifers. In this study, the two-phase flow process of H₂–brine in a pore medium is simulated. The solubility of H₂ in brine is only 0.002 mole fraction at 15 MPa, 323 K reservoir conditions; thus, the dissolution of H₂ in brine is neglected in the simulation [35,36]. H₂ is highly compressible under reservoir conditions [19,37]. Therefore, this study seeks to solve the two-phase compressible, immiscible fluid flow process. The governing equations used for the pore-scale simulations of H₂–brine two-phase flow include the mass conservation and momentum conservation equations, and in the presence of gravity as well as a source term, the momentum equation for the immiscible multiphase system is expressed as follows [38,39]:

$$\nabla ·\rho \mathit{U}+{\displaystyle \frac{\partial \rho }{\partial t}}=0,$$

(1)

$${\displaystyle \frac{\partial \rho \mathit{U}}{\partial t}}+\nabla ·\left(\rho \mathit{U}\mathit{U}\right)=-\nabla p+\rho \mathit{g}+\nabla ·\mathit{\tau }+{\mathit{F}}_{\mathit{\sigma }},$$

(2)

where $\rho$ refers to fluid density, $t$ is time, $\mathit{U}$ is the velocity vector, $p$ is the flow pressure, $\mathit{\tau }$ is the stress tensor, ${\mathit{F}}_{\mathit{\sigma }}$ is the interfacial force, and $\mathit{g}$ is the gravitational acceleration.

H₂ and brine are Newtonian fluids; hence, the stress tensor $\mathit{\tau }$ can be further expressed as Equation (3) [40]:

$$\nabla ·\mathit{\tau }=\nabla ·(\upsilon (\nabla \mathit{U}+\nabla {\mathit{U}}^{\mathit{T}})).$$

(3)

This study uses a VOF (Volume of Fluid) model to simulate the interfacial flow of a two-phase immiscible fluid [41,42], whereas in the VOF method, a variable $\alpha$ is defined to represent the phase volume fraction of the fluid, and the overall fluid properties of each grid block are defined through the following equation:

$$\rho ={\alpha }_{H_2}{\rho }_{H_2}+{\alpha }_{brine}{\rho }_{brine},$$

(4)

$$\upsilon ={\alpha }_{H_2}{\upsilon }_{H_2}+{\alpha }_{brine}{\upsilon }_{brine},$$

(5)

where the subscripts denote H₂ and brine, and $\rho$ and $\upsilon$ stand for density and viscosity, respectively.

We use the flux-corrected transport (FCT) algorithm-based Volume of Fluid (VOF) method to characterize the movement of H₂–brine interfaces. The VOF method derives the interfacial force ${\mathit{F}}_{\mathit{\sigma }}$ as a function of the phase fraction gradient; see Equation (6):

$${\mathit{F}}_{\mathit{\sigma }}=\sigma \kappa \nabla \alpha ,$$

(6)

where $\sigma$ is the interfacial tension coefficient, $\kappa$ is the curvature radius at the interface, and $\alpha$ is the phase fraction indicating the volumetric percentage of a fluid.

The phase transport equation in the compressible VOF model is expressed as follows:

$${\displaystyle \frac{\partial {{\rho }_{H_2}\alpha }_{H_2}}{\partial t}}+\nabla ·\left({\rho }_{H_2}{\alpha }_{H_2}\mathit{U}\right)=0,$$

(7)

$${\displaystyle \frac{\partial {{\rho }_{brine}\alpha }_{brine}}{\partial t}}+\nabla ·({\rho }_{brine}{\alpha }_{brine}\mathit{U})=0,$$

(8)

OpenFOAM uses the Finite Volume Method (FVM) to discretize Equations (1)–(3) in order to solve the numerical solutions of these equations [40]. The Finite Volume Method (FVM) is a commonly used numerical computation method that discretizes partial differential equations by dividing the computational region into a number of control volumes and integrating over these volumes. The governing equations are solved through discretization, linearization, prediction, correction, and output steps. First, the Navier–Stokes and continuity equations are discretized to a finite volume computational domain, and time is also discretized to capture the flow dynamics. The nonlinear convection term is linearized by the value of the last iteration and simplified to an easily solvable form. The velocity field is then predicted from the linearized equations and the pressure gradient from the previous time step, and the pressure is corrected to satisfy the continuity equation by solving the Poisson equation. This correction process is repeated twice to ensure stability and accuracy, and the solution is ultimately determined and used for the next time step.

The compressibleInterFoam solver in OpenFOAM^® is used in this study to solve multiphase flow problems. The compressibleInterFoam solver is an improved version of the interFoam solver [40,43]. This solver is based on the Volume of Fluid (VOF) method for multiphase flow modeling, which enables an accurate simulation of the interaction and flow behavior between different phases by tracking the location of the interface between the two phases.

To normalize the flow rate for multiscale analysis, the dimensionless capillary number (Ca) is defined, which represents the ratio between the capillary and viscous forces. The Bo number is used to characterize the relative importance between gravity and capillary forces under reservoir conditions where gravity is considered. The two uncaused numbers ($Ca$ and $Bo$) together provide a quantitative means of analyzing the balance between viscous, capillary, and gravitational forces:

$$Ca={\displaystyle \frac{\mu u}{\gamma cos\theta }},$$

(9)

$$Bo={\displaystyle \frac{∆\rho gh{l}_t}{\gamma cos\theta }},$$

(10)

where $\mu$ represents the viscosity; $u$ is the injection velocity; $\gamma$, $\theta$, and $∆\rho$ are the interfacial tension coefficients, the contact angle, and the difference in the density of the two phases, respectively; $g$ is the gravitational acceleration; $h$ denotes the total pore height; and $l_t$ is the pore throat size.

2.3. Initial and Boundary Conditions

We use initial and boundary conditions, mesh scenarios, and solver configurations for our simulations that have been shown to be reliable in the literature. The left boundary of the heterogeneous model is the inlet, and the right boundary is the outlet. The inlet velocity boundary is set using fixedvalue, and the outlet velocity boundary is set to pressureInletOutletVelocity, which is applicable when pressure is specified on the velocity boundary. The internal pore wall boundary is set to the noSlip solid wall condition. The outlet pressure boundary is set to fixed pressure using totalPressure, and the inlet and wall boundaries use the fixedFluxPressure condition, which automatically adjusts the pressure gradient to make the boundary fluxes conform to the velocity boundary condition. The initial phase field conditions are defined using the setFields dictionary file, making the porous medium saturated with brine before the simulation. The phase field boundary is set to the gas-phase inlet, and the pore walls use the constantAlphaContactAngle condition to set the contact angle between the rock and the two-phase fluid. The specific boundary condition settings are shown in

Table 3. The boundary conditions used in this study.

Parameter	Inlet Boundary Condition	Outlet Boundary Condition	Solid Walls
p	calculated	calculated	calculated
p_rgh	fixedFluxPressure	totalPressure	fixedFluxPressure
U	fixedValue	pressureInletOutletVelocity	noSlip
α	fixedValue	zeroGradient	constantAlphaContactAngle
T	fixedValue	fixedValue	fixedValue

.

In this study, the two-phase H₂–brine flow in the pore space is simulated, and the process of H₂ driving out the brine is the target to be investigated. The physical property parameters are set with reference to the ranges given in the literature [27,44,45,46]. Four heterogeneous structures and three wettability simulations are set up with contact angles of 80°, 55°, and 30°. The gas–water physical parameters at 10 MPa and 35 °C are selected for the simulation, and the inlet velocity is set to 0.01 m/s. The Ca and Bo numbers are calculated based on the parameters, and the specific parameter values are shown in

Table 4. Fluid physical parameters in simulation.

Parameter		Value
Pressure (MPa)		10
Temperature (℃)		35
Density (kg/m3)	H2	7.43
	Brine	1400
Viscosity (Pa·s)	H2	9.22 × 10−6
	Brine	7.19 × 10−4
Contact angle (°)	H2	100	125	150
	Brine	80	55	30
Interfacial tension (N/m)		0.077
Inlet velocity (m/s)		0.01
Ca		6.90 × 10−6	2.09 × 10−6	1.38 × 10−6
Bo		0.235	0.071	0.047
Bo·Ca−1		3.41 × 104	3.40 × 104	3.41 × 104

.

3. Results and Discussion

3.1. The Influence of Heterogeneous Structures

Pore space is the basic storage space for fluid storage in tight reservoirs, and the throat is an important channel for controlling fluid seepage. The shape and connectivity of pore throat structures control the physical interaction between rock and fluid and the storage capacity. The heterogeneity of the pore structure significantly affects the pore-scale filling mechanism and the development of preferential flow paths, so we constructed four different heterogeneous structures to study the effect of heterogeneity on the two-phase flow patterns of H₂ and brine. Figure 3 illustrates the flow distributions of the four heterogeneous structures at different times when the contact angle of the brine is 80°. The differences in the flow distributions mainly result from the influence of the pore structure on the flow paths, and even though the pore throat sizes remain consistent in the three basic geometries, they still exhibit significant flow differences due to the different arrangements and distributions of the pore sizes.

Observing the flow distributions of the four nonhomogeneous structures in Figure 3, it can be found that the flow distributions of the HML and LMH structures are more similar at different times. In both structures, the flow in the medium-pore–throat ratio region is more favorable, and the H₂ flow in the high-pore–throat ratio region is restricted. This similarity is due to the consistent adjacency between the heterogeneous pores of these two structures; the intermediate regions both have a medium pore–throat ratio. The difference between these two structures is mainly due to the effect of gravity. There is a density difference between the liquid and the gas, and the gas is more likely to migrate upwards due to buoyancy. This causes H₂ to invade the upper pores more easily. The flow distance of the low-pore–throat ratio structure is longer when flowing in the upper space compared to that in the lower, and the high-pore–throat ratio structure can break through the first throat into the second pore space when in the upper. From the results of these two heterogeneous structures, it can be seen that when the pore heterogeneous structures are arranged according to high to low or low to high pore–throat ratios, the structural part with high pore–throat ratios is unfavorable for the flow of H₂, and more H₂ is retained in the larger pore spaces. Comparing the results of the flow distributions in the HML, HLM, and MHL structures, obvious differences can be found, where the flow in the pores in the low-pore–throat ratio regions of the HLM and MHL structures is significantly restricted. It is interesting to note that the flow in the region of the medium pore–throat ratio is favorable in all four heterogeneous structures, for which this phenomenon is analyzed.

Among the three pore–throat ratio pore structures designed in this study, the flow topology of the low-pore–throat ratio structure is more complex with more flow bifurcation and focusing routes, and the low-pore–throat ratio structure needs to pass through a more narrow throat when H₂ flows the same distance. The fluid will generate larger pressure loss when passing through a smaller throat, so the pressure in the region of the low pore–throat ratio is lower than that in the regions of the medium and high pore–throat ratios when H₂ flows to the same distance in the lateral direction. As can be seen from Figure 3a₁,b₁, the flow distribution in the three pore–throat ratio regions is relatively uniform in both the HML and HLM structures at t < 0.02 s. When t increases to 0.04 s due to the difference in flow resistance, the distance of H₂ breakthrough in the pores with low pore–throat ratios is much smaller than that in the pores with medium and high pore–throat ratios on both sides. Hydrogen flows from the pores with medium and high pore–throat ratios on both sides into the pores connected with the central low-pore–throat ratio region, resulting in a localized pressure increase in the central region. This will further exacerbate the uneven pressure in the central region. The low-pore–throat ratio region is more susceptible to the influence of adjacent regions, and the flow obstruction is more obvious when adjacent to the high-pore–throat ratio region than when adjacent to the medium-pore–throat region. This phenomenon can be understood as a localized obstruction of flow in small-pore regions due to the complex pressure distribution in heterogeneous pores, where hydrogen cannot continue to flow, and the breakthrough process is blocked or stalled.

The three pore–throat ratio structures have the same throat sizes and large differences in pore sizes. The high-pore–throat ratio structures have large differences in pore and throat sizes and can be seen as many converging–diverging capillaries. Previous studies have shown that the wetting behavior of converging–diverging capillaries is much slower than that of straight capillaries. Moreover, when the diameter of the downstream capillaries is significantly larger than that of the upstream ones, the local capillary pressure is low, and the viscous force is very high [47,48]. This is due to the fact that in real porous regions, the velocity of the interface changes abruptly as it passes through a wide pore and adjacent narrower throat, and the interfacial velocity at the local pore scale may increase by several orders of magnitude, leading to non-negligible viscous flow effects that change the apparent contact angle at the interface. The reversal of interfacial curvature also occurs in converging–diverging capillaries [49,50], and this reversal of interfacial curvature can be clearly seen in Figure 3 in the region of the high pore–throat ratio. The reason for the coexistence of concave–convex interfaces is that the interface morphology becomes more and more dependent on the pore angle (the angle at which the pores converge or diverge) as the wettability changes from strongly to moderately wetting. In this study, the actual interfacial curvature, as well as the flow state in the pores, is the result of a complex interaction between the pore geometry, the angle, and the effective contact angle. In this study, the simulation results show that the pore–throat ratio change from 6.35 (low pore–throat ratio) to 12.12 (medium pore–throat ratio) promotes H₂ flow, while the change from 12.12 (medium pore–throat ratio) to 23.67 (high pore–throat ratio) negatively affects H₂ flow. Even with the same throat size, the pore size, arrangement, and heterogeneity still lead to different flow characteristics as a result of the combined effect of capillary, viscous, and inertial forces in the heterogeneous pore structure.

Figure 4 shows the pore pressure distribution curves in the above four sets of simulations. The existence of the throat makes the pressure in the pore space show obvious nonlinear distribution characteristics, the pressure drop mainly occurs in the throat, and the larger the pore–throat ratio, the larger the pressure drop in the throat. At the initial moment, the liquid phase fills the pore space, so the liquid flow will dominate the pressure distribution within the pore space. In Figure 4, the curve display for t = 0.005 s reflects the pore pressure characteristics at the initial stage of flow. As the gas enters the pore space, the pressure distribution will become more complicated due to the existence of the gas–liquid interface. It can be found that the pressure curve at t = 0.015 s increases in the region of 0–1000 μm compared to that at t = 0.005 s, which is due to the presence of the gas–liquid interface that significantly increases the flow resistance. At t = 0.03 s, the heterogeneous expansion of the gas–liquid interface causes the pressure curve to show several fluctuations. At t = 0.045 s, the pressure curve of Figure 4b has lower values in the range of 0–3000 μm compared to the other three curves. This is due to the restricted H₂ flow in the pores in the low-pore–throat ratio region of the HLM pore structure, which leads to a steep drop in the pressure curve at the gas–liquid interface at a distance of 1000 μm. At t = 0.06 s, H₂ in the four heterogeneous structures breaks through the outlet, the gas forms a stable and dominant channel in the system, and the effect of capillary force is weakened.

From the pore-scale simulation results, it can be seen that parameters such as heterogeneity, pore shape, and pore–throat ratio are crucial for clarifying the H₂ flow characteristics in underground pores. Therefore, the influence of these factors should be fully considered in the core-scale or reservoir-scale studies of underground hydrogen storage to establish more accurate reservoir models and flow simulations.

3.2. The Influence of Wettability

In the study of gas–liquid two-phase flow, wettability is a key factor affecting the behavior of the gas–liquid interface and flow characteristics. Wettability has an important effect on the preferential flow path and saturation of H₂. Normally, reservoir rocks are hydrophilic, which means that the rock surface is more easily wetted by brine [51]. In this study, flow in four heterogeneous structures is simulated by setting three different contact angles between rock and brine (80°, 55°, and 30°). The effects of wettability on gas breakthrough, flow paths, and system pressure distribution are analyzed.

Figure 5 presents the variation curves of H₂ saturation in the four heterogeneous structures. In the case of smaller contact angles, the capillary resistance of H₂ through the pores is larger, which makes it difficult to expand rapidly, thus leading to lower gas saturation in the pores. Therefore, usually, the saturation of H₂ in the pores is the lowest when the contact angle is 30°, while the gas breakthrough is easier and the saturation in the pores is the highest when the contact angle is 80°. From Figure 5, it can be seen that in the three heterogeneous structures, HML, LMH, and MHL, the final saturation of H₂ decreases gradually with the decrease in contact angle. However, in the HLM structure, the situation is different. The saturation of H₂ is the highest when the contact angle is 55° and the lowest when the contact angle is 80°. By analyzing the phase distribution diagrams, it can be speculated that this phenomenon may be related to the longitudinal breakthrough of H₂ occurring in the pores. Specifically, in the HLM structure, the transverse breakthrough in the central low-porosity ratio region becomes more difficult as the contact angle decreases. The heterogeneity of the gas–liquid interfacial distribution makes the heterogeneity of the longitudinal pressure distribution increase, leading to the longitudinal breakthrough of H₂ from the medium- and high-pore–throat ratio regions on both sides to the central low-pore–throat ratio region.

The variation in the values of the final H₂ saturation for different heterogeneous structures (HML, HLM, MHL, LMH) at contact angles of 80°, 55°, and 30° is illustrated in Figure 6. Specifically, the H₂ saturation is generally higher for heterogeneous structures when the contact angle is 80°, reaching about 0.5 for the HML, MHL, and LMH structures. In contrast, the HLM structure has the lowest value of only about 0.41. As the contact angle decreases to 55°, the H₂ saturation of the HML, LMH, and MHL structures decreases to a value in the range of 0.38~0.46. As the contact angle is further reduced to 30°, H₂ saturation decreases in all heterogeneous structures. Overall, H₂ saturation usually increases with increasing contact angle. However, there are differences in the effects of different pore structures on H₂ saturation, especially the HLM structure, where H₂ saturation increases and then decreases with increasing contact angle. It can be seen that the H2 saturation of the MHL structure is higher at all three contact angles, and the magnitude of the change in H2 saturation is smaller with the decrease in the contact angle. This is because the main flow channels in the MHL structure are concentrated in the regions of medium and high pore–throat ratios, where the number of throat channels is small, and the pore space is relatively large. Therefore, the change in contact angle has a relatively small effect on these regions. In contrast, the effect of wettability change is more significant in the low-pore–throat ratio region due to the higher number of throat channels and the greater influence of the capillary force. There is also a big difference in the heterogeneous model of the maximum and minimum H2 saturation at the end of drainage under the three contact angle conditions. This suggests that the heterogeneity of the pore structure and the combined effect of the contact angle are important factors affecting H2 saturation at the end of drainage, and their complex interaction determines the final distribution state of the gas in the reservoir. This also suggests that the heterogeneity and wettability characteristics of the reservoir need to be considered comprehensively in practical engineering to optimize gas injection and drainage strategies and thus improve the efficiency of gas storage and utilization.

Overall, the results show significant differences in the effect of wettability in different heterogeneous structures. In heterogeneous structures, the effect of wettability does not always follow a simple pattern. This means that although normally, a decrease in wettability enhances capillary resistance and reduces gas saturation, in some complex heterogeneous structures, localized regions of pore structures and pressure distribution can change this trend and lead to different flow behaviors than expected. Therefore, the study of the interaction between heterogeneity and wettability is crucial for understanding H₂–saline two-phase flow in complex porous media.

3.3. Gravitational Influence

For two-phase H₂–brine flow in pores, in addition to capillary forces, there is also the effect of gravity. Under gravity, the liquid will tend to gather downward due to its higher density, while the gas is more likely to migrate upward due to the buoyancy effect. This effect leads to a clear upward and downward divergence in the distribution of two-phase fluids in porous media. We simulated the MHL structure without considering the effect of gravity, and only the results for a contact angle of 80° show differences for the three wettability conditions. This is due to the fact that the effect of capillary forces becomes more significant as the contact angle decreases, thus weakening the effect of gravity on the liquid flow path. In this case, the capillary forces dominate the flow process, and the difference in the gas–liquid distribution mainly stems from the effects of pore structure and wettability rather than gravity. Figure 7 shows the simulation results with and without considering gravity at a contact angle of 80°. The difference in the flow distribution in the region of the high pore–throat ratio occurs in both cases, where the gas is affected by an upward force due to the presence of density difference when gravity is considered, and in the upper pores, the gas in the pore space is unable to pass the throat that extends downward due to this force acting upward. However, in the simulation, this gravity effect on gas–liquid flow has a small range of influence, mainly concentrated in the local pore space region. However, in actual reservoirs, especially in large-scale reservoirs, the effect of gravity on the two-phase flow paths and phase distribution cannot be ignored. Gravity not only changes the flow paths of gas and liquid in porous media but also affects the breakthrough mode of gas, the gathering location of the liquid phase, and the overall flow behavior of the system.

4. Conclusions

To better understand the pore-scale behavior of hydrogen in the underground heterogeneous porous structure, this study explores the effects of wettability, pore–throat ratio, and pore structure heterogeneity on the two-phase flow behavior of H₂–brine by constructing heterogeneous pore models with different pore–throat ratios and arrangements. Based on the OpenFOAM platform, the two-phase flow process of H₂–brine in four heterogeneous pore media is simulated by the Volume of Fluid (VOF) method. According to the simulation results, the differences in H₂–brine flow distribution under different heterogeneous structures and wettability conditions are analyzed, as well as the effects of different pore–throat ratios and pore arrangement structures on the H₂ flow behavior. The results show the following:

(1)

There are obvious differences in the H₂–brine flow distributions in the four heterogeneous structures, which are due to the combined effects of capillary, viscous, and inertial forces in the heterogeneous pore structures.

(2)

H₂ passing through the region of the low pore–throat ratio suffers from greater capillary resistance, and the breakthrough distance is shorter compared to the neighboring regions. The low-pore–throat ratio region is more easily affected by the neighboring pore region, and the obstruction effect is more obvious when it is adjacent to the high-pore–throat ratio region than when it is adjacent to the medium-pore–throat region. This is a result of the complex pressure distribution in heterogeneous pores.

(3)

In high-pore–throat ratio structures, an abrupt change in interfacial velocity is associated with the difference between wide pores and adjacent narrow throats. The interfacial velocity may increase by several orders at the local pore scale, leading to non-negligible viscous flow effects.

(4)

In this study, the increase in the pore–throat ratio from 6.35 (low pore–throat ratio) to 12.12 (medium pore–throat ratio) promotes H₂ flow. In contrast, the increase in the pore–throat ratio from 12.12 (medium pore–throat ratio) to 23.67 (high pore–throat ratio) negatively affects H₂ flow. This suggests that the complex interaction between contact angle and pore geometry plays a crucial role in controlling the pattern of displacement.

(5)

In the HLM structure, the saturation of H₂ is the highest at a contact angle of 55° and the lowest at a contact angle of 80°. This phenomenon suggests that although normally, a decrease in wettability enhances the capillary resistance and reduces the saturation of the gas, the pore structure and pressure distribution in localized regions in complex heterogeneous structures may change this trend.

(6)

The influence of gravity on the two-phase H₂–brine flow should not be neglected, and the gravity effect will change the gathering location of the liquid phase as well as the flow path of H₂.

This study reveals the influence of heterogeneous pore structure, wettability, and pore–throat ratio on the flow behavior of hydrogen in underground porous media. Insights are provided for understanding the role of the heterogeneity of pore structures in H₂–brine two-phase flow during underground hydrogen storage.