Ostwald ripening has been studied extensively in bulk fluids, where bubbles are always spherical with curvature κ inversely proportional to radius R [9][10][11][12][13]. When comprised of only one component, mass transfer between bubbles is purely curvature-driven. A hallmark of such ripening is that all bubbles coarsen into a single large bubble, because any configuration with more than one bubble has further room to minimize surface energy by reducing bubble count. Extensions to multicomponent bubbles in bulk fluids have also been made, where inter-bubble mass transfer is more complex [14][15][16][17]. By contrast, studies of ripening in porous media are a recent development and have focused almost exclusively on single-component bubbles.

Early microfluidic experiments [18] and pore-network modeling [19,20] of single-component ripening revealed that equilibration is far more complex within the confinement of a porous medium. A bubble population in disequilibrium no longer evolves toward a single bubble, but toward another population with different bubble volumes V b and uniform curvature κ. This is because confinement exerts first-order control on κ, rendering it a non-monotonic function of V b -not so for spherical bubbles. Subsequent advances in pore-scale modeling [21][22][23][24], microfluidic experiments [25,26], and in-situ X-ray µCT observations [27][28][29] have provided further evidence of ripening in complex microstructures, shed insight into controlling factors (e.g., pore-size distribution, wettability, initial condition), and enabled higher-fidelity simulations. Theories developed for single-component ripening fall into two categories: pore-scale and continuum. Continuum theories [26,30,31] describe how trapped-bubble saturation evolves due to gradients in κ or capillary pressure P c at large scales. They have been used to assess the timescale and trajectory of bubble redistribution due to gravity and geothermal gradients [32][33][34][35]. Crucially, continuum theories assume local equilibrium within each representative elementary volume (REV). Pore-scale theories, by contrast, can capture the kinetics prior to local equilibrium and predict population statistics of bubbles, providing far more detail than a single aggregate quantity like saturation. The present work focuses on developing a pore-scale theory.

Pore-scale theories build upon the seminal works by Lifshitz & Slyozov [9] and Wagner [10] (LSW) for single-component bubbles in bulk fluids. The idea is to introduce a bubble number-density function g(R)

that describes the statistical distribution of bubble radii R and evolve it with a population balance equation:

Here, dR/dt is the rate of change of radius for bubbles with radius R, which we refer to as the phase velocity (i.e., velocity in the 1D statistical space of R). It is computed by assuming each bubble interacts with a mean field that represents the rest of the population. For bubbles with radius R, the phase velocity reads:

where x is the dissolved concentration next to the R-bubble and x mf is the mean-field concentration. Since x is proportional to κ (or P c ) by Henry’s law, and κ ∝ 1/R by Young-Laplace for spherical bubbles, the second proportionality in Eq.2 follows. Note x mf ∝ 1/R c , where R c is the critical radius below which bubbles shrink and above which they grow. R c is a property of the mean field and equals the average radius of all bubbles in the population, obtained by imposing mass conservation on the total mass of bubbles. The factor 1/R can be viewed as the “distance” of the mean field from an R-bubble. LSW further showed that g(R) approaches a self-similar regime described by g(R) = R -4 c f (R/R c ), where f (•) is a complicated function. The first extension of LSW theory to porous media was proposed by Yu et al. [36], under the assumption that trapped bubbles remain spherical throughout ripening. Therefore, the only difference from bulk fluids is that the pore network topology imposes a fixed distance across which bubbles interact. When this distance, l, is set to some multiple of the throat length l t , Eq.2 becomes:

although l continuously grows with time as bubbles dissolve and disappear. Yu and co-workers derived similar self-similar solutions for g(R) and showed that R c is now the harmonic average of all R in the population (again obtained by imposing mass conservation). Note that since bubble deformation by porous confinement is neglected, pore sizes R p do not factor into this theory, making it best suited for homogeneous pore networks.

Bueno et al. [37] formulated the first theory accounting for bubble deformation, pore-size heterogeneity, and spatial correlation. The extension was non-trivial because it required defining the bubble number-density g(s) in a 2D statistical space, where each bubble state s = (R p , S b ) has two coordinates: size of the occupied pore R p and volume fraction of the pore occupied S b (bubble saturation). This required expressing Eq.1 as:

where dR p /dt = 0 was assumed, entailing bubbles, once trapped in a pore, remain there forever or dissolvethey do not mobilize due to flow or capillarity. The phase velocity dS b s /dt for a bubble at state s is:

where κ s denotes the curvature of the s-bubble and κ mf s is that of its mean field. Notice each s has its own mean-field curvature, and that bubble radius R no longer enters the formulation because it is ill-defined for a deformed bubble. Instead, κ is used, which is a U-shaped function of S b (bubbles are spherical at low S b and deformed at high S b ). To understand Eq.5, note the left-hand side is proportional to dV b s /dt, where V b s denotes bubble volume and V b s ∝ R 3 p,s S b s . To account for spatial correlations in R p , the κ mf s was defined by:1

involving conditional probabilities:

The details of how these were derived are not important here. What matters is p(R ′ p |R p ; t), the probability that a bubble residing at a pore of size R p interacts with one at a pore of size R ′ p , is proportional to its initial condition at t = 0. Bueno and co-workers assumed that an s-bubble interacts only with its nearest neighbors, thereby neglecting interactions with distant bubbles. Hence, if two bubble groups reside in mutually distant pore sizes R p and R ′ p such that p(R ′ p |R p ; 0) = 0, then the clusters evolve independently, unaware of each other’s presence. This shortcoming was highlighted as a challenge in [37] for future theories to overcome. We also note that this theory did not conserve total bubble mass for networks with correlated pore size.

All above theories are for single-component bubbles. The only theoretical work on multicomponent bubbles trapped in porous media is that of Bueno et al. [38], who successfully described g(s) at equilibrium, now defined in a 3D statistical space s = (R p , S b , y) where y is the mole fraction of bubbles. The kinetics, or timescale to reach equilibrium, could not be estimated. This theory built on an earlier equilibrium theory for single component bubbles by Mehmani & Xu [20], where the statistical space was 2D with s = (R p , S b ).

In this work, we generalize the single-component kinetic theory of [37] to multicomponent bubble populations trapped in porous media with heterogeneous and spatially correlated pore sizes. The proposed theory also addresses the challenge of remote-bubble interactions, ensures total bubble-mass conservation through provable theorems, and features a simpler formulation, making it superior to [37]. We systematically validate our theory against the pore-network model (PNM) of [38] across 19 cases spanning homogeneous, heterogeneous, correlated, and uncorrelated networks with different bubble populations relevant to UHS. Good agreement is observed throughout. Pending challenges and limitations are discussed subsequently.

The paper is organized as follows: Section 2 describes the problem we aim to solve, its conceptualization, and assumptions. Section 3 briefly reviews the PNM of [38] used to validate our theory. Section 4 details the theory’s formulation. Section 5 discusses the validation test cases, and Section 6 presents the results. We discuss key limitations and implications for UHS in Section 7 and summarize main findings in Section 8.

We consider a porous medium with a rigid microstructure, as shown in Fig. 1. The void space is conceptualized as a pore network, a computational graph consisting of nodes (or pores) and links (or throats).

Pores provide storage volume for fluids, while throats impose mass-transfer resistance between them. Pore and throat sizes may be homogeneous or heterogeneous, and spatially correlated or uncorrelated. We assume the void space is occupied by a perfectly wetting (zero contact angle) liquid, and some pores contain a compressible gaseous bubble. The gas consists of two components that are partially soluble in the wetting phase. For concreteness, and without loss of generality, we consider bubbles to be binary mixtures of H 2 and CO 2 immersed in water as the wetting phase. The gas is assumed to be ideal and devoid of H 2 O, and the dissolved fractions of H 2 and CO 2 in water to be dilute. Bubbles exchange mass through Ostwald ripening, where concentration gradients of dissolved species drive Fickian diffusion. Such gradients exist because of differences in the species concentrations surrounding each bubble, which are dictated by the composition and curvature of each bubble. We assume each bubble occupies only one pore, and that pore shapes (but not sizes) are identical throughout the network. Lastly, we assume the interfacial tension between liquid and gas is not altered by bubble composition or dissolved concentrations. Our goal is to develop a theory that predicts the temporal evolution of bubble statistics, including composition, volume, and occupied-pore size.

Figure 1: The porous microstructure is composed of a solid matrix and an interconnected void space. We conceptualize the void as a pore network, a computational graph that consists of nodes (pores) and links (throats). All pores contain a wetting phase (H 2 O), and some are occupied by partially miscible, two-component (CO 2 and H 2 ) gaseous bubbles. Bubbles are confined to one pore and exchange mass via Ostwald ripening, driven by concentration gradients of dissolved species around each bubble.

We briefly review the pore-network model (PNM) of Bueno et al. [38], which we use to validate our theory for multicomponent Ostwald ripening. The two share many of the same equations. The governing equations of the PNM consist of a mass balance imposed at each pore i and for each species α:

The variables n α,i , J ij , and z i denote the number of moles of species α in pore i, the diffusive flux (mol/s) from neighboring pore j to pore i, and the number of pores connected to pore i by a throat, respectively.

To the right of the second equality, D α is the molecular diffusion coefficient, ρ w the molar density of water, and x α,i the mole fraction of species α dissolved in the wetting phase within pore i. The A ij and L ij are the cross-sectional area and length of the throat connecting pores i and j, through which species diffuse.

To transform Eq.8 into a system of ordinary differential equations (ODEs) in terms of the unknown x α,i , we need constitutive relations expressing n α,i in terms of x α,i . The first such relation is

, provided pore i is occupied by a gaseous bubble. The first term accounts for moles of α in the bubble, and the second term for dissolved moles in the wetting phase. Here, V p,i , V b i , ρ b i , and y α,i are the pore volume, bubble volume, gas molar density, and mole fraction of α in the gas, respectively. If pore i is vacant, then V b i = 0 and the expression simplifies to n α,i = x α,i ρ w V p,i . Recall bubbles occupy only one pore. The second constitutive relation expresses y α,i in terms of x α,i by assuming each bubble is in local equilibrium with dissolved species in its pore. Under the dilute-solution assumption (Section 2), equilibrium takes the form of Henry’s law: x α,i H α = y α,i P b i , where H α is Henry’s constant and P b i the bubble pressure. P b i is related to the wetting-phase pressure P w through Young-Laplace: P b i = P w + σκ i , where σ is interfacial tension and κ i the bubble curvature in pore i. Since the wetting phase is stagnant (Section 2), P w is constant.

The third relation expresses ρ b i in terms of P b i using the ideal-gas law P b i = ρ b i RT , per our assumption in Section 2. Here, R and T are the universal gas constant and absolute temperature, respectively. The final constitutive equation we need relates κ i to V b i . For a semi-cubic pore i, defined as one whose maximal inscribed radius is R p,i and volume satisfies V p,i = (2R p,i ) 3 , this relationship takes the form [21]:

where V c,i = 4πR 3 p,i /3 is the volume of the maximal inscribed sphere, and κ min i = 2/R p,i the corresponding curvature. Eq.9 is non-monotonic, with κ i attaining its minimum

, the bubble is small and spherical. For V b i ≥ V c,i , it is large and deformed by the pore’s confinement. The parameters a and b control this deformation and depend on pore shape (here a = 1 and b = 0.01). Following [38], we refer to V c,i as the critical volume, bubbles with V b i < V c,i as sub-critical, and those with V b i > V c,i as super-critical. The species balance Eq.8 combined with the above constitutive relations yields a system of nonlinear ODEs, which we solve using a modified Newton scheme with adaptive time stepping [20]. The modification avoids instabilities and divergence of Newton iterations near bubble disappearance events or when V b i crosses the critical volume V c,i , where the derivative dκ i /dV b i is discontinuous. We also refer the reader to [38] for a lower bound on V b i , below which bubbles must be marked as “vanished” to prevent numerical instabilities.

We formulate a kinetic theory that predicts the temporal evolution of the statistics of a bubble population driven by two-component Ostwald ripening. The PNM of Section 3 models this evolution spatially, meaning each bubble is tracked individually and assigned a spatial position, while interacting with neighboring bubbles at other positions. In our theory, we do not track individual bubbles, but groups of bubbles sharing the same statistical state, s. We let this state s = (R p , S b , y 1 ) consist of three coordinates in a three-dimensional statistical space, referred to as the phase space. The first coordinate is the pore size R p that bubbles at state s occupy. The second is bubble saturation S b = V b /V p , which denotes the volume fraction of the host pore each bubble in the group occupies. The third coordinate is the mole fraction of H 2 for each bubble in the group, denoted by y 1 . Since bubbles are binary mixtures, the mole fraction of CO 2 follows from y 2 = 1 -y 1 .

To quantify the number of bubbles in each state s at any given time t, we introduce the number-density function (NDF) g(s; t). Concretely, the number of bubbles in an infinitesimal volume of phase space dΩ = dR p dS b dy 1 centered at s equals g(s; t) dΩ. Figs.2-3 illustrate this conceptualization. Fig. 2a depicts the 3D phase space with g(s; t) visualized as a cloud, along with its projections (or marginals) onto the coordinate planes. Fig. 3 provides a simple example of how spatial information in the PNM is translated into statistical information in the theory. The left panel is a schematic of pure H 2 (cyan; y 1 = 1) and CO 2 (red;

bubbles residing in pores of different sizes (R p ) with different saturations (S b ). The remaining three panels visualize 2D projections of the corresponding g(s; t) onto the coordinate planes. Bubbles from the spatial picture are overlain on the phase-space plots to highlight which part of g(s; t) they represent. Pure-species bubbles (y 1 = 0 or 1) are used for simplicity, as bubbles are generally mixtures with arbitrary mole fraction

). The shadows on each plot represent 2D marginals of g(s; t), and the pink curves along the axes represent 1D marginals. The goal of our theory is to predict the temporal evolution of g(s; t).

We describe the evolution of g(s; t) in time with the following population balance equation:

where u = (u Rp , u S b , u y1 ) = ( dRp dt , dS b dt , dy1 dt ) is the velocity that advects g(s; t) through phase space. This velocity is determined purely by the physics of Ostwald ripening for a population of gaseous bubbles composed of two components (here H 2 and CO 2 ). We assume u Rp = 0, implying that bubbles cannot change pore sizes phase space (Fig. 2). Bubbles from the first panel are annotated on these phase-space plots to emphasize which part of g(s; t) they represent. Shadows represent 2D marginals of g(s; t), and pink curves along the right/bottom axes the 1D marginals.

by, for example, mobilizing from one pore to the next due to background flow or capillary action. 2 This is consistent with our assumption in Section 2 that each bubble occupies only one pore and the wetting phase is stagnant (i.e., no flow). Computing u = (0, u S b , u y1 ) closes the theory, which we address in Section 4.2.

Given some initial condition (IC), g(s; t = 0), solving Eq.10 requires boundary conditions (BCs). Since we have assumed u Rp = 0, no BCs are needed at R p = 0 or R p = R max p , where R max p is the maximum pore size in the network. At y 1 = 0 and y 1 = 1, the normal velocity u y1 always points toward the interior of the phase space, preventing bubbles from ever crossing either boundary. Otherwise, mole fractions would fall outside the valid interval [0, 1]. To also prevent new bubbles from being introduced (or “injected”), we set g = 0, or equivalently u y1 = 0, at y 1 = 0 and y 1 = 1. Following similar arguments, we set g = 0, or equivalently u S b = 0, at S b = 1. The only tricky boundary is S b = 0, because shrinking bubbles during ripening can indeed vanish and leave the phase space, in which case u S b < 0 holds. Since Eq.10 is first-order hyperbolic, no BC is needed at such outflow portions of the boundary-numerically, the upwind value of g gets used. If, however, u S b > 0 holds, then we set g = 0 at such inflow portions to prevent introducing new bubbles into the population.

To determine the phase-space velocity u, we must compute dS b /dt and dy 1 /dt. We adopt a mean-field approximation that reduces the problem of a bubble interacting with multiple neighboring bubbles to the simpler problem of interaction with a single entity: the mean field. Consider a bubble at state s = (R p , S b , y 1 ), as shown in Fig. 2b, that exchanges mass with surrounding bubbles. On average, its nearest neighbors lie a distance l s away. We therefore idealize the region r < l s as bubble-free, which is a disk/sphere of radius r = l s in 2D/3D. All bubbles in r > l s define the mean field, to which average properties are assigned, such as dissolved concentration x mf α,s , bubble pressure (P b ) mf s , and gas composition y mf α,s . Mean-field properties are computed from g(s; t), as shown later. Since bubbles can vanish during ripening, l s grows with time. For tractability, we assume the bubble-spacing parameter l s depends primarily on the pore-size coordinate R p , with S b and y 1 playing secondary roles. This is reasonable because spatial correlations in R p dictate where bubbles are initially trapped in the network (e.g., by cyclic injections of H 2 ) and how curvature-driven mass transfer evolves. We therefore write l R hereafter to make this dependence explicit (Fig. 2). Note the mean-field approximation requires a sufficiently large bubble population for ensemble statistics to have converged.

We now derive expressions for dS b s /dt and dy 1,s /dt. To compute mass exchange between the s-bubble and the mean field (Fig. 2b), we idealize transport in the wetting phase as radial diffusion toward a spherical (3D) or cylindrical (2D) surface at r = l 0 , representing the bubble’s interface (specified later). The annular region from r = l 0 to l R is bubble-free, with mass transfer governed by Fick’s second law in radial coordinates:

Eq.11 assumes steady-state, radially symmetric transport, implying the wetting phase acts merely as a solute transport medium, not storage. Here, d is the network dimension (2 or 3), and D e α = D α ϕ a /τ the restricted diffusion coefficient, which depends on molecular diffusivity D α , tortuosity τ , and areal porosity ϕ a . Appendix A shows ϕ a = CA t /l t , with C = z/(2h) in 2D and C = z/(2l t ) in 3D. Here, z is the mean coordination number (neighbor count) of pores, h the network’s out-of-plane thickness in 2D, and A t and l t are the mean crosssectional area and length of throats. Solving Eq.11 with BCs x α (l 0 ) = x α,s and x α (l R ) = x mf α,s yields:

Recall x α,s is the dissolved concentration at the interface of the s-bubble.

Given Eq.12, we write a mole balance for the bubble at state s:

where j is the molar flux of dissolved solute toward/away from the s-bubble, and A the surface area of the cylinder/sphere through which transport occurs. Note the product jA does not depend on r. Differentiating Eq.12 with respect to r and substituting into Eq.13, along with the expression for ϕ a above, yields:

We next relate dn α,s /dt to dS b s /dt and dy 1,s /dt. Following Section 3, the number of moles in the s-bubble

, where all variables retain their original meaning, except the pore index i is now replaced with the bubble state coordinate s. Invoking the same constitutive equations in Section 3, it is straightforward to see n α,s is fully determined by s = (R p , S b , y 1 ). Recall these constitutive equations include the ideal-gas law ρ b s RT = P b s , the Young-Laplace equation P b s = P w + σκ s (with P w assumed constant), and the curvature-volume relation κ s (V b s , V p,s ) in Eq.9 for the pore size R p occupied by the s-bubble. As a result, the time derivative of n α,s = n α,s (R p , S b , y 1 ) can be expanded using the chain rule:

where the first term in Eq.15a is zero since dR p,s /dt = 0 was assumed in Section 4.1. Eq.15b follows from substituting the above constitutive relations into the remaining two terms in Eq.15a and simplifying.

We now substitute Eq.15 into the left-hand side (LHS) of Eq.14, which results in two equations corresponding to α ∈ {1, 2}. To express this 2×2 ODE system in terms of only S b s and y 1,s as our primary variables, we use Henry’s law x α,s = y α,s P b s /H α to replace x α,s on the right-hand side (RHS) of Eq.14. This yields:

where

In Eq.16a, we have used dy 2 /dt = -dy 1 /dt, which holds for binary mixtures.

Right-multiplying both sides of Eq.16a by the inverse of the mass matrix on the LHS, we obtain:

which are explicit expressions for dS b s /dt and dy 1,s /dt as intended. The Kronecker delta δ 1α assumes a value of one for α = 1 and zero for α = 2. We note (y α P b ) mf s ̸ = y mf α,s (P b ) mf s and that (y α P b ) mf s = x mf α,s /H α holds by Henry’s law. The latter implies that the solution of Eq.17 at large t asymptotes as x α,s → x mf α,s for all α and s. Unfortunately, the mathematical form of Eq.17 is numerically problematic when the population balance Eq.10 is solved with an Eulerian (grid-based) method, as done here with finite volumes: the curves dS b s /dt = 0 and dy 1,s /dt = 0 form a narrowing corridor in phase space that traps bubbles before they reach equilibrium (see Appendix B for details). An additional problem is that each state s interacts with a different mean-field

x mf α,s , and there is nothing in the expressions provided later that forces these mean-fields to converge to the same state at t → ∞. To address these challenges, we alter the driving force in Eq.17b as follows:

Eq.18 constitutes our final expressions for dS b s /dt and dy 1,s /dt. Notice the RHS of Eq.18b is now driven by the difference between the s-bubble composition, y 1,s , and the mean field, y mf 1 . Crucially, y mf 1 does not depend on s, ensuring all bubbles approach the same composition at t → ∞. We prove in Section 4.4 that because our theory conserves mole exchange between bubbles, y mf 1 → y mix 1 at large t, where y mix 1 is the bulk (or mixture) composition of all bubbles in the population. Therefore, bubbles are guaranteed to converge to a physically meaningful equilibrium composition. Notice Eq.18 drives x α,s and y 1,s towards x mf α,s and y mf 1 for all s. To derive Eq.18b, we have assumed P w ≫ σκ s for all s, which is valid for high-pressure subsurface environments (e.g., UHS). This allows setting P b ≈ P w and pulling it out of the summation in Eq.17b. We caution that applying the same simplification to Eq.17a as well would be fatal, as the ODE system would then be driven solely by differences in bubble composition (i.e., y mf 1 -y 1,s ), not curvature.

We close this section by returning to the inner radius l 0 in Eq.12 and deriving an expression for it. Notice the smallest l R can be is a single throat length l t , corresponding to when all neighboring pores of the s-bubble are occupied by bubbles. Since there are z such neighbors, we can write analogously to Eq.8 in the PNM:

The first equality follows from Eq.8, where by radial symmetry the dissolved concentration within all neighboring pores is x α (l t ). The second equality follows from substituting for x α (l t ) the expression from Eq.12 with r = l R = l t . Comparing Eq.19 to Eq.14, we arrive at the following formula for l 0 :

This, in turn, yields the expression below for the geometric factor G s from Eq.16b:

Note G s = z when l R = l t , which holds for a densely bubble-occupied network, and it decreases with l R . In summary, the closure equations from this section are Eqs.18 and 21. To complete the theory, we now need expressions that allow us to compute (y α P b ) mf s and y mf 1 at any given time t. This is addressed next.

In this section, we postulate explicit expressions for (y α P b ) mf s and y mf 1 . We begin with (y α P b ) mf s and propose the following, which accounts for potential spatial correlations in pore size R p in the network:

The p(s ′ , l|s) is the conditional probability density function (PDF) of finding a bubble at state s ′ a distance l from the bubble at state s (see Fig. 2b). Eq.22 is essentially a weighted average of (y α P b ) s ′ , where the weights p(s ′ , l|s) / l decline with distance from the s-bubble according to spatial correlations in bubble states.

To compute p(s ′ , l|s), we first apply Bayes’ theorem:

and set p(l|s ′ , s) = δ(l -l ss ′ ), where δ(•) is the Dirac delta function. This expression for p(l|s ′ , s) is consistent with our theory’s conceptualization in Fig. 2b, meaning that once s and s ′ are chosen, the interaction distance between them is fixed at l ss ′ . Recall the space between bubbles s and s ′ is free of other bubbles. To estimate p(s ′ |s) and l ss ′ , we assume much of the spatial correlation in bubble states is dictated pore size (i.e., the R p coordinate of s). This justifies renaming l ss ′ to l RR ′ (Fig. 3b) and approximating p(s ′ |s) as follows:

where the second equality follows from expanding s ′ in terms of its coordinates and applying Bayes’ theorem.

The first PDF on the RHS of Eq.24 can be computed exactly from the NDF g(s ′ ; t) as follows:

where t is suppressed for brevity. The second PDF on the RHS of Eq.24 is approximated as:

where f (R ′ p , R p ; l) is the joint PDF of pore sizes R ′ p and R p separated by a distance l. This PDF describes the pore network itself (not bubble occupancy) and is assumed to be known; it can be computed from Xray images of porous samples. The smaller the distance l, the larger the value of f (R ′ p , R p ; l). The marginal

) for all l. The factor F R ′ represents the fraction of bubble-occupied pores in the network that have sizes between R ′ p and R ′ p +δR ′ p . Its multiplication by f (R ′ p |R p ; l RR ′ ) accounts for bubble occupancy, since p(R ′ p |R p ) represents the probability of two bubble-occupied pores interacting (not just any two pores). The denominator in Eq.26 ensures p(R ′ p |R p ) integrates to one. F R ′ can be computed exactly from g(s ′ ; t) and f (R ′ p ) as follows:

where n p denotes the total number of pores in the network.

Substituting Eqs.23-24 into Eq.22, and relabeling the variables l ss ′ and l s as l RR ′ and l R to explicitly signal our assumption that they depend solely on correlations in R p (not S b or y 1 ), we obtain:

as the final expression for (y α P b ) mf s , with p(S b ′ , y ′ 1 |R ′ p ) and p(R ′ p |R p ) evaluated from Eqs.25-27. The bubble spacing l RR ′ in Eq.28 grows as bubbles dissolve during ripening. To compute it, we propose:

where d RR ′ is the average distance of pores with size R p from their nearest pore with size R ′ p . Thus, d RR ′ is a property of the network itself (not bubble occupancy). For most networks, d RR ′ ≈ l t is a good approximation, but more accurate values can be computed from an X-ray image. Recall d is the network dimension and F R ′ is given by Eq.27. While empirical, Eq.29 has a few desirable properties. In homogeneous networks, where R p assumes a single value (R ′ p = R p ), Eq.29 reduces to l RR = L/n is the spacing between bubbles if spread uniformly across the network. In heterogeneous networks, when all pores of size R ′ p are bubble-occupied, l RR ′ = d RR ′ as expected. As bubbles dissolve due to ripening, F R ′ decreases from one to zero, causing l RR ′ to increase.

In Appendix C, we compare l RR ′ computed from Eq.29 with direct PNM simulations and find acceptable agreement in heterogeneous networks with both correlated and uncorrelated pore sizes. 3 Having formulated (y α P b ) mf s , we now postulate the following expression for y mf 1 :

As a bubble population described by g(s; t) evolves due to Ostwald ripening, the number of moles of each species must be conserved. Specifically, moles leaving one bubble must enter another, keeping the total moles of each species constant. Recall the wetting phase in our theory has no storage capacity for dissolved solute.

In Appendix D, we prove the theorem below, establishing a condition equivalent to mass conservation:

Theorem 1. The following two expressions are equivalent, where Ω is the entire phase space:

In words, the left expression states that the total number of moles of species α is constant for all t.

Setting α = 1 and substituting Eq.15 for dn 1,s /dt into the rightmost expression in Theorem 1, with dS b s /dt and dy 1,s /dt given by Eq.18, we obtain:

  g(s; t) ds = 0 (32) 3 We make two technical remarks about l RR ′ : (1) Eq.29 ignores any dependence on F R , the bubble occupancy of pores with size Rp. Hence, it lacks the necessary symmetry to be considered a “distance metric” for bubble spacing; (2) For all choices of Rp and R ′ p , the inequality l RR ′ ≥ lt > l 0 holds as can be verified from Eqs.20 and 29 and the fact that d RR ′ ≥ lt.

Unfortunately, the (y α P b ) mf s and y mf 1 postulated by Eqs.28 and 30 in the previous section do not guarantee that Eq.32 is satisfied. Therefore, a correction to these mean-field properties is needed to conserve bubble mass. We propose the following multiplicative corrections:

where the asterisk indicates a corrected mean-field property.

To compute the correction multiplier λ, we write an equation similar to Eq.32 for α = 2, add it to Eq.32, use the fact that y 1,s + y 2,s = 1 and y mf 1 + y mf 2 = 1, and substitute the corrected mean-field properties in Eq.33 into the resulting expression. After some algebra, this yields:

Given λ from Eq.34, we can compute the second correction multiplier γ by substituting the corrected meanfield properties in Eq.33 into Eq.32 and solving for γ. This yields:

To summarize, we first compute the uncorrected mean-field properties (y α P b ) mf s and y mf 1 via Eqs.28 and 30, then correct them via Eqs.33-35. These corrected mean-field properties are used in Eq.18 to compute the phase-space velocities dS b s /dt and dy 1,s /dt needed in the population balance equation Eq.10. We close this section with a corollary to Theorem 1, which guarantees that bubbles converge to a physically meaningful equilibrium composition. The proof is provided in Appendix E.

Corollary 1. The bulk (or mixture) composition of all bubbles in the population, defined as:

remains constant for all t. Moreover, as t → ∞, the mean-field mole-fraction satisfies y mf, *

This result further supports our choice of the s-independent y mf 1 in Eq.30 and the reformulation of dy 1,s /dt in Eq.18. By contrast, the precursor ODE in Eq.17 does not guarantee convergence to y mix

.

In Appendix E, we also show γ → 1 as t → ∞, implying y mf 1 → y mix

. Finally, we offer the following remark about λ:

Remark 1. While it is obvious that γ corrects the mean-field composition y mf 1 , the use of a single multiplier λ for both species α ∈ {1, 2} is equivalent to correcting the mean-field bubble curvature κ mf s . To see this, note the mean-field operator (•) mf s is linear, so by Henry’s law and the Young-Laplace equation

s . This is why a single multiplier λ suffices for both species.

We validate the theory in Section 4 against PNM simulations for both one-and two-component bubble populations. We consider pore networks with a 2D lattice topology but 3D pore and throat shapes. Pores are semi-cubic of size R p , defined in Section 3, and throats are prisms with cross-sectional area A t and length We validate our theory in two stages. In the first stage, we consider single-component bubbles composed of pure H 2 . The goal is to demonstrate that our theory not only generalizes the formulation of [37], but also removes key limitations that enhance predictive accuracy. In the second stage, we consider two-component bubbles composed of H 2 and CO 2 . Table 1 lists the cases considered in each stage, which we detail next. We In stage 1 of our validation, we consider cases in the top segment of Table 1, all taken from [37]. The first three correspond to pure H 2 bubbles (y 1 = 1) occupying all pores in the network at fixed local saturation S b . We label these Hom-One-Full, Unc-One-Full, and Cor-One-Full, depending on whether the network is Hom, Unc, or Cor. In Unc-One-Full and Cor-One-Full, local saturation S b is uniform across all pores (thus

), whereas in Hom-One-Full, S b is randomly drawn from [ 1, 5 ]%. These three cases correspond to the more difficult problems in [37], which we use to demonstrate that our theory reproduces those results.

The fourth case, Cor-One-Part, corresponds to a scenario presented as a challenge in [37] (Section 7.2), whose ripening dynamics that theory could not capture due to large initial bubble spacing. We show our theory overcomes this limitation, thereby generalizing [37]. The setup places sub-critical bubbles in the largest 20% of pore sizes, with local saturation sampled uniformly from S b ∈ [ 5, 35 ]%, and super-critical bubbles in the smallest 20%, with S b sampled uniformly from [ 65, 95 ]%. Here too, all bubbles are composed of pure H 2 .

In stage 2 of the validation, we consider cases in segments 2-6 of Table 1. These comprise different combinations of network type (Hom, Unc, Cor ), bubble population (Pop1, Pop2, Pop3 ), and spatial distribution of y 1 (yRand, ySmall ). Since network types were described above, we focus here on the three bubble populations: In Pop1, all pores contain bubbles at local saturation S b = 60%, with roughly half by volume being pure H 2 and the other half pure CO 2 . Pop2 is similar, but the pure CO 2 bubbles are replaced by bubbles that are an equimolar H 2 -CO 2 mixture. In Pop3, roughly half of the bubbles by volume are pure H 2 at S b = 60%, and the other half are pure CO 2 at S b = 30%. For each population, the labels yRand and ySmall specify how bubbles are assigned to pores. In yRand, pores are chosen at random, implying no correlation between pore size R p and bubble composition y 1 . In ySmall, the smallest pores are filled with H 2 first, then the remaining pores are filled with other bubbles (pure CO 2 in Pop1 and Pop3, equimolar mixture in Pop2 ).

In the Cor network, ySmall produces spatial correlation in y 1 , not just R p . Hereafter, we refer to each of the resulting 15 cases by the code names listed in Table 1 (e.g., Cor-Pop2-yRand ).

The above initial conditions for two-component ripening are motivated by underground hydrogen storage.

Cyclic injections of H 2 are often preceded by a cushion gas like CO 2 , which can mix with subsequent H 2 injections. In addition, cyclic injections can cause capillary trapping of bubbles within preferential flow paths (larger pores), diverting subsequent H 2 into smaller pores. This motivates the equimolar mixture in Pop2 and the ySmall placement. We note that placing H 2 bubbles in larger pores would be a less challenging test of our theory, since early-stage mass transfer is exclusively from CO 2 to H 2 due to differences in solubility and diffusion coefficient, causing more bubble deformation when H 2 resides in smaller pores.

In all cases, the temperature T is set to 333 K, wetting-phase pressure P w to 10 8 dynes/cm 2 , surface tension σ to 32 dynes/cm, molecular diffusion coefficients to D 1 = 5.13 × 10 -5 cm 2 /s for H 2 and D 2 = 1.6 × 10 -5 cm 2 /s for CO 2 , and Henry’s constants to H 1 = 7.1 × 10 10 dynes/cm 2 for H 2 and H 2 = 0.16 × 10 10 dynes/cm 2 for CO 2 . For each case in Table 1, the temporal evolution of the bubble population due to Ostwald ripening is predicted by the theory and compared to PNM simulation. The population balance Eq.10 is discretized over the phase space using the finite-volume method with a first-order upwind scheme. Although Eq.10 must be solved numerically, it is far more efficient than PNM because there is no limit to the number of bubbles that can be simulated. Moreover, networks with large correlation lengths in pore size require prohibitively large representative elementary volumes (REVs) in the PNM, which is not a limitation for the theory.

We begin by comparing the theory to PNM for the one-component ripening cases listed in Table 1.

Focusing on Hom-One-Full, Fig. 5 shows the theory successfully predicts the temporal evolution of g(s; t).

Since the network is homogeneous with uniform R p , and bubbles are composed of one component with y 1 = 1, the phase space is one-dimensional, consisting of only the S b axis. Focusing next on Unc-One-Full and Cor-One-Full, Fig. 6 compares the evolution of g(s; t) from the theory to the PNM, where good agreement is again observed. The PNM is depited by a heat map, and the theory by white contour lines. Here, the phase space is two-dimensional with axes S b and R p , because y 1 = 1. At t = 0, g(s; t) is a vertical stripe since all bubbles have the same saturation. Over time, the higher curvature in smaller pores drives mass flux toward larger pores, causing g(s; t) to shift left (toward S b = 0) at small R p and right at large R p . At equilibrium, curvature is uniform across all bubbles, producing the oblique line in Fig. 6. The results in Figs. 56for Hom-One-Full, Unc-One-Full, and Cor-One-Full were also successfully predicted by the theory of [37].

Figure 5: Comparison of the bubble number density g(s; t) predicted by theory (white contours) and PNM (heat map). Since bubbles are single-component (y 1 = 1) and network is homogeneous, the phase space is one-dimensional with S b axis. Snapshots are shown at t = 0, 0.1, 10 h, and equilibrium for the Hom-One-Full case in Table 1. The theory reproduces results from [37].

Figure 6: Comparison of the bubble number density g(s; t) predicted by theory (white contours) and PNM (heat map). Since bubbles are single-component (y 1 = 1), the phase space is two-dimensional with Rp and S b axes. Snapshots are shown at t = 0, 10, 100 h, and equilibrium for the Unc-One-Full and Cor-One-Full cases in Table 1. The theory reproduces results from [37].

We now turn to Cor-One-Part, which was presented as a challenge in [37]. Fig. 7 shows the complex initial 1. The proposed theory (solid lines) shows good agreement with PNM (symbols). Predictions from [37] (dashed lines) are slightly worse for Hom-One-Full, Unc-One-Full, and Cor-One-Full, but completely off for Cor-One-Part. condition of g(s; t) and its subsequent evolution in time. The top row corresponds to the theory of [37], which fails to predict this evolution, resulting in an isolated cluster of bubbles in the bottom-right of the phase space that remains disconnected or “unaware” of other bubbles in the population. This failure is due to the fundamental inability of that theory to capture interactions between bubbles that are multiple throats apart.

Our theory in the bottom row of Fig. 7 removes this limitation, and is therefore superior to that of [37]. Fig. 8 shows the corresponding bulk properties of the above populations versus time. They consist of the average bubble volume ⟨V b ⟩, mean curvature ⟨κ⟩, and survived fraction of bubbles F s (i.e., current number of bubbles divided by the number of bubbles at t = 0). These properties are readily computed from g(s; t):

Good agreement between our theory and PNM is observed in all cases. Fig. 8 also includes predictions from [37], which are slightly worse for the first three cases (Hom-One-Full, Unc-One-Full, Cor-One-Full ) but completely off for Cor-One-Part. In summary, Figs.5-8 demonstrate that our theory not only successfully reproduces results from the theory of [37] for single-component ripening but also generalizes it by removing key limitations. The next section takes this further by capturing two-component ripening dynamics.

6.2. Two-component ripening

We begin our comparison of the theory against PNM for two-component ripening by considering cases in the second segment of Table 1, namely, Hom-Pop1-yRand, Hom-Pop2-yRand, and Hom-Pop3-yRand. These correspond to the homogeneous network with random placement of bubbles at t = 0, i.e., no correlation between y 1 and R p . Fig. 9 shows the evolution of bubble number density g(s; t) versus time, where we see good agreement between theory and PNM. Since R p is uniform, the phase space is two-dimensional with y 1 and S b axes. At t = 0, H is driven by curvature differences in a stage called coevolution [38], causing the two groups to approach each other along the constant-y 1 line. At equilibrium, only one cluster remains, with all bubbles having the same composition and curvature (thus saturation, since R p is uniform). In Fig. 9, the difference in S b at equilibrium between the theory and PNM is 4-5%, and the difference in y 1 is 0.5% for Pop1 and Pop2, but 14% for Pop3.

The larger y 1 discrepancy in Pop3 arises because its CO 2 bubbles are sub-critical (S b = 30%) and reside in large pores, leaving substantial water volume that absorbs H 2 during partitioning. The PNM captures this, lowering the equilibrium y 1 predicted, whereas the theory does not as it assigns no storage capacity to water. Fig. 10 shows bulk properties ⟨V b ⟩, ⟨κ⟩, F s , and mean composition ⟨y 1 ⟩ versus time, where ⟨y 1 ⟩ is computed analogously to Eq.37. Note this is identical to y mf 1 in Eq.30. The agreement is decent except for early stages of ⟨V b ⟩, and late stages of F s and ⟨y 1 ⟩ for Pop3. Focusing on ⟨V b ⟩, we see a temporary sag in PNM that is present but much less pronounced in the theory. This sag in PNM has two causes: (1) The wetting phase can store dissolved solute, which introduces a lag in mass transfer between H 2 and CO 2 /mixture bubbles on account of differences in solubility and diffusion coefficient. This effect is not captured by the theory, as mass exchange between bubbles is instantaneous; (2) The rapid initial influx of mass into H 2 bubbles causes them to grow and become deformed, resulting in the observed peak in ⟨κ⟩. Since bubbles are gaseous, their compressibility causes a temporary decrease in mean volume, before deformed bubbles relax back to lower curvatures. This effect is captured by our theory, which is why the spike in ⟨κ⟩ is predicted accurately.

Focusing next on F s and ⟨y 1 ⟩, we see good agreement for Pop1 and Pop2, but significant error for Pop3.

The ⟨y 1 ⟩ deviation in Pop3 occurs early (∼ 10 -4 s), consistent with the water storage effect discussed above.

The F s deviation occurs later (∼ 10 0 s), after coevolution is complete, and reflects a different mechanism: in Pop3, CO 2 bubbles are sub-critical (spherical) from the start and remain so as they shrink during partitioning.

The subsequent curvature-driven ripening of such spherical bubbles is unstable [38], with survival depending on local neighborhood configuration. In the theory, sub-critical bubbles interact only with the mean field, through which they tend to lose mass to super-critical (deformed) bubbles. In PNM, sub-critical bubbles also exchange mass with neighboring sub-critical bubbles, some of which gain mass at the expense of others.

This allows more bubbles to survive in PNM than in the theory, explaining why the theory underpredicts F s .

We now focus on cases in the third segment of Table 1, namely, Unc-Pop1-yRand, Unc-Pop2-yRand, and Unc-Pop3-yRand. These correspond to the heterogeneous and uncorrelated network with random placement of bubbles at t = 0, i.e., no correlation between y 1 and R p . Now, g(s; t) lives in a three-dimensional phase space. To visualize it, we consider three projections (or marginals) of g(s; t) onto the S b -y 1 , R p -S b , and R p -y 1 planes. The S b -y 1 projection is shown in Fig. 11 at t = 0, 1 s, 10 s, and equilibrium. The evolution of g(s; t) is very similar to Fig. 9, where H 2 and CO 2 /mixture bubbles first approach a line of constant y 1 , then approach each other whilst staying close to that line. One difference compared to the homogeneous network in the previous section is that the evolution here is slower. Another is that at equilibrium, bubbles attain different S b to ensure uniform curvature due to R p variability. Overall, the theory is in good agreement with PNM. The largest discrepancy corresponds to Pop3 for reasons already discussed in Section 6.2.1.

The corresponding R p -S b projection is shown in Fig. 12. At t = 0, g(s; t) appears as a vertical stripe in Pop1 and Pop2 because S b = 60% is uniform across all bubbles, but as two stripes in Pop3 because S b = 60% holds for H 2 bubbles and S b = 30% for CO 2 bubbles. As time progresses, g(s; t) splits into two stripes corresponding to the CO 2 /mixture (left) and H 2 (right) bubbles, since each pore size is occupied by both bubble types at t = 0. Once bubble compositions have homogenized, the two stripes approach each other and merge into an oblique line similar to that in Fig. 6. This is because for the same curvature, smaller pores 1.

Network is heterogeneous and uncorrelated with random bubble placement. Snapshots are at t = 0, 1, 10 s, and equilibrium.

require less bubble deformation, thus lower S b . Here too, the agreement between theory and PNM is good.

The observations and takeaways remain similar for the R p -y 1 projection shown in Fig. 13.

Finally, Fig. 14 shows how bulk properties ⟨V b ⟩, ⟨κ⟩, F s , and ⟨y 1 ⟩ evolve with time. Overall, theory and PNM are in good agreement, with larger deviations observed in Pop3. Similar to Fig. 10, the ⟨y 1 ⟩ deviation in Pop3 arises because PNM accounts for solute storage in water, whereas the theory does not. The F s deviation reflects instability of sub-critical bubble ripening. Namely, in Pop3, CO 2 bubbles are spherical from the start and remain so during partitioning, after which their survival depends on local neighborhood configuration rather than interaction with the mean field. That said, F s and ⟨y 1 ⟩ are in slightly better agreement here than in the homogeneous network, since the variability in R p likely provides some stabilization.

The results for the cases in the fourth segment of Table 1 are very similar to those from the previous section. These cases are Cor-Pop1-yRand, Cor-Pop2-yRand, and Cor-Pop3-yRand, corresponding to the 1.

Network is heterogeneous and uncorrelated with random bubble placement. Snapshots are at t = 0, 1, 10 s, and equilibrium.

heterogeneous and spatially correlated network with random bubble placement, i.e., no correlation between y 1 and R p . Due to this similarity with the uncorrelated network, we only present the bulk properties in Fig. 15, which is nearly indistinguishable from Fig. 14. The only difference is that late-time ripening dynamics are slightly slower in the correlated network (see F s ). The likely explanation is that spatial correlations in R p play little role in early ripening, which is driven primarily by gradients in bubble composition. What matters is which bubbles (H 2 or CO 2 /mixture) are initially neighbors-or equivalently, how y 1 is arranged spatially. This arrangement is random both here and in the previous section, hence the small difference in results. At late times, however, when ripening is driven mainly by gradients in curvature, larger correlation entails smaller gradients in κ (since R p controls κ through Eq.9), thus slower mass transfer between bubbles.

We now turn to cases in the fifth segment of Table 1, namely, Unc-Pop1-ySmall, Unc-Pop2-ySmall, and Unc-Pop3-ySmall. These correspond to the heterogeneous uncorrelated network with H 2 bubbles placed in 1.

Network is heterogeneous and uncorrelated with random bubble placement. Snapshots are at t = 0, 1, 10 s, and equilibrium. illustrate the remaining two projections, where we see qualitatively different evolution patterns compared to prior sections when the initial arrangement of y 1 was spatially random. The agreement between theory and PNM is good, with errors larger in Pop3 due to reasons already discussed in Section 6.2.1.

Bulk properties are plotted in Fig. 18. The main difference compared to Fig. 14 concerns ⟨y 1 ⟩: while roughly constant in Fig. 14, it now decreases with time. Note that ⟨y 1 ⟩, which equals y mf 1 by Eq.30, is not the same as y mix

. By Corollary 1, y mix 1 is constant over time in the theory, and only as t → ∞ does y mf

The values of y mix 1 for each case are listed in Table 1. One can verify that the equilibrium values of ⟨y 1 ⟩ in Fig. 18 match the y mix 1 values in Table 1. The discrepancy in equilibrium ⟨y 1 ⟩ for Pop3 arises because some H 2 is absorbed by water in PNM-an effect not captured by the theory-resulting in lower ⟨y 1 ⟩ predicted by PNM. The deviation in F s reflects the instability of sub-critical ripening, as discussed in Section 6.2.1.

We now examine the final segment of Table 1: Cor-Pop1-ySmall, Cor-Pop2-ySmall, and Cor-Pop3-ySmall, which correspond to the heterogeneous correlated network with H 2 bubbles placed in small pores and CO 2 /mixture bubbles in large pores. These cases are the most challenging because both R p and y 1 are spatially correlated. As in the previous section, we omit the S b -y 1 projection of g(s; t) because it is very similar to Fig. 11. The remaining two projections, R p -S b and R p -y 1 , are shown in Figs. 19-20. We see satisfactory agreement between theory and PNM, but errors are significantly larger than in previous sections. This is especially apparent in the equilibrium configurations of Pop1 and Pop3, and in the theory overestimating the speed of ripening. Fig. 21 shows bulk properties versus time, with consistent observations.

The most likely explanation relates to two fundamental assumptions in the theory. First, bubbles at state s interact with a mean-field having composition y mf 1 and curvature κ mf s -equivalently (y 1 P b ) mf s . Only κ mf s accounts for spatial correlations through Eq.28, hence the subscript s. The mean-field composition y mf 1 , computed via Eq.30, is s-independent and informed by all bubbles in the population, whether proximate or 4 From Appendix E, γ → 1 as t → ∞. Therefore, both y mf, * 1. Network is heterogeneous and uncorrelated with H 2 bubbles placed in small pores. Snapshots are at t = 0, 1, 10 s, and equilibrium.

distant relative to the s-bubble. This causes overestimation of the phase velocities through Eq.18b, directly impacting dy 1,s /dt and indirectly impacting dS b s /dt through coupling terms, thus accelerating ripening. Second, spatial correlations accounted for in κ mf s are of pore size R p , not composition y 1 . This assumption was necessary to render the computation of mean-field properties tractable via approximations made between Eqs.22 and 28. The errors in Figs. 19-20 are the price we pay for these simplifications. The reason Pop2 predictions are more accurate is that chemical-potential gradients between H 2 and mixture bubbles are less than gradients between H 2 and CO 2 bubbles in Pop1 and Pop3, resulting in weaker disequilibrium.

This work presents the first theory for multicomponent ripening of bubbles in porous media, generalizing that of [37] for single-component ripening. It also simplifies the formulation and removes key limitations of 1. Network is heterogeneous and uncorrelated with H 2 bubbles placed in small pores. Snapshots are at t = 0, 1, 10 s, and equilibrium. that theory, as demonstrated in Section 6.1. While we focused on H 2 and CO 2 bubbles, other gases relevant to underground hydrogen storage (e.g., N 2 , CH 4 ) or otherwise could be used without loss of generality. Below, we discuss pending limitations of the theory to guide future research. They are divided into methodological, 1. Network is heterogeneous and correlated with H 2 bubbles placed in small pores. Snapshots are at t = 0, 1, 10 s, and equilibrium.

which exist within the assumptions of the theory and may be addressed by reformulating certain parts, and foundational, which stem from core assumptions that may require a completely new formulation.

Even under the assumptions of the theory, Section 6.2.5 showed that predicting ripening of bubble populations faces three remaining challenges: (1) capturing both kinetics and equilibrium state is difficult when composition is spatially correlated (Cor-ySmall ); (2) late-time predictions of bubble composition ⟨y 1 ⟩ are overestimated when the wetting phase is preferentially undersaturated in one component at t = 0 (Pop3 ); (3) early-time predictions of bubble volume ⟨V b ⟩ exhibit a temporary sag whose magnitude is underestimated.

Addressing the first limitation requires alternative approximations to those in Eqs.22-28 and Eq.30. The former postulates dominant control of pore-size correlations over other spatial correlations, while the latter removes state-dependence from y mf 1 and replaces the closure Eq.17 with Eq.18 to avoid numerical difficulties and ensure convergence to equilibrium. The analysis in Appendix B shows this reformulation is necessary 1. Network is heterogeneous and correlated with H 2 bubbles placed in small pores. Snapshots are at t = 0, 1, 10 s, and equilibrium. for Eulerian (grid-based) approaches used to discretize and solve the population balance Eq.10 (e.g., finite volumes herein). A Lagrangian (particle-based) method could obviate such reformulation and use the more accurate closure Eq.17 directly. This in turn would enable a state-dependent y mf 1,s that accounts for spatial correlations in composition. The second and third limitations arise because the wetting phase has no storage capacity for dissolved solute, rendering mass transfer between bubbles effectively instantaneous. This causes overprediction of ⟨y 1 ⟩ and the inability to capture the sag in ⟨V b ⟩. One remedy is to augment Eq.18 with an auxiliary balance equation for dissolved solute in the wetting phase. These remedies are left as future work.

While the theory was tested for two-component bubbles, extension to multiple components is straightforward: each new component adds another dimension to the phase space and an extra phase velocity to Eq.18.

Finally, contact angle and pore shape can be generalized beyond those assumed here (zero contact angle, semi-cubic pores), as they enter the theory solely through the curvature-volume relation in Eq.9.

Foundational limitations question the very choice of phase-space axes we have made for the theory, namely, R p , S b , and y 1 . The most significant is the assumption that each bubble occupies no more than one pore, which renders the use of pore size R p as a state coordinate inappropriate. However, a clear alternative is not self-evident because even if a bubble occupies only two pores, it begs the question: which pores?-since their sizes matter. Let the bubble occupy three or more pores and the question grows with combinatorial complexity. Another limitation of the same type is when bubbles, even if confined to one pore, occupy pores with different shape or wettability (i.e., contact angle). Put differently, there is heterogeneity in the curvature-volume relation Eq.9, which the theory assumes to be uniform across all pores. Introducing a separate axis for each shape or contact angle would make the theory intractable. A viable option might be to parameterize the curvature-volume relation with contact angle and shape-controlling constants and add those as extra phase-space axes. Either way, an increase in dimensionality of the phase space seems unavoidable.

Since the population balance Eq.10 must be discretized and solved numerically (here using the finitevolume method), readers might wonder if there are true gains in computational cost over solving the problem spatially using PNM. The answer is an emphatic yes for two reasons: (1) there is no limit to the number of bubbles that can be simulated in the theory, whereas the cost of PNM scales with the number of bubbles;

(2) when the correlation length l c of the porous medium is large, truly enormous pore networks are needed to ensure ripening predictions are representative of an REV, which is computationally prohibitive. The theory’s cost is independent of l c , enabling predictions out of reach with PNM or any other pore-scale model.

For the results in Section 6, the theory was faster by a factor of 192 for homogeneous networks and 24 for heterogeneous networks. These speedups were achieved despite using a very fine mesh for the phase space, which from experience we have found unnecessary, meaning even higher acceleration is possible.

The proposed theory has practical applications in underground hydrogen storage, where major concerns include: (1) inability to recover injected H 2 due to capillary trapping; (2) loss of purity of recovered H 2 ; and (3) maintaining high storage capacity and injectivity (i.e., low flow resistance) across multiple injectionwithdrawal cycles. The theory predicts which pore sizes are occupied by bubbles, their volumes, and compositions at all times. When coupled to a permeability model like Carman-Kozeny [39,40], where the porous medium is represented by a bundle of capillary tubes of radii R p , predictions of g(s; t) can be translated into estimates of relative permeability evolution. This would enable prediction of storage capacity and injectivity versus time, and could be used to maximize H 2 retrieval and minimize purity loss in cyclic injections.

The ySmall cases in Section 5 were motivated by the scenario where cushion gas (e.g., CO

We have presented the first kinetic theory for multicomponent Ostwald ripening of bubbles in porous media. The formulation describes a trapped bubble population with a number-density function g(s; t) in a three-dimensional statistical space of bubble states s. The coordinates of this space consist of the size of the occupied pore R p , occupied fraction of the pore volume S b (bubble saturation), and bubble composition y 1 .

The theory predicts the evolution of g(s; t) in time through a population balance equation. Closure is achieved through mean-field approximations that embed the physics of ripening, account for spatial correlations in pore size, and ensure convergence to a physically meaningful equilibrium with provable mass conservation.

The theory generalizes that of [37] for single-component ripening, removing key limitations such as the inability to capture interactions between distant bubbles. Systematic validation against pore-network simulations across 19 cases-spanning homogeneous, heterogeneous, correlated, and uncorrelated networks with various bubble populations-shows good agreement with PNM in both the evolution of g(s; t) and bulk properties of the population, without any adjustable parameters. The most challenging cases involve spatially correlated bubble composition, where errors arise due to methodological limitations of the theory discussed in Section 7.1.1. Unlike PNM, the theory allows for arbitrarily large bubble populations and arbitrary correlation lengths in pore size, enabling predictions beyond PNM’s reach. This makes the theory a practical tool for predicting Ostwald ripening in a variety of applications, including underground hydrogen storage.

At equilibrium, s * -bubbles will satisfy dS b s * /dt = 0 and dy 1,s * /dt = 0, with s * corresponding to the meanfield state s mf . This point, s mf , is the intersection of two curves: one defined by dS b s /dt = 0 and the other by dy 1,s /dt = 0-s here is a free index. Fig. B.1 visualizes this in the slightly altered phase space with coordinates (κ/κ mf , y 1 ), instead of (S b , y 1 ). We pause to note that the majority of bubbles must be super-critical (nonspherical) because otherwise they would not be in equilibrium (contradiction); ripening of spherical bubbles is unstable to perturbations [20]. This makes the coordinate mapping S b → κ/κ mf well-posed, otherwise the map would not be injective (i.e., spherical and deformed bubbles with different S b could map onto the same s /dt = 0 and dy 1,s /dt = 0 correspond to curves, use them to set the LHS of Eq.17 to zero:

Substituting y 2,s = 1 -y 1,s and P b s = P w + σκ s , and recognizing (y 1 P b ) mf ≈ y mf 1 (P b ) mf due to the population being in near equilibrium, the conclusion follows. Notice we have dropped the subscript s from (•) mf s because the mean field here does not depend on s. Now, we claim that regardless of where the s * -bubbles are located in the (κ/κ mf , y 1 ) plane, they will move towards s mf depicted by the yellow dot in s /dt > 0 in the region to its right. An identical argument proves that dy 1,s /dt > 0 must hold left of the black curve, and dy 1,s /dt < 0 to its right. However, this time the arbitrary point is chosen on the horizontal s /dt = 0 and dy 1,s /dt = 0 for Eq.17 in phase space for a homogeneous network. The phase space is 2D because pore size (Rp) is uniform, with equivalent coordinates (κ/κ mf , y 1 ) instead of (S b , y 1 ).

The two curves form an ever narrowing corridor that bubbles must traverse to reach equilibrium (yellow dot). At some point, the corridor’s width is smaller than the grid used to discretize the population balance Eq.10, halting further progress towards the yellow dot. (Right) The reformulation in Eq.18 removes this drawback by repositioning the black line to widen the corridor. line κ s ′ = κ mf , implying P b s ′ = P b,mf . Substituting this into Eq.17b yields:

where C 2 > 0, and the sign of each region is determined.

What we have shown above is that the vertical component of the phase velocity dS b s /dt, which is proportional to dκ s /dt by Eq.9 and the chain rule, points down in the region left of the red curve and points up to the right of it. Similarly, the horizontal component of the phase velocity dy 1,s /dt points right in the region left of the black curve and points left to the right of it. This means once s * -bubbles fall into the wedge between the red and black curves, they cannot escape it and the numerical problem above is inevitable.

The reformulation of Eq.17 into Eq.18 solves this problem by repositioning the black curve to the vertical line shown in the right panel of Fig. B.1. This broadens the wedge and prevents the gap between the red and black curves from becoming smaller than the grid size unless s * is very close to the yellow dot (equilibrium). Empirically, we observed the same problem in heterogeneous networks, which are far less tractable because the phase space is 3D. As discussed in Section 7.1.1, the reformulation of Eq.17, which comes with approximation errors, may be avoided if a Lagrangian (particle-based) method is used to solve the population balance Eq.10.

We validate the l RR ′ computed from Eq.29 against direct computations from the PNM in Section 3. We consider heterogeneous correlated and uncorrelated networks corresponding to the SC-Cor and SC-Unc cases in Section 5, respectively. Since the exact position of all bubbles are known in the PNM, we can compute l RR ′ directly using the following formula corresponding to the variable’s definition:

Namely, l RR ′ is the average nearest distance from bubbles residing in pores of size R p to bubbles in pores of size R ′ p . In Eq.C.1, x R and x R ′ denote the positions of these bubbles. The minimum is taken over all bubble-occupied R ′ p and the average ⟨•⟩ over all bubble-occupied R p . Incidentally, Eq.C.1 is also used to compute d RR ′ in Eq.29 with the crucial difference that the minimum and average are taken over all pores in the network regardless of bubble occupancy. We thus focus on the L ∞ norm of n α,s dS b s /dt over Γ ϵ in Eq.D.4. As ϵ → 0, it is evident S b s → 0. In this limit, bubbles are circular/spherical in shape whose curvature scales like κ s ∼ (S b s ) -1/d . We use the symbol a ∼ b to mean: there exist constants C for the term before the summation in Eq.D.5. For the summation itself, we can write:

2 and CO 2 bubbles in Pop1 have local saturation S b = 60%. Pop2 is the same, but instead of pure CO 2 bubbles, we have equimolar H 2 -CO 2 bubbles. In Pop3, H 2 bubbles have S b = 60%, while CO 2 bubbles have S b = 30%. In all three cases, as bubbles evolve, g(s; t) approaches a horizontal line of constant y 1

2

We have simplified this and following expressions by assuming uniform throat size across the network. See Eq.17 in[37].

One consequence of u Rp = 0 is that Eq.10 reduces to a series of 2D PDEs in the variables y 1 and S b for each fixed Rp.

Because dκ s ′ /dS b s ′ > 0 in Eq.17a on account of bubbles being super-critical.