A level set approach to Ostwald ripening of trapped gas bubbles in porous media

Ostwald ripening of gas bubbles is a thermodynamic process for mass transfer, which is important for both foam enhanced oil recovery and geological CO2 storage. We present a methodology for simulating Ostwald ripening of gas ganglia surrounded by liquid in arbitrary pore geometries. The method couples a conservative level set model for capillary-controlled displacement and a ghost-bubble technique that calculates mass transfer based on difference in chemical potentials. The methodology is implemented in a software framework for parallel computations. As a validation of the model, we show that simulations of bubble ripening in a pore throat connecting two pore bodies are consistent with previously reported trends in similar geometries. Then we investigate the impact of gas type, compressibility factor, and local capillary pressure on gas-bubble ripening in various water-wet pore geometries. The results confirm that gas solubility and compressibility factor are proportional to the rate of mass transfer. Our simulations suggest that Ostwald ripening has largest impact in heterogeneous or fractured porous structures where differences in gas-bubble potentials are high. However, if the liquid separating the gas bubbles is also a disconnected phase, which can happen in intermediate-wet porous media, the resulting local capillary pressure can limit the coarsening and stabilise smaller bubbles. Finally, we simulated Ostwald ripening on a 3-D pore-space image of sandstone containing a residual gas/water configuration after imbibition. Characterization of gas-bubble morphology during the coarsening shows that large ganglia get more ramified at the expense of small spherical ganglia that cease to exist. A level set approach to Ostwald ripening that handles arbitrary pore geometries and gas-ganglia configurations. The role of a disconnected liquid phase in Ostwald ripening is to stabilize more bubbles and prolong bubble lifetime. During coarsening in sandstone spherical bubbles dissolve while large ganglia grow and develop a more ramified structure. A level set approach to Ostwald ripening that handles arbitrary pore geometries and gas-ganglia configurations. The role of a disconnected liquid phase in Ostwald ripening is to stabilize more bubbles and prolong bubble lifetime. During coarsening in sandstone spherical bubbles dissolve while large ganglia grow and develop a more ramified structure.

may also render the surrounding liquid a discontinuous phase which will affect the bubble pressures and bubble interactions.
Recent studies have explored implications of Ostwald ripening on the permanency of CO 2 ganglia in residual (or, capillary) trapping as a mechanism for CO 2 storage Li et al. 2020Li et al. , 2022. Such disconnected ganglia can take many irregular shapes and span multiple pores. If Ostwald ripening causes large ganglia to grow further, they can get more susceptible to capillary instabilities (e.g., coalescence and splitting) and remobilization by viscous pressure differences or gravitational/buoyant effects, which ultimately may lead to less CO 2 stored by residual trapping. Capillary pressures for disconnected gas ganglia, obtained from calculation of interface curvatures on segmented microtomography images of two-phase gas/water configurations in porous rocks, generally show a wide distribution Herring et al. 2017;Li et al. 2018), providing prime condition for Ostwald ripening. At larger scales and in sandstones with significant laminations and sorting difference, the capillary pressure distribution could be even broader, and thus the potential for remobilisation by Ostwald ripening could be even more significant . Li et al. (2020) proposed a mathematical model to compare mass transfer through diffusion and mass transfer through Ostwald ripening. They concluded that even though gas redistribution by Ostwald ripening is a prolonged process at reservoir scale, at millimetre scales, the phenomena can occur at a timescale of days to years.
Only a few pore-scale models have been developed to explore Ostwald ripening of twophase fluid systems in porous media. de  and Mehmani and Xu (2022) presented pore-network models that approximate the pore space as a network of interconnected pores with idealized geometry. de  constructed an algorithm based on identified diffusion paths between interacting bubbles, while Mehmani and Xu (2022) solved a set of mass balance equations for every pore and validated their model with micromodel experiments (Xu et al. 2017). Apart from the idealization of pore geometry, these studies make a number of assumptions: First, each gas bubble occupies only one pore. While this assumption allows calculation of pressure-volume relations for the gas bubbles in the pertinent pore geometry as input to the simulations, it excludes possible capillary instabilities arising from further bubble growth. Second, the liquid phase surrounding the gas bubbles is assumed to be continuous and assigned a fixed, uniform pressure, and hence all interfaces of a gas bubble will have the same capillary pressure. This assumption eradicates pressure interactions with neighbouring phases from the calculations of gas bubble pressures. Such interactions are important when both the gas and liquid form disconnected ganglia with different pressures: Mass transfer between two gas ganglia is an outcome of their pressure difference, which in turn depend on the pressures of surrounding liquid ganglia.
In this work, we propose a new methodology to study Ostwald ripening which couples a previously developed multiphase level set (MLS) model for direct simulation of quasi-static, capillary-controlled displacement on arbitrary pore geometries (Helland et al. 2019;Jettestuen et al. 2021) with a "ghost-bubble" (GB) method that describes mass transfer between gas bubbles based on differences in chemical potential. We refer to the coupled approach as the GB-MLS model. Advantages with level set methods are their natural ability to handle interface redistribution and topological changes due to capillary instabilities. Hence, we do not impose any constraints on either the pore geometry or the gas/water configuration and allow scenarios where gas ganglia occupy multiple pores. Further, unlike standard level set methods, the MLS model has the option to conserve volumes of disconnected ganglia of either one or both phases and predict their pressures (Jettestuen et al. 2021). As such, pressure-volume relations for gas bubbles are calculated naturally during interface evolution, accounting for scenarios where the presence of surrounding liquid ganglia with different pressures impacts the gas bubble pressures. Finally, we have implemented three equation of states (EoS) to assess the effect of gas type and gas compressibility factor on coarsening in porous media at reservoir conditions. Based on gas bubble pressures, the model computes the gas compressibility factor of each isolated gas bubble during evolution.
The paper is organized as follows: Sect. 2 provides a brief description of fundamentals required to study Ostwald ripening, and Sect. 3 describes the GB-MLS method and algorithm. Section 4 applies the GB-MLS model to different porous geometries and fluid configurations. First, we show that the model captures the expected behavior for bubble coarsening in bulk liquids in the absence of porous media. Second, we show that ripening simulations in a pore throat connecting two pore bodies occupied by gas bubbles are consistent with the results de  obtained in a similar geometry. Third, we use our model to investigate the impact of EoS, gas type, and gas compressibility factor, on simulations of ripening in a 2D micromodel with a pore-scale heterogeneity. Then, we perform simulations on a gas-bubble configuration in a sinusoidal pore channel to demonstrate how coarsening behavior changes when the surrounding liquid is modelled as a disconnected phase rather than as a connected phase. Finally, we apply the GB-MLS model directly on a segmented microtomography image of water-wet sandstone to explore the significance of Ostwald ripening on a residual gas/water configuration after imbibition. Section 5 gives the conclusions from our work.

Theory
In bulk gas-liquid systems, Ostwald ripening is a diffusive interaction mechanism that leads to the growth of larger clusters at the expense of smaller clusters based on their pressure difference. The rate of the diffusive mass transfer depends on the properties of both gas and liquid. For example, the diffusive transfer rate is higher for a lighter gas molecule than for a heavier gas molecule (as follows from Graham's law). Alternatively, higher solubility can lead to faster mass transfer by providing a higher concentration gradient for a similar capillary pressure difference. The (partial) pressure in a gas bubble, P g [Pa], affects the mass fraction X of gas species dissolved in the continuous phase next to the bubble, according to Henry's law: where H [Pa m 3 mol −1 ] is Henry's constant.
Capillary pressure P c , which is the pressure difference between gas (g) and liquid (l), is related to interface curvature by Young-Laplace's equation: Here, P g and P l are the phase pressures, gl is the gas/liquid interfacial tension, r gl,1 and r gl,2 are the principal radii of curvature, and gl is the interface curvature. Thus, for known P l and P c we can calculate gas bubble pressure P g .
Instead of being regarded as a direct bubble-to-bubble interaction, the mass transfer process can be considered as a two-step process (Lemlich 1978): First, gas diffuses from the bubbles into the surrounding liquid, followed by the second step of mass transfer from the (2) P c = P g − P l = gl 1 r gl,1 + 1 r gl,2 = gl gl .
liquid to the gas. The mass transfer rate is then derived from a general rate equation (Lemlich 1978): where n is the number of moles of gas in a bubble, t (s) is time, J (mol s kg −1 m −1 ) is the effective permeability for the mass transfer, A (m 2 ) is the gas/liquid interfacial area through which the mass transfer occurs, and ΔP (Pa) is the pressure difference between the gas bubble and a fictitious gas bubble (or, ghost bubble) representing the surrounding liquid region. The ghost bubble pressure is related to gas concentration in liquid through Henry's law. The pressure difference ΔP can be derived from the Young-Laplace equation and is given as where is the ghost bubble radius and r is the gas bubble radius. In literature, Ostwald ripening for gas bubble systems is generally studied as a process driven by the bubble pressure difference (Lemlich 1978;de Chalendar et al. 2017;Xu et al. 2017;Li et al. 2020;Mehmani and Xu 2022). It is described as an effort by the system to decrease its internal energy (Voorhees 1992). In the absence of elastic stress, the total interface area must be reduced spontaneously to reach thermodynamic equilibrium. At constant pressure P and temperature T, spontaneous changes to reach thermodynamic equilibrium reduce Gibbs free energy (Castellan 1983). The change in Gibbs free energy, G, per unit mole of substance i, n i , defines chemical potential, i : Chemical potential is an intensive property of a system. In case of potential difference between two gas bubbles, the spontaneous process will lead to the transfer of matter from higher potential to lower potential (Castellan 1983). The chemical potential in a system containing only one chemical species is given by the molar Gibbs free energy, For real gases, the chemical potential is calculated as Castellan (1983) where μ std is the chemical potential of the pure gas at standard pressure (typically P std ∼ 101, 325 Pa), R (8.314 J mol −1 K −1 ) is the universal gas constant, and f is the fugacity. The methodology for calculating fugacity using Soave-Redlich-Kwong (SRK) and Van der Waals' (VdW) equation of state is described in Supplementary information. In the case of ideal gases, fugacity is replaced by pressure P (Castellan 1983): (6) = G n .
(7) = std + RT ln f , de Vries (1958) carried out experimental investigations on foam to study the effect of gas diffusion on foam stability. Lemlich (1978) used this data to formulate an Ostwald ripening model using the concept of an imaginary ghost bubble to represent the liquid in its interaction with the gas bubbles. Here, we adopt a similar approach that we refer to as the "ghost bubble" (GB) method. We describe the GB method for a real gas; however, the method can be used for an ideal gas by replacing fugacity with pressure. The GB approach quantifies the relative net mass transfer between gas bubbles at different chemical potentials toward a state of thermodynamic equilibrium. The purpose of a ghost bubble is to provide a relative measure of net mass transfer from a bubble with high chemical potential to one with lower chemical potential through a liquid region. The ghost bubble (denoted by subscript 0) is defined so that its pressure P 0 and chemical potential μ 0 represent the stable state for the gas bubbles connected to the liquid region at a particular time step. The GB method makes the following assumptions: (i) The interface readjustment is instantaneous in comparison to the diffusive transfer. This is similar to the assumption made by de Chalendar et al. (2017). (ii) If a liquid droplet (or region) contacts only one gas bubble, they are in mutual equilibrium with equal chemical potentials. In case of any changes in the chemical potentials of neighbouring ganglia, the chemical potential of the liquid droplet changes spontaneously to adapt to the new chemical potential equilibrium (Castellan 1983). (iii) The mass transfer occurs due to chemical potential difference across the liquid region boundary. Here, we consider mass transfer only through bulk liquid regions. (iv) Two gas bubbles interact with each other only if a path can be found between them that does not include another gas bubble ). (v) There is no mass accumulation or loss (dissolution) from gas bubbles to the liquid during Ostwald ripening, instead the liquid region acts as a path for mass transfer between two or more bubbles. (vi) Reservoir temperature is above the critical temperature of the gas. This assumption ensures that the chemical species is in one phase, i.e., gas or supercritical fluid, and prevents phase change calculations.

Ghost bubble method
Regarding assumption (v), experimental studies on CO 2 solubility in water over temperatures ranging from 4.5 to 175 • C and pressures ranging from 1.5 to 18 MPa suggest that the binary fluid system reaches solubility equilibrium in a timescale of orders of hours (Ma et al. 2017;Hou et al. 2013;Prutton and Savage 1945;Wiebe and Gaddy 1939). Similarly, studies with CH 4 and water also reported a solubility equilibrium timescale of hours to weeks (Serra et al. 2006;Blount and Price 1982). Generally, the contact area to volume ratio at the pore scale is higher than in the aforementioned experiments. So, a smaller time interval should be expected at the pore scale in the subsurface environment for the water surrounded by gas bubbles to get saturated with the gas phase and for assumption (v) to be valid.
We begin by writing the general rate equation of Lemlich (1978) (Eq. (3)) in terms of chemical potential difference between a gas bubble and ghost bubble, Δ = i − 0 . The molar flow rate from neighbouring gas bubbles to bubble i through a liquid region l is where k l is the proportionality parameter, J ′ the effective phase permeability, and A is the gas/liquid interfacial area through which mass transfer takes place. The negative sign in Eq. (9) indicates mass transfer direction from higher to lower potential. Now, in case of only two bubbles, this mass transfer rate is equal to the mass loss rate from one bubble j and mass gain rate by the other bubble i: As there is no mass accumulation in the liquid, n i,l = −n j,l is the mass transferred from bubble j to bubble i through the liquid region l. Using Eq. (7), the potential difference between gas bubble i and the ghost bubble is Combining Eqs. (10) and (11) gives For an intermediate liquid region in direct contact with N gas bubbles, the total mass transfer at any time instance t ′ is zero. Hence, The above equation leads to the following expressions for the ghost bubble parameters (see derivations in "Appendix 1"). Ghost-bubble fugacity (in case of real gas) is whereas ghost-bubble pressure (in case of ideal gas) is In the case of mass accumulation in the liquid, the right-hand side of Eq. (13) will not be zero. Instead, it will equal the average rate at which the gas concentration in the liquid increases due to gas dissolution from bubbles surrounding it. The ghost bubble fugacity (or pressure) predicted from Eq. (14) (or Eq. (15)) will need to be reduced to accommodate the additional dissolution of gas in the liquid phase. The effective phase permeability, J ′ , sets the timescale of our simulations, and it can be determined from experiments. We have derived a value through dimensional analysis for our model. Micromodel experiments from Xu et al. (2017) suggest that the timescale of ripening at the pore scale is in hours. The effective phase permeability coefficient used in this work is (see dimensional analysis in "Appendix 2"): where s (mol m −3 ) is the solubility of the gas in the liquid phase, H (Pa m 3 mol −1 ) is Henry's constant, and M (kg mol −1 ) is the molar mass. For real gases, Eqs. (12), (14), and (16) describe the mass transfer at each time step. In the case of ideal gas, we replace fugacity with pressure in Eq. (12) and use Eq. (15) instead of Eq. (14). The change in gas/water configuration caused by the mass transfer, including phase pressures of disconnected ganglia, is calculated using the MLS model (Helland et al. 2019;Jettestuen et al. 2021), which we describe in the next section.

Multiphase level set method
The level set method is a subgrid-scale interface tracking method in which the zero contour of a signed distance function describes the fluid/fluid interfaces implicitly (Osher and Fedkiw 2003). The method has previously been used to simulate curvature-driven flows like two-phase capillary-controlled displacement in porous media at different wetting states (Prodanović and Bryant 2006;Jettestuen et al. 2013). The multiphase level set (MLS) method (Helland et al. 2019;Jettestuen et al. 2021) extends this approach further to handle interactions in porous media between two or more fluids, using one level set function for each fluid phase . Another static level set function describes the pore/solid geometry, so that = 0 represents the pore walls. Then, the pore space occupied by the fluids (in our case, a two-phase gas/liquid system) is defined as: where Ω is the computational domain. While it is sufficient to use a single level set function to describe interfaces in two-phase systems, imposing local volume conservation on isolated ganglia of both phases still requires a phase-by-phase approach (Jettestuen et al. 2021). Therefore, in this work we proceed with the MLS method for Ostwald ripening investigations in the gas/liquid system. This also allows for further development to investigate ripening in the presence of a third phase which is a planned future work.
In the MLS method, each fluid phase has its surface tension , phase pressure P , and solid/fluid interaction angle (see Fig. 1). The surface tensions are defined from the interfacial tension as The effect of rock wettability is handled by the solid/fluid interaction angles and the surface tensions , which Young's equation relates to the contact angle gl (measured through the liquid phase) as follows Helland et al. (2019): Here, �� ⃗ l = ∇ l ∕|∇ l | and ��⃗ s = ∇ ∕|∇ | are unit normal vectors of the liquid and pore/ solid level sets, respectively, see Fig. 1.
The evolution equations for the MLS method are derived by energy minimization and then parameterising the solutions in a fictitious iteration time (Helland et al. 2019): In Eq. (20), a sharp Heaviside function H separates the velocity for capillary-controlled displacement in the pore space from the velocity which forms the contact angle on the pore walls and in solid. Further, = ∇ ⋅ ∇ ∕|∇ | is a scalar curvature field for the level sets , Δx is grid spacing, and S( ) is a sign function given by Note that the level set Eq. (20) for the gas and liquid phases are uncoupled. Hence, during evolution, the zero contours of g and l may develop overlap or void regions. We solve this problem by performing the projection step from Losasso et al. (2006) at the end of each iteration step Δ . In our two-phase system, the projection step simply updates the level set functions everywhere as ̄g = ( g − l )∕2 and ̄l = −̄g . This sews the level sets together (see Fig. 1) and ensures that the steady-state solution for capillary equilibrium satisfies Young-Laplace's Eq. (2) (where gl = g = − l ) and Young's Eq. (19) for the contact angle (Helland et al. 2019).
A drawback with level set methods is their need for explicit conservation, which is especially important when dealing with fluid flow in porous media. Jettestuen et al. (2021) included volume conservation of isolated fluid ganglia (accounting for splitting and merging) in the MLS method by updating their pressures in Eq. (20) to prevent volume changes. The pressure of an isolated ganglion i of phase at iteration time , subjected to local volume conservation, is calculated as Saye and Sethian (2011): Here, V ,i (0) is the initial or original ganglion volume, while V ( ) ,i and A ( ) ,i are the calculated ganglion volume and surface area at iteration time . The method has the option to enforce (20) Fig. 1 Illustration of the MLS method for a fluid/fluid interface in a pore geometry, before (left image) and after (right image) executing the projection step of Losasso et al. (2006) 1 3 conservation on either one or both phases. It can also distinguish between connected and disconnected parts of the phases by specifying a pressure to the connected part (which contacts inlet or outlet boundaries of the computational domain and is not conserved), while Eq. (22) calculates the pressures in the conserved disconnected parts. The volume and surface area of the fluid phases are and respectively. Here, H is a smoothed Heaviside function, is a smoothed delta function, and = 1.5 × Δx is the smoothening parameter. In two-phase systems, the interfacial area between gas and liquid is given by Eq. (24), or alternatively, by A gl = (A g + A l )∕2.
The method uses non-dimensional parameters for pressure, surface tension and length (P, , and L), which are derived from the input parameters ( P actual , actual , and L actual ), and a set of characteristic parameters ( P * , * , and L * ) for the physical problem at hand. Thus, During evolution, the level set functions get distorted and require repeated reinitialisation to maintain their signed distance nature. The reinitialisation is carried out using the following equation (Osher and Fedkiw 2003): where S( ) is given by Eq. (21).

Combined GB-MLS method: Algorithm
This section describes the complete algorithm of the GB-MLS method for investigating mass transfer and evolution of gas ganglia driven by Ostwald ripening. Figure 2 shows a flowchart of the complete process. Following de Chalendar et al. (2017), the method assumes the interface adjustment is instantaneous, and hence only the mass transfer contributes to the time development. The different steps in the GB-MLS method are as follows (with reference to Fig. 2): Step 1: The method takes the following as input: (i) Temperature T and pressure P that describe the reservoir conditions, and gas properties required to solve the EoS (molar mass M, acentric factor , critical pressure P crit and critical temperature T crit ); (ii) interfacial tension gl and contact angle gl from which we determine and , = g, l , using Eqs. (18) and (19); and (iii) binary image data of the pore geometry and fluid locations.
Step 2: We calculate the signed distance functions , g , and l , using Eq. (26). In addition, pressures P g and P l are assigned for any connected parts of the phases. Then, following Jettestuen et al. (2021), we identify any isolated ganglia of the conserved phase, calculate their initial pressure, volume, and surface area, using Eqs. (22), (23), and (24).
Step 3: Use Eq. (20) to evolve the level set functions forward one iteration-time step Δ and enforce the projection step of Losasso et al. (2006). Periodically reinitialise to

Fig. 2 GB-MLS methodology
for determining gas/liquid interface evolution due to Ostwald ripening. The workflow shows GB method steps (green), MLS method steps (orange), decision boxes (red), and numbers in circles that provide a reference to each step signed distance functions using Eq. (26). In Eq. (20), we approximate normal and advective velocity terms using a weighted essentially non-oscillatory (WENO) scheme with appropriate upwinding techniques (Osher and Fedkiw 2003), while we approximate area and curvature terms with central differences. The iteration-time discretisation uses a thirdorder Runge-Kutta method with time step determined from a standard Courant-Friedreichs-Levy (CFL) condition (Osher and Fedkiw 2003).
Step 4: Redistribute volumes of isolated domains and allocate pressures to them (see Jettestuen et al. (2021)).
Step D1: Check if the fluid system has reached an equilibrium state based on convergence criteria for (Jettestuen et al. 2021): where we use c = 0.001 and = 1.5x , and n and m are the two last iteration steps where reinitialisation of level set functions occurred.
Step 5: Based on Eqs. (22), (23), and (24), provide the pressure, volume, and interfacial area, respectively, of all isolated gas and liquid domains for the capillary-equilibrium state achieved with the MLS method. The GB method exploits these data further.
Step 6: For each liquid region (which we represent by a ghost bubble region), generate a list of all neighbouring gas ganglia. The neighbourhood identification finds the location of interfaces shared by the different gas ganglia and the ghost bubble region, from which we calculate the area of each interface separately, which is required by Eq. (14) (or, Eq. (15)) (see also "Appendix 1").
Step 7: For real gases, calculate the fugacities (and compressibility factors) of all gas ganglia based on their pressures, using VdW or SRK EoS (see Supplementary information). Then, use these fugacities together with the neighbourhood lists to calculate the fugacity of ghost-bubble regions with Eq. (14). In case of ideal gas, we calculate the pressures of ghost-bubble regions with Eq. (15) using the pressures of neighbouring gas ganglia.
Step 8: Calculate the mass transfer Δn i,l between each gas bubble i and the surrounding liquid region l for a particular timestep Δt by converting Eq. (12) to a finite difference equation: In case of ideal gas, we replace fugacity with pressure in Eq. (28). This net mass transfer can be positive or negative to indicate mass loss or mass gain for a particular gas bubble through the liquid. If a particular gas bubble i is in contact with N l different liquid ganglia, we calculate the mass evolution of the bubble by where n i is the mass of the bubble, superscript b and a represent the next and current iterative time, respectively, and ∑ N l l=1 Δn i,l is the mass transfer from the gas bubble in timestep Δt . We calculate Δt by setting a threshold for the maximum gas bubble volume change per time step equal to 2 to 10 grid cells for the 2D geometries and 10 grid cells for the 3D geometry.
Step 9: Update the original volume of each isolated gas bubble, V (0) g,i in Eq. (22), before progressing the multiphase system toward the next capillary equilibrium state with the MLS method. As chemical potential of a gas is pressure-dependent, we assume in the numerical scheme that the pressure is constant for small changes in volume.
Step D2: After each timestep Δt of the GB method, check if the Ostwald ripening process has reached a steady state based on the difference between highest and lowest dimensionless gas-bubble pressures in the system: where is an error threshold. In this work, we use = 0.001 for 2D geometries and = 0.05 for the 3D geometry.
As for the MLS method (Helland et al. 2019), we have implemented the combined GB-MLS code within SAMRAI, which is a software framework that enables parallel simulations using patch-based data structures (Hornung and Kohn 2002;Hornung et al. 2006;Anderson et al. 2013).
A challenge with the GB-MLS method is that the bubble distribution evolution can have a cyclic mass transfer during simulation. When a discrete amount of mass transfer from one bubble to another reverses the direction of mass transfer for the next step, we refer to it as "cyclic mass transfer". It is also possible that more than two bubbles participate in cyclic mass transfer that leads to gas bubble pressure fluctuations in the vicinity of the stable state. Error tolerance values, computational errors, and sharp changes in pore geometry or interface curvatures can cause cyclic mass transfers. Since our method is developed to work on any pore geometry, we cannot have an a priori method to control cyclic mass transfer. So in our model, we characterize a particular bubble distribution by the minimum, maximum and average bubble pressures, and the number of bubbles. If we find a similar bubble distribution among the last 10 GB-MLS iterations, we consider the mass transfer to be cyclic, and terminate the simulation. Cyclic mass transfer is linked to the maximum volume-change threshold set to determine Δt , and a very high threshold can significantly impact the results. Conversely, a higher in Eq. (30) lowers the possibility of cyclic mass transfer at the cost of accuracy.

Results and Discussion
We apply the GB-MLS method to investigate the impact of Ostwald ripening on gas-bubble configurations in different porous medium geometries for different physical conditions. After investigating the behavior through several 2D numerical examples, we apply the method to evaluate the significance of Ostwald ripening on a residual gas/water configuration in a 3D segmented image of sandstone.

Bubble interactions in bulk liquid
We begin by simulating the interaction of gas bubbles in the absence of porous media. For this purpose, we consider an initial configuration of five CO 2 bubbles with radii 7, 8, 9, 10, and 12 μ m separated by CO 2 -saturated water in a 2D spatial domain of size 1.2L * ×1.1L * m 2 , see Fig. 3. The characteristic parameters are set to L * = 10 −4 m and * = 10 −1 N m −1 , which gives P * = 1 kPa. We assume reservoir conditions with P = 15 MPa and T = 323.15 K, for which CO 2 is a supercritical fluid, the CO 2 /water interfacial tension is gl = 32.6 × 10 −3 N m −1 , and Henry's constant for CO 2 is H = 7.57 × 10 −5 mol m − 3 Pa −1 ). In the MLS method, grid spacing is Δx = 0.01 . Further, we use reflective boundary conditions and reinitialise the level set functions , = g, l , to signed distance functions every 10 iterations. In the ghost bubble method, we use SRK EoS to solve for CO 2 fugacity and compressibility factor at every timestep. Figure 3 depicts the bubble evolution toward equilibrium. As expected, growth of larger bubbles occurs at the expense of smaller bubbles until only the largest bubble (with lowest capillary pressure) exists in the system. Note also that the location of the bubble centres remain stationary during the growth.
The evolution of bubble radii over time, as shown in Fig. 4, highlights some important points. First, the bubbles disappear in the order of their sizes (Clark and Blackman 1948). Second, mass transfer rate increases with bubble-size differences. We observe this from the radius evolution when the system consists of only two bubbles: The radius of the smaller bubble decreases at an increasing rate, leading to a convex trend in the curve for the shrinking bubble. Conversely, the rate of coarsening decreases when the bubble sizes come closer, which has also been observed for bubble-size distributions in foam (Lemlich 1978). Third, the second largest bubble does not decrease in size right from the start. Instead, it grows as the mass gain from smaller bubbles is larger than the mass loss to the bigger bubble. With time this imbalance turns and the size of the bubble starts decreasing. This instantaneous nature of growth and reduction is similar to those reported by Lifshitz and Slyozov (1961) for grain growth in solid solutions. This reversal shows that the bubble cluster responds to its immediate environment without preexisting knowledge of the future stable state.

Bubble-pair interactions in a 2D pore-throat/pore-body geometry
We proceed by investigating different scenarios for Ostwald ripening that can occur when two gas bubbles confined in a pore geometry are in capillary equilibrium with different pressures (or, interface radii) initially. For this purpose, we will simulate the scenarios that de Chalendar et al. (2017) explored with a typical pore-network model geometry composed Fig. 3 Evolution of five CO 2 bubbles in bulk water. The bubbles are coloured orange, cyan, red, green, and black in order of increasing size, and water is shown in dark-blue colour of a 3D converging/diverging pore throat (filled with water) connected by two pore bodies containing CO 2 bubbles. They identified four scenarios: Mass transfer from small to large bubble (case 1); mass transfer from large to small bubble (case 2); mass transfer leading to increased interface radii of both bubbles (case 3); and mass transfer leading to decreased interface radii of both bubbles (case 4). We explore these cases with the GB-MLS method in a similar 2D pore geometry, as shown in Fig. 5. The size of the spatial domain is 2.1L * ×0.9L * m 2 . We use brine/CO 2 contact angle of gl = 7.28 • , while all other input and simulation parameters, as well as EoS, are the same as in the previous subsection. Figure 6 shows the initial and final equilibrium states from GB-MLS simulations of the four cases, as well as the corresponding interface radius evolution of the two gas bubbles. The results are consistent with those reported by de . In case 1 the bubbles interact through bulk liquid without contacting the pore walls, and hence they exhibit the same behavior as described in the previous subsection. Case 2 describes a scenario where the pore geometry confines both bubbles, so that the large bubble has lower interface radius (or, higher pressure) than the small bubble. Thus, Ostwald ripening causes the larger bubble to retract from the pore throat while the smaller bubble ingresses until their interface radii are equal, leading to bubble "anti-coarsening" Xu et al. 2017).
In case 3 and 4, the pore geometry confines the large bubble while water completely surrounds the small bubble. Case 3 is another example of "anti-coarsening" as the small bubble has highest interface radius so that mass transfers from the large bubble to the small bubble. This causes the large bubble to retract from the pore throat and the small bubble to expand freely, accompanied by a substantial increase of the large-bubble interface radius and a modest increase of the small-bubble interface radius. In case 4, the situation is reversed as the small bubble completely surrounded by water has lowest interface radius. This leads to a mass transfer from the small bubble, which forces the large bubble deeper into the pore throat. The interface radii of both bubbles decrease with time until they become equal. Note that the mass transfer is least in case 4 since the initial difference between the bubble radii is minimal. Cases 2, 3, and 4, all emphasize the potential impact of porous media on gas-bubble ripening, as both bubbles co-exist at the stable state with equal interface radii. In all cases, the rate of evolution reduces with time as the bubble radii come closer. This leads to a convex trend in the radius-evolution curves from the GB-MLS model.

Bubble interactions in a 2D micromodel with permeable channels
We construct a 2D micromodel with two high-permeable channels to investigate the impact of pore-space heterogeneity, gas type and EoS, on gas-bubble interactions at reservoir conditions in porous media. The size of the micromodel is 3.84L * ×2.80L * m 2 and it consists of square grains separated by pore space, see Fig. 7. The width of the regular pore channels is 8 μ m, while the width of the two vertical, high-permeable channels in the middle of the micromodel is 16 μ m. The micromodel reflects a conducive environment for gas trapping. We place 15 gas bubbles with volumes (area in 2D) ranging from 501.88 to 971.22 μm 2 inside the regular pore space and one large bubble with volume 1530.65 μm 2 in the middle of the two high-permeable channels, see Fig. 7. Brine saturated with gas species occupies the surrounding pore space. In the GB-MLS method, we use the same input and simulation parameters as before, except that we now consider a slightly different reservoir condition ( P = 15 MPa and T = 333.15 K) and investigate the impact of gas type (CO 2 , CH 4 , and N 2 ) and EoS. Interfacial tensions for these fluid systems are as follows: gl = 32.6 × 10 −3 N m −1 (brine/CO 2 ) (de Chalendar et al. 2017), gl = 60 × 10 −3 N m −1 (brine/CH 4 ) (Shariat 2014), and gl = 60.636 × 10 −3 N m −1 (brine/N 2 ) (Chow et al. 2016). These data were measured at slightly different temperatures where the variations of interfacial tensions are small and negligible (Chiquet et al. 2007;Shariat 2014;Chow et al. 2016). Henry's constants for CO 2 , CH 4 , and N 2 at T = 333.15 K were calculated as H = 1.1417 × 10 −4 , 4.048 × 10 −6 , and 7.9667 × 10 −6 mol m -3 Pa −1 , respectively, using relations from Sander (2015). We set the contact angle to gl = 7.28 • in all cases.
The large bubble present in the high-permeable channels will grow during coarsening of this system as it has the lowest gas-bubble pressure (chemical potential) among all bubbles. Therefore, we end the simulations when the large bubble spans the entire channel length. Figure 7 depicts the gas-bubble evolution for the brine/CO 2 system, using SRK EoS: Image (i) shows the initial capillary equilibrium state for the trapped gas Illustration of the pore geometry used to explore interaction between two bubbles confined by pore walls (pore space in cyan, and solid space in orange). Circles of equal radius represent the pore bodies where gas bubbles are present, while four equal truncated triangles form a symmetrical pore throat where water is present 1 3 bubbles, while images (ii)-(vi) show the growth of the large bubble and the corresponding depletion of the other bubbles, leaving 13 tiny bubbles left in the pore space (image (vi)). The bubble interactions occur similarly for the other fluid systems and EoS, however, the mass transfer rates are different. Fig. 6 Interaction of two gas bubbles (green) surrounded by brine (blue) in a pore geometry, as simulated with GB-MLS model: a Initial configuration, b evolution of normalized radius with respect to normalized time for bubble in left (red) and right (blue) pore body, and c final configuration Figure 8a shows volume evolution of the large bubble over time for the different gases and EoS. The evolution is nearly linear and misses the typical convex shape observed in the other numerical examples. This is because the pressure difference between the large bubble and the smaller bubbles in this micromodel geometry remains nearly similar throughout the simulation. Capillary pressure fluctuations due to interface movement between narrow pore throats and wider pore regions occur for different bubbles at different timesteps. Yet, the cumulative effect of these variations on the mass transfer rate is small, as seen by the nearly linear growth of large bubble volume in Fig. 8.
As expected, the ideal gas law overpredicts the rate of coarsening for CH 4 and CO 2 , whereas SRK and VdW EoS predict longer coarsening times. The different predictions of coarsening rate can be attributed to different gas compressibility factors Z calculated with the three EoS. For the initial gas-bubble distribution (see Table 1) the range of calculated Z deviates significantly from unity for CO 2 because system temperature and pressure are close to the critical point of CO 2 . When comparing Table 1 with Fig. 8a, we observe that the EoS which yields lower compressibility factor predicts a longer coarsening time for the large bubble. This is intuitive since gas compressibility factor is a measure of deviation from ideal gas behaviour: When Z < 1 , attractive forces dominate, which corresponds to lower molar volume and longer time to transfer a given volume. On the other hand, when Z > 1 , the calculated rate of volume change is higher for real gas than ideal gas (as observed for N 2 with SRK EoS). Figure 8a also shows that the growth rate is highest for CO 2 and lowest for N 2 . Zeng et al. (2016) made similar observations while experimenting with CO 2 , CH 4 , N 2 , and their binary mixture foams. They observed that foam strength depends on gas solubility in the liquid phase and found that the N 2 -foam was strongest and the CO 2 -foam weakest. In our model the solubility (at T = 15 MPa) of CO 2 , CH 4 , and N 2 , is 2124.87 mol m −3 , 119.50 mol m −3 , and 60.72 mol m −3 , respectively. Figure 8a shows that the growth rates for the gases are similarly separated, with the rate for CH 4 being closer to N 2 than CO 2 .  In Fig. 8b, we transform Fig. 8a by non-dimensionalizing the time by the total time required for the large bubble to span the channel in each simulation. Figure 8b shows that, when time is rescaled in this manner, the bubble evolution follows a unique curve for all the cases (deviations are within 2.5%), irrespective of gas type and compressibility factor.

Bubble interactions in a pore channel with disconnected liquid ganglia
Thus far we have investigated Ostwald ripening of gas bubbles in porous media by assuming the surrounding liquid is a continuous wetting phase with constant pressure. This can be a reasonable assumption under strongly wetting states where the liquid maintains connectivity across gas bubbles through wetting films along the pore walls (Blunt 2017). However, at less wetting states, or if the wetting films rupture, both gas and liquid form disconnected ganglia with different pressures. During Ostwald ripening, any growth or decay of gas bubbles will pressurize or depressurize the surrounding liquid ganglia. Hence, the interfaces between a gas bubble and different liquid ganglia will generally have different curvatures or capillary pressures.
We explore the impact of a disconnected liquid phase on bubble coarsening using a 2D sinusoidal pore-channel geometry of size 2.65L * ×0.7L * μm 2 , see Fig. 9. In dimensionless variables, the pore space Ω p is defined as: where Initially, the pore channel contains three gas bubbles separated by saturated brine. As before, we model the brine/CO 2 fluid system at the reservoir condition P = 15 MPa and T = 323.15 K, this time combined with VdW EoS. All other input and simulation parameters are the same as before. We perform one simulation with conservation of gas only, which describes gas bubbles surrounded by connected brine, and another simulation with conservation of both gas and brine, which describes gas bubbles separated by disconnected brine. Figures 9a and 10 show evolutions of the fluid configurations, gas-bubble interface radii, and gas-bubble volumes, for the simulation with connected brine. The results show that the middle bubble grows while its interface radius changes non-monotonically over time due to the confining pore geometry. First, the middle bubble grows at the expense of the right bubble. This growth leads to a continuous increase of the middle bubble interface radius over several capillary equilibrium states. To incorporate these larger interface radii while maintaining the contact angle, the bubble slightly shifts its contact location with solid surface and moves gradually into the wide pore region on the right where it contacts the pore wall on the opposite side (image (ii) of Fig. 9a). Subsequently, the middle bubble grows with decreasing interface radius while the left bubble shrinks and the right bubble dissolves completely (image (iii)). The system reaches the stable state when the middle-bubble interface radius has started to increase again and (31) y 1 (x) = 1 100 52 + 5 sin 25 9 x + 8 sin 50 9 x , y 2 (x) = 1 100 17 − 5 sin 25 9 x − 8 sin 50 9 x . the left bubble ceases to exist (image (iv)). The left and right bubbles contract freely during the evolution, and hence both their interface radii and volumes exhibit convex trends over time, as shown in Fig. 10. Figures 9b and 11 show corresponding results for the simulation with disconnected brine, in which both fluids are trapped inside the pore space. As opposed to the case with connected brine, the liquid ganglia will generally have different pressures, implying that a gas bubble can have interfaces with different radii. By numbering the water and gas domains as illustrated in image (i) of Fig. 9b, we label the interfaces from left to right as W1G1, G1W2, W2G2, … , W3G3, G3W4 . Initially, the radii of interfaces G1W2, W2G2 and G2W3 are equal, implying equal gas pressures in bubbles G1 and G2 and equal water pressures in domains W2 and W3. Gas bubble G1 is almost static as its volume and the radii of interfaces W1G1 and G1W2 are about constant in time, see Fig. 11. The radius of interface W2G2 is also constant in time, so that any pressure change in bubble G2 leads to similar pressure changes in regions W2, G1, and W1. Thus, in this simulation the mass transfer occurs mainly from the right bubble G3 to the middle bubble G2. Due to the conservation of neighbouring water domains, the middle bubble G2 evolves differently in this simulation compared to the connected-brine case. To maintain conservation of water domain W2 and W3, interface W2G2 is static while interfaces G2W3 and W3G3 move to the right over time with respectively increasing and decreasing radius (see image (ii) of Fig. 9b). Eventually, water domains W3 and W4 merge when the right bubble G3 loses contact with the pore walls (image (iii)). In the stable state (image (iv)), the pressures in the two remaining bubbles are equal as the radii of G1W2 and W2G2 are equal. However, the radii of interfaces W1G1 and G2W3 differ significantly because the remaining water domains have different pressures, see Fig. 11a. Figure 12 compares the mass transfer rate in the two simulations. During evolution, our model has managed to maintain a constant total mass, with a total mass loss of 0.3% over time. The presence of a disconnected brine shortens the time to reach a stable state, although it prolongs the life of smaller bubbles. The time required to reach stability in the simulation with disconnected brine is around 54% of the time required to reach stability with connected brine.

Coarsening of a residual gas configuration in sandstone
Finally, we apply the GB-MLS method to explore the significance of Ostwald ripening on a residual configuration with trapped CO 2 ganglia surrounded by brine in sandstone. For this purpose, we have extracted a 3D pore-space image containing 200 3 voxels from a segmented micro-tomography data set of Castlegate sandstone, which is publicly available (Sheppard and Schroeder-Turk 2015). Porosity of the extracted image is 22.7% and voxel length is 5.6 μ m. We model a brine/CO 2 system at a reservoir condition described by P = 15 MPa and T = 338.15 K, and use SRK EoS to calculate fugacity. Henry's constant is H = 7.57 × 10 −5 mol m −3 Pa −1 , interfacial tension is gl = 32.6 × 10 −3 N m −1 , and contact angle is = 20 • . The characteristic length is set equal to voxel length, L * = 5.6 × 10 −6 m, while * = 10 −2 N m −1 . For the level set calculations, we use grid spacing Δx = 1 and perform reinitialization every five iteration steps. The simulations utilize 320 computing cores.
We create the residual gas configuration by simulating a saturation history composed of primary drainage and imbibition, using level set methods for capillary-controlled displacement (Helland et al. 2017;Jettestuen et al. 2021). In primary drainage, gas enters the water-filled rock sample through bottom boundary where we place a layer of gas-filled pore space at a given constant pressure. We set the values at the boundary by extrapolation from bulk values. The other faces of the rock use reflective boundary conditions. Gas invasion occurs by increasing the capillary pressure stepwise after each equilibrium state, without conserving any of the phases (Helland et al. 2017). The following imbibition adds a water-filled, water-wet porous plate at the top boundary to allow water invasion, while gas can only escape through bottom boundary. Following Jettestuen et al. (2021), we enforce local volume conservation to disconnected gas ganglia and simulate water imbibition by decreasing stepwise the capillary pressure between connected phases, starting from an initial water saturation of 5% after primary drainage. In the final state after imbibition, the residual gas saturation is 24% and all of the gas is capillary trapped in ganglia.
Then we simulate Ostwald ripening on the residual configuration (containing 76 gas bubbles), using GB-MLS model with gas conservation and reflective boundary conditions  Figure 13 compares the gas bubble configurations before and after Ostwald ripening. It shows that some bubbles disappear while others extend into, or retract from, the nearby pore spaces. Figure 14 shows that the difference in bubble pressures generally decreases with time. However, a temporary incremental difference in bubble pressures is a precursor to bubble loss. This is evident from Fig. 15, in which the loss of a bubble follows a sharp rise in pressure deviation. During evolution, this pressure deviation slowly approaches zero. Figure 15 also shows the progression of the number of bubbles in the system over time. The number of bubbles is 76 initially and 48 at the stable state, which yields a total bubble loss of 37%. We note that 50% of this loss occurs within the first 0.2 time fraction, while 75% of the loss occurs within a time fraction of 0.6. This trend is similar to our previous observations, where the mass transfer rate decreases with time as the pressure and the chemical potential of the gas bubbles come closer.
We calculate the surface area-volume relationship for each CO 2 ganglion to evaluate impacts of coarsening on the morphology of the residual gas configuration. Previous works have shown that ganglia distributions imaged in 3-D porous media exhibit a power-law  (Iglauer et al. 2013;Carroll et al. 2015). Here, p = 2∕3 describes the situation where ganglia preserve their shape upon volume changes, which is the case for spherical bubbles that grow or shrink in response to Ostwald ripening. However, due to the confined pore space, larger ganglia get more ramified inside porous media, and thus ganglia distributions of the nonwetting phase typically exhibit p-exponents between 0.7 and 0.8, where p ≈ 0.75 is a common value (Iglauer et al. 2013;Carroll et al. 2015). Figure 16a shows the area-volume relationship from our simulation. By fitting the power law to the data, we obtain p = 0.734 for the initial CO 2 -bubble configuration and p = 0.774 for the final configuration after Ostwald ripening. The slight increase in p-exponent is a result of the coarsening, in which large ganglia get more ramified as they grow while small, almost spherical, bubbles dissolve completely. In Fig. 16a, this leads to relocation of the data toward higher ganglion areas and volumes over time. However, the total interfacial area between gas and water decreases from 6.37 to 5.705 mm 2 during the coarsening, which leads to a reduction in interfacial energy of the system. Wang et al. (2021) suggested a linear relationship between surface energy and volume for multi-pore bubbles in a regular 2D pore geometry consisting of circular solid grains. Figure 16b shows the energy-volume relationship from our simulation. Here, surface energy includes both fluid-fluid and fluid-solid contributions. It can be seen that the trends are similar to those in Fig. 16a. In our 3D Castlegate sandstone model, most bubbles are smaller than 10 pore body spaces. The smaller size, irregular pore geometry, and an additional spatial dimension might have been the cause of this sublinear relationship.

Conclusions
In this work, we have proposed a new methodology for investigating Ostwald ripening of gas bubbles in porous media. The approach (called GB-MLS method) couples a conservative multiphase level set method for capillary-controlled displacement with Fig. 15 Time evolution of the number of CO 2 bubbles, and the standard deviation of CO 2 -bubble pressures, from the simulation of Ostwald ripening in Castlegate sandstone a ghost-bubble method that describes gas-bubble interactions based on differences in chemical potential. We implemented the method in a software framework for parallel computations. The GB-MLS method handles arbitrary pore geometries and arbitrary gas/water configurations at various reservoir conditions. Specifically, the method accounts for the following challenging scenarios: (i) Gas bubbles that occupy multiple pores; (ii) capillary instabilities induced by bubble growth or shrinkage; and (iii) a liquid phase composed of ganglia with different pressures so that their interfaces with a single gas bubble generally have different interface curvatures. GB-MLS simulations confirm that gas-bubble coarsening in bulk liquid (in absence of porous media) behaves like Ostwald ripening in solid solutions, whereas in porous media both coarsening and anti-coarsening occur, which generally leads to stable states with multiple bubbles confined by the pore geometry. The interface evolution from simulations of four different cases of mass transfer between two bubbles placed in a pore-body/pore-throat geometry are consistent with previously reported results in a similar geometry. The conclusions are as follows: • Simulations on a micromodel with permeable channels suggest that the impact of Ostwald ripening is more significant in heterogeneous or micro-fractured porous media due to larger difference in chemical potentials (or pressures) among the residing gas bubbles. Consequently, growth of bubbles in a micro-fracture can eventually mobilize trapped gas and reduce residual gas saturation. • The mass transfer rate depends strongly on gas type and equation of state. Gas solubility and compressibility factor are proportional to the mass transfer rate. • Simulations of bubble coarsening in the presence of a disconnected liquid phase suggest that the role of disconnected liquid ganglia is to stabilize more bubbles and to prolong the lifetime of bubbles that eventually cease to exist. We expect this scenario to be most relevant for intermediate-wet media. As such, Ostwald ripening is more significant for strongly-wet states where the liquid could be treated as a connected phase with constant pressure. • Simulations on a residual gas configuration in sandstone show that 37% of the gas bubbles disappear due to Ostwald ripening. Analysis of the surface area-volume relation for gas ganglia during the coarsening reveals that larger ganglia get more ramified as they grow, while small spherical bubbles dissolve completely.
GB-MLS method provides a methodology to assess the impact of Ostwald ripening on stable real gas configurations in complex pore geometries at reservoir pressure. This is important to assess the permanency of residual trapping for CO 2 and gas storage, as well as to estimate gas/water interfacial area for other mass transfer processes, like dissolution trapping. A natural extension of this work is to investigate how the presence of residual oil ganglia impacts gas-bubble interactions caused by Ostwald ripening, which is relevant for gas storage in depleted oil reservoirs.

Appendix 1: Ghost bubble pressure and fugacity
For an intermediate liquid region, the total mass transfer at any time instance should be zero. This ensures that there is no increase in the concentration of gas in liquid beyond its saturation level. The ghost bubble represents this intermediate liquid region (Fig. 17).
The following section derives ghost-bubble fugacities for real gases. The same methodology is used to derive ghost-bubble pressures for ideal gases. At any time, t ′ , for a region surrounded by N bubbles, Using Eq. (10) the above equation can be rewritten as If the proportionality parameter k l is proportional to the y th power of the interface area, A y , and B is a scalar quantity, then From Eqs. (34) and (35): Since B is a constant, Eqs. (36) and (12) yield Thus, at time step t ′ we obtain For ideal gases the pressures replace the fugacities. Thus, the ghost-bubble pressure for ideal gas is

Appendix 2: Proportionality parameter
The proportionality parameter, k l , links molar rate of mass transfer to chemical potential difference. It helps in quantifying the relative mass transfer among the gas bubbles and will have unit dimensions of mol 2 s kg −1 m −2 . For simulations, a reasonable value of the constant of proportionality will be required. This value can be estimated from dimensional analysis. We know that the diffusive transfer of gas molecules through a liquid phase is proportional to diffusion constant D, Henry's constant H, the solubility of the gas in the liquid s, and interface area A i . By Graham's law of gas mixture, the mass transfer rate is also inversely proportional to the molar mass of gas M. Verdes et al. (2004) have reported that interfacial tension can also affect the mass transfer. Thus, Eq. (10) can be rewritten as where B is a dimensionless factor. As we are calculating net mass transfer during iterative time, we can assume B to be 1. Hence,  From dimensional analysis we can determine possible values for the exponents by solving the following set of equations: Table 2 presents the possible solutions. We choose parameter combination no. 3 in Table 2, as it is reasonable to expect that the mass transfer through an interface is directly proportional to the interface area. This choice also maintains similarity with Eq. (3) used by Lemlich (1978). Thus, for the rate Eq. (9) that we use in this work, we obtain an effective phase permeability coefficient J ′ given by (49)

Conflict of interest
The authors declare no conflict of interest.

Code availability Available on request
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.