Ostwald Ripening in Underground Gas Storage

TL;DR

This study investigates Ostwald ripening in underground gas storage, revealing a two-stage process through microfluidic experiments and a validated continuum model, with implications for long-term storage stability.

Contribution

It provides the first detailed characterization of Ostwald ripening in confined porous media and develops a predictive model linking microscale processes to field-scale storage timescales.

Findings

Ripening occurs faster than convective dissolution in CO2 storage.

A two-stage dynamic of local equilibration followed by slow global depletion was observed.

The model accurately predicts gas saturation evolution across conditions.

Abstract

Underground gas storage is central to both climate mitigation and energy transition strategies, supporting both long-term carbon sequestration and seasonal hydrogen storage. A key mechanism governing the fate of injected gases is Ostwald ripening, the curvature-driven mass transfer between trapped gas ganglia in porous media. While ripening is well understood in open systems, its behavior in geometrically confined porous structures remains poorly characterized, especially over long timescales relevant to subsurface operations. Here, we present ultra-high-resolution microfluidic experiments that capture the evolution of residually trapped hydrogen over weeks in realistic, heterogeneous pore geometries under well-defined boundary conditions. We observe a distinct two-stage dynamic: a rapid local equilibration among neighboring bubbles, followed by slow global depletion driven by long-range diffusion toward low-chemical potential boundaries. Building on these insights, we develop and validate a continuum model that couples microscale capillary pressure-saturation (Pc-s) relationships, extracted via the pore-morphology method, with macroscopic diffusive transport. The model accurately predicts gas saturation evolution without fitting parameters and collapses experimental results across a range of experimental conditions. Extending the model to reservoir scales, we estimate equilibration timescales for CO2 and H2 in homogeneous sandstone aquifers. We find that ripening occurs much faster than convective dissolution in CO2 sequestration, and on timescales comparable to seasonal H2 storage operations. These findings establish a quantitative framework linking pore-scale heterogeneity to field-scale gas redistribution, with implications for the design and longevity of subsurface storage strategies.

Equations

$$\left\{\left[\frac{H}{RT{\bar{\rho }}_w}-1\right]\left[H\left(\frac{\text{d}{S}_g}{\text{d}{P}_c}\right)\chi +S_g\right]+1\right\}\frac{\partial \chi }{\partial t}=\frac{D_w}{\tau }{\nabla }^2\chi ,$$

$$\frac{\partial \chi }{\partial t}=D_w\varphi {\left[H\frac{{\bar{\rho }}_g}{{\bar{\rho }}_w}\frac{\text{d}{S}_g}{\text{d}{P}_c}+1\right]}^{-1}\frac{{\partial }^2\chi }{\partial x^2},$$

$$\frac{\partial S_g}{\partial t}=D_{\text{eff}}\frac{{\partial }^2{S}_g}{\partial x^2},$$

Full text

I Introduction

Underground gas storage is essential to addressing contemporary energy and environmental challenges [1]. Geologic sequestration of carbon dioxide (CO2) provides a long-term solution for reducing atmospheric greenhouse gas concentrations by permanently storing CO2 in deep geological formations [2, 3]. Concurrently, seasonal underground hydrogen (H2) storage offers a scalable means of buffering the intermittency of renewable energy sources, thereby facilitating the reliable integration of low-carbon energy into the broader energy infrastructure [4].

As gases such as CO2 and H2 are injected into the subsurface, they displace the wetting ambient brine in a drainage process. Upon subsequent reinvasion of the brine—referred to as imbibition—portions of the injected gas become immobilized as discrete ganglia surrounded by the wetting phase, in a process known as residual trapping. In geological carbon sequestration, residual trapping primarily occurs as buoyant CO2 migrates upward [5]. In underground hydrogen storage, it arises from the cyclical drainage and imbibition of the pore space during repeated injection and withdrawal cycles [6].

Once trapped, these isolated gas ganglia partially dissolve into the surrounding aqueous phase. Over time, molecular diffusion drives mass exchange between bubbles of differing sizes and curvatures—a process known as Ostwald ripening [7, 8, 9]. This redistribution of mass causes some bubbles to shrink and eventually vanish, while others to grow sufficiently to overcome capillary constraints, regain buoyancy, and migrate upward to rejoin the mobile gas plume. Such dynamics may prove advantageous for underground hydrogen storage by reducing residual hydrogen losses; however, they pose a significant concern in geological carbon sequestration, as bubble mobilization can increase the risk of CO2 leakage. In addition, ripening-induced redistribution of bubbles modifies local flow and chemical conditions, thereby influencing geochemical reactions such as rock dissolution and mineral precipitation [10].

Ostwald ripening, first described in the late 19th century, has been extensively studied in open systems, where gas bubbles grow unconstrained by geometry. Its investigation within porous media, however, is a more recent development. Pore-scale observations have revealed that geometric confinement fundamentally controls the ripening dynamics. In micromodels with homogeneous pore structures, an initially polydisperse population of gas bubbles evolves toward a monodisperse state, with bubbles attaining uniform size [11]. This behavior changes markedly in the presence of pore-scale heterogeneity. Pore-scale modeling studies predict that under such heterogeneity, ripening leads to a stable equilibrium distribution of bubble sizes, governed by the interplay between initial conditions and the underlying pore geometry [12, 13, 14]. Recent X-ray computed tomography (CT) experiments of hydrogen storage in rock cores reveal that ripening behavior is dictated by ganglion size: medium-sized ganglia fragment, larger ganglia grow, and smaller ones dissolve completely [15, 16, 17]. These CT observations, however, are typically limited to small subdomains and short timeframes ($\sim$1 day), providing only a partial—and potentially unrepresentative—snapshot of the ripening process. Concurrently, continuum-scale models have been developed that couple capillary pressure–saturation ($P_{c}$–$S$) relationships with Henry’s law to describe the macroscopic evolution of Ostwald ripening in porous media [18, 19]. However, the appropriate choice of $P_{c}$–$S$ relationships for modeling ripening remains unclear, and the predictive accuracy of these models has yet to be tested against experimental data. A promising approach for extracting $P_{c}$–$S$ relationships from microscale geometry was recently proposed by Mehmani and Xu [19], and is validated in this work.

In this study, we perform microfluidic experiments to investigate Ostwald ripening of hydrogen in realistic pore geometry under well-defined boundary conditions. By achieving ultra-high-resolution imaging (0.3 $\mu$m/pixel) across the entire microfluidic domain, we track the collective evolution of a large bubble population while simultaneously resolving the interface curvature—and thus the capillary pressure—of individual bubbles with high accuracy. Each experiment spans over 15 days, an order of magnitude longer than in previous studies, allowing us to capture long-term dynamics. We observe that neighboring bubbles rapidly reach local equilibrium, as pore-scale heterogeneities provide a multitude of metastable configurations that facilitate equilibration through Ostwald ripening. In contrast, equilibration with the boundaries—which are much larger than a typical pore throat—occurs significantly more slowly, governed by the timescale of long-range diffusive mass transfer.

Crucially, we demonstrate for the first time that a continuum model incorporating microscale pore-geometry characteristics can accurately predict the evolution of gas saturation without the use of fitting parameters. When applied to representative sandstone formations, the model yields characteristic equilibration times for CO2, H2, and CH4, showing that Ostwald ripening proceeds significantly faster than convective dissolution in carbon sequestration scenarios and occurs on timescales comparable to seasonal hydrogen storage cycles. These findings establish the first validated, predictive framework that links pore-scale heterogeneity to field-scale gas redistribution.

Microfluidic Experiments

We conduct microfluidic experiments to investigate Ostwald ripening among residually trapped hydrogen bubbles in porous media. The micromodel is fabricated from a two-dimensional X-ray CT image of a Canadian sandstone, preserving the pore-scale heterogeneities characteristic of natural porous media in the subsurface (Fig. 1). The resulting structure exhibits a median pore size of 42 $\mu$m, a median pore-throat size of 10.6 $\mu$m, and a porosity of 0.3 (SI Appendix, Fig. S1). The micromodel is laterally bounded by inlet and outlet channels with widths of 100 $\mu$m and 170 $\mu$m, respectively. The device is fabricated by Interface Fluidics (Canada) using deep reactive ion etching (DRIE) to a depth of 7.8 $\mu$m on a silicon wafer, followed by anodic bonding to borosilicate glass. This approach provides high mechanical integrity and hermetic sealing, effectively preventing hydrogen leakage during long-term experiments. A heating element is placed beneath the micromodel to maintain uniform and controlled temperature conditions throughout the experiments.

To initiate each experiment, the micromodel is first saturated with deionized water (MilliporeSigma). Research-grade hydrogen gas (Air Liquide) is then introduced at a constant rate for a total of 50 pore volumes. Subsequently, water is reinjected for a total of 100 pore volumes, resulting in the residual trapping of hydrogen [20] with an initial gas saturation $S_{g}^{\text{init}}$ ranging from 0.49 to 0.65. To establish well-defined boundary conditions, hydrogen bubbles are deliberately retained in both the inlet and outlet channels, which are significantly larger than the pore throats in the porous matrix. These larger bubbles exhibit lower chemical potential, thereby driving diffusive mass transfer from the porous matrix toward the boundaries. The system is allowed to equilibrate until high-precision single-crystal quartz resonant pressure sensors (XTALX LUQS1, Phase Sensors) at both the inlet and outlet return to atmospheric pressure. Once pressure equilibrium is reached, the inlet and outlet valves are closed, thereby isolating the micromodel and initiating the Ostwald ripening process. Experiments are conducted at two distinct temperatures, 40∘C and 80∘C, with each experiment lasting 24 days and 15 days, respectively, except in one case where the experiment was terminated prematurely due to an overnight power outage. The micromodel is mounted on the motorized stage of an upright microscope (BX51, Olympus), and time-lapse imaging is performed by scanning the domain in a serpentine pattern. This imaging protocol yields an ultra-high resolution of 0.3 $\mu$m/pixel, with each full-frame image composed of 120 stitched sub-images (Materials and Methods).

Results

Pore-scale Observations. The sizes of trapped hydrogen bubbles span more than two orders of magnitude, ranging from single-pore bubbles to ganglia that extend across multiple pores. Due to spatial confinement, large ganglia possess multiple hydrogen–water interfaces. The ultra-high resolution of our experiments enables accurate measurement of in-plane interface curvature via least-squares circle fitting. We find that all interfaces associated with a given bubble exhibit the same in-plane curvature, although neighboring bubbles display distinct values (Fig. 2A). This uniformity within individual bubbles indicates local capillary equilibrium, consistent with 3D tomography observations of residually trapped CO2 bubbles [21, 22]. A survey of the entire bubble population in the micromodel shows that initial in-plane curvatures vary by an order of magnitude, with no apparent correlation to bubble size (SI Appendix, Fig. S2).

We observe water droplets within some trapped bubbles at the onset of the experiment (Fig. 2B). These droplets arise from condensation triggered by a sudden increase in system pressure during water reinjection, which is carried out to displace hydrogen immediately prior to the experiment. The abrupt pressure rise increases the total pressure inside the hydrogen bubble, but the partial pressure of water vapor cannot adjust rapidly enough via diffusion. As a result, the vapor becomes supersaturated and condenses, preferentially nucleating on surface imperfections within the micromodel. Once formed, these condensed water droplets remain largely unchanged throughout the experiment (SI Appendix, Fig. S3) and have a negligible impact on the ripening process, as demonstrated later.

We examine a 1.5 mm$\times$1.5 mm region near the center of the micromodel in the experiment conducted at 80∘C with an initial gas saturation of $S_{g}^{\text{init}}=0.49$. Within this domain, some bubbles grow, others shrink, and some remain unchanged, reaching apparent equilibrium within approximately 5 hours (Fig. 2C). This rapid re-equilibration is facilitated by pore-scale heterogeneity, which offers a wide range of metastable configurations. As a result, individual ganglia require only minimal volume adjustments to achieve capillary equilibrium with their neighbors. At this stage, boundary effects have not yet reached the central region. Over a longer timescale of 1 to 15 days, however, we observe that all bubbles gradually shrink due to diffusive mass loss driven by Ostwald ripening toward the left and right boundaries (Fig. 2D). A high-resolution time-lapse video of the experiment is provided in SI Appendix, Movie S1.

Two-stage Ripening. Pore-scale observations of individual bubbles reveal a distinct two-stage evolution during the ripening process. To quantify this behavior, we track the residual gas saturation ($S_{g}$), the average in-plane radius of curvature ($\bar{r}$), and the standard deviation of radius of curvature ($\sigma_{r}$) for a population of bubbles within a 1.5 mm-long window that spans the full width of the micromodel. This analysis window is centered within the device, equidistant from the left and right boundaries. Figure 3A presents these measurements for an experiment conducted at 80∘C. We find that $\sigma_{r}$ decreases rapidly at the start of the experiment, reaching a minimum value near zero within five hours, while $S_{g}$ and $\bar{r}$ remain relatively constant over the same period. The near-zero value of $\sigma_{r}$ indicates that bubbles in the central region have reached capillary equilibrium, each attaining the same interface curvature. We refer to this rapid equilibration phase as local ripening (Stage 1).

As the ripening process progresses, bubbles in the central region begin to experience the influence of the boundaries, which exhibit lower chemical potentials due to the larger aperture sizes. This chemical potential gradient drives hydrogen diffusion from the center toward the boundaries, leading to a gradual decrease in $S_{g}$ as bubbles shrink. This depletion phase is accompanied by a monotonic increase in $\bar{r}$. Meanwhile, $\sigma_{r}$ rises to nonzero values, as diffusive mass loss disrupts the capillary equilibrium established in Stage 1. Shrinking bubbles struggle to maintain metastable configurations with their neighbors, resulting in a broader distribution of curvatures. We refer to this slow, boundary-driven depletion as global ripening (Stage 2) (Fig. 3A). Plotting the evolution of $S_{g}$ across the entire micromodel reveals the same two-stage ripening process across all experimental conditions, though the distinction between Stage 1 and Stage 2 becomes less pronounced, as bubbles near the boundaries begin to lose mass as soon as the experiment starts (Fig. 4A).

To elucidate the thermodynamic driving force underlying the ripening process, we examine the evolution of interfacial energy as a function of gas saturation for the same population of bubbles located at the center of the micromodel. We define the reference energy as that of a fully water-saturated micromodel and assign it a value of zero. Owing to the quasi-two-dimensional geometry of the micromodel, the in-plane component of interfacial energy scales linearly with gas saturation and thus offers limited insight into ripening dynamics. We therefore focus our analysis on the out-of-plane component of interfacial energy given by $F_{\text{o-o-p}}=(\gamma_{\textrm{gs}}-\gamma_{\textrm{ls}})bl_{\textrm{gs}}+\gamma b\pi l_{\textrm{gl}}$, where $\gamma_{\textrm{gs}}$, $\gamma_{\textrm{ls}}$, and $\gamma$ are the interfacial tensions of the gas-solid, liquid-solid, and gas-liquid interfaces, respectively; $l_{\textrm{gs}}$ and $l_{\textrm{gl}}$ are the corresponding interfacial lengths; and $b$ is the micromodel gap thickness. Assuming that water perfectly wets the solid, Young’s equation yields $\gamma_{\textrm{gs}}-\gamma_{\textrm{ls}}=\gamma$, simplifying the expression to $F_{\text{o-o-p}}=\gamma bl_{\textrm{gs}}+\gamma b\pi l_{\textrm{gl}}$. This formulation represents the total out-of-plane interfacial energy as the sum of the gas–solid and gas–liquid contributions and serves as a thermodynamic indicator of the system’s progression toward equilibrium during the ripening process.

As expected, $F_{\text{o-o-p}}$ decreases monotonically in time as the system evolves toward a lower-energy configuration (Fig. 3B). The steepest decline occurs during Stage 1, as large bubbles retract from narrow pore throats (Fig. 2C), leading to a rapid decline in gas–solid interfacial energy. This decline, however, is partially offset by a concurrent increase in gas–liquid interfacial energy. These trends persist into Stage 2, albeit at a much slower rate, as diffusive mass loss to the boundaries drives bubbles to further reduce their contact with the solid surface while transitioning toward lower-curvature interfaces (Fig. 2D). The thermodynamic description of multiphase flow in porous media, which couples capillary pressure, saturation, and interfacial area, has been extensively explored [23, 24, 25]. However, its application to Ostwald ripening remains very limited [26] and represents a promising direction for future investigation.

Continuum Model of Ripening. To capture the macroscopic dynamics of Ostwald ripening in porous media with pore-scale heterogeneity, we develop a continuum model based on the pore-network formulation of Bueno et al. [27]. The model describes the spatiotemporal evolution of the molar fraction of dissolved gas in the aqueous phase, denoted $\chi$, and is given by

[TABLE]

where $H$ is Henry’s constant in volatility form (with units of pressure per molar fraction), $R$ is the universal gas constant, $T$ is the absolute temperature, $\bar{\rho}_{w}$ is the molar density of the aqueous phase, $D_{w}$ is the molecular diffusion coefficient of the gas in the aqueous phase, and $\tau$ is the tortuosity of the porous medium. Assuming both $\chi$ and $S_{g}$ vary only along the longitudinal direction, we solve Equation 1 in one dimension (1D). Following the Bruggeman relation [28], we further assume $\tau=\phi^{-1}$, where $\phi$ is the porosity. The term $\textrm{d}S_{g}/\textrm{d}P_{c}$ captures the pore-scale coupling between interfacial curvature and gas redistribution, linking saturation changes to capillary pressure gradients. Rather than imposing an empirical capillary pressure–saturation relationship, we compute $\textrm{d}S_{g}/\textrm{d}P_{c}$ directly using a modified pore-morphology method (PMM) adapted for ripening, with the micromodel pore geometry as direct input (Materials and Methods). This approach allows us to explicitly encode the influence of pore-scale geometry on macroscopic ripening dynamics [19].

We impose Dirichlet boundary conditions on the molar concentration at the left and right inlet channels of the domain, specified by Henry’s law as $\chi=(P_{w}+P_{c}-P_{v})/H$, where $P_{w}$ is the water pressure (equal to atmospheric pressure in our experiments) and $P_{v}$ is the saturated vapor pressure of water. Assuming the liquid phase is perfectly wetting to the micromodel, we calculate the capillary pressure $P_{c}$ at the boundaries using the Young-Laplace equation, $P_{c}=2\gamma({d^{-1}+b^{-1}})$, where $\gamma$ is the gas-water interfacial tension, $d$ is the constant width of the boundary channels and $b$ is the gap thickness of the micromodel. To relate changes in gas saturation to variations in dissolved concentration, we solve for the capillary pressure as $P_{c}=H\left(\chi-\chi_{0}\right)+P_{v}$, where $\chi_{0}$ is the equilibrium molar concentration of dissolved gas under ambient temperature and pressure, assuming zero gas–liquid interface curvature. We then invert the capillary pressure–saturation relationship to obtain $S_{g}$. Remarkably, the model described by Equation 1 captures the experimentally observed evolution of gas saturation without the need for any fitting parameters (Fig. 4A). A full derivation of the model and underlying assumptions is provided in the SI Appendix.

Noting that the dimensionless prefactor $H/(RT\bar{\rho}_{w})\gg 1$, due to the low solubility of gases in water, and that $H\left(\frac{\textrm{d}S_{g}}{\textrm{d}P_{c}}\right)\chi=\left(\frac{\textrm{d}S_{g}}{\textrm{d}P_{c}}\right)P_{g}\gg S_{g}$, Equation 1 in 1D simplifies to a diffusion equation:

[TABLE]

where $\bar{\rho}_{g}$ is the partial molar density of the gas phase. We define the characteristic timescale as the diffusion timescale, $T_{c}=L_{c}^{2}/D_{\text{eff}}$, where the characteristic length $L_{c}$ is set by the length of the micromodel, and the effective diffusivity is given by $D_{\text{eff}}={D_{w}}{\phi}\left[H\frac{\bar{\rho}_{g}}{\bar{\rho}_{w}}\frac{\textrm{d}S_{g}}{\textrm{d}P_{c}}+1\right]^{-1}$. Rescaling time $t$ by $T_{c}$, and residual gas saturation $S_{g}$ by its initial value, perfectly collapses the saturation evolution curves across different experimental conditions (Fig. 4B).

Large-scale implications. Saline aquifers are widely regarded as a promising geological setting for large-scale underground gas storage due to their widespread availability and substantial storage capacity [29, 3, 4]. Among these, macroscopically homogeneous saline aquifers—particularly those composed of consolidated sandstones—are especially well-suited owing to their enhanced predictability, containment security, and operational reliability. Notable examples include the Mount Simon Sandstone in the Illinois Basin and the Basal Cambrian Sandstone in Western Canada. Both formations have been extensively characterized and exhibit a high degree of homogeneity over meter to tens-of-meter scales, with relatively uniform porosity and permeability, despite underlying heterogeneity at the pore scale [30, 31].

Residual trapping in such formations has been investigated through core flooding experiments conducted under reservoir-relevant conditions. These experiments typically involve gas injection into a brine-saturated core (drainage), followed by brine reinjection (imbibition), thereby mimicking the primary processes responsible for residual trapping. When coupled with X-ray computed tomography (CT) imaging, these tests enable direct visualization of the spatial distribution of trapped gas along the core [15, 32]. Analysis of four saturation profiles reveals that the residual gas saturation varies on the centimeter scale and can be reasonably approximated by a sinusoidal function: $S_{g}^{\textrm{init}}=\hat{S}_{g}\sin(2\pi/\lambda~x)+\bar{S}_{g}$, where $x$ is the longitudinal distance, $\lambda\sim{0.11}-0.63$ cm is the wavelength, $\hat{S}_{g}\sim{0.03}-0.1$ is the amplitude of saturation variation, and $\bar{S}_{g}\sim{0.26}-0.46$ is the mean saturation (SI Appendix). Furthermore, we apply our modified PMM code to extract $P_{c}$-$S_{g}$ relationships from publicly available 3D X-ray CT scans of homogeneous sandstone cores (SI Appendix, Fig. S5). These curves are more appropriate for predicting ripening, as they reflect equilibrium configurations of disconnected ganglia and avoid the viscous effects or hysteresis commonly associated with experimental measurements (Materials and Methods). Noting that $\textrm{d}S_{g}/\textrm{d}P_{c}$ is approximately constant over the relevant saturation range ($S_{g}\in{0.2}-0.6$), Equation 2 can be expressed as an evolution equation for $S_{g}$:

[TABLE]

which admits a well-known analytical solution in the absence of a concentration sink (e.g., a fracture). Specifically, the initial variations in ${S}_{g}$ decay exponentially with a characteristic equilibriation time $t_{\text{eq}}={\lambda^{2}}/({4\pi^{2}D_{\textrm{eff}}})$.

We estimate $t_{\text{eq}}$ for three gases commonly considered for subsurface storage: H2, CO2, and natural gas (CH4). To account for depth-dependent effects, we develop constitutive models that capture the nonlinear variation of $D_{w}$, $H$, and $\bar{\rho}_{g}$ with depth, incorporating geothermal and hydrostatic gradients while neglecting the influence of brine salinity. These constitutive models are validated against published experimental and field data (SI Appendix). Across all three gases, $t_{\text{eq}}$ increases with depth, but the change becomes marginal beyond approximately 1 km, where increases in $D_{w}$ are largely offset by rising $H$, and variations in $\bar{\rho}_{g}$ are minor. Across all depths, CO2 exhibits the shortest equilibriation time, followed by H2 and then CH4. The equilibriation time also varies substantially with $\lambda$. For a macroscopically homogeneous saline aquifer at a depth of 1 km, with $\phi=0.2$ and $\lambda\sim 0.11-0.63$ cm as inferred above, the estimated $t_{\text{eq}}$ for CO2 ranges from approximately 20 to 500 days, with a median value of 300 days (Fig. 5). This is substantially shorter than the onset time of convective dissolution in mobile CO2 plumes, typically 10–100 years in similar saline aquifers [33, 34]. In the absence of sinks such as fractures, ripening-induced equilibration is not expected to compromise residual trapping, as only minimal volume adjustments are required to reach local capillary equilibrium given the broad range of metastable bubble configurations at the pore scale (Fig. 2C). Meanwhile, the estimated $t_{\text{eq}}$ for H2 ranges from approximately 50 to 1,500 days, with a median of 700 days—comparable to the 6–12 month injection–withdrawal cycles characteristic of seasonal underground hydrogen storage, yet orders of magnitude shorter than the million-year timescales required to dissolve an entire column of residually trapped gas via ripening [35]. Depth-dependent ripening arises from the increase in hydrostatic pressure with depth, which elevates the equilibrium concentration of dissolved gas, thereby creating a vertical concentration gradient and sustaining an upward diffusive flux. Our analysis does not account for ripening in the vertical direction.

The results above pertain to local ripening (Stage 1). A second implication arises from global ripening (Stage 2), also captured by our experimentally validated continuum model (Eq. 2). Mehmani and Xu [19] solved a linearized form of this model to examine bubble depletion in an infinite, microscopically heterogeneous porous medium, where large bubbles in cavities or fractures act as constant-concentration sinks. Their analysis showed that the mass flux into such sinks may decrease, remain constant, or increase depending on whether ripening proceeds in linear, cylindrical, or radial geometries. The latter two geometries are especially relevant to annular gaps behind wellbore casings and their cement sheath, where ripening can sustain or even accelerate leakage. Because these outcomes are direct mathematical consequences of Equation 2, now corroborated by experiments, they underscore wellbore integrity as a critical vulnerability in underground gas storage—an issue evidenced by recent explosions and uncontrolled releases linked to wellbore corrosion in underground gas storage facilities [36].

Conclusions

Ostwald ripening plays a critical role in governing the long-term fate of residually trapped gas in subsurface gas storage. Through high-resolution microfluidic experiments performed under well-defined boundary conditions, we present the first long-term investigation of Ostwald ripening dynamics within heterogeneous porous media. Our results provide direct experimental evidence that ripening proceeds in two distinct stages: a rapid local equilibration driven by pore-scale heterogeneity, followed by a slower global depletion governed by long-range diffusive mass transfer toward boundary regions with lower chemical potentials (Figs. 2–3).

Leveraging these insights, we develop and validate a continuum-scale model that couples pore-scale capillary pressure-saturation ($P_{c}$-$S$) relationships to diffusive mass transport. Crucially, this study is the first to demonstrate that $P_{c}$-$S$ curves extracted using a simplified pore morphology method adapted for ripening can accurately predict the evolution of gas saturation driven by Ostwald ripening, without the need for fitting parameters (Fig. 4). Applying the validated model to macroscopically homogeneous sandstone aquifers, we estimate characteristic equilibration times ($t_{\text{eq}}$) for CO2, H2, and CH4 under reservoir-relevant conditions. We find that $t_{\text{eq}}$ increases with depth but plateaus beyond approximately 1 km due to compensating thermophysical trends. Among these gases, CO2 equilibrates most rapidly, followed by H2 and CH4 (Fig. 5).

Our findings show that Ostwald ripening occurs significantly faster than the onset of convective dissolution in mobile CO2 plumes during geologic carbon sequestration. In contrast, the characteristic timescale of ripening is comparable to the seasonal injection–withdrawal cycles typical of underground hydrogen storage. Notably, our analysis reveals that the equilibration timescale is sensitive to the spatial wavelength of saturation heterogeneity—a parameter we quantify using existing core-scale experiments. Future studies should focus on how pore geometry, wettability, and flow conditions influence the initial saturation heterogeneity.

Altogether, this study provides a quantitative framework for understanding and predicting Ostwald ripening in subsurface porous media. By directly linking pore-scale geometry to continuum-scale gas transport, our approach establishes a physics-based foundation for assessing the long-term security and efficiency of geological gas storage systems.

II Materials and Methods

Imaging

We visualize the ripening process using an upright brightfield microscope (BX51, Olympus) equipped with a 10$\times$ objective lens and a scientific CMOS camera (BUC3D-1000C, BestScope) featuring a 10-megapixel sensor. To capture the entire micromodel at ultra-high resolution, we employ a motorized XY stage (EK 14 mot, Marzhauser Wetzlar), which is programmed to traverse the field of view in a serpentine pattern with 35% overlap between adjacent frames. Each scan acquires 15 columns and 8 rows of images, yielding a total of 120 individual frames, with each frame spanning 1.14 mm $\times$ 0.85 mm. These images are stitched using the Grid/Collection Stitching plugin in Fiji [37], producing a composite image with a spatial resolution of 0.3 $\mu$m/pixel that enables detailed observation of pore-scale ripening dynamics across the entire micromodel.

To characterize the trapped hydrogen bubbles, we developed a custom image analysis pipeline in MATLAB. First, a dry image of the micromodel is processed using intensity thresholding and morphological operations to extract the pore structure. This step also generates a binary mask that identifies the solid matrix, which is subsequently applied to all time-lapse frames to exclude static grain regions. Gas bubbles are then identified via adaptive thresholding and mask subtraction. Each identified bubble is uniquely labeled and analyzed to determine its centroid position, projected area, pore occupancy number, gas–water and gas–solid interfacial lengths, and Euler characteristic. This enables robust tracking and morphological analysis of individual bubbles throughout the ripening process.

Pore morphology model

We parameterize the continuum model using a capillary pressure–saturation ($P_{c}$–$S_{g}$) relationship derived directly from the binary image of the micromodel pore geometry. Our approach follows the method proposed by Mehmani and Xu [19], a variant of the classical pore-morphology method (PMM) [38] adapted to capture Ostwald ripening dynamics. Let $I$ denote the binary image, where solid and void regions are represented by 0 and 1, respectively. Our PMM involves two successive morphological operations—erosion followed by dilation—using a circular structuring element $B_{r}$ of radius $r$. The erosion operation, $I\ominus B_{r}$, identifies all pixels for which the entire structuring element lies within the void space. The dilation operation, $I\oplus B_{r}$, subsequently expands these pixels outward by a distance $r$. In this framework, $r$ corresponds to the in-plane radius of curvature of gas–liquid interfaces, and thus to a specific capillary pressure via the Young–Laplace equation. To compute the saturation $S$ associated with $P_{c}=2\gamma/r$, we evaluate $(I\ominus B_{r})\oplus B_{r}$ and calculate the ratio of nonzero pixels in the resulting image to the total number of void pixels in $I$. The discrete $P_{c}$–$S$ data obtained from this modified PMM are then fit to a smooth, differentiable curve [19], enabling computation of $\mathrm{d}S_{g}/\mathrm{d}P_{c}$ in the continuum model (Eq. 2).

The key distinction between the classical PMM [38] and the variant of Mehmani and Xu [19] lies in their treatment of disconnected ganglia. Classical PMM excludes isolated clusters to ensure that the nonwetting phase remains connected to an inlet, making it well suited for modeling invasion percolation. In contrast, the ripening-specific variant used here retains disconnected ganglia, consistent with the absence of advective flow and the dominant role of local curvature in establishing capillary equilibrium. As a result, the $P_{c}$–$S$ relationship for ripening is effectively isotropic and independent of flow direction. In invasion-driven processes, by contrast, the $P_{c}$–$S$ curve is highly sensitive to inlet geometry: an inlet comprising a single narrow pore can produce a nearly flat $P_{c}$–$S$ curve, while broader or more connected inlets yield steeper relationships. In some cases, the initial distribution of trapped gas among pores of varying sizes can have a non-trivial influence on the resulting $P_{c}$–$S$ curve [13]. In such cases, PMM must be applied selectively to the initially occupied pores. For the pore geometry considered in this study, however, this sensitivity was found to be minor and is therefore neglected.

Acknowledgements.

Acknowledgements

This research is jointly supported by the Natural Sciences and Engineering Research Council of Canada (NSERC) Alliance Grants under ALLRP 592525-23 and ALLRP 567631-24 (for MS, TL, BZ) and the National Science Foundation, United States under Grant No. CBET-2348723 (for NB, ZL, YM).