Pore-scale modelling and sensitivity analyses of hydrogen-brine multiphase flow in geological porous media

Underground hydrogen storage (UHS) in initially brine-saturated deep porous rocks is a promising large-scale energy storage technology, due to hydrogen’s high specific energy capacity and the high volumetric capacity of aquifers. Appropriate selection of a feasible and safe storage site vitally depends on understanding hydrogen transport characteristics in the subsurface. Unfortunately there exist no robust experimental analyses in the literature to properly characterise this complex process. As such, in this work, we present a systematic pore-scale modelling study to quantify the crucial reservoir-scale functions of relative permeability and capillary pressure and their dependencies on fluid and reservoir rock conditions. To conduct a conclusive study, in the absence of sufficient experimental data, a rigorous sensitivity analysis has been performed to quantify the impacts of uncertain fluid and rock properties on these upscaled functions. The parameters are varied around a base-case, which is obtained through matching to the existing experimental study. Moreover, cyclic hysteretic multiphase flow is also studied, which is a relevant aspect for cyclic hydrogen-brine energy storage projects. The present study applies pore-scale analysis to predict the flow of hydrogen in storage formations, and to quantify the sensitivity to the micro-scale characteristics of contact angle (i.e., wettability) and porous rock structure.

www.nature.com/scientificreports/ behaviour through cyclic transport, and whether empirical functions in the literature are valid. Note that the supplementary materials provide complementary data sets that were left out of the main manuscript, and a list of all UHS projects.

Methods
Multiphase flow properties in subsurface reservoirs can be predicted by several methods including laboratory measurements and numerical simulations [31][32][33][34][35][36] . Experimental approaches are costly and time consuming, and are found for a limited samples and conditions 31 . Numerical modelling and simulations are therefore crucial to complement the laboratory studies to allow for wider range of studies and sensitivity analyses.
In this study, the most computationally-efficient pore scale approach, which is called quasi-static pore network modelling "PNM", was used to simulate fluid flow of UHS at the pore-scale and to estimate the macro-scale properties, i.e., capillary pressure and relative permeability. Quasi-static PNM is an appropriate method for the capillary-dominated flow regime where the capillary number is less than 10 −437 . This is a common range for immiscible two-phase flow in many subsurface applications. All the simulations in this paper were implemented by using an open-source software from Imperial College London, called "pnflow" 38 , which is validated with experimental results for hydrocarbon reservoirs 37,39 .
Below, the key components of pore-network modelling are first presented. Then the specific developments relevant for H 2 -brine transport are discussed. Detailed information about the well-developed pore-network modelling approach can be found in the literature 31 . Here, for the sake of brevity, we focus on new developments relevant for H 2 -brine systems.
• Pore-Network Modelling (PNM) description PNM uses a topologically and geometrically equivalent network of a rock sample to predict capillary pressure and relative permeability by simulating fluid flow through the elements of the network. Processing 3D micro-CT images of a core plug provides information about the network and its elements including radius, volume, length, coordination, clay volume, etc. 40,41 . The elements with larger volume represent "pores" and the elements which connect the pores are called "throats" 31,42 . Figure 1 shows an illustration of a 3D image of a sandstone and its extracted pore network 37 .
All the elements are uniform ducts with various cross-sectional shapes which are classified into circle, square and triangle shape factors G, which is defined as A/P 2 , where A is the cross-sectional area and P is the perimeter length 31,43 . Figure 1 presents an illustration of pore shapes. Most of the earliest network models assumed circular cross-sections for simplicity 44 . One major disadvantage of a circular shape is that it cannot accommodate more than one fluid in a stable configuration in a single pore. Therefore, it does not allow for any films or layers of additional phases to be formed during displacement. In contrast, it can be clearly seen from thin section images of rocks that real pore shapes are very irregular and have many corners. Experimental studies have shown that in pores with angular cross sections, the wetting phase can occupy the corners with the non-wetting phase in the centre. These corner wetting layers provide additional phase connectivity and hence have an impact on trapping. Thus, networks that have pore elements with angular cross sections provide the opportunity to predict experiments where corner flow is crucially important.
In a quasi-static network model the capillary pressure is the driving force that determines the saturation evolution in the model. This quantity is changed gradually and, for each capillary pressure value, the equilibrium position of the fluid-fluid interfaces is determined. The displacements or the changes of fluid configurations occur sequentially, according to the entry capillary pressure criteria and the connectivity of each fluid to the inlet and outlet of the network. Since the network is initially saturated with the wetting phase, displacement process starts with primary drainage followed by secondary imbibition and finally secondary drainage to simulate cyclic injection and withdrawal. The details of the fluid-flow simulation procedure and governing equations have been explained in the literature 37 .
• Fluid and rock parameters  Table 2.
In this work, we define the base-case for the H 2 -brine system as the one obtained based on the experiment described previously 2,45 . Since the experimental data for hydrogen-brine is limited to only the primary drainage cycle (initial injection of H 2 ), to estimate advancing contact angles in all models, the Morrow relationship, as shown in Fig. 2, was used. However, for receding angles more than approximately 12 • , if one follows the Morrow curves, the corresponding advancing contact angle becomes bigger than 90 • . This shows the rock becomes hydrogen-wet, which is very unlikely to be realistic. To resolve this challenge, in this study, the maximum advancing contact angle is, therefore, set to 85 • , as shown in Fig. 2. This modification guarantee that hydrogen is the non-wetting phase in all cases. Additionally, a uniform distribution of contact angle (a constant value in all pores) is employed for all the simulations in this paper. Also, as for the rock, relatively homogeneous Berea sandstone is considered, unless otherwise stated. The extracted network data is illustrated in Fig. 3 47 .
The base case scenario for our analyses is defined as follows. The base rock is considered to be the Berea sandstone with the characteristics presented above. The receding contact angle and the interfacial tension for H 2 -brine in the base case are defined following the available experimental data at the core scale 2,45 . The receding contact angle, assumed to be constant across the micro-scale sample, was found by fitting to the measured capillary pressure. As presented in the literature, this results in θ r = 21.56 • and σ = 51 [mN/m], as shown in Table 2. Based on the modified-Morrow hysteresis curve, Fig. 2, the advancing contact angle for the base case is set as θ a = 85 • . Using these data sets, the relative permeability and capillary pressure for the base case is presented in Fig. 4. Note that the base case properties correspond to the Shallow formation reported in Table 2. Moreover, the results corresponding to the fluid properties of the Deep formation are also presented in Fig. 4. Consistent with the reported experiments, our pore-network model results also confirm insensitivity towards different pressure and temperature values of Shallow and Deep formations. This is because we only have a small change in interfacial tension, which affects the capillary pressure, and no change in advancing contact angle, leading to identical relative permeabilities in the two cases.
What comes next is a systematic analysis and quantification of uncertainty for cyclic transport and within the parameter range of the available experimental data 2, 45 . In addition, we investigate whether empirical Corey     Table 3. Note that to differentiate the effect of different rock parameters such as clay volume and coordination number (i.e., the number of throats connected to a pore), statistical models were generated 47 .

Results
Hydrogen-brine transport for energy storage applications is unique in some aspects. Firstly, the upscaled transport functions need to be benchmarked against empirical functions employed in the literature for other gas-brine systems, especially Corey functions. Our first test case addresses this important aspect. Then, the cyclic nature of the storage application, requires the simulation of cycles of injection and production. This leads to adding the secondary drainage into the pore-network modelling framework. The second test cases deals with this important aspect, and quantifies the hysteretic nature of the upscaled functions around the base case scenario. Lastly, uncertainty in the fluid and rock properties needs to be investigated, to find their impact on the upscaled multiphase flow functions. The last two test cases address the uncertainty in the fluid and rock properties, respectively. With these studies we aim to highlight the key aspects of hydrogen-brine transport for energy storage. A summary of the test cases is given in the list below.
Test Case 1 Benchmarking pore-network upscaled functions with those obtained with the widely-used Corey equation.
Test Case 2 Hysteretic upscaled transport functions for cyclic transport Test Case 3 Impact of fluid properties on the hysteretic upscaled functions Test Case 4 Impact of rock properties on the hysteretic upscaled functions Test Case 1: Benchamrking pore-network modelling (PNM) and Corey functions. Figure 5 shows the relative permeabilities for hydrogen injection into brine-saturated rock. Since in the literature no hysteresis effects were considered 30 , the PNM results are also plotted only for primary drainage. The best fit to the PNM results, constrained with Corey equation parameters 48 , resulted in different values for the exponents of hydrogen and brine, 1.31 and 4.36 respectively. However, in the literature the constant exponent of 2.5 for both hydrogen and brine has been used 30 . As such, the PNM-based studies indicate that different exponents for hydrogen and brine need to be considered. This makes physical sense, since the exponents are related to pore structure and wettability and are higher for the wetting (brine) phase as it fills the smaller pores, as opposed to lower exponents (higher relative permeabilities) for the non-wetting phase (hydrogen) that preferentially occupies the larger regions of the pore space. www.nature.com/scientificreports/ Test Case 2: Hysteretic upscaled transport functions for cyclic transport. In real-field hydrogen storage projects, there are repeated cycles of injection and withdrawal. As such, both primary and secondary drainage processes will occur in addition to secondary imbibition. To investigate this, Fig. 6 shows the impact of cyclic injection and withdrawal of hydrogen on capillary pressure and relative permeability. Note that, as discussed in the Methods Section, the modified Morrow values for advancing and receding contact angles have been used. Compared with primary drainage, secondary drainage shows lower relative permeability values at the same saturation. This is due to disconnection of the hydrogen phase after secondary imbibition, as well as the presence of trapped wetting-phase (brine) in some pores. Also, the remaining hydrogen phase, after secondary imbibition, leads to a decrease of capillary pressure values for the secondary drainage stage. Additional cycles of injection and production of hydrogen into the network of Berea sandstone are also investigated and provided in the supplementary material. Our observations indicate that the hysteresis effect stays the same for subsequent drainage and imbibition displacement, after the first injection-production cycle.

Test Case 3: Impact of fluid properties on the hysteretic upscaled functions. Among the fluid
properties, wettability is found to have major impacts on the upscaled functions. The findings are discussed in two separate parts as follows.

Effect of wettability
The range of simulated contact angles for H 2 -brine is given in Table 4. Simulation results are shown in Fig. 7. Note that increasing contact angles (which are less than 90 • ) changed the wettability of system and made it less water-wet. Since the receding contact angles varied between 5 • to 30 o , there is no significant effect on the results for primary drainage. However, during secondary imbibition and secondary drainage, increasing contact angles resulted in a smaller amount of trapped hydrogen (lower capillary pressure). Also, the maximum relative permeability of the brine phase is shown for the strongly water-wet system with θ i = 51 • . Relative permeabilities of hydrogen during secondary drainage are mostly similar. However, during secondary imbibition, increasing contact angles towards neutral wettability (i.e., when intrinsic contact Figure 6. The impact of cyclic transport on capillary pressure and relative permeabilities for primary drainage (PD), secondary imbibition (SI), and secondary drainage (SD). Table 4. Fluid and rock properties used for the wettability sensitivity analysis, by changing advancing ( θ a ) and receding ( θ r ) contact angles. *Indicates the base-case from the literature 2 . www.nature.com/scientificreports/ angles are close to 90 • ) increases the relative permeability of hydrogen. This is due to the reduced amount of trapped hydrogen. 2. Effect of the difference between advancing and receding contact angles, i.e., �θ = θ a − θ r , with θ r = 21.56 • In this set of tests, all the parameters of the base-case remained constant except the advancing contact angle. Figure 8 shows the impact of the difference between receding and advancing contact angles on capillary pressure and relative permeabilities. The maximum trapped hydrogen was observed for the case with equal advancing and receding contact angles. For the higher advancing contact angles rock becomes less water-wet. Therefore, relative permeability of hydrogen increases. There is also less trapping of hydrogen.
Test Case 4: Impact of rock properties on the hysteretic upscaled functions. Rock structure has an important role in the transport behaviour of fluids. To study changes in the saturation-dependent functions for various rock types, five extracted pore network models from images of sandstones and carbonates were used 47 . The properties of these samples and the system of the fluids (H 2 -brine) are given in Table 5. Moreover, detailed rock structures for all samples are provided in the supplementary materials. The outcome of the simulations are shown in Fig. 9. Among the three sandstone models, small Berea and A1 showed almost zero water saturation trapped after injecting hydrogen which implies that they contain no clay. Similar patterns are www.nature.com/scientificreports/ observed for flow properties during primary drainage. Sample A1 which had the highest porosity and absolute permeability had the smallest residual hydrogen saturation and, as expected the highest values for the relative permeability of brine in secondary imbibition. Comparison of the two models from carbonate reservoirs indicates the impact of a complex pore structure. Sample C2 has the smallest pores and restricted connectivity with therefore the highest capillary pressures during all displacement cycles and the highest residual (trapped) hydrogen saturation after secondary imbibition. We define the maximum trapped hydrogen saturation as the one found after production corresponding to a capillary pressure P c = −100 kPa. Figure 10 shows that many injection-production cycles do not change the maximum amount of produced hydrogen significantly. However, the maximum trapped hydrogen is affected by the rock structure. Characterisation of these models (as provided fully in the supplementary material) indicates that lower connectivities and the smaller pores in model C2 is indeed the reason behind its the highest residual  www.nature.com/scientificreports/ hydrogen saturation, among all samples. On the other hand, model A1 with the maximum average coordination number and high permeability allows for the highest hydrogen production (i.e., minimum residual hydrogen saturation). In addition, the two models of Berea sandstone, with similar characteristics but different sample sizes, have almost equal storage efficiency. To differentiate the effect of different rock parameters such as clay volume and coordination number, some statistical models were generated using an open source software 47 . The percentage of clay volume in the network directly changes the trapped water saturation after primary drainage, but the patterns of capillary pressures and relative permeability remain similar. Increasing the coordination number reduces capillary entry pressures and trapped hydrogen saturation during secondary imbibition. Moreover, its effect on relative permeabilities becomes significant during displacement of hydrogen by water: higher coordination number facilitates flow and suppresses trapping. Full details are presented in the supplementary material.

Conclusions
H 2 -brine transport properties are quantified at the continuum scale through capillary pressure and relative permeability. These functions were predicted based on pore network modelling (PNM) which simulates the porescale displacement of fluids. Through several systematic studies we first benchmarked the PNM fluid parameters with existing experimental data. This allowed us define a meaningful base case configuration. The brine remains as the wetting phase. When the relative permeabilities are fit to a power-law type empirical model, the exponents