Analytical and Numerical Study of the Impact of Methanogenic Bacteria on Gas Composition in Underground Hydrogen Storages

Unlike natural gas, hydrogen gas mixtures stored in underground reservoirs undergo active chemical transformations under the influence of methanogenic microorganisms inhabiting in porous reservoirs. They lead to reduction of hydrogen and carbon dioxide concentrations and increase methane concentration. This chemical activity coupled with bacterial dynamics and gas/water flow through porous medium causes the phenomenon of self-organization such as the occurrence of autowave spatial structures, whose dynamics is characterized by a multiplicity of scenarios and bifurcations between them. In this paper we continue to develop the qualitative theory of self-organization in underground hydrogen storage, for more complicated cases that include the mechanism of chemotaxis, which is one of the main types of bacterial movement, and takes into account the flow of both phases. The analysis of scenarios is based on the model of two-phase compositional flow coupled with population dynamics.


Introduction
Renewable energy sources such as hydrogen (H 2 ), solar energy, geothermal energy are considered as the energy of the future [1][2], due to low greenhouse-gas emissions during their production, absence of toxic byproducts when they are utilized and less negative affection to the environment. Among these energy sources, hydrogen energy is considered as the cleanest one with promising application at large scale in the future. The global energy cycle related to H 2 includes [11]: the renewable energy production from windmills and solar cells, the conversion of the exceedingly produced electricity to hydrogen through electrolysis, the hydrogen storage to regulate the difference between intermittent production and permanent gas consumption, and application in fuel cells to produce electricity and for vehicles as a fuel.
In case of producing large amounts of hydrogenous gas, the most efficient and most inexpensive method of storing large amounts of hydrogen is to inject it in geological formations like aquifers, depleted gas reservoirs, or salt caverns [3-4, 11, 15-19].
Several underground storages of hydrogen or town gas exist in the world, for instance, Teeside in the UK, Kiel in Germany, Lobodice in Czech, Beynesan ex-storage in France and storages in Texas state and Russia.
As described in [7,8,10,11] the behaviour of hydrogen in natural rocks is very different from that of the natural gas, as H 2 is chemically active in presence of anaerobic bacteria, which initiate the following global chemical reaction between injected H 2 and CO 2 : The in-situ and laboratory observations really have revealed the increase of CH 4 in the stored gas contents and the decrease of CO 2 and H 2 . Along with this other effects have been observed as the creation of a spatial alternation of the areas saturated preferably by hydrogen or methane. The explanation to these effects have been done in [10] where it was shown that the coupling between hydrogen and CO 2 transport in the reservoir and bio-chemical reactions leads to the appearance of auto-waves equivalent to the observed spatial alternations.
In contrast to papers [9][10] where the fluid was considered as a single-phase gas with residual immobile water, and only the diffusion was retained as the mechanism of bacterial motion, we introduce the following new elements in the present paper: the flow is two-phase; water is mobile; -the appearance of additional form of bacterial colony: the neuston, which is a thin film living at the interface between water and gas; -the chemotaxis as the second important form of bacterial motion, which is the main mechanism of neuston formation. The chemotaxis is the mechanisms of bacterial motion to the direction of nutrients.
The first version of the developed mathematical model was published in [12], in which the impact of the chemotaxis was not considered.

Model of population dynamics
Let us assume that a mixture of H 2 and CO 2 with large domination of hydrogen is injected in an aquifer, which contains water, gas and an initial population of methanogen bacteria. We introduce the chemotaxis, which is the mechanisms of bacterial motion to the direction of nutrients (to the neuston).
The two-phase system in porous medium represents a fine dispersed alternation of gas bubbles or channels with water channels of droplets. At the macrsocale such a system is considered as two interpenetrating continua coexisting at each space point. The water-gas interfaces which are observed on the pore scale disappear in macroscopic description. At any point two phases are identified by saturation of water S.
Both phases can consist of several chemical components: The gas phase essentially consists of H 2 and CO 2 , while liquid consists mainly of H 2 O with low concentration of CO 2 , H 2 , and CH 4 (the injected gas contains low concentration of CO 2 , and hydrogen is low soluble in water). This determines the specific situation when bacteria live in water but the major part of nutriments is concentrated in gas phase.
We consider two kinds of bacteria: -bacteria present in water: they can be plankton or biofilms attached to pore walls wetted by water; -the neuston: a biofilm situated just at the interface between water and gas.
Bacteria living in water consume dissolved H 2 and CO 2 . Bacteria from neuston consume H 2 and CO 2 directly from the gas phase. On the macroscopic scale (Darcy's scale) both phases contain both kinds of bacteria which can be found at any spatial points. Despite the fact that CO 2 is highly soluble in water, it is low present in the injected gas, while hydrogen is very low soluble in water. Therefore, we should assume that the concentrations of both these components in water are of the same order.
In gas we have an abundant resource of H 2 and a sufficiently low resource of CO 2 . Then the eating rate of bacteria in neuston is controlled only by the concentration of CO 2 .
Bacterial population can grow due to replication of species and can decay due to natural or forced death. As usually, we assume that the population grow rate is proportional to the eating rate. For the kinetics of bacterial growth we accept the nonlinear form suggested in [10], which represents a combination of the Turing's and Monod's kinetics. The kinetics of decay is assumed to be constant.
Bacteria also can move. We distinguish several types of their motion: -the plankton can move chaotically similar to Brownian motion (bacterial diffusion); -they can move due to chemotaxis; -bacteria living in water can be transported by water flow (single-phase bacterial advection); -bacteria living in neuston can be transported simultaneously with the movement of the water-gas interface (two-phase bacterial advection) We assume that bacteria in neuston are not transported by chemotaxis but can diffuse. We keep diffusion as it is the mechanism which stabilizes the mathematical properties of the solution, which is considered in the paper [12].
The disappearance of gas-water interfaces in macroscopic equations imposes some difficulties in describing the neuston which represents a pore-scale object. This means that the movement of neuston in macroscopic equations can be obtained by homogenization of its pore-scale motion.

Asymptotic model for low gas saturation
Let's consider the asymptotic model for low gas saturation as S → 1. In this case the neuston is neglected, and bacteria living in water dominate far from the interface. Consequently, the chemotaxis, which determines the neuston formation, should be taken into account. Since the reaction kinetics depends on the concentrations of both CO 2 and H 2 , the model of the process resulting from [12] and consists of three equations in this case (for k = 1.2): Assuming that water density, and Henry coefficients are constant, and neglecting variations of saturation S which is assumed to be close to 1, we obtain the following system of three equations: is the mole fraction of chemical component k in phase i; N is number of bacteria ; S is volume water saturation; i = w, g is the phase (wwater, ggas); D b is the coefficient of bacterium diffusion in bulk water; is the diffusion coefficient of chemical component k in water phase; t e,w is characteristic time of eating at vanishing resource; t d is the time of decay; ρ i is the molar density of phase i (mpl/m 3 ); ϕ is the porosity; G inj is the molar rate of gas injection (mol/s); Ω is the total volume of the reservoir; c (k), inj is the injection concentration of component k in the injected gas (constant value); D ch is the maximal chemotaxis rate; λ ch is empirical parameter; H (k) is the Henry coefficient of component k, which is a given function of pressure; η w is the rendering coefficient (the coefficient of proportionality between the eating rate and growth rate).
Moreover, when the concentration of one of the components is very low, we obtain the model which may be analyzed without simplifications. Let us assume that water contains very low concentration of hydrogen, that is, . Then concentration may be considered as variable with small change. From (2) the following expression is obtained: which is the Turing model [13], if the chemotaxis term is neglected.
For the case without chemotaxis, three cases were analysed in [10]: (i) The case of spatially invariant solutions (no special gradient, only the time derivatives are retained) can be reduced to a dynamic system of the second order: the existence of the limit cycle was revealed for some intervals of parameter variation (α 1 = α 3 = β = 1, 0.90032 < q 1 < 1.0), which means the system oscillates in time and the oscillations tend to be periodic at the infinite time; (ii) The case of time independent, stationary, solutions (no time derivatives, only the special gradients are retained) also can be reduced to another dynamic system of the second order: the appearance of the focus was shown, which means the periodic solutions in space.
(iii) The non-stationary and non space-invariant case (all the derivatives have been retained) has been analyzed numerically: it was shown that an instability appeared in the system, which evolved next to a stable solution stabilizing in time and having the form of periodic spatial waves.

Reduced model of the process, limit stationary spatial waves
To analyse the impact of the chemotaxis, we have to keep the spatial gradients. In contrast we can try to neglect the time variation and analyze only the stationary solutions, which represent the limit of the where ∆ -Laplace operator non-stationary solutions at t → ∞. In the 1D case our model problem becomes: ; . We add the natural boundary conditions: x N which means that the right-hand boundary is impermeable, while a bacterial concentration and H 2 concentration are specified (injected) at the left-hand boundary.
Similar to [10] we will neglect the bacterial diffusion: D b = 0.
The case of this system without chemotaxis was analyzed in [10]. To analyse the impact of the chemotaxis we apply the method of perturbation, by assuming that the chemotaxis parameter, D ch is small. Then the solution of (4)-(5) may be searched in the following form:

Analysis of the first approximation: impact of the chemotaxis
From the second equation of system (8) we can calculate N 1 if the function C 1 is known: This equation is no longer autonomous (due to the explicit presence of x), but is linear with non-constant periodic coefficients. Its analytical solution in terms of the known special functions does not exist [14]. But we can analyze this equation qualitatively.
The periodicity of the coefficient χ(x) and the right-hand part f(x) does not mean the periodicity of the solution of Eq. (10). In contrast, for the homogeneous equation: , the periodic solution can exist only if coefficient χ(x) <0, which imposes the inequality on the parameters: β -2α 3 C 0 (x) N 0 (x)<0. Otherwise, i.e., when β -2α 3 C 0 (x)N 0 (x)>0, the solutions of the homogeneous equation (10) are exponential, and contain both the exponentially decreasing and increasing components. In all these cases, the right-hand part, which is responsible for the chemotaxis, plays the role of an external force that progressively attenuates or amplifies the behavior of the homogeneous solution. Therefore, the periodic solutions of non-homogeneous equation (10) can not exist.
Consequently, we expect to have two scenarios for (10): (i) the solution is non-oscillating with almost exponential growth in space; (ii) the solution is oscillating but non-periodic in space. The total solution (6) is expected to be oscillating (due to the zero-order term), but non-periodic in space, due to the first-order term.
We see that the chemotaxis perturbes the periodic solutions and involves some disordered behaviour, such that the concentrations becomes non-periodic, but still oscillating functions.
This is explained by the fact that the chemotaxis means the oriented motion of bacteria to the direction of nutriments. If a zone with an exceeding amount of nutriments arises, then bacteria tend to move towards this zone. This oriented motion perturbs the system periodicity.
The results of numerical simulations presented in following section confirm this qualitative analysis.

Numerical Study of the Impact of the Chemotaxis
We also analyzed the total non-stationary behaviour of the system (3) in general case. This analysis can be performed only numerically.

Problem formulation and parameters
We analyze system (3) of gas injection in two-dimensional case with constant initial conditions and Neumann boundary conditions which correspond to impermeable boundaries: where n is the normal direction to the boundary. Table 1 shows the data used in the calculations.
The initial values are located within the zone of attraction of the limit cycle, so that the solution of this problem is space-invariant and oscillating in time, if the chemotaxis would be neglected. The flow rate q 1 in (3) represents the hydrogen injection into the reservoir. The numerical simulations show that the solution invariable in space and varies only in time, as expected. We have obtained the oscillations in time.
After this, since a fixed moment, this space-invariant solution was perturbed in the form of an instantaneous non-zero concentration gradient applied to the small vicinity of the origin. The evolution of the perturbation is shown in Figs. 1 and 2.

The evolution of the solution in time
After perturbation, the irregular waves travelling throughout the overall domain were observed. Their evolution was very fast establishing to the structure presented by regular periodic waves invariable in time. The Figs. 1 and 2 represents the results of numerical calculation of the evolution of the hydrogen concentration, the variation in the number of bacteria at t = 50 ÷ 500.
As seen, sufficiently regular ring waves are developed with excess and deficiency of hydrogen and bacteria in some areas, which alternate with each other. In areas with high bacterial concentrations, where the reaction CO 2 + 4H 2 = CH 4 + 2H 2 O is rapid, the alternation with the ring excess and deficiency of bacteria appear, whereby the methanogenic bacteria generates methane. In case of taking into account the chemotaxis of bacteria, the bacteria forms neuston.
In this work an attempt has been carried out to qualitatively analyze the impact of methanogenic bacteria on the dynamics of the generation of methane in underground hydrogen storage. The occurrence of undamped oscillations during the time, which tends asymptotically to periodic waves, means that the system undergoes self-organization of new structures in the form of methane accumulations. It should be noted that in several cases the damping oscillations were observed in space caused by the chemotaxis. In the limit of computational time, the steady-state spatial pattern of frozen waves is observed. Following the results presented in Fig.  2, the effect of a natural in situ separation of hydrogen gases was observed, which corresponds to the observations in underground storage Lobodice [7,8].
According to these simulations, the solution is periodic in the case without chemotaxis (the lefthand figures). For the case when the chemotaxis is the dominating mechanism of bacterial motion (the right-hand figures), the solution remains non-periodic for any moment of time.

Conclusions
In papers [7] and [8] it was proved that an underground storage of hydrogen can work as a natural chemical reactor producing methane from hydrogen and carbon dioxide. The reaction between H 2 and CO 2 (1) is catalyzed by methanogenic bacteria and happens in the form of the bacterial metabolism.
In paper [10] the first mathematical model of the process was developed. It was based on the single-phase flow model coupled with population dynamics equation. The bacterial population was considered in the average; various forms of its existence were reflected in special nonlinear kinetic function of population growth.
In paper [12] we suggested the two-phase flow model coupled with the dynamics of two bacterial populations: bacteria living in water and neustona thin biofilm situated at the interfaces between water and gas. We also suggested the mathematical model of chemotaxis in two-phase fluid, which is the main mechanism of neuston formation: bacteria living in water feel the presence of nutriments concentrated in the injected gas and move to the direction of the interfaces water-gas.
In the present paper, we studied the above mentioned mathematical model [12] analytically and numerically. To simplify the analysis, several simplifications have been accepted, as, for instance, the low gas saturation. For this case, we obtained the asymptotic model (3). The main objective of the analysis was to study the impact of the chemotaxis on the storage behaviour.
The analytical study was performed under the assumption that the chemotaxis was weak, so that we could apply the method of perturbation. We analyzed the behaviour of stationary solutions, which correspond to the infinite times. We obtained the non-linear autonomous dynamic system, which depends on space coordinates, (8). In the zero approximation (without chemotaxis), we obtained the system analyzed previously in [10]. The influence of the chemotaxis appears through the first approximation, which leads to a linear ordinary differential equation with periodic coefficient. Its solution combines the properties of an oscillating function and an exponential function, therefore it is not periodic. This means that the chemotaxis perturbs the system periodicity, which corresponds to the physical meaning of the oriented motion of bacteria caused by the chemotaxis.