The Role of Diffusion on Reservoir Performance in Underground Hydrogen Storage

The European Union has set a goal to achieve net zero emissions by 2050 (European Union 2021). The target is to reduce net greenhouse gas emissions by expanding the use of renewable energy sources. To achieve this goal, the International Energy Agency predicts a large-scale development in renewable energy business in Europe (IEA 2021). Hydrogen is emerging as a key element of the sustainable energy system of the future (European Commission 2019). The usage of UHS provides flexibility in short- and long-term energy supply with a climate-neutral energy carrier (e.g., Mouli-Castillo et al. 2021; Heinemann et al. 2021a; Clemens and Clemens 2022).

Various types of reservoirs can be utilized for UHS, including salt caverns, aquifers, lined hard rock caverns, and depleted gas reservoirs (e.g., Heinemann et al. 2021b, Małachowska et al. 2022; Delshad et al. 2022; Ali et al. 2022). Currently, depleted natural underground gas reservoirs account for a majority (64%) of the gas storage capacity in Europe (Cihlar and Mavins 2021). Many of these reservoirs have been used for natural gas storage for several decades, but with the shift toward cleaner energy, existing natural gas storage facilities or depleted gas fields can be repurposed for hydrogen storage. This approach offers several advantages, including the use of existing infrastructure, the ability to store large hydrogen volumes, proximity to end customer, lower subsurface risk as the reservoirs were previously used for gas storage, and cost saving for cushion gas (e.g., Olabi et al. 2021; Zivar et al. 2021).

According to Cihlar and Mavins (2021), cushion gas typically comprises 50–60% of the gas initially in place in underground gas storage fields. The purpose of the cushion gas is to maintain reservoir pressure, thereby ensuring that requested production rates can be achieved. When storing hydrogen in depleted gas reservoirs, the remaining hydrocarbon gas can serve as cushion gas. This approach can significantly reduce the capital expenditure necessary for UHS projects.

The injection of hydrogen into natural gas reservoirs can lead to the contamination of the backproduced hydrogen due to mixing (e.g., Feldmann et al. 2016; Pfeiffer et al. 2017; Wang et al. 2022; Maniglio et al. 2022). To understand the extent and process of mixing, it is essential to consider the mass transfer mechanisms at play. When hydrogen is injected into a porous medium, two mass transfer mechanisms occur: advection and diffusion. Advective mass transport refers to the bulk flow of the fluid which is governed by its viscosity. Diffusion, on the other hand, is associated with molecular motion (Bird et al. 2006).

In conventional petroleum engineering for oil fields, advective flow considerations are dominant and diffusion effects are often considered negligible. The reason is that molecular diffusion between gas, oil, and water phases in oil fields is very small and has often little influence on operational parameters and is therefore disregarded. In gas fields, there is typically no component concentration gradient, and diffusive mixing of components does not occur. However, when hydrogen is introduced into the system during storage operations, concentration differences develop, and molecular diffusion becomes significant. As the diffusion of components in gases is several orders of magnitude larger than the diffusion in liquids, a thorough understanding of gas mixing in porous media is crucial for designing efficient UHS systems.

Molecular diffusion has been extensively studied across various industries and the theories described (e.g., Bird et al. 2006; Cussler 2009). Although there is significant literature on gas diffusion (e.g., Gotoh et al. 1974; Chu et al. 1974; Sigmund 1976), there remains a lack of hydrogen-hydrocarbon diffusion studies. In subsurface systems, diffusion processes have been shown to be important in various contexts, including fractured reservoirs, (e.g., Hoteit and Firoozabadi 2009; Yanze and Clemens 2012; Shojaei and Jessen 2014; Ghasemi et al. 2017), gas condensate reservoirs (e.g., Paraschiv et al. 2012), heavy oil (e.g., Ratnakar and Dindoruk 2015; Marciales and Babadagli 2016), shale reservoirs (e.g., Wu et al. 2016; Tian et al. 2021; Bellal and Assady 2022), coalbed methane (e.g., Yi et al. 2006; Liu and Emami-Meybodi 2022), and wettability modification (e.g., Stoll et al. 2008). When modeling flow in a porous medium, the understanding of advective and diffusive fluxes is important as well as dispersivity effects due to local velocity variations and heterogeneity (Ratnakar and Dindoruk 2022). Although several studies have investigated the effects of heterogeneity in UHS (e.g., Lysyy et al. 2021; Mahdi et al. 2021), the coverage of the impact of diffusion on UHS is limited in the literature. Notably, the influence of pressure on diffusion and the reservoir setting has yet to be thoroughly explored.

In this study, we analyze the impact of diffusion on the mixing of hydrogen and methane during cyclic operations in a UHS system. We investigate how diffusion affects the distribution of hydrogen within the reservoir and the composition of the backproduced gas. By conducting a sensitivity analysis of various reservoir properties, we identify the conditions under which diffusion cannot be neglected, providing insights for the design and optimization of UHS systems.

This paper is structured as follows: we first present the methodology used in the study and provide the base-case model parameters. We then describe the impact of numerical dispersion on the simulation results. Next, we discuss the simulation outcomes and demonstrate the significance of diffusive mass transport in heterogeneous reservoirs. Finally, we perform a sensitivity analysis of various parameters on gas mixing, including layer permeability contrast, Kv/Kh ratio, reservoir heterogeneity, binary diffusion coefficient, and pressure-dependent diffusion coefficient. The models are compared based on the hydrogen recovery factor and methane-to-hydrogen ratio during the first two production cycles.

The influence of diffusion on UHS in a depleted gas reservoir was investigated using the fully implicit numerical reservoir simulator Intersect 2022.1 (Schlumberger 2022). To ensure the accuracy of the simulations, we investigated the effect of numerical dispersion by comparing the results of multiple 2D models with varying grid sizes (ranging from 0.1 to 100 m in length). A base-case simulation was conducted to establish a benchmark for comparing the model performance with and without diffusion. To incorporate diffusion into the simulator, we enabled diffusive mass transport and specified the diffusion coefficients for the simulated components. We then performed a sensitivity analysis using a series of numerical models to examine the influence of diffusion on UHS with respect to various input parameters.

In the next section, the calculation of diffusion as used here is described. Afterward, the base-case numerical model is covered.

Pressure-volume-temperature experiments were performed for a mixture of a gas in an existing field and for hydrogen. An equation of state was generated which was used in the simulations reported here. For details of the measurements and equation of state, see Nielsen et al. (2023).

In the simulator, diffusive mass transport is calculated by applying Eq. 1 (Schlumberger 2022):

where:

$Jα,i-$ diffusive molar flux for component $i$ in phase $α$⁠;

$Sα-$ phase saturation of phase $α$⁠;

$Dαi-$ effective diffusion coefficient of component $i$ in phase $α$ (intersect formulation);

$φ-$ porosity;

$bα-$ phase molar density of phase $α$⁠;

$xα,i-$ mole fraction of component $i$ in phase $α$⁠.

As it can be seen from the equation, the diffusive mass flux is associated with composition gradient and scales with the diffusion coefficient.

To simulate molecular diffusion in the simulator, effective diffusion coefficients need to be provided (Ratnakar and Dindoruk 2022). Generally, effective diffusion considers both porosity and tortuosity of the rock. However, in the simulator Intersect, porosity is included in the diffusive mass transport equation (see Eq. 1). Therefore, the value of the effective diffusion coefficient provided to the simulator needs to include tortuosity only (see Eq. 2):

where:

$τ-$ tortuosity;

$D-$bulk diffusion coefficient.

The hydrogen-methane diffusion coefficients used here were estimated with the Chapman-Enskog equation that is based on the kinetic gas theory. Furthermore, it was cross-validated with existing laboratory-measured data reported in the literature (Arekhov et al. 2023; Chu et al. 1974; Fejes and Czarán 1961). Fig. 1  shows the comparison of measured and calculated bulk diffusion coefficients of the hydrogen-methane system. As can be seen from Fig. 1 , the binary diffusion coefficient is inversely proportional to the pressure and directly proportional to the temperature. The pressure and temperature dependency of the diffusion coefficient is associated with the mean free path of molecules and can be explained by the kinetic gas theory (Bird et al. 2006).

To transfer the measurements shown in Fig. 1  to reservoir conditions studied here, the diffusion coefficients were calculated using the Chapman-Enskog correlation. This correlation provided a good fit for the measured data at 25°C, as well as to the literature data. This correlation was also used to estimate the temperature effect on the diffusion coefficient. Fig. 2  shows the calculated bulk diffusion coefficients as a function of temperature at a constant pressure of 50 bar.

As can be seen in Fig. 2 , the diffusion coefficient does not show a strong dependency on temperature. The magnitude of the change of the diffusion coefficient with temperature is small in particular as temperature effects are expected to occur in the near-wellbore area which is advection dominated. Owing to the limited effect, an isothermal model was used here.

For simplicity, the base-case numerical model does not include the pressure dependency of the diffusion coefficient. The pressure dependency of the diffusion coefficient is covered in a sensitivity analysis.

Another important aspect that needs to be considered to ensure comprehensive understanding of the subsurface gas mixing is multicomponent diffusion. The composition of the underground hydrocarbon fluid can be very different depending on the reservoir and injected gas composition. In this work, we consider UHS in a depleted dry gas field, containing 96% methane and minor amounts of impurities. To understand the possible effects on the diffusion process caused by the presence of impurities, the DaSilva method [see da Silva and Belery (1989) for the details of the method] was used to calculate multicomponent diffusion coefficients. It was decided to perform a sensitivity analysis on the ternary system to reduce complexity, namely, hydrogen diffusion into the methane-ethane gas mixture. In Fig. 3 , the molar fraction of hydrocarbon components was changed and the diffusion coefficients were calculated.

As can be seen in Fig. 3 , adding 5 molar percent of ethane in the mixture does not influence diffusion coefficients significantly. In this case, the hydrogen diffusion coefficient is reduced by 2%. Thus, to reduce model complexity, it was decided to proceed with hydrogen-methane binary diffusion.

At the same time, it is crucial to draw attention to the importance of multicomponent diffusion in the case of hydrogen storage in wet gas or gas condensate reservoirs. As can be seen from Fig. 3 , a larger concentration of heavier hydrocarbon molecules can significantly alter diffusion coefficients of all components.

In the end, using the Chapman-Enskog correlation, the bulk diffusion coefficient at reservoir conditions (41.5°C) was estimated for the case reported here as 0.133 m²/d. However, before implementation in the simulator, this value must be adjusted to account for the tortuosity of the rock. The tortuosity of rock can vary widely depending on its internal structure ranging from 1.2 to 4 (e.g., Donaldson et al. 1976; Arekhov et al. 2023). In this study, we used a tortuosity of 1.7 for the sandstone being studied. This results in an effective diffusion coefficient of $Dαi=0.08m2/d$ which is used in the base-case numerical model. Geomechanical effects such as changing porosity and permeability with pressure might have an effect on the effective diffusion coefficient and could be the subject of additional studies.

Here, we are using binary diffusion; therefore, the same diffusion coefficient was introduced for methane and hydrogen as suggested for binary mixtures by Morel et al. (1993).

A small-scale 2D model was generated to study the influence of diffusion on the UHS storage process. The model parameters are described in Table 1 . The grid is shown in Fig. 4 .

To model UHS, no-flow boundary conditions were applied to all edges of the model. The injection and production rates were selected to simulate reservoir pressures ranging from 10 to 50 bar. We took care that well productivity did not impact the results of any of the scenarios. To accomplish this, we ensured that all target rates were reached. Fig. 5  shows the reservoir and well bottomhole pressure response of the base-case model.

To evaluate the molecular diffusion effects for more complex reservoir simulation models, a full-field reservoir simulation model was created. The model incorporates fluid flow in all three special dimensions, the stratigraphy, and the heterogeneous distribution of reservoir parameters. Structurally, the reservoir is a four-way dip closure. The modeled reservoir consists of two distinct reservoir zones: a thin high-permeability layer on the top as well as a thick low-permeability zone below. Fig. 6  displays the structure of the reservoir including zonation (left) and permeability distribution (right).

The dynamic reservoir model was history-matched for the hydrocarbon production period, but the details are not given in the paper. The main reservoir properties are listed in Table 2 . As the objective of this study is the investigation of diffusion effects within the reservoir, the caprock was not modeled and a no-flow boundary is applied on the top of the high-permeability layer.

In this section, first, numerical dispersion effects are described. Then, the results of a base-case model are given followed by a parametric study and the 3D model results.

A uniform permeability of 100 md in the x-direction was used to study numerical dispersion. Kv/Kh was set as 0.01. The grid cell size in the x-direction varied from 0.1 to 100 m. No molecular diffusion was implemented for the cases described here to focus on numerical dispersion. Fig. 7  represents the hydrogen concentration profile after 3 months of injection for various grid sizes.

In the theoretical case, a sharp front between hydrogen and methane would exist. However, a significant numerical dispersion can be observed even in the case of using a small gridblock size of 0.1 m. Numerical dispersion results in smearing out of the concentration front. Numerical effects become more prominent with increased grid size. Refining the grid size leads to decreasing numerical dispersion but a significant increase in computation time.

In this paper, we do not discuss the effect of reservoir dispersivity vs. numerical dispersion as this is outside the scope of the paper and requires an analysis of reservoir dispersivity and numerical dispersion as described in Gelhar et al. (1992) and Appelo and Postma (1991). The grid size that should be used needs to account for the internal geological structure within a numerical gridblock, local dispersivity, and numerical dispersion and will be case-specific (e.g., Arya et al. 1988; Coats et al. 2009; Garmeh and Johns 2010). Here, we use a gridblock size of 0.25 m and 40 grid blocks in the vertical direction to limit numerical dispersion but to achieve acceptable simulation times.

The numerical model parameters are described in the “Methodology and Model Setup” section. Two simulation runs of hydrogen injection into a methane-containing reservoir were performed. One case covers the advective transport without diffusion. The second case includes the simulation of diffusion. Fig. 8  shows the hydrogen and methane production profiles during cyclic operations. Fig. 9  displays the comparison of the hydrogen concentration profile at the end of each cycle.

The importance of diffusive mass transport becomes apparent when analyzing the propagation of hydrogen in a layered system. With diffusion, hydrogen penetrates deeper into the reservoir (Fig. 9 , left) than for the case without diffusion (Fig. 9 , right), resulting in a larger reservoir volume being utilized. The mass flux in the vertical direction into the low-permeability layer is significant reaching the top of the reservoir. The propagation of hydrogen into the low-permeability layer is diffusion dominated. The diffusion front is also smoothened when diffusion is included.

When diffusing is neglected, most of the hydrogen enters the high-permeability layer (Fig. 9 , right). This causes an increase in pressure in the lower zone, which leads to a hydrogen flux toward the upper zone. Due to the low permeability of the upper zone, only a small portion of the hydrogen can reach the farthest parts of the reservoir, where methane is effectively trapped.

During the production cycle, contrary effects occur. Methane becomes more mobile if diffusive flux is included compared with the simulations without diffusion (Fig. 9). Thus, the methane breakthrough in the production well occurs much faster compared to the case without diffusion. Some of the hydrogen remains in the low-permeability zone at the end of the production cycle because of the slow advective process. When diffusion is activated, the methane-to-hydrogen production ratio is 13.5% and 4.8% for the first and the second production cycle, respectively.

During the injection cycle, the pressure distribution governs the behavior of the production cycle in the absence of diffusive flux. The main well inflow occurs from the high-permeability zone, which is saturated with hydrogen. Therefore, the initial production stream has only trace quantities of methane when diffusion is not activated. The pressure in the lower zone decreases, which initiates the flow of methane from the low-permeability zone. The methane breakthrough occurs later in time than for the case if diffusion is included. However, methane production rises much faster as soon as methane reaches the production well in the high-permeability layer. For the case without diffusion, the methane-to-hydrogen production ratio is significantly lower than in the absence of diffusion. The methane-to-hydrogen ratio without activating diffusion is 5.5% and 3.7% for the first and the second production cycles accordingly.

The results of the base case indicate a significant impact of diffusion on the distribution of hydrogen within the reservoir. This, in turn, affects the purity of the backproduced hydrogen during UHS operations in depleted gas reservoirs. Ignoring molecular diffusion effects might lead to the erroneous prediction of the backproduced hydrogen quality as well as the onset of contamination of the backproduced hydrogen with methane. The contamination occurs in the case without diffusion after 3 months, whereas in the case with diffusion, the backproduced gas is almost immediately contaminated with methane.

In the following sections, the various parameters that have an influence on the mixing are investigated.

The model was initialized with constant permeability of 1,000 md in the y-direction. Kv/Kh is 0.01. In the case of uniform permeability, the diffusion does not influence the hydrogen production profile and slightly changes hydrogen propagation and distribution (see Figs. 10 and 11 ). Diffusion in this case only leads to countergravity flux and stabilizes the vertical front. When the permeability of the reservoir is high, advection becomes dominant over diffusion. Diffusion effects can be neglected in such conditions.

Smearing of the front mainly occurs due to the numerical dispersion. Due to large advective flux, the diffusion does not contribute significantly to the mass transport. Thus, diffusion plays a major role only in the direction of low permeability (vertical and horizontal).

The layered model contains two layers with a large permeability contrast. The lower layer always has a permeability of 1,000 md. The permeability of the upper layer varies between 1 and 100 md. Fig. 12  shows the hydrogen recovery and methane to hydrogen production ratio during the first two injection/production cycles with respect to the permeability contrast.

Fig. 12  shows that with increasing permeability contrast, the effect of diffusion becomes more pronounced. When diffusion is not considered, there is a small difference in hydrogen recovery factor as well as methane-to-hydrogen production ratio during the two cycles. This occurs because most of the injected hydrogen propagates into the high-permeability zone. In such a case, only a part of the reservoir is efficiently used for storage. There is only a minor advective flux between the two layers.

In contrast, when diffusion is taken into account, the hydrogen propagates into the low-permeability zone owing to the concentration differences leading to diffusion fluxes. During production, some of the injected hydrogen stays in the low-permeability zone and is not produced back. At the same time, methane flux toward the high-permeability layer is increased for the case with diffusion compared to neglecting diffusion and higher methane production is observed. However, during the second cycle, there is significantly less mobile methane in the reservoir, which reduces methane production and leads to a higher hydrogen recovery factor.

This model setup is the same as the base-case scenario with variable Kv/Kh ratio ranging from 0.005 to 0.5. Interestingly, the anisotropy does not affect the hydrogen recovery factor or hydrogen purity (see Fig. 13 ). It also has a slight influence on hydrogen distribution in the reservoir (see Fig. 14 ).

The reason is that, in the case of a lower vertical permeability, the longitudinal propagation of hydrogen is further, especially in the low-permeability zone than for the cases with higher vertical permeability. The simulation cases used a constant injection rate. No influence of the well performance is present. This is the reason for the small sensitivity of anisotropy on the hydrogen recovery factor and methane-to-hydrogen production ratio. The same amount of gas is injected/produced in all simulation runs. In the case of a smaller vertical permeability, higher mass flux is observed in the longitudinal direction. It also leads to higher well bottomhole injection pressure as the well deliverability is reduced. Anisotropy does not influence the production as well as mixing of gasses in this specific setup. However, it is important to keep the anisotropy as a parameter in mind, as this parameter might be important in full-field simulations.

To investigate the aspects of reservoir heterogeneity, a multilayer reservoir system was modeled. The parameters of the models are the same as in the base case with the only difference in the number of layers. Low- and high-permeability layers are alternating sequentially with a total of two, four, eight, and ten layers.

As can be seen in Fig. 15 , the heterogeneity of the reservoir significantly influences hydrogen production and hydrogen purity for both cases with and without diffusion. Also, the hydrogen distribution in the reservoir is significantly affected by heterogeneity (see Fig. 16 ).

We can observe a large difference in hydrogen distribution in the reservoir when the number of layers is increased. When diffusion mass transport is active, increasing the number of layers leads to an almost homogeneous reservoir behavior. First, hydrogen enters the high-permeability zones. After that, hydrogen migrates into low-permeability layers owing to diffusion. In the 10-layer system, the area of the low- and high-permeability zones interface is significantly higher. This increases the component interchange into the adjoined zones. Consequently, methane is being more efficiently displaced from the near-wellbore region, which leads to an improvement in hydrogen purity and recovery factor.

For the case without diffusion, the simulation outcome is the opposite to the cases with diffusion. With an increasing number of layers, the distance that methane migrates into high-permeability layers is reduced. Thus, the amount of mobile methane near the wellbore is increased, which leads to higher methane-to-hydrogen production ratio as well as a decrease in the hydrogen recovery factor.

These cases illustrate the importance of diffusion mass transport when gas mixing effects are simulated in heterogeneous reservoirs.

These simulations were performed to study how the binary diffusion coefficient affects hydrogen storage operations.

The molecular diffusion coefficient is an uncertain parameter. It is possible to measure binary diffusion in the laboratory accurately (e.g., Arekhov et al. 2023). However, when applied to reservoir conditions, the effective diffusion coefficient might vary with porosity, tortuosity, pressure, water saturation, and gas composition. To understand the sensitivity of the diffusion coefficient on UHS, a set of dynamic models were simulated with diffusion coefficient ranging from 0 to 0.8 m²/day. Fig. 17  shows the hydrogen recovery factor and methane-to-hydrogen production ratio for the base-case model.

As can be noticed, the response on hydrogen recovery is nonlinear. When diffusion is equal to zero, mainly the high-permeability zone is exposed to cyclic operations. Hydrogen is injected and produced back without a large inflow of methane, which stays in the low-permeability zone.

When the diffusion coefficient is below about 0.1 m²/day, the recovery factor is reduced with increasing values of the diffusion coefficient. Although more hydrogen propagates into the low-permeability zone, the diffusion process is still slow. Hydrogen does not displace methane far from the well, which allows methane flux during production. The production is governed by advection in the high-permeability layer.

At a diffusion coefficient of about 0.12 m²/d, the diffusive flux becomes dominant. Hydrogen becomes mobile enough to displace methane from the low-permeability layer during injection. Therefore, during the production cycle, a smaller amount of methane is recovered than for smaller diffusion coefficients.

So far, we were assigning a single value of the diffusion coefficient to all simulation grid blocks. However, as it can be seen from the experimental results (Arekhov et al. 2023), the diffusion coefficient is a pressure-dependent parameter (see Fig. 1 ). As the simulations cover a pressure range from 10 to 50 bar, one can expect a higher diffusion mass flow in the reservoir at lower pressures. The simulator that was used allows us to introduce pressure-dependent diffusion coefficients. Fig. 18  displays the component production rate comparison between three cases: without diffusion, constant diffusion, and pressure-dependent diffusion.

Fig. 19  shows the diffusion coefficients in the grid blocks for the pressure-dependent simulations. Using pressure-dependent diffusion coefficients slightly changes the gas mixing process. Taking pressure dependency into account alters hydrogen distribution in the reservoir as well as the production behavior. In the beginning of the production cycle, the diffusion coefficients are small due to higher pressures. Therefore, the mobility of methane is poorer, and, as a result, it is produced at lower quantities compared with the case with a constant diffusion coefficient. On the other hand, the lower the reservoir pressure, the more prominent the effects of diffusion. The mobility of methane is enhanced at low pressures. Thus, the production of methane is accelerated in the end of the cycle compared to the case with the constant diffusion coefficient.

This section details the simulation outcome from the 3D simulation model described in the “Methodology and Model Setup” section of this paper. Two dynamic simulation runs, with and without diffusion effects, were performed and compared. To simplify the interpretation, a constant diffusion coefficient for hydrogen-methane binary diffusion was used. Due to the heterogeneous property distribution in the 3D model, the production profiles differ from the results of 2D simulation shown earlier; however, the main observations still hold. Introduction of diffusive mass flux leads to increased mixing in the subsurface and, consequently, higher methane composition in backproduced stream. Fig. 20  shows the component production profile of hydrogen obtained from the 3D simulation.

According to the simulation results, the methane concentration at the end of the first production cycle differs by 2.5% between cases with and without simulating diffusion. During the second cycle, the methane concentration difference is reduced to 1.3%. It is important to point out that the composition of backproduced stream depends on the applied development strategy. In this simulation, we produced half of the injected volume to ensure that no wellbore limitations were reached during the production period. The effect of diffusion might be different if another well operating mode is selected. The effects of diffusion might be more pronounced for other reservoir settings as shown in the parameterization study.

Besides composition of backproduced hydrogen, diffusion is also affecting hydrogen distribution in the reservoir and the contacted pore volume (Fig. 21 ). Fig. 21  illustrates the reservoir volume saturated with hydrogen (filter 0.01 mole fraction of hydrogen) at the end of the second production cycle. One can see that hydrogen propagates vertically toward the low-permeability layer if diffusion is included, whereas in the case when no diffusion is considered, hydrogen spreads more horizontally.

UHS into depleted gas fields can make use of the remaining hydrocarbon gas as cushion gas. However, in particular, in the early cycles, the backproduced hydrogen will be contaminated with hydrocarbon gas. Simulations including diffusion and comparing them with simulation results neglecting diffusion show the following:

• Numerical dispersion might lead to significant simulation errors when predicting gas mixing. Optimum gridblock size shall be determined taking the reservoir dispersivity into account.

• Molecular diffusion plays an important role in the process of hydrogen storage in depleted gas reservoirs. Including diffusive mass transport, the gas contained in low-permeability zones is becoming more mobile. This leads to larger accessible pore volume for hydrogen, but, at the same time, higher content of methane in the backproduced gas stream.

• Diffusion affects the hydrogen distribution in the reservoir. It smears the concentration front out allowing further propagation of hydrogen into low permeability zones.

• Diffusion effects are pronounced in case of large variation in lateral and vertical permeability.

• A sensitivity analysis showed a strong dependency of the back-produced hydrogen purity on reservoir permeability and heterogeneity and on the value of the diffusion coefficient. The anisotropy did not influence the production profile significantly. However, anisotropy can be an important parameter when simulating large-scale 3D geological models.

• The diffusion coefficient is a pressure-dependent parameter. When simulating UHS, it might be important to account for changes of diffusion coefficients due to reservoir pressures during injection/production cycles dependent on the reservoir setting and operating parameters.

• The differences in the amount of backproduced hydrogen can amount to 9% if diffusion is neglected in the simulation.

• The onset of methane contamination in the backproduced hydrogen can be underestimated by almost half a cycle if diffusion is neglected.

• 3D full-field reservoir simulation results reveal that molecular diffusion is affecting hydrogen distribution in all spatial dimensions and is important to consider while modeling composition of the back-produced fluid stream.

 
• $bα$

  phase molar density of phase $α$

 
• C

  concentration

 
• $C1$

  highest concentration

 
• $C2$

  lowest concentration

 
• D

  bulk diffusion coefficient

 
• $Dαi$

  effective diffusion coefficient of component in phase $α$

 
• $Jα,i$

  diffusive molar flux for component in phase

 
• $Sα$

  the phase saturation of phase

 
• t

  time

 
• x

  distance

 
• $xαi$

  mole fraction of component in phase

 
• $φ$

  porosity

 
• $τ$

  tortuosity

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

This paper (SPE 214435) was accepted for presentation at the SPE 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 6 February 2023. Revised manuscript received for review 21 June 2023. Paper peer approved 26 June 2023.