Measurement of Effective Hydrogen-Methane Gas Diffusion Coefficients in Reservoir Rocks

According to the 2030 climate target plan of the European Union, greenhouse gas emissions should be reduced by at least 55% by 2030 (European Commission 2019). This is the first step in achieving the European Union’s net zero emission strategy by 2050 (European Union 2021). Reaching these ambitious targets is only possible with a step change in the energy sector, namely, substitution of fossil fuels by renewable energy sources (IEA 2021). One of the challenges with renewable energy is the intermittent character of energy generation (e.g., wind and solar) and the difference in seasonal energy production levels (Becherif et al. 2015). Therefore, energy storage becomes exceptionally vital for balancing energy supply and ensuring energy security (e.g., Crotogino et al. 2010; Juez-Larré et al. 2019; Scafidi et al. 2021). While daily energy consumption fluctuations can be covered by electricity storage in batteries, seasonal energy storage requires technologies providing significantly larger storage capacities (Parra et al. 2018).

Hydrogen is believed to be a suitable climate-neutral energy carrier, which can be stored underground in large volumes and with reasonable discharge times (e.g., Zivar et al. 2021; Mouli-Castillo et al. 2021; Clemens and Clemens 2022; Delshad et al. 2021). Different types of reservoirs are being investigated for the purpose of hydrogen storage, such as aquifers (Sainz-Garcia et al. 2017; Alhotan et al. 2023), salt caverns (Malachowska et al. 2022; Michalski et al. 2017), and deleted gas reservoirs (Ali et al. 2022; Tarkowski 2019). When considering storage in depleted gas reservoirs, some of the already existing infrastructure can be repurposed for hydrogen storage (McKinsey 2021). Therefore, the use of depleted hydrocarbon fields as UHS sites is widely considered (e.g., Muhammed et al. 2022; Delshad et al. 2022).

When hydrogen is stored in depleted hydrocarbon reservoirs, the mixing of hydrogen and in-situ hydrocarbons occurs. This leads to alteration in the composition of the backproduced gas stream. It is crucial to understand the extent of mixing to properly design surface facilities as well as to meet customer requirements (e.g., Leicher et al. 2022). One of the factors affecting the distribution of hydrogen in the subsurface is reservoir heterogeneity (Alhotan et al. 2023; Shojaei and Jessen 2012). In addition, molecular diffusion is an important process, which affects the magnitude of gas mixing (Arekhov et al. 2023). Proper accounting of molecular diffusion is necessary to plan successful UHS projects.

Molecular diffusion is a fundamental physical process, which is well known and was studied by researchers in many industries (e.g., Cussler 1976; Barr 1997; Thambynayagam 2011; Graue 1965; Hu and Kobayashi 1970). In the petroleum industry, diffusion was widely investigated for hydrocarbon gasses, carbon dioxide, and nitrogen (Sigmund 1976; Ghasemi et al. 2017; Ratnakar and Dindoruk 2022; Dawson et al. 1970; Hayduk and Cheng 1971; Lo 1974; Reid et al. 1987; Shinji and Hiroji 1966; Woessner et al. 1969). Methane and hydrogen diffusion into water was also widely studied (e.g., Coats and Richardson 1967; Hemme and van Berk 2018). However, laboratory work concerning the measurement of binary molecular diffusion of methane-hydrogen systems is still limited (e.g., Chu et al. 1974; Keshavarz et al. 2022). When dealing with hydrogen storage in a porous medium, the effective diffusion coefficients must be taken into account. Effective diffusion scales with porosity as well as tortuosity of the rock. Measurements of effective diffusion coefficients were performed on water-saturated samples (e.g., Chen 1973; Grathwohl 1998; Lou et al. 2021; Eltahan et al. 2020; Shen and Chen 2007; Jacops et al. 2013). Nevertheless, there is a lack of diffusion measurements for hydrogen-methane systems in gas systems at various elevated pressures including reservoir rocks.

Numerous approaches exist to quantify diffusion effects and estimate molecular diffusion coefficients (Riazi and Whitson 1989; Da Silva and Belery 1989; Ghorayeb and Firoozabadi 2000; Yanze and Clemens 2012). The diffusion of hydrogen into the caprock was simulated by Heinemann et al. (2021); here, we investigate diffusion within heterogeneous reservoir rocks. The reason is that diffusion fluxes of gas in gas in heterogeneous reservoirs are significantly higher than diffusion into water or very low permeability caprocks and can significantly contribute to gas mixing in UHS.

In this paper, we describe the procedure of how to measure molecular diffusion of hydrogen and hydrocarbon gas through porous rock samples. We measured binary diffusion coefficients first without porous medium and then with three different rock types to investigate the effective diffusion coefficients of hydrogen in reservoir rocks. The result of the paper can be used to address diffusion-related mixing of hydrogen with hydrocarbon gas in reservoir rocks including pressure dependency of the effective diffusion coefficient.

The paper has the following content. In the beginning, a summary of molecular diffusion theory is given with a description of available models for calculating diffusion. Then, the experimental setup is described, including information on materials used, the procedures applied, the calibration process, as well as data acquisition and interpretation. The following section contains the experimental results of binary diffusion coefficients measured through the capillary. Afterward, the effective diffusion measurement results are given.

Molecular diffusion is a fundamental mass transfer phenomenon that can be expressed in two different ways. The first approach links the diffusive flux and driving force via a constant called diffusion coefficient. This approach, for instance, is applied in Fick’s law of diffusion, which has several formulations. In this work, diffusive transport is expressed by average molar flux, where the driving force is the concentration gradient (see Eq. 1; Ahmadi et al. 2021). Thus, the diffusive flux occurs from the high-concentration zone toward the low-concentration zone. Other driving forces which can be described similarly are pressure and temperature gradients, electrical potential, and others:

where J_i is the diffusive flux (mol/m²·s), D_ij is the diffusion coefficient of the binary mixture (m²/s), c is the molar concentration of component i, and x_i is the mole fraction of a gas component i.

The second approach takes the Brownian motion of the diffusing molecules into account. In this case, diffusion is associated with spreading out of molecules in space. Diffusion is defined as the mean square displacement of particles (Einstein 1905).

In this paper, we focus on the first approach to quantify the molecular diffusion of gasses in terms of the molecular diffusion coefficient, which can be further used for numerical forecasting of the mass transport due to diffusion.

As far as diffusion in the porous media is concerned, the molecular diffusion coefficient needs to account for the amount of accessible space in the rock (porosity) and the shape of the pore throats where the molecules move (tortuosity). Thus, the diffusion coefficients in the porous medium can be estimated by means of an effective diffusion coefficient (see Eq. 2; Grathwohl 1998):

where $τ$ is the tortuosity, $φ$ is the porosity, and D is the diffusion coefficient (without porous medium).

The magnitude of diffusive flux in the pore space depends on the relative size of the molecules and pore throats. The type of diffusion can be defined based on Knudsen number, which is the ratio of mean free path to pore diameter (see Eq. 3). It indicates whether the molecules rather collide with the pore wall or with another molecule (e.g., Ratnakar and Dindoruk 2022):

where Kn is the Knudsen number, $λ$ is the mean free path, $dp$ is the pore diameter, $kB$ is the Boltzmann constant, T is the temperature, P is the pressure, and $σ$ is the molecular size (radius).

If the Knudsen number is larger than 1, Knudsen diffusion occurs. It appears mainly in the porous medium with small pore sizes (<0.2 µm). The diffusion coefficient is reduced with decreasing pore diameter. On the other hand, when the pore size is large (>0.2 µm), regular bulk diffusion occurs. For the gas-gas system, bulk diffusion can be the dominant mass transfer phenomenon. The bulk diffusion is independent of the pore effects.

When considering underground gas storage in a porous medium, the typical pore size of reservoir sandstones is large enough for bulk diffusion to occur in the reservoir. Thus, in this paper, we estimate bulk diffusion by measuring binary diffusion coefficients in gas-gas systems.

The binary bulk diffusion coefficient can be calculated analytically from the kinetic gas theory. In this paper, we propose the following approach for cross-validation of measured diffusion coefficients (see Eqs. 4 and 5):

Chapman-Enskog formula (Bird et al. 2006):

where $Dij$ is the binary diffusion coefficient (cm²/s), $P$ is pressure (bar), $σij$ is the mean molecular size $σij=(σi+σj)2$ (Å), $Ωij$ is the collision integral obtained from Lennard-Jones parameters (Lennard-Jones 1931), $T$ is temperature (°C), and $Mw$ is molecular weight (g/mol).

Stiel and Thodos (1962) proposed the following form of the Chapman-Enskog correlation (see Eq. 5):

where $ρm0Dij0$ is the zero pressure limit of the density-diffusivity product (cm²/s); $R$ is the gas constant, which here is 83.14 due to unit conversion [cm³·bar/(K·mol)]; $σAB$ is the mean molecular size $σij=(σi+σj)2$ (Å); $ΩAB$ is the collision integral obtained from Lennard-Jones parameters (Lennard-Jones 1931), $T$ is temperature (K); and $Mw$ is molecular weight (g/mol).

Detailed description on how the parameters in Eqs. 4 and 5 are derived can be found in the corresponding literature (e.g., Lennard-Jones 1931; Hirschfelder et al. 1949; Neufeld et al. 1972). It is important to note that the described approach requires an estimate of an average mixture composition corresponding to the final equilibrium state. For the binary system, we have used 0.5–0.5 mole fraction composition as suggested by DaSilva for two-component fluids (Da Silva and Belery 1989).

Although the focus of this paper is diffusion in the binary hydrogen-methane system, this method can also be used for multicomponent mixtures. Da Silva and Belery (1989) described the method of diffusion coefficients estimation for each component in the mixture. Even if the effect of minor amounts of heavier components in dry gas reservoirs might be negligible, this approach can be important when dealing with molecular mixing effects during hydrogen storage in wet gas, condensate, or oil reservoirs.

For the ideal gas mixtures with fixed composition, the product $ρm0Dij0$ is constant (Balenovic et al. 1969). However, the experimental data (this work and literature) show that the density-diffusion product given by Eq. 5 does not remain valid for higher pressure ranges as encountered in UHS.

Sigmund (1976) proposed a correction term to the ideal behavior by fitting a third-order polynomial to a set of experimental data. By plotting $ρmDij/ρm0Dij0$ vs. reduced molar density $ρmr$⁠, he determined the correlation equation (see Eq. 6) (Sigmund 1976):

Sigmund’s correction model:

DaSilva et al. subsequently extended Sigmund’s model to avoid negative correction coefficients at high reduced molar densities. This might be applicable to high-pressure, high-temperature systems (Da Silva and Belery 1989). They proposed to introduce an exponential relationship of density-diffusivity product correction at reduced molar densities larger than 3:

Sigmund’s extended correction model:

This approach makes binary and self-diffusion coefficients of various fluid systems, pressures, temperatures, and concentrations comparable. Hence, more literature data have been found and transferred to the plot initiated by Sigmund, as shown in Fig. 1 . Although the experimental data are not perfectly represented by the polynomial, the overall trend of the data is honored.

Molecular density, which is necessary for the $ρmDijρm0Dij0$ calculation, was calculated using the REFPROP (V10) default equation of state for all literature data, where the measured density was not reported. The reduced molar density $ρmr$ is defined as

where $ρmc$ is the mixture critical molar density, which is expressed as a function of molar fraction $xi$ and critical molar volume $Vci$ of each component (see Eq. 10):

In principle, if all the parameters can be calculated (see Eqs. 5, 9, and 10), rearranging Eq. 6 will give an approximation of a binary mixture diffusion coefficient with known initial volumes of each component, system pressure, and temperature:

To assess the validity of the polynomial proposed by Sigmund for the prediction of effective diffusion coefficients in different gas storage reservoirs, bulk diffusion of a binary methane-hydrogen mixture has been measured and compared with the general trend as seen in Fig. 1 . Additionally, to ensure that the measurements of this work are representative, diffusion coefficients of the methane-nitrogen mixture have been identified and compared with the literature and Sigmund’s correction for the kinetic theory.

Although Sigmund included a variety of diffusion data for different gas mixtures including CH₄, N₂, C₂H₆, Ar, nC₄H₁₀, CF₄, SF₆, CO_2, and C₃H₈, he did not consider hydrogen-hydrocarbon gas mixtures. Moreover, there are no reported investigations of hydrogen-hydrocarbon diffusion coefficients at high pressures (above 100 bar). The results of this work, however, show a rapid deviation from the kinetic gas theory at earlier reduced density values. This is also observed from the data reported by Islam and Stryland (1969) for argon-methane diffusion at pressures ranging from 20 to 143 bar. This indicates that the correlation (see Eq. 6) is not representative of all the mixtures and that modifications to the polynomial fit might be required. However, for pressures below 100 bar, the experimental results of hydrogen-methane binary diffusion are well within the kinetic gas theory for ideal gases, where $ρmDijρm0Dij0=1$⁠.

Nitrogen-methane diffusion coefficients reported in this work are in a good alignment with the literature and show a good agreement with Sigmund’s correlation for “dense” mixtures. Upon the validation of the bulk diffusion results with the theoretical calculations, the results of the hydrogen-methane diffusion experiments were carried forward to the evaluation of the diffusion through porous media experiments.

In this section, the experimental setup of the diffusion measurement, the materials, procedures, and data acquisition and analysis are described.

The experimental method used in this study involves a diffusion apparatus that was developed for measuring concentration changes of dissolved gases. Various techniques are being successfully used for detecting concentrations of gas components. In diffusion-related works, the conventional gas chromatography method is primarily used (Sigmund 1976; Jacops et al. 2013). However, optical methods show significant advantages over conventional methods by offering real-time measurements and eliminating the need for sampling. In this work, Fourier transform infrared spectroscopy was put into practice. Coupled together with a broadband light source, this combination offers an effective analysis of multicomponent gas systems with high resolution and sensitivity (Lacome et al. 1984; Chambers 2005; Shemshad et al. 2012).

This technique is fundamentally based on the Beer-Lambert law (see Eq. 12), which correlates light absorbance phenomena by different gas species to their concentrations in the system (Chambers 2005; Shemshad et al. 2012; Chen et al. 2020):

where $It(i)Io(i)$ is the ratio of the incident and detected light intensity, $ε(i)$ is the molar absorption coefficient, $x$ is the gas concentration, and $L$ is the path length.

To interpret spectroscopy analysis and convert it to gas concentration change, an empirical correlation is required. The view cell is used to link the output light intensity received by a detector to the absorption characteristics of the test gas. In other words, the reduction of the recorded power output at a specific wavelength can be directly related to the increase in concentration (partial pressure) of a gas component $i$⁠. To achieve an accurate correlation model, the view cell is pressurized to the expected partial pressure range of the component in a binary or multicomponent system. Measured pressures are converted to concentrations using REFPROP’s (v10) default equation of state.

In the apparatus, two diffusing gases are placed on both sides of a shut-in/shutoff valve (V9) in Fig. 2 . As the diffusion phenomena evolve, the molecules of the test gas initially placed on the left-hand side of the valve (V9) travel through a diffusion chamber. The optical transmission through the diffusant is then measured continuously through the view cell (CT-1) using a Fourier transform optical spectrum analyzer. The response of the light transmission can be then converted to concentrations using the predetermined correlation model.

Different methods have been proposed to measure diffusion coefficients, including but not limited to closed tube/Loschmidt cell, Taylor dispersion, two-bulb apparatus, gas chromatography, and others (Marrero and Mason 1972; Baranski et al. 2003; Jacops et al. 2013). One of the approaches was introduced by Loschmidt in 1870 (Loschmidt 1870). Since the first application, this method has been modified for higher pressures. In this approach, two chambers are separated by a sliding element. Each chamber contains a different gas sample. By opening the sliding element, gases are connected and start diffusing into each other.

The system used in this work (Fig. 2) is a modified version of the conventional Loschmidt cell. These modifications, for example, allow for facilitation of measurements of effective diffusion coefficients through rock samples at confining stress. The setup can be subdivided into gas supply, pressure equilibration, and gas concentration measurement systems.

Two high-pressure piston accumulators (PA-1 and PA-2) containing diffusant gases are connected to a high-pressure metering pump (P-1). The highly accurate pump ensures identical and precise pressures in both piston accumulators. The outlets of both cylinders are connected to the pressure equilibration system that further ensures no differential pressure (ΔP = 0) between the two gas lines.

The gas measurement system includes a high-pressure view cell, broadband infrared (IR) globar (IR light source) light collimation and the Fourier transform optical spectrum analyzer. The gas concentration within the view cell is identified based on the light absorption by the gas molecules theory described earlier. The wavelengths at which the light is absorbed are specific to various gas molecules; thus, the optical spectrometer is ideal for this application.

The diffusion chamber is adaptable to experimental requirements. Fig. 2  illustrates a core holder (CH-1) containing a core sample as the diffusion chamber. Initially, this was replaced by capillary (L = 0.1 m, ID = 1×10^–3 m) to measure bulk diffusion coefficients of methane-hydrogen and methane-nitrogen binary systems and validate the accuracy of the setup.

All the lines, connections, and sealings were specially selected for high-pressure applications involving hydrogen gas.

High priority has been given to a thorough approximation of a probable error in the measurements. As the experiments were performed at no-flow conditions, primary sources that may lead to significant uncertainties include pressure-temperature measurements and transmitted light detection. The maximum expanded measurement uncertainty analysis was conducted to estimate the accuracy of the measured data (Bell 2001). The analysis proceeds in two stages, which are uncertainty analysis of the calibration and uncertainty analysis of the recorded data.

The uncertainty of the correlation equation was evaluated using the expression below:

where $u1(C)$ is the uncertainty of the gas concentration calculated from measured pressure, temperature and known gas component, and $u2(E)$ is the uncertainty of the power measurement due to transmitted light.

Concentration of real gases can be calculated using the equation $C=Pz⋅R⋅T$⁠. Thus, the combined uncertainty of concentration can be calculated using the following equation:

In the equation above, the uncertainty of concentration is subdivided into the root sum of squared uncertainties of pressure and temperature measurements. This is further separated into Eqs. 15 and 16:

Eqs. 14 and 15 express the total uncertainties of thermodynamic conditions as the sum of the standard estimated uncertainties of the measurement devices $u(P)$⁠, $u(T)$ and the concentration gradients due to the pressure-temperature variations $∂C∂P$⁠, $∂C∂T$⁠. The standard uncertainties of the measurement devices reported by the manufacturer were ±0.03 bar and ±0.1 K, respectively. The concentration gradients can be calculated using the following equations:

Gas compressibility factors $z$ required for concentration estimations at upper and lower limits of pressure-temperature conditions were calculated using REFPROP’s (v10) default equation of state.

The second term of Eq. 13 is calculated using the following expression:

The reported uncertainty of the detector is ±3 mV. The gradient due to measurement fluctuations was calculated using Eq. 20:

where $a$ and $b$ are calibration constants obtained from fitting an exponential trend curve to plotted data points (calculated concentrations vs. light detected).

The measurement uncertainty is primarily dominated by the uncertainty of the light detector used during the experiments. Thus, the total uncertainty introduced during the data recording process can be expressed as follows:

where $u(E)$ can be calculated using Eq. 19.

The expanded total uncertainty can be found by multiplying the result of the combined standard uncertainty by a coverage factor $k$⁠. Multiplying it by a factor of 2 to give the confidence level of 95% gives us a maximum expanded uncertainty of 1.05%.

The uncertainty calculated is small for diffusion experiments (1%, k = 2). However, it is important to note that the actual measurements show oscillations in concentrations that are considerably higher (see Figs. 3 and 4 ). One explanation could be that the error margin of one or multiple hardware components is underestimated. Nevertheless, a good fit is obtained to determine diffusion coefficients.

Before starting an experiment, the entire apparatus was pressurized to the maximum measurement pressure with hydrogen to make certain that the system was leakage free. A sufficient period of equilibration time was given for pressure and temperature stabilization. When the system was found to be leakage free, the lines and connections were flushed thoroughly with tested gases before they were brought to the required pressure. The apparatus was placed into an isolated oven with a temperature accuracy of ±0.1°C.

To ensure that concentration measurements, and subsequently, diffusion coefficients were representative, tests were performed through a capillary to measure bulk diffusion coefficients for a quality check. As a benchmark, the literature contains some data on experimental measurements of diffusion coefficients in methane-nitrogen and methane-hydrogen binary systems. Moreover, these can also be calculated analytically using various models. Once the setup passed the quality check, the measurements were performed through different core samples of known properties.

All the core samples were flushed and dried before experiments to eliminate any residual moisture content. Once the sample was ready for the experiment, it was placed into a core holder and pressurized to a certain confining pressure by pumping water at the annulus between the core holder and the sample. When the lines were connected, the gases were brought to the required pressure and temperature. An accurate pump provides minimal pressure difference between the two gas channels. The pressure equilibration system eliminates any residual pressure difference, ensuring an accurate measurement.

Once the gases were equilibrated, the experiment started by opening the valve (V9) and establishing an interface between the two gasses. The light transmitted through the test gas was constantly recorded and saved on a high-speed computer.

The apparatus described earlier measures the concentration change of a component $i$ with time. To determine diffusion coefficients, the experimental data are interpreted by computing forward simulation of a 3D gas transport model. In this model, the diffusive transport through porous media is represented by Fick’s law (see Eq. 1).

The diffusion coefficient $Dij$ of a binary mixture can be calculated analytically or numerically. Analytical solutions of methane-hydrogen and methane-nitrogen diffusion coefficients were discussed in the previous section. However, to interpret molar flux obtained from experiments, this problem needs to be solved numerically.

A 3D numerical model was constructed in COMSOL Multiphysics (v6.0) that represents the geometry of the experimental setup (COMSOL 2022). COMSOL Multiphysics is a finite element method (FEM) based software that has many advantages in solving fluid dynamics problems with complex geometries and fluid mixtures. The accuracy of the FEM calculations is influenced by mesh size (mesh density). Because FEM calculations are performed based on the model variables in each element, better results are generally achieved with finer mesh size. Additionally, timestep has a significant effect on the outcome of the numerical calculations. Therefore, element size and timesteps were analyzed until further refinement had no impact on the results. As such, the timestep was reduced to t = 0.02 hours, with an extremely fine grid size for all simulations.

Fick’s law is implemented in the transport of diluted species physics (tds) of the diffusive transport model. This allows us to incorporate the initial gas concentrations in each domain at corresponding pressure and temperature conditions. Figs. 5 and 6  illustrate the primary part of the model, where the majority of the migration takes place. The gas components interface domain is represented by the shut-in/shutoff valve (V9), where methane is on the left side and hydrogen/nitrogen are on the right. The rest of the geometry in the model is conditioned to have no flux boundaries. As the simulation is initiated, methane molecules will migrate through the diffusion capillary (CT-1 in Fig. 5 ) or rock sample (CH-1 in Fig. 6 ) to the view cell (VC-1). The average concentration profile of methane inside the view cell is then calculated and compared with the experimental results until a match between computed and experimental results is achieved.

The rock sample was implemented in the model using a porous medium subdomain, where the gas transport only occurs in the porous volume. Porosity and other model parameters for the simulations are summarized in Table 1 .

Several tests through diffusion capillary and rock samples were performed according to the procedure described in the Experimental Setup Description section. Fitting the numerical calculations and the experimental profiles led to identifying the diffusion coefficients. The results are presented in Tables 2 through 8 . Some examples of the fitting match are shown in Figs. 3 and 4 . The figures show the best-fit match for methane concentration profiles from the experiments with Berea and Thann rock samples with porosities $φ=20.94%$ and $φ=25.01%$⁠, respectively.

Tables 2 and 3  summarize the results of bulk diffusion through the capillary tube for nitrogen-methane and hydrogen-methane binary systems. To validate the experimental results of the binary gas systems, the diffusion coefficients are also calculated using the equations described in the Molecular Diffusion Summary section. Figs. 7 and 8  depict the comparison between analytical and experimental data points for methane-nitrogen and methane-hydrogen binary diffusion through the capillary, respectively.

The results show that the measured diffusion coefficients using the capillary tube are in agreement with the data in the literature. Deviations in the hydrogen-methane data at higher pressures related to real gas behavior are discussed earlier in the Molecular Diffusion Summary section. The diffusion coefficient is inversely proportional to the pressure. This is expected from the kinetic theory of gases when the mean free path of molecule movement is smaller than the pore (or here capillary size) (Ratnakar and Dindoruk 2022). For hydrogen in the capillary tube at 27°C, the diffusion coefficient is about 20 times larger for low pressures than for high pressures.

The setup of the diffusion cell is described above. A core holder (CH-1) for the core samples is introduced and the concentration increase is measured (diffusion chamber). Each sample was thoroughly preflushed using the gas of the respective experiment to remove any contaminations.

Four cores were used for the experiments (Table 4). Tables 5 through 8  summarize the experimental data of binary methane-hydrogen gas system diffusion through the core samples. The data summary includes the effective diffusion coefficients (⁠$Deff$⁠), which represent the diffusivity through core samples and bulk diffusion coefficients through capillary (⁠$Di$⁠) at the same pressure and temperature conditions. The tortuosity was calculated using Eq. 2.

Fig. 9  summarizes the results of effective diffusion through cores and bulk diffusion coefficients through the capillary tube.

The results (Fig. 9) show that similar to the capillary tube experiments, the effective diffusion coefficient is inversely proportional to the pressure. For the pressure range of 5–40 bar, the effective diffusion coefficient for Berea was seven times higher for the low than for the high pressure.

The effective diffusion coefficients of the measurements, including the rocks, are lower than the diffusion coefficients measured in the capillary tube. For the reservoir rocks used in this study, the diffusion is characterized as bulk diffusion, and most of the molecules collide with each other rather than with the wall (Ratnakar and Dindoruk 2022). The diffusive flux in porous media can be characterized by the tortuosity. The tortuosity is the average length of the true flow paths to the length of the system in the direction of the macroscopic flux (Koponen et al. 1996; Shen and Chen 2007). The movement of species through a porous medium depends on many factors such as porosity, tortuosity, grain shape, and size distribution (Matyka et al. 2008). The tortuosity and the relationship between tortuosity and porosity were evaluated by various authors (e.g., Allaire et al. 2008; Peng et al. 2012). Several theoretical studies were performed to determine the porosity-permeability relations of rocks. The studies included image analysis and deep learning (Graczyk and Matyka 2020), percolation theory (e.g., Ghanbarian et al. 2013), and fractal models (e.g., Luo et al. 2014). The tortuosity as determined in the experiments did not substantially change with the pressures (see Tables 5 through 8 ). This indicates that the mean free path was still short compared to the pore size and the assumption of bulk diffusion is still valid. Fig. 10  shows the measured relation $Deff/Di$ as a function of porosity in comparison with different equations given in the literature (Eqs. 22 through 24) (Grathwohl 1998):

The measured data are similar to the calculated data using Eq. 21 for the lower porosities. The other equations lead to significantly higher or lower values for $Deff/Di$ than the measured data. The reason for the deviation for the high porosity might be a different pore structure; the value of $Deff/Di$ depends on the rock type and pore structure as well (Grathwohl 1998). Here, we neglect the effect of adsorption; however, for rocks containing a larger amount of clay minerals or in case that carbon dioxide is considered as cushion gas, adsorption might need to be considered.

The diffusion distance can be estimated by $L=Defft$⁠. Injection of hydrogen into a depleted gas field for half a year results in a diffusion distance of larger than 3 m from high-permeability layers into low-permeability layers. Dependent on the permeability contrast and spatial configuration of the low-permeability and high-permeability layers, diffusion can be an important factor (Arekhov et al. 2023).

The effects are summarized in Fig. 11 . During the hydrogen injection part of the cycle, hydrogen is moving by advection into higher permeability areas of the reservoir. Through the large interface of high- to low-permeability layers, hydrogen is diffusing into low-permeability layers. The pressures are constantly increasing and the diffusion coefficients decreasing accordingly. During production cycles, hydrocarbon gas and hydrogen are diffusing from the low-permeability layer into the high-permeability layer. Due to decreasing pressures, the diffusion coefficients and fluxes are increasing. Dependent on the reservoir conditions, the diffusive fluxes can be substantial and can influence the quality of the backproduced gas.

Depleted gas reservoirs or gas fields used in underground gas storage might be good candidate fields for UHS operations. The hydrocarbon gas still present in the underground gas storage field could be used as a cushion gas to reduce costs. However, the composition gas backproduced in UHS will be influenced by the mixing of the hydrogen with hydrocarbon gas. The mixing of the gas depends on dispersivity, permeability contrast, and geological setting. The results of laboratory experiments showed the following conclusions:

• Hydrogen and methane gas diffusion can be accurately measured with and without porous medium at various pressures.

• The effective diffusion coefficient increases with increasing porosity and can be described by correlations. The correlations are dependent on the rock type.

• For the investigated reservoir rocks, the effective diffusion coefficient was about a factor 10 lower than the diffusion coefficient without rock.

• The tortuosity determined by the laboratory experiments was only slightly influenced by the reservoir pressure.

• The difference in diffusion coefficients at high and low pressures can be more than a factor of 20.

• The difference in diffusion coefficients with pressure might influence the diffusive fluxes into and from low-permeability into high-permeability zones of the reservoir.

• Diffusion might have an influence on the backproduced gas and hence should be evaluated for potential UHS projects in depleted or underground gas storage fields as surface facilities need to be designed accordingly.

 
• a,b

  apparatus constants

 
• c

  molar concentration

 
• d_p

  pore diameter

 
• D

  diffusion coefficient, m²/s

 
• D_eff

  effective diffusion coefficient, m²/s

 
• D_i

  diffusion coefficient through capillary, m²/s

 
• E

  transmitted light energy, mV

 
• $It(i)Io(i)$

  ratio of the incident and detected light intensity

 
• J_i

  diffusive flux, mol/m²·s

 
• k_B

  Boltzmann constant

 
• Kn

  Knudsen number

 
• L

  path length

 
• M_w

  molecular weight, g/mol

 
• P

  pressure, Pa

 
• R

  gas constant, [m³·Pa/(K·mol)]

 
• t

  time, seconds

 
• T

  temperature, K

 
• u₁(C)

  uncertainty of the gas concentration calculated from measured pressure, temperature and known gas component

 
• u₂(E)

  uncertainty of the power measurement due to transmitted light

 
• x_i

  mole fraction of component i

 
• z

  compressibility factor

 
• ε_AB

  Lennard-Jones energy

 
• $ε$(i)

  molar absorption coefficient

 
• $λ$

  mean free path

 
• ρ_m

  molar density, g-mol/cm³

 
• $ρmc$

  critical molar density

 
• ρ_mr

  reduced molar density

 
• σ

  molecular size (radius), Å

 
• ^$(∑$_A$V^)$

  diffusion volume of Gas A

 
• τ

  tortuosity

 
• φ

  porosity

 
• $Ωi,j$

  collision integral obtained from Lennard-Jones parameters (Lennard and Jones 1931)

 
• $∇C$

  spatial gradient of concentration

Thanks to OMV Energy for the permission to publish the paper.

This paper (SPE 214451) was accepted for presentation at the SPE EuropEC - Europe Energy Conference featured at the 84th EAGE Annual Conference & Exhibition, Vienna, Austria, 5–8 June 2023, and revised for publication. Original manuscript received for review 13 March 2023. Revised manuscript received for review 5 May 2023. Paper peer approved 22 May 2023.