Dynamics of a compressible gas injected into a confined porous layer

TL;DR

This work develops a long-wave, two-phase model for compressible gas injection into a confined brine-filled porous layer, deriving coupled evolution equations for gas pressure and interface height in a long domain. By nondimensionalizing and performing asymptotic analyses, it reveals distinct regimes where compressibility either transiently or persistently modulates spreading, including an inner incompressible-like region and an outer compressible thin film. The model predicts how gas compression can raise bubble pressure, slow spreading, and alter breakthrough times, with explicit scalings for pressure rise, breakthrough time, and gas mass delivery. Numerical results validate the reduced models and demonstrate the interplay between compression, buoyancy, and viscous dissipation across parameter space relevant to hydrogen storage. Practically, the findings inform injection strategies and reservoir design to optimize storage while mitigating pressure buildup and unwanted spreading in subsurface systems.

Abstract

Underground gas storage is a critical technology in global efforts to mitigate climate change. In particular, hydrogen storage offers a promising solution for integrating renewable energy into the power grid. When injected into the subsurface, hydrogen's low viscosity compared to the resident brine causes a bubble of hydrogen trapped beneath caprock to spread rapidly into an aquifer through release of a thin gas layer above the brine, complicating recovery. In long aquifers, the large viscous pressure drop between source and outlet induces significant pressure variations, potentially leading to substantial density changes in the injected gas. To examine the role of gas compressibility in the spreading dynamics, we use long-wave theory to derive coupled nonlinear evolution equations for the gas pressure and gas/liquid interface height, focusing on the limit of long domains, weak gas compressibility and low gas/liquid viscosity ratio. Simulations are supplemented with a comprehensive asymptotic analysis of parameter regimes. Unlike the near-incompressible limit, in which gas spreading rates are dictated by the source strength and viscosity ratio, and compressive effects are transient, we show how compression of the main gas bubble can generate dynamic pressure changes that are coupled to those in the thin gas layer that spreads over the liquid, with compressive effects having a sustained influence along the layer. This coupling allows compressibility to reduce spreading rates and gas pressures. We characterise this behaviour via a set of low-order models that reveal dominant scalings, highlighting the role of compressibility in mediating the evolution of the gas layer.

Figures (12)

• Figure 1: Schematic showing the displacement of an ambient liquid in a porous medium (regions II$_w$ and III) due to the injection of a compressible gas (occupying regions I and II$_g$) from a line source at $\boldsymbol{x}^* = \boldsymbol{x}_0^*$. The interface separating the fluids is sharp and located at ${y}^* = {F}^*({x}^*, {t}^*)$. The pressure in the ambient brine is hydrostatic at the outlet at $x^*=L$.
• Figure 2: (a) Isothermal equations of state for hydrogen, methane, nitrogen, and carbon dioxide at 333K, sourced from NIST. (b) The dimensionless function $\mathcal{P}$\ref{['eqState1']} for hydrogen, computed using the reference pressure $p_{g0} = 23$ MPa and reference density $\rho_{g0} = 15$$kg.m^{-3}$. The dashed line represents the tangent at the reference state, with slope $\mathcal{P}^{\prime}(1)$.
• Figure 3: A schematic map of $(\zeta, \mathcal{M})$-parameter space. The full problem (\ref{['eq:put2']}) is derived for $\mathcal{M}\sim \mathcal{L} \sim \zeta\sim 1$. The map specialises to the case $\mathcal{L}\gg 1$, focusing on high gas mobility and weak compressibility ($\mathcal{M}\ll 1$, $\zeta \ll 1$); features of the model in the shaded regions in the map become increasingly distinct as $\mathcal{L}\rightarrow \infty$. The distinguished limit $\mathcal{M}\sim \mathcal{L} \zeta \sim 1$ (problem $\Pi_a$) is given by (\ref{['eq:put2t']}) below. An inner/outer structure given by (\ref{['eq:put2ti']}, \ref{['eq:put2to']}) emerges from this system along $\mathcal{L}^{-1}\ll\mathcal{L} \zeta\sim \mathcal{M} \ll 1$ (problem $\Pi_{12}$). The outer problem simplifies to (\ref{['eq:put2toa']}) for $\max(\mathcal{L} \zeta,\zeta^{1/2})\ll \mathcal{M} \ll 1$ (shaded pink, region $\Pi_1$) and (\ref{['eq:put2toc']}) for $\mathcal{L}^{-1}\ll \mathcal{M}\ll \min(1,\mathcal{L} \zeta)$ (shaded blue, region $\Pi_2$). Buoyancy effects influence the inner region for $\mathcal{L} \mathcal{M}^2\lesssim 1$, allowing $X_l$ to recede. Problem $\Pi_b$ arises at $\mathcal{L}\mathcal{M}^2\sim 1$ and $\mathcal{L}\zeta\sim\mathcal{M}$. Ultra-low viscosity effects emerge via (\ref{['eq:put21']}) (problem $\Pi_c$) at $\mathcal{M}\sim \mathcal{L}^{-1}$ and $\zeta\sim \mathcal{L}^{-2}$, from which emerge sub-limits shown in green (region $\Pi_3$; $\mathcal{LM}\ll 1$, $\mathcal{L}^{-2}\ll\zeta\ll 1$) and orange (region $\Pi_4$; $\mathcal{M}\ll \zeta^{1/2}$, $\mathcal{L}^{2}\zeta\ll 1$). Scales for the pressure at the origin $P$, breakthrough time $t_b$ and pressure rise time $t_r$ are shown in blue, green and magenta respectively in the coloured regions; $t_r\ll t_b$ above the line $\mathcal{L} \mathcal{M} \sim 1$ for $\mathcal{L}\zeta \gtrsim 1$. Relative locations within parameter space of results shown in figures \ref{['fig:M_0.1_zeta_10-4']}--\ref{['fig:M_10-3_zeta_10-4']} below are indicated with symbols.
• Figure 4: A simulation of (\ref{['eq:put2']}) for $\zeta = 10^{-4}$, $\mathcal{M} = 0.1$, and $\mathcal{L} = 100$, within the pink region of figure \ref{['fig:overview']}. (a) and (b) show the pressure field $P(x,t)$ and the interface height $F(x,t)$, respectively, at $t = \{0, 0.1, 1, 4, 8, 10.4\}$, following the light-to-dark color gradient. Black dots in (a) show samples of the pressure at the upper contact line $P(X_u,t)$. (c) shows the pressure evolution at the origin and at the upper contact line; the early-time approximation given by \ref{['eq:Boyles']} (dashed orange) captures the transient rise in pressure due to gas compressibility. (d) shows contact-line locations; $X_u$ approaches the late-time solution (\ref{['eq:Peg']}) (dashed black line).
• Figure 5: A simulation of (\ref{['eq:put2']}) for $\zeta = 10^{-3}$, $\mathcal{M} = 0.1$, and $\mathcal{L} = 100$. (a) and (b) show $P(x,t)$ and $F(x,t)$ at $t = \{0, 0.2, 0.8, 4.4, 12, 15\}$. Solutions in (a-d) are plotted using the format described in figure \ref{['fig:M_0.1_zeta_10-4']}.