Bioreaction coupled flow simulations: Impacts of methanogenesis on seasonal underground hydrogen storage

Assessing microbial risks is key to feasible hydrogen storage in geological formations. This work quantitatively analyses the impacts of bio-methanation on hydrogen storage performance. Fine-scale flow simulations, coupled with the bio-methanation reaction, are presented to analyse its impact on the storage performance. Based on the reported rates in literature, methanogenesis may slightly degrade the recovery performance of hydrogen but is considered minor compared with the issue of gas mixing. The impacts of methanogenesis on a time scale of months (330 days) becomes observable in the system configured here, when the methanation rate is above 1746 nano molality per hour. The assumed methanation rate is two times greater than the rate reported from the Olla filed. Validated scaling theory generalises findings for gravity-dominated scenarios. But viscous-dominated flows see complications from property variations due to pressure changes at high rates. This study provides definitions of “ target properties ” (e.g., acceptable methanogenesis rates) for screening hydrogen storage projects.


Introduction
H 2 can play a vital role in the future energy mix, as it may provide a viable solution to address the intermittency challenge associated with renewable energy sources [1,2].Excess renewable energy can be used to produce H 2 through electrolysis, and temporarily stored [3].When energy demand increases, the stored H 2 is withdrawn and combusted to produce electricity or heating [4].Currently, the value chain of H 2 used as an energy vector has not been fully established [5].The scale-up of H 2 energy to advance the Net-Zero target still faces several technical challenges.For instance, H 2 is the lightest gas and necessitates temperatures below − 252.87 • C to be liquified [6].This means the volume of H 2 involved in the future energy mix would be enormous.Therefore, the need for a safe and efficient storage approach becomes critical.
Subsurface porous media potentially provides low-cost storage capacity at scale (gigawatts to terawatts) for H 2 storage, particularly for accommodating seasonal variations in energy demand [7,8].However, the reliability of subsurface H 2 storage in porous media is not yet proven [9,10].One of the central concerns is the potential microbial risks and how they could negatively impact operational performance.This is because H 2 is a superb electron donor with the potential to trigger various microbial reactions [11,12].Microbial communities are ubiquitous in subsurface rocks when conditions including temperature, pressure, salinity, and nutrient availability are favourable [13].For example, H 2 can be consumed by sulfate-reducing microorganisms, which may lead to variations in pH and the generation of hydrogen sulfide [14,15].In addition, H 2 can be converted to methane (CH 4 ) via the metabolism processes of methanogen [16,17].These metabolic processes can impact fluid compositions, reservoir geochemistry and flow properties, potentially leading to operational issues.Detailed investigation into the impacts (and mitigation) of microbial risks is therefore required.
Although CH 4 is currently still a major energy source, biomethanation is increasingly seen as an unfavourable process for the purpose of underground H 2 storage when considering the ultimate goal of achieving net-zero emissions.There is an increasing number of research studies on bio-methanation process in the context of underground hydrogen storage [6,18,19].Ebigbo et al. [20] investigated the governing mechanisms of conversion of H 2 to CH 4 based on a coupled pore-scale model.Farajzadeh et al. assessed the exergetic efficiency and carbon footprint in the process of subsurface bio-methanation [21].As seen in Equation (1), methanogenesis or bio-methanation entails the generation of CH 4 from substrates of H 2 and carbon dioxide (CO 2 ) during the metabolic process carried out by methanogens [16,22,23].This is particularly relevant when storing H 2 in depleted hydrocarbon reservoirs or when CO 2 is used as cushion gas to provide pressure support and prevent water breakthrough during the H 2 withdrawal process [24].Depending on the end users of H 2 , the requirement for the H 2 purity can be strict [25].According to the standardization report published by ISO [26], the hydrogen purity for combustion is requited be above 98 %.Therefore, there is great necessity to understand the extent of methanogenesis, which may lead to H 2 loss and the impurity of CH 4 into the storage system.
In fact, there is very limited data available in the literature on the reaction rate of bio-methanation at the scale of a potential H 2 storage reservoir [27,28].Thaysen et al. [29] reported an estimated rate of 0.02-1205 nano molarity per hour (nm h − 1 ) of H 2 consumption at reservoir scale from their review of laboratory results.Tyne et al. [30] used isotope methods to estimate the methanogenesis in the Olla oil field, where CO 2 was injected for incremental oil recovery purposes.The hydrogen source was from the hydrocarbon biodegradation.They claimed that at least 13-19 % of the CO 2 retained in the reservoir had been converted to CH 4 via methanation reactions over 29 years.Note that the methanation rate reported in the latter study is one to two orders greater than that reported in the former one, as correlated in calculations of this work.
It is not surprising to see significantly different forecasts of the methanogenesis rate in the subsurface, as there is considerable uncertainty over the various factors which impact the rate of the biomethanation process [31].To address such uncertainties, instead of applying a single value of methanation rate into flow simulations, the impacts of methanation process were systematically assessed with a range of correlated reaction rates.The objective was to gain insights into the conditions under which the methanation process becomes significant and to identify instances where its impact can be considered negligible.
Note that the focus of this study is not to investigate the controlling factors of the methanation reactions, such as the reservoir temperature and aquifer salinity.Instead, the significant uncertainty over the fate of microorganisms is acknowledged and therefore a range of possible scenarios using numerical simulations are tested in this work.The Arrhenius formula, which is an empirical relationship based on a rate constant and widely used for chemical reactions, is incorporated in flow simulations [32].This means a set of reaction rates are tested to determine if the methanations processes become important or are negligible in the systems configured here.To the best of knowledge, there are very few quantitative studies that address the impact of biological reactions coupled with flow dynamics during hydrogen storage in subsurface porous media.We acknowledge that utilising Arrhenius equations is a simplification that does not fully capture the complexity of microbial activity.However, this simplification reveals key interactions between hydrodynamic flow and bio-methanation processes, and their influence on overall storage and recovery performance.In another of our papers [33], we explored bio-reaction kinetics using a dual Monod model to more properly reflect the biochemistry.We are currently working towards developing coupled simulations including hydrodynamic flow behaviour, microbial activity and geochemical reactions.Our vision is to leverage both simplified representations and rigorous mechanistic modelling to develop a systematic understanding of the complex dynamics governing subsurface H 2 system.
In the previous work [34], the focus was to analyse the hydrodynamic behaviour in various flow regimes and recognize the key factors influencing hydrogen recovery.Here, that study is extended to quantitively estimate the impacts of bio-methanation process during injection and withdrawal of H 2 .To do so, the methanation reaction rates are at first correlated with data in the literaturementioned above based on 1D models [29,30].Then, the interactions between flow behaviour (viscous instability, gravity segregation, fluid mixing etc.) and methanogenesis are analysed in the 2D cross-sectional heterogeneous system.With the use of very-fine scale simulations, the impacts of the methanation process on the recovery performance (H 2 purity, water-to-gas ratio, H 2 recovery etc.) are tracked and explained in detail.
To generalise the findings for other systems of different sizes and dimensions, the applicability of scaling theory is successfully validated in the context of H 2 storage with methanation process considered.This is for the first time that scaling theory has been applied in this context.The advantages and the limitations of the scaling theory applied in this work are also briefly discussed.This numerical analysis has identified situations where the impact of methanation process become relevant on engineering timescales (months -years).The outcome of this work provides useful insights for the experimentalists and the industry in screening projects of hydrogen storage in subsurface porous media.

Hydrodynamic flow modelling
Fig. 1 shows the 2D vertical heterogeneous permeability field which was generated to trigger the flow features of interest (gravity segregation, fluid mixing, viscous fingering etc.) based on Petrel [35].Since the process of H 2 storage involves multiple fluid components (H 2 , CO 2 and CH 4 ), full compositional simulations were conducted using CMG/GEM  G. Wang et al. [36].The input for the hydrodynamic flow functions, including the Equation of State (EoS), two-phase relative permeabilities, capillary pressure and gas solubilities, are inherited from previous work [34,37].
Here, the main modelling methods and the data sources used in this work are only summarised here.See Table 1.

Modelling of methanogenesis-arrhenius reactions
The Arrhenius equation (Equation ( 2)) is applied in this work to model the methanation process [36].
where r is the reaction rate (molality/sec); F is the frequency factor (1/ (sec ⋅ molality)); Ea is the activation energy for the reaction; R is gas constant; T andT * are reservoir temperature and laboratory reference temperature, respectively; c1 and c2 are molality of H 2 and CO 2 , respectively.To simplify the system, it is assumed that there is no activation energy required for the reaction and E a is set to 0. F is varied to correlate with the reaction rates in the literature.Since the brine salinity is not included here and water density throughout the system is close to 1000 kg/m 3 , the unit of molality and molarity is interchangeable in this work.See more discussion on the correlation of methanation rates in the section of Results and Discussion.

Single-cell numerical analysis
A single-cell model was generated to represent a closed vessel for biomethanation reactions.Without any fluid injections or productions, this closed system was initialized with a gas zone consisting of H 2 and CO 2 at irreducible water saturations.By varying the frequency factor (F) in Equation ( 1), the reaction rates in single-cell simulations were matched with the data published by Tyne et al. [30] and Thaysen et al. [29].The correlated reaction rates were then used for the following 2D numerical analysis.A new synthesis, which incorporates complex multiphase flow behaviour and the process of methanogenesis is therefore developed here.See further discussions of these issues in the result section.

2D vertical cross-sectional simulations
Table 2 presents a simple seasonal operational strategy of H 2 storage which is utilised here to provide a testbed for analysing the joint impacts of the hydrodynamic flow behaviour and the methanation reaction on the storage/withdrawal performance of H 2 .The whole operation process can be split into 4 stages: CO 2 injection (the "cushion gas"), H 2 injection, a shut-in (or "no flow") period and finally H 2 back-production.The volume ratio of CO 2 (cushion gas) to H 2 (working gas) is approximately 2:1 under reservoir conditions, considering the difference in gas solubility between CO 2 and H 2 .In the base case (Case A), a gravity-dominated scenario (330 days in total) is designed to mimic the seasonal operations of H 2 storage and back-production.Also, the rates of injection and back-production are increased by 10 times to trigger the viscous instability in Case B, with the aim of reflecting the practical concerns over the fluctuations in the gas supply and energy demand.Both the injector and the producer were horizontally perforated across the top and controlled by the flow rate under reservoir conditions.The length of the shut-in period is 100 days in both cases.

Results and discussion
A number of simulations were conducted to investigate the interactions between the hydrodynamic flow behaviour and the methanation reactions, in terms of their joint effect on H 2 production, purity etc.To clarify the purpose of each step of the analysis and the logic of this research, the results are presented in the following four steps (#1 to #4), as explained in outline in Table 3.

Step #1: correlation of reaction rates using single-cell simulations
As mentioned above, Tyne et al. [30] reported that 13%-19 % of the retained CO 2 in Olla field was converted into CH 4 via methanogenesis in 29 years (10585 days).Therefore, the component ratio of CO 2 to H 2 in moles is set to 1.86:1 in the single-cell system.This is to ensure that up to 13 % of CO 2 can be converted to CH 4 , if H 2 is completely consumed within a time length of 29 years.Note that the isotope method utilised in the study by Tyne et al. [30] can estimate the bio-methanation rates quantitively at the reservoir scale.Further results are still required to indicate if or when the conversion was actually completed.In other words, the methanation process could be still an ongoing process or have been completed well before the final field measurement (at t = 29 years).To address such uncertainties, a set of reaction rates using various frequency factors (outlined in Table 4) were tested in this work.Fig. 2 shows the gas mole fraction of CH 4 throughout 29 years (10585 days) in Cases 1-5 with rates in a descending order.As seen in Fig. 2, the CH 4 mole fraction increases sharply at the early stage and its subsequent production flattens out over time as shown in Cases 1-4.This means that H 2 has been completely consumed in these cases and approximately 13 % of CO 2 has been converted by the end of 10585 days.Obviously, the greater the value of the input conversion rate, the quicker the conversion is completed.Note that Case 4 (F = 5e-8) is the lowest possible input rate to complete the target conversion (13 % CO 2 ) in this system.With an even lower input of reaction rate in Case 5, the profile of CH 4 mole fraction increases almost linearly with time but in a much more gradual slope compared with that at the early stage of Cases 1-4.The methanation process is still ongoing in Case 5 and there is H 2 remaining at the end of simulation, which means less 13 % CO 2 is converted (9 %).
On the other hand, Thaysen et al. [29] reported the rates of 0.02-1205 nm h − 1 of H 2 consumption at reservoir scale.Since the methanation process entails 4 units of H 2 in a single reaction, the reaction rate is therefore ¼ of the conversion rate of H 2 , viz.0.005-301.25nm h − 1 .With a trial-and-error method, the frequency factor is set to Henry's law [43] Note that the CH 4 solubility is ignored for simplicity, as the CH 4 dissolution in water is very minor and does not trigger any further reactions.Tests with the solubility of CH 4 considered have been conducted, and no evident changes in flow performance or reactive behaviour was observed.1.4e-9 to achieve a methanation reaction rate of 296 nm h − 1 (Rate 6 in Table 4).The time profile of CH 4 mole fraction of Case 6 is plotted together with other cases in Fig. 2, with the aim of the showing the possible range of methanation rates at field scale published in literature.

Step #2: baseline of hydrogen recovery without methanogenesis in 2D cross-sectional simulations
Before the correlated reaction rates are incorporated in the 2D crosssectional simulations, it is important to first recognize the key flow behaviour occurring in the system configured here.That is, some-zero methanation simulations are generated as baseline cases with which to compare the later cases with the methanation process occurring.Fig. 3 is plotted to show the gas saturations (top) and gas mole fraction of H 2 (bottom) at three time points: the end of H 2 injection, the end of the shut-in, and the end of the back production in Case A (left) and Case B (right), respectively.As indicated by the white lines in Fig. 3, the gas displacing front (G/W) in Case B advances further than that in Case A at the end of gas injection.Since the injection rate of Case B is 10 times of that in Case A, the viscous instability becomes evident and leads to the fingering flow of gas (CO 2 ) into water.The flow features of viscous fingering have been well recogised in the context of CO 2 storage and the hydrocarbon industry [44,45].In the meantime, the less viscous H 2 also infiltrates the underlying cushion gas of CO 2 in case B at the end of H 2 injection.This has led to a much expanded mixing zone at the end of the shut-in period, where the fluid redistribution is driven by gravity force.In contrast, H 2 almost uniformly accumulates at the top of the reservoir at the end of H 2 injection in Case A. The size of the mixing zone is restricted and much smaller than that of Case B at the end of shut-in period (indicated by the black bracket).
Fig. 4 is plotted to compare the water-to-gas ratio under reservoir conditions and H 2 purity of the produced gas throughout the back production period.As seen in Fig. 4, the water-to-gas ratio remains negligible in Case A throughout the process.However, significant water breathrough occurs at the later stage of withdrawal in Case B. This is because the viscous fingering of the cushion gas (CO 2 ) occurring in Case B has led to more extensive fluid spreading and thus the increased CO 2 dissolution into water.In the meantime, the upward flow of gas driven by gravity leads to more trapped gas in the bottom area of Case B. As a result, the total volume of the mobile cushion gas is much reduced.This phenomenon compromises the function of the cushion gas, which is to prevent water breakthrough during back production.In the system configured here, the amount of cushion gas is adequate for the gravitydominated scenario (Case A) but not for Case B, where viscous force is also important.At the same time, the H 2 purity in the produced gas also decreases more rapidly in Case B than that in Case A. This is because a much greater degree of mixing between H 2 and CO 2 occurs at the end of shut-in period in Case B due to the flow of viscous fingering.For this reason, the back produced gas in Case B experiences a much quicker reduction in H 2 purity than that in Case A.
Fig. 5 is plotted to compare the H 2 recovery against the pore volume production (left).With the same pore volume of production (0.075 PV), the H 2 recovery in Case B is 10 % lower than that in Case A, due to the     joint impacts of water breakthrough and fluid mixing.According to the International Organization for standardization, the H 2 purity required for combustion is 98 % [26].And the general purifying capacity of a processing plant (based on ammonia) can normally tolerate up to 10 % of CO 2 content [46].Therefore, the H 2 recovery against these two threshold values (98 % and 90 %) are plotted in the right side of Fig. 5. Due to the significant impact of fluid mixing, the H 2 recovery with the purity above 98 % and 90 % in Case B is 17 % and 18 % lower than those in Case A.

Step #3: H 2 recovery performance including methanogenesis
In the previous section, the significance of the flow dynamics on the fluid distributions (H 2 , CO 2 and water) has been recognised during the processes of H 2 injection and production.Depending on the operation rates, gravity segregation and viscous instability are two driving factors for H 2 recovery performance.In this section, the work is extended by incorporating the methanation process to investigate the interactions between flow dynamics and the process of methanogenesis.To do so, six of the correlated reaction rates outlined in Table 4 were tested in a descending order using 2D flow simulations.This is to detect the end members of the reaction rate at which the methanogenesis process becomes significant or negligible.
Here the cases with the Reaction Rate 1 (denoted as met1 in 2D simulations) is taken as an example to demonstrate the variations in flow behaviour caused by the methanation reaction.The methanation rate selected here is approximately twenty times higher than the reported rate (lower boundary) from Olla field [30].This intentional selection represents an extreme case, utilising an exaggerated reaction rate, to effectively illustrate the interplay between the bio-methanation process and the flow dynamics.Fig. 6 shows that the size of the mixing zone (as indicated by the black brackets) of both cases (Case A_met1 and Case B_met1) are both much smaller than their corresponding cases without the methanation process activated (Case A and Case B).This is because the methanation reaction consumes the total of 5 units of gases (H 2 and CO 2 ) but produces 1 unit of CH 4 and 2 units of water (Equation ( 1)).For this reason, the total volume of the mobile gas decreases with the methanation process.As seen in Fig. 7 (left), the water breakthrough, which does not happen in Case A, occurs at the late stage of withdrawal in CaseA_met1, where the methanation reaction is included.Similarly  for Case B, the methanation process further aggravates the issue of water breakthrough in Case B_met1 and the increase in the water-to-gas ratio in Case B_met1 becomes even more severe (see left side of Fig. 8).This is because the newly generated water from methanation process increases the local saturation of water and its relative water permeability and thus the water mobility.In addition, the H 2 purity also decreases much more dramatically in both Case A_met1 and Case B_met1 than in the cases without methanation included (right side of Figs. 7 and 8), since the H 2 is continuously consumed by the methanation process.Note that the relative permeability between CO 2 and water is used here.This is a valid assupmtion because CO 2 is injected in advance as cushion gas and therefore the two-phase flow is predominantly CO 2 and water, not H 2 and water.However, recent work measuring relative permeability between H 2 and water shows large differences in flow functions compared to CO 2 /water [47][48][49].This concern becomes important when methanogenesis is significant as water is generated during methanogenesis.Therefore, the "correct" choice of relative permeability during back production should be addressed with further analysis when methanogenesis plays a major role.
Another finding is on the rate of H 2 conversion in the entire system throughout the operation period, which is plotted in Fig. 9.In the 1D cases, the methanation reaction is simply a time-dependent process and the amount of CH 4 is incremental with time before H 2 is completely consumed.In other words, the longer time length of the process, the greater the quantity of H 2 and CO 2 that can be converted to CH 4 in the single-cell system, given sufficient reactants.However, as seen in Fig. 9, the total percentage of H 2 conversion in Case B_met1 is 4 % higher than that in Case A_met1 at the end, although the time length of operation in Case B_met1 is much shorter.In particular, the H 2 conversion during the shut-in period is 14 % higher in Case B_met1 than that in Case A_met1.This clearly occurs as a result of the flow dynamics, as both cases have the same time length of shut-in period (100 days).The higher the flow  rate, the more pronounced the viscous instabilities become.As explained previously, the viscous instability has led to the evident H 2 infiltration into the cushion gas of CO 2 in Case B_met1.This in turn amplifies the size of the mixing zone and thus increases the degree of fluid mixing between H 2 and CO 2. Therefore, the extent of the methanogenesis is also much amplified in Case B_met1 compared to Case A_met1 during the shut-in period.
As mentioned above, one of the main objectives of this work is to identify the boundary of the reaction rates at which the methanogenesis process becomes important.Therefore, the H 2 recovery in cases with various reaction rates are plotted in Fig. 10.Based on the reported rate from Olla field (met4), H 2 recovery almost overlaps with the case without including methanogenesis (Case A).The impact of methanogenesis is negligible (<1 %).On the other hand, the methanation process at an average rate of 1746 nm h − 1 (met3) starts to become observable in both Case A and Case B. This rate is twice the reported rate from Olla field [30].This suggests the methanation process at any reaction rate above that may have a noticeable impact on the flow behaviour and the H 2 recovery performance in the system configured here.The extent of the methanation impact on the H 2 recovery becomes clear now and the lower boundary (met3) for the necessity of considering methanogenesis has been identified for this particular system.However, this does not tell us what an acceptable rate of methanation is for a larger or smaller system, and this is dealt with in the next section using scaling theory.

Step #4: validation of scaling theory
In Steps #1 to #3, the impacts of methanation process on the fluid behaviour and thus the H 2 recovery has been identified.To generalise the findings of this work for systems of different sizes and dimensions, scaling theory for subsurface flow in porous media is also tested here [50,51].In the previous study [34], the applicability of scaling theory for flow behaviour during hydrogen storage has been successfully validated, but without considering the methanation process.In this study, the work is extended to conduct some preliminary tests on rescaling the flow simulations with the methanation process considered.
According to Shook et al. [52], the flow behaviour in two different-size systems should be equal if the dimensionless groups (viscous-gravity number, gravity-capillary number etc.) are matched.Taking simulations of this work as an example, if the model of Case A_met1 is expanded by 3 times (or 2 times) in both horizontal and vertical directions, a range of other parameters are adjusted accordingly to achieve the same values of the similarity group and thus the flow behaviour.The full list of adjusted parameters is outlined in Table 5.A full explanation on the application of the system rescaling method is presented in the previous work [34].
Note that the method of modelling the methanation process in this work is based on the Arrhenius equation (Equation ( 2)), which is an empirical formula and based on a rate constant.The reaction will not decrease significantly when the reactants (H 2 or CO 2 ) are still abundant, under certain reservoir conditions (temperature and pressure).Since the total length of Case A_met1 is only 330 days, there is abundant H 2 and CO 2 for the reactions during the operation.This means the modelled reaction performance in a local cell is still at early stage (without running out of H 2 or CO 2 ) and the number of reactions should increase linearly with time.Therefore, with the same rationale of rescaling the capillary pressure, the methanation reaction rate is decreased by 3 times for the expanded case (3 × 3 times larger).This is because the total time length in the expanded models is prolonged by 3 times.Another case with 2 times larger horizontally and vertically is also conducted as a  comparison.Fig. 11 shows the gas saturation and H 2 gas mole fraction between the original case (Case A_met1) and the two enlarged ones.Fig. 12 shows the dynamic parameters including H 2 recovery, H 2 purity, water to gas ratio, and the gas mole fraction of CH 4 throughout the back production.There is negligible difference in the fluid distribution and the operational performance.The rescaling theory works very well in the gravity-dominated Case A_met1, indeed it is almost exact for this case.
In effect, this says: (i) suppose the system is rescaled to a larger one from the "base case" but which is exactly scaled in terms of the viscous/ gravity force balance, (ii) then the maximum acceptable methanation rate constant, r max , is determined for this base case (i.e.gives a negligible effect on H 2 production purity, H 2 loss, etc.); (iii) in an n x n size scaled system, then the corresponding acceptable maximum methanation rate is, r max /n.The rule is the same if a smaller scaled system than the base case is required (where n < 1).
System rescaling tests were also conducted for Case B_met1, where the viscous force and the gravity force are both important.The results in Figs. 13 and 14 show that the overall fluid features and operational performance in the original-size model are reasonably well produced in the rescaled models, but the agreement is not exact.The differences are small and are not obvious by eye in the fluid distributions at each stage  in Fig. 13.However, the H 2 purity and the CH 4 profiles are a little different as shown in Fig. 14.In addition, the water breakthrough occurring in Case B_met1 (see Fig. 14) does not happen in the expanded models.This discrepancy occurs because of differences in the absolute pressures in the base case and larger rescaled models since absolute pressure affects the gas PVT behaviour.The pressure drops during the back production are greater in the larger models than that in the original model.The larger volume of gas expansion has reduced the total quantity of produced gas and therefore postponed the water breakthrough.This issue becomes increasingly evident with the system size and should be treated carefully when applying the scaling theory in the context of geological H 2 storage in porous media [34].However, although scaling theory is not exact when absolute pressures are different (as in the more viscous dominated systems here), the rule that the corresponding acceptable maximum methanation rate is, r max /n (as stated above) is still sufficiently correct for all practical purposes.

Conclusions
In the previous work [34], the range of possible flow behaviour which can be observed during H 2 storage in subsurface porous media has been demonstrated.That work is focused on the balance of forces (viscous/gravity and viscous/capillary), stable and unstable (fingering) flow, the effects of heterogeneity etc.In this paper, the work is extended to incorporate another important mechanism (generally deleterious)the methanogenesis process, which may occur in subsurface H 2 storage.With the use of very-fine scale flow simulations coupled with the bio-methanation reaction, the impacts of the methanogenesis on the fluid behaviour and the operational performance are described and explained in some detail.Scaling theory is again applied for flow simulations coupled with the methanogenesis in the system configured here.The main conclusions from this study are as follows.
1.The methanogenesis process may degrade the operational performance of H 2 storage if the system rate constant for methanation is greater that some notional r max .When the rate constant is above this limit, then the H 2 consumption and CH 4 production are such that the H 2 purity during the back production and the total H 2 recovery at the end of simulation are much decreased.Also, there is a reduction in the quantity of the total gas remaining in the system and the product of water during the methanation process may also accelerate water breakthrough, even in a gravity-dominated scenario.Of course, the corollary is clearly that, if the methanation rate constant can be shown to be significantly « r max , then this issue can be neglected in the H 2 storage assessment simulations.2. Following on from the first conclusion, it is recommended that, in a given field assessment, then an acceptable target r max is established by numerical simulation (using the approach in this paper) and then field and laboratory efforts are made to establish the actual value of this number.In the base case system configured here, the lower boundary of the average methanation rate is 1746 nm h − 1 , above which the impacts of methanogenesis on a time scale of months (330 days) should not be ignored in calculations of the H 2 recovery.This rate is approximately twice as high as the rate reported from the Olla filed.3. Simulations of H 2 storage with CO 2 cushion gas were carried out both when the force balance is known to be gravity dominated (i.e.gravity stable) and when both viscous and gravity forces are significant.The gravity dominated case is generally preferable in terms of hydrogen production and purity, even when methanation is not included in the calculations.When methanation is included, then this general situation continues to be the case.When the viscous forces become significant (higher rate H 2 placement), then the viscous instability which results may lead to H 2 infiltration into the underlying cushion gas of CO 2 during the injection period.This in turn amplifies the size of the mixing zone between CO 2 and H 2 during the shut-in period.As a result, the methanogenesis is amplified, which further degrades the H 2 purity during back production and accelerates water breakthrough.4. Scaling theory has again been applied in this work to rescale the calculations for the "base case" simulation to larger (or smaller) models which are "similar", i.e. they are identical in terms of the balance of forces (viscous/gravity ratio etc.).This has been done for both gravity dominated cases and where both viscous and gravity forces are significant, including methanation reactions in each of these cases.It is shown that the gravity dominated case with the methanation reaction can be scaled almost exactly.Furthermore, if the maximum acceptable methanation rate constant is r max in the

Fig. 1 .
Fig. 1.Correlated random permeability field (correlation range = 1 m in both horizontal and vertical directions) with a high level of heterogeneity: V dp = 0.6 [34].

Fig. 3 .
Fig. 3. Gas saturations (top) and mole fractions of H 2 (bottom) at the end of each cycle of Case A (left) and Case B (operational rates × 10, right).The interface between gas and water (G/W) is indicated by the white line.The sizes of the mixing zone between H 2 and CO 2 are indicated by the black brackets.

Fig. 4 .Fig. 5 .
Fig. 4. Water to gas ratio of under reservoir conditions (left) and H 2 purity (right) throughout the back production of Case A (blue) and Case B (orange), respectively.(For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Fig. 6 .
Fig. 6.Gas saturations (top) and mole fractions of H 2 (bottom) at the end of each cycle of Case A_met1 (left) and Case B_met1 (right).

Fig. 7 .
Fig. 7. Water to gas ratio of under reservoir conditions (left) and H 2 purity (right) throughout the back production of Case A_met1 (blue) compared with Case A (black), respectively.(For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Fig. 8 .Fig. 9 .
Fig. 8. Water to gas ratio of under reservoir conditions (left) and H 2 purity (right) throughout the back production of Case B_met1 (orange) compared with Case B (black), respectively.(For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Fig. 10 .
Fig. 10.Comparison of H 2 recovery versus pore volume productions of cases with various methanation rates.

Fig. 12 .
Fig. 12. Dynamic parameters including H 2 recovery (top left), H 2 purity (top right), water to gas ratio (bottom left) and gas mole fraction of CH 4 (bottom right) against pore volume production throughout the back production in Case A_met1 and the rescaled cases.

Fig. 14 .
Fig. 14.Dynamic parameters including H 2 recovery (top left), H 2 purity (top right), water to gas ratio (bottom left) and gas mole fraction of CH 4 (bottom right) against pore volume production throughout the back production in Case B_met1 and the rescaled ones.

Table 1
Modelling methods and data sources used in this work.

Table 2
Operational strategy of Case A (top) and Case B with 10 times of the operation rates (bottom).

Table 3
Structured outline of numerical tests and objectives of each stage for this study.

Table 4
Correlation of reaction rates based on single-cell simulations.

Table 5
Parameters adjusted to achieve same values of similarity groups.