The dynamics of gas-bubble formation at saturated conditions in porous media flow

We investigate the stability of gas bubbles formed at saturated (bubble-point) conditions during two-component (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {CO}_2$$\end{document}CO2, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {H}_2$$\end{document}H2O), two-phase (gas, liquid) flow by developing and analyzing a \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2\times 2$$\end{document}2×2 dynamical system describing flow through a single pore to study the dynamics of gas bubble formation and evolution. Our analysis indicates that three regimes occur at conditions pertinent to petroleum reservoirs. These regimes correspond to a critical point changing type from an unstable node to an unstable spiral and then to a stable spiral as flow rates increase. In the stable spiral case gas bubbles will achieve a steady-state finite size only if they form within the attractor region of the stable spiral. Otherwise, all gas bubbles that form undergo, possibly oscillatory, growth and then dissolve completely. Under steady flow conditions, this formation and dissolution repeats cyclically.

Compositional flow involving a dissolved gas is of importance in many areas, including oil reservoir production 3,14,22 , pipeline transport 12,19,25,28 , CO 2 sequestration 6,27 , and the disposal of radioactive waste 4,5 . Such flow involves the inherent possibility of creation of a gas phase and its subsequent transport. Excluding specific tertiary recovery practices such as CO 2 foam flooding 24 , keeping potential gas components dissolved in fluid phases is important for efficient extraction in reservoirs; the presence of gas bubbles and the resultant fluid-gas menisci complicates flow and can compete with fluid movement. The ability to prevent bubble formation through control of formation pressure or flow rates is therefore important for extraction efficiency.
Once formed in a porous medium, the gas phase is typically non-wetting. During imbibition (displacement by a wetting phase) a non-wetting phase can either exit pore space completely under piston-like displacement, or a fraction of it may become trapped in the form of one or more bubbles by a process known as snap-off. The snapoff process is strongly dependent on pore-geometry, wettability, viscosity and interfacial-tension conditions 26 . While snap-off is an important process in trapping the non-wetting phase, the process requires a starting condition in which both phases-non-wetting and wetting-are present.
In this study, we focus instead on the process by which a gas phase comes out of a saturated solution-forming bubbles. We therefore directly address the transition between single-and two-phase flow and the initial dynamics of those gas bubbles that do form. We are consequently studying flow conditions that would occur before sufficient (local) volumes of non-wetting phase have formed to begin any process of snap-off.
A challenge to the numerical simulation of compositional flow in porous media is the change in the system of equations that accompanies the appearance or disappearance of the gas phase. This difficulty has been addressed in several computational approaches 1, 2, 6-8, 10, 21 . In our computations of two-phase, two-component ( H 2 O , CO 2 ) flow in a 3D pore network 6 , we noted the periodic appearance and dissolution of the gas phase in certain pores. Intensive evaluation of our algorithms led us to conclude that this observed cyclic phenomenon was not numerical in origin. In reviewing the literature on gas transport in porous media at reservoir scales 13 , in micromodel studies 23 , specific studies on gas bubble formation 17 , and mathematical studies of gas phase disappearance in water-hydrogen systems 15 , we have been unable to find any mention of this periodic phenomenon. We have therefore pursued a mathematical investigation. In this article, we extract a 2 × 2 dynamical system from the mathematical model upon which our computations in Chang and Lindquist 6 were based in order to study the dynamics of gas bubble formation and evolution and the mechanics of this phase-cycling phenomenon. We summarize the mathematical model and derive the dynamical system. We analyze the direction fields for this the mathematical model. Under single-phase flow conditions (flow remains under-saturated with dissolved CO 2 ) the fluid transport is governed by the equations where C Cl and p l are, respectively, the concentration of CO 2 and the pressure of the liquid phase in the pore and l is the conductivity of the liquid phase from the pore to the outlet. The liquid phase is assumed incompressible with H 2 O remaining the dominant species; therefore the concentration of water in the liquid phase is assumed to be the constant value C W = 1/18 mol cm −3 . Under these assumptions, the equations in (1) simplify to which provide the single-phase liquid pressure of the water and the CO 2 concentration in the pore.
A gas phase is generated in the pore when the concentration of CO 2 exceeds the solubility (see equation (17) of Chang and Lindquist 6 ) of dissolved CO 2 in the liquid phase at the pressure p sp l , and the resulting two-phase flow follows the system below (these are, respectively, Eqs. (1), (9), (3), (4), (14) and (15) of Chang and Lindquist 6 applied to our single pore geometry): In this system, s l and s g are the saturations of the liquid and gas phases; p l and p g are the pressures in the liquid and gas phase in the pore; C Hg is the gas phase H 2 O concentration; C Cl and C Cg are the liquid and gas phase CO 2 concentrations; m H , m C are the pore averaged concentrations of H 2 O and CO 2 ,

Figure 1.
A cross section of the axially symmetric 3D geometry. The arrow indicating flow direction is placed on the axis of symmetry.
Scientific RepoRtS | (2020) 10:13175 | https://doi.org/10.1038/s41598-020-69506-w www.nature.com/scientificreports/ and Q l , Q g are, respectively, the volumetric flow rates of the liquid and gas phases exiting the pore, In (7), f (s g ) and g(s g ) are the relative permeabilities of the liquid and gas phases respectively. As in Chang and Lindquist 6 , we adopt the models 11 Equations (4) and (8) are written to reflect the assumption that the outlet reservoir is sufficiently large that it can be assumed to remain pure water (i.e. any dissolved CO 2 in the outlet reservoir mixes to negligible concentration and exiting gas bubbles rise due to buoyancy and cannot flow back into the pore). Therefore, under the backflow condition, p l < p , (4) and (8) reflect the fact that pure water flows into the pore from the outlet. The variables l and g are the channel conductivities of the liquid and gas phase between the pore and the outlet reservoir. Conductivity values are modeled from the Hagen-Poiseuille equation, � l = πa 4 /8ν l L and � g = πa 4 /8ν g L . The viscosities, ν l and ν g , of liquid and gas phases are given in units of Pa s by Lide and Kehiaian 18 , Equations (5) express the equality (between the liquid and gas phases at equilibrium) of the chemical potentials of H 2 O and CO 2 . Inherent in our analysis is the assumption that phase changes occur rapidly compared to fluid transport so that equilibrium conditions are essentially maintained. µ ⊖ βl and µ ⊖ βg are the chemical potentials of the species β = H 2 O or CO 2 in the liquid and gas phase at standard conditions. p * H and p * C are standard pressures which we take to be 1 bar. Equations (5)  Assuming the liquid phase is perfectly wetting, the capillary pressure in (3b) is modeled by the Young-Laplace equation, where γ is the surface tension and r is the radius of the gas bubble. We assume that the gas phase forms as a single bubble in the pore and therefore estimate r = s 1/3 g R pore . The surface tension, γ , is evaluated by the Eötvös rule, where T c = 647 • K is the critical temperature for water and κ = 2.1 × 10 −7 JK −1 mol −2/3 . We also assume the gas phase is ideal; therefore the total pressure of the gas phase is the dynamical system. It will be convenient to divide concentration variables by C W , pressure variables by C W RT , flow rates by V, and denote the resultant variables using "over-bar" notation (e.g. As a result, all variables become dimensionless except for the volumetric flow rates, Q , Q l , Q g , which have dimension of time −1 . The natural time variable (through which to introduce dimensionless flow rates) is the quantity Q −1 . However, the flow rate Q will be a critical variable in our analysis. We therefore retain the variable Q in our equations and consequently keep an explicit non-dimensional time variable in the dynamical system.
We identify the fundamental variables, w, y, of the dynamical system as (6) m H = s l C W + s g C Hg , m C = s l C Cl + s g C Cg ; (7) Q l = � l f (s g )(p l − p), Q g = � g g(s g )(p g − p).
ln(ν l ) = −10.4349 − 507.881 (149.39 − T) , Scientific RepoRtS | (2020) 10:13175 | https://doi.org/10.1038/s41598-020-69506-w www.nature.com/scientificreports/ Since the inlet solution is an aqueous phase, we require C C < 1 . Therefore w ≤ C Cg (C C ) < K C . From (9), the concentrations C Hg and C Cl can be expressed in terms of w: Note that the definition of C Cg (C C ) in (12) and the relation for C Cl in (13) gives the identity From (10), (11) and (3b), the dimensionless pore phase pressures are expressed in terms of w and y as where we have defined s = 2γ /(R pore C W RT) . From (14) we note p l < p g . Equations (4) become We can obtain alternate expressions to (15) for H and C . Starting from (6), we have Taking the time derivative of (16) gives the alternate expressions, where for notational brevity, we have introduced the functions G 1 , G 2 , F 1 , F 2 . Solving (17) for dw/dt and dy/dt gives the dynamical system Using (13) to evaluate the concentrations and concentration derivatives in (17) and noting that K H < 2 × 10 −5 and K H /K C < O(10 −3 ) over the range of temperatures of interest, we have, to good approximation, The statement K H /K C ≪ 1 has a physical consequence, namely that the H 2 O component in the gas phase remains small in comparison to the CO 2 component. Approximations (19) imply the simplifications: of (14) to of (15) to  (18) to The second and third equations of (14′) hold only for w ≫ K H = O(10 −5 ) . This approximation will ultimately effect computation of phase space trajectories for values of w 10 −4 . Our results will demonstrate that the region of phase space governing bubble formation occurs well away from "small w" values, and we proceed with the approximation (14′). In (18′) we have used a standard notation, F(w, y) and G(w, y), for the right-hand side of this two-by-two system of ODEs. The dynamical system we consider is (18′) with: H , C defined by (15 ′ ); F 1 , F 2 defined by (17b); Q l , Q g given by (7) and (8) appropriately normalized, where � α = � α C W RT/V , α = l, g ; and the phase pressures given by (14′).

Analysis of the dynamical system
It contains two important curves. The first is the curve Ŵ l0 defined by the condition Q l (w, y) = p l (w, y) − p = 0 which, from (14′), can be written either as The second important curve in the phase space, Ŵ g0 , is defined by i.e., Ŵ g0 is the vertical line w = p.
The curves Ŵ l0 and Ŵ g0 divide the phase space into three regions. Let: Ŵ ++ denote the region of the phase space where p g (w) > p l (w, y) > p ; Ŵ −+ denote the region where p g (w) > p > p l (w, y) ; and Ŵ −− denote the region where p > p g (w) > p l (w, y) . Figure 2 provides sketches of Ŵ l0 , Ŵ g0 and indicates the regions Ŵ ++ , Ŵ −+ and Ŵ −− . The direction field dy/dt = G(w, y). The direction field dy/dt is analyzed in Appendix A of the supplementary material by considering the isoclines of G(w, y). All isoclines in the region Ŵ −− ∪ Ŵ g0 ∪ Ŵ −+ ∪ Ŵ l0 are of the monotonic form (A.9) given in Appendix A and are negative-valued. The curve Ŵ g0 is not an isocline; explicit evaluation using (18′), (15′a), (7′a), (3b), (21) and (14′) gives There exists a maximum value of Q such that, for flow strengths above this value, the isocline G(w, y) = 0 does not (and hence no positive-valued isoclines can) appear in the physical phase space and no gas bubble growth is possible. Computation for this maximum flow strength Q max = α max � l s , is illustrated in Fig. A.1(b) of Appendix A. With the condition Q < Q max satisfied, the form of the isocline G(w, y) = 0 is shown in Fig. 3 where, following the notation introduced in Appendix A, it is referred to as W G (y) . W G (y) intersects the phase space boundary w = C Cg (C C ) at the values y G1 and y G2 .
The direction field dw/dt = F(w, y). The direction field F(w, y) is more complex and we have relied on estimation procedures to determine its behavior. To proceed, we express F(w, y) in (18′) as Since the denominator in (22) is always positive, we note that F and F have the same zero-points and signs. Using (15 ′ ), and (7 ′ ) we obtain the following expressions for F :  (17b) and (13), F 2 (w) = w − w/(K C − w) . Note F 2 (w) < 0 providing K C < 1 , which is indeed the case in our numerical computations (see Table 1). However, F 2 (w) + C C > 0 in the domain of interest. The expressions for F(w, y) in (23) are continuous across the curves Ŵ l0 and Ŵ g0 . We introduce four restrictions on the range of Q to make our estimations tractable. Evaluating (23b) at the point (p + s, 1) , we have Restricting Q < (p + s)� g s/ F 2 (p + s) + C C will guarantee that F(p + s, 1) < 0 . Evaluating (23b) at the point (C Cg (C C ), y l0 (C Cg (C C ))) , we have Restricting Q > � g s s/(C Cg (C C ) − p) 8 will guarantee that F(C Cg (C C ), y l0 (C Cg (C C ))) > 0 . Evaluating (23d) at the point (p, 1) we have Restricting Q < −F 2 (p)� l s/(F 2 (p) + C C ) will guarantee that F(p, 1) < 0 . We also note from the previous section that the restriction Q < Q max guarantees that there is a region G(w, y) > 0 within Ŵ ++ . We combine these restrictions on Q into the single statement, The behavior of F(w, y) in the phase space is discussed next.
F(w, y) on Ŵ l0 . In Appendix B (supplementary material) we show F(w, y) has the following behavior on Ŵ l0 : The existence of the point (w * , y * = y l0 (w * )) within the physical domain is guaranteed by assumption (24). The point (w * , y * ) is illustrated in Fig. 3.
F(w, y) on Ŵ g0 . Using (23d) we see that the derivative ∂ F(p, y))/∂y = −F 2 (p)(K C − p)� l sy −4/3 /3 > 0 . Hence F(p, y) is a strictly increasing function of y. From restriction (24) we have F(p, 1) < 0 . Therefore, F(p, y) < 0 for all y ∈ (0, 1]. F(w, y) in Ŵ ++ . F(w, y)| Ŵ ++ is analyzed in Appendix C (supplementary material). There we prove the existence of a single curve segment, W 1 (y) , in Ŵ ++ on which F(w, y)| Ŵ ++ = 0 . W 1 (y) connects the point (w * , F(w, y) in Ŵ −+ . F(w, y)| Ŵ −+ is analyzed in Appendix D (supplementary material). There we show Ŵ −+ contains a curve, y max (w) , on which ∂ F(w, y)/∂y| Ŵ −+ = 0 . y max (w) joins a unique point ( w 1 , 1) , where p < w 1 < p + s , to a unique point (w Ŵ , y Ŵ ) on Ŵ l0 . In addition, F(w, y)| Ŵ −+ = 0 in Ŵ −+ on a single curve segment W 2 (y) joining the point (w * , y * ) to (C Cg (C C ), y F2 ) , where y F2 is a unique value satisfying 0 < y F2 < y l0 (C Cg (C C )) . The curve W 2 (y) can have two possible forms, depending on the relative size of w Ŵ and w * (cases DI and DII in Appendix D). The curves y max (w) and W 2 (y) are sketched in Fig. 3  The curve W F is complicated to analyze. However, examining (15 ′ b) in Ŵ ++ we recognize that the curve W C defined by the equation also crosses the curve W G at exactly the same critical points. Equation (25) can be explicitly solved for y| W C , We are only interested in the form of W C within Ŵ ++ .
By examining the behavior of the curves W G and W C with respect to Q , it is straightforward to determine that there exist values Q 1 and Q 2 such that, for Q ∈ (Q 1 , Q 2 ) , W G and W C cross at two critical points c 1 = (w c1 , y c1 ) and c 2 = (w c2 , y c2 ) in Ŵ ++ . At each of the values Q = Q 1 and Q = Q 2 , W G and W C touch tangentially at a single critical point. And for Q < Q 1 or Q 2 < Q there are no critical points. A sketch of this behavior is given in Fig. 5. Critical points c 1 and c 2 are also shown in Fig. 3, which is drawn for the case Q 1 < Q < Q 2 .

Critical point classification.
A critical point is classified by analyzing the solution of the linearized form of (18′) in the local neighborhood of equilibrium 16,20 . With u = w − w c , v = y − y c denoting the components of a small perturbation from the critical point, the linearized system is  Solution trajectories. Single-phase flow trajectory. In the single-phase flow regime, the dynamical system variable y satisfies y(t) = 0 . From (13) we have the relation w = K C C Cl /(1 + C Cl ) . We use this relation to extend the definition of the variable w into the single-phase regime. Then (2) provides the time development of w(t) during single phase flow.
Transition to two-phase flow; initial condition for gas bubble formation. An initial condition (w 0 , y 0 ) is needed to solve the two-phase flow system (18′). This is set by the initial size of the gas bubble, y 0 , and its initial CO 2 concentration, w 0 . Initial bubble size is determined by microscopic, non-linear dynamics at a nucleation site 9 that are beyond the scope of this paper. In order to perform numerical computation using a transition from the single-phase to the two-phase region we determine initial bubble formation time, size and CO 2 concentration using approximations outlined in Appendix F. We recognize that conditions (F.1)-(F.3) provide an approximation to actual bubble formation. While we use these conditions in some of our numerical solutions, we also explore other bubble trajectories under the recognition that our initial conditions may not be sufficiently accurate. We note that unless ( w 0 , y 0 ) lies on a trajectory that either starts in, or enters, the region G(w, y) > 0 , the bubble size will monotonically decrease to 0 (the bubble dissolves).
Solution trajectories in Ŵ ++ and Ŵ −+ . The sign of the direction field F(w, y) changes across the curves W 1 (y) , W 2 (y) , W 3 (y) ; the sign of G(w, y) changes across the curve W G (y) . The critical points c 1 and c 2 play a decisive role in trajectory behaviors. In our numerical computations we show that, as Q increases, c 1 changes type: from an unstable node, to an unstable spiral point and then to a stable spiral point; while c 2 always remains a saddle point.
For purposes of illustration of solution trajectories we consider the case when c 1 is a stable spiral point and c 2 is a saddle point. Let σ (t) = (w(t), y(t)) , w(0) = w 0 , y(0) = y 0 , be the trajectory of a solution of (18′), where the initial condition w(0), y(0) is assumed to lie in the region Ŵ ++ . There are four general patterns for the behavior of σ (t) . These are illustrated in Fig. 6 which is drawn for W 2 (y) having the form of case DII in Appendix D. The first behavior, illustrated by the trajectory labeled σ 1 (t) , is that the trajectory remains in the region of Ŵ ++ and spirals into the critical point c 1 . The remaining patterns shown involve those trajectories that cross the curve Ŵ l0 . The second trajectory type, σ 2 (t) , does not escape the zone of attraction of the spiral point; it enters the region F + G − in Ŵ −+ , re-enters Ŵ ++ , and spirals into the critical point c 1 . The third, σ 3 (t) , is outside of the zone of attraction; after it enters the region F + G − in Ŵ −+ , it never leaves Ŵ −+ , ultimately entering the region F − G − and subsequently reaches the w-axis. The last behavior, illustrated by σ 4 (t) , enters and remains in the region F − G − after crossing Ŵ l0 ; subsequently reaching the w-axis.
For trajectories σ 1 (t) and σ 2 (t) , the gas bubble in the pore approaches a steady-state size. For trajectories σ 3 (t) and σ 4 (t) the bubble first grows to a maximum size before decreasing and ultimately dissolving completely. (See the section below that discusses how trajectories approach the w-axis.) Solution trajectories Ŵ −− . From the signs of the direction fields on each side of W 3 (y) (Fig. 7) we deduce the following behaviors. No trajectory in the F − G − region to the right of W 3 (y) can enter the F + G − region to the left of W 3 (y) , and any trajectory that starts in the F + G − region must exit. In particular, the trajectory labelled σ d (t) Figure 6. Illustration of possible trajectory behaviors in the vicinity of ( w * , y * ) when the critical point c 1 is a spiral point and c 2 is a saddle. There are four trajectory behaviors, σ k (t) , k = 1, . . . , 4 , illustrated for W 2 (y) having the form in case DII.
Solution trajectories approaching the w-axis. Finally, we address the slope of trajectories approaching the w− axis. From (18′) and (22) it is straightforward to show that, in the regions Ŵ −− and Ŵ −+ , Thus the slope, dy/dw, of a trajectory approaches a finite, w-dependent value as y → 0 . Consequently, solution trajectories like σ 3 (t) and σ 4 (t) in Fig. 6 will evolve toward the w-axis. For such a trajectory, σ (t) = (w(t), y(t)) , there will therefore be a time, t d , such that y(t d ) = 0 (i.e. the bubble is completely dissolved) and w(t d ) is undersaturated. Since C Cg (C C ) is greater than w(t d ) , w(t) will then increase (under single-phase flow conditions) for t > t d until w(t) reaches a critical saturation value (Appendix F of the supplementary material), a new bubble will form, and a new trajectory cycle will commence.

numerical computations
Using Maple, numerical computations were performed to verify and illustrate our analytic results. The numerical values of all pore, physical and fluid parameters used in the computations are summarized in Table 1. With these parameter values, (24) sets a flow range limit of Q lo = 7.19379 × 10 −5 ≤ Q ≤ Q hi = 3.41059 × 10 3 . However, the discussion in the section on critical point existence states that critical points can only exist for flows between Q 1 and Q 2 (where Q 2 = Q max ). From the data in Table 1, we have Q 1 = 1.6 × 10 −3 and Q 2 = 9.61673 × 10 3 . Therefore we have the relative sizes Q lo < Q 1 < Q hi < Q 2 = Q max . Note that two critical points always exist for the flow range Q 1 < Q hi (see Fig. 5b). We perform our numerical computations in this range.
Critical point changes and flow regimes. We first examine the Q dependence of the character of the critical point using the linearized model (27). We find three flow regimes: R1: Q 1 ≤ Q < 14.6431 ; R2: 14.6431 ≤ Q < 70.7865 and R3: 70.7865 ≤ Q ≤ Q hi . In R1 the critical point c 1 is an unstable node; in R2 it is an unstable spiral; and in R3 it is a stable spiral. In all regimes c 2 is a saddle point. Figure 8 shows the movement of the point w * , y * and the critical points c 1 and c 2 as a function of Q . We make the following observations. (O1) The huge difference in scales between the w-and y-axes distorts perception of angles. Fig. G (supplementary material) shows a properly scaled view of angles in the vicinity of c 2 for Q = 494.   www.nature.com/scientificreports/ (O2) The curve W G closely approaches Ŵ l0 , both in the vicinity of c 2 and as Q decreases. This is evident in Fig. 8 and Fig. G (Appendix G).
(O3) The critical point c 2 moves very slowly with Q , remaining in the close vicinity of the point (C Cg (C C ) , y l0 (C C )) with w c2 ≤ C Cg (C C ) , y c2 > y l0 (C C ) . At the value Q = Q max , the critical points c 1 and c 2 coalesce on the w = C Cg (C C ) boundary.

Gas bubble dynamics.
To detail the evolution of gas bubbles, we performed computations using the nonlinear system (18′) for a value of Q from each regime. The first computation is for regime R3 using the value Q = 475. Figure 9a shows a portion of a phase-space trajectory ( σ 5 ). For this value of Q , using the bubble-point equations (F.1)-(F.3) the pore reaches critical saturation at t = 8.3296 × 10 −3 . A gas bubble forms in Ŵ ++ with (w 0 , y 0 ) = (6.9806 × 10 −3 , 8 × 10 −3 ) . The gas bubble growth trajectory crosses Ŵ l0 at a time t 1 = 8.6718 × 10 −3 at the point (w l0 , y l0 ) = (6.733216 × 10 −3 , 0.484318) where y l0 ≥ y * . At this point, dw/dt < 0 and dy/dt < 0 , σ 5 enters Ŵ −+ , and evolves (Fig. 9b) to the y-axis with the complete dissolution of the bubble at a finite time t 2 = 8.7116 × 10 −3 . For t > t 2 , σ 5 follows the w-axis under single-phase flow conditions with w(t) increasing   Fig. 9b shows the complete trajectory cycle for σ 5 . Figure 10a,b plot w(t) and y(t) over two cycles of σ 5 . The dissolved CO 2 concentration in the liquid phase, w(t), accumulates in the pore until the bubble point is reached. Once the bubble point is reached, time scales shorten dramatically as the bubble expands and subsequently dissolves very rapidly (Fig. 10b). Over this period, the CO 2 concentration in the liquid phase drops precipitously. Figure 9a also shows the trajectory σ 6 for a bubble that is assumed (i.e. not using equations (F.1)-(F.3)) to form with (w 0 , y 0 ) = (6.8437 × 10 −3 , 8 × 10 −3 ) . In this case the trajectory lies within the attractor region of c 1 and spirals into the critical point. The values w(t) and y(t) for σ 6 are plotted in Fig. 10c,d, documenting how the dissolved CO 2 concentration in the liquid phase and the gas phase saturation approach steady state values. Note there is a rather sharp oscillation of the bubble prior to reaching steady state.
The second computations explore regime R2 using the flow rate Q = 68 . Only cyclic bubble formation occurs in this flow regime. If a bubble forms sufficiently close to the unstable spiral point, c 1 , the cycle may involve an initial period of instability in the bubble growth (Fig. 11a). If the bubble initializes on a more "remote" trajectory, the cycle of bubble formation and dissolution is more regular (Fig. 11b). Fig. 11c,d plot w(t) and y(t) over a time range covering two cycles of the trajectory σ 7 in Fig. 11a.
The third computations explored regime R1 using the flow rate Q = 13 . Only cyclic bubble formation occurs in this flow regime. Trajectories are similar to that shown for σ 8 in Fig. 11b.

Discussion
The mathematical analysis presented in this paper, motivated and supported by our earlier network flow computations 6 , leads us to a somewhat unexpected conclusion: namely that when a gas component "bubblesout" from a saturated solution under steady state flow in porous media, the phase space conditions under which such bubbles achieve a stable size is limited. Our results indicate that, for sufficiently low flow rates (regimes R1 and R2), any gas bubble that forms will re-dissolve, possibly (regime R2) after undergoing oscillatory behavior. In the (much wider) range of faster flow rates (regime R3), whether a gas bubble remains stable or decays depends on the fine-scale details of bubble initiation. In all cases involving bubble dissolution, under steady flow conditions the behavior of bubble formation and re-dissolution repeats cyclically.
Our numerical computations show that these regimes cover a range of conditions pertinent to reservoir flow. Therefore, with respect to flow management in reservoir operations (specifically trying to maintain multicomponent, saturated, single phase flow), our results imply that, once pressure conditions allow bubble formation, there is a small "window of opportunity" at low flow rates (regimes R1 and R2) to restrict bubble formation to at Figure 10. Graphs of (a) w(t) and (b) y(t) for two cycles of σ 5 , and of (c) w(t) and (d) y(t) for σ 6 of Fig. 9.
Scientific RepoRtS | (2020) 10:13175 | https://doi.org/10.1038/s41598-020-69506-w www.nature.com/scientificreports/ most cyclic behavior. However once flow rates exceed the rather low threshold delineating regime R3, whether bubbles form stably or cyclically depends on microscopic processes beyond reservoir management control. Our investigation is not exhaustive; we have analyzed the dynamical system (18′) by employing a number of assumptions. These assumptions are summarized in Appendix H (supplementary material) and fall into two general categories, physical (e.g. liquid phase is perfectly wetting) and technical. The technical assumptions rely on having parameter values similar to petroleum reservoirs (Table 1). We have not investigated removing these assumptions, except to make a comment regarding the flow-rate restriction (24). Two of the inequalities in (24) ensure that the point ( w * , y * ) on Ŵ l0 lies within the phase space. This restriction can be lifted to some extent, allowing ( w * , y * ) to move somewhat outside of the phase space. Reducing this restriction will extend the analysis to a larger range of flow values but will produce little change in the qualitative aspects of our analysis.
The y-dependence in the dynamical system is a critical determinant in the resulting behavior. This dependence comes from three sources: the capillary pressure term (10) where we assume r = y 1/3 R pore , the distribution of water and CO 2 mass (16) where the y-dependence is linear, and from the fractional flow functions f(y) and g(y) which introduce complicated nonlinear dependence. It is possible that different forms for the fractional flow functions could significantly alter the gas bubble behavior of this dynamical system.
One possible critique of our analysis is that bubble dissolution, which is governed by the evolution of trajectories in Ŵ −− and Ŵ −+ (i.e. under back-flow conditions), is solely due to our assumption that back-flow involves pure water invading the pore from the outlet reservoir. Such back-flow reduces C Cl concentration in the pore, potentially forcing dissolution of the gas phase. In actual porous media flow situations, back-flow would involve water containing CO 2 (at concentrations less than C Cl ) infiltrating the pore. Back-flow under these more realistic conditions would affect the shape of trajectories in the regions Ŵ −− and Ŵ −+ , opening the possibility that trajectories that enter Ŵ −+ from Ŵ ++ may re-enter Ŵ ++ . However, changing back-flow conditions should not affect the analysis in Ŵ ++ , and specifically not affect the type of critical points. Thus, at low flow conditions there would still be no attractor node to produce stable bubbles. At higher flow rates a stable spiral node exists, but it has a finite-sized region of attraction (as can be inferred from trajectory segments in the vicinity of the saddle c 2 in Fig. 6). Thus not all trajectories would lead to the growth of stable bubbles. Finally, we note that our computations in Chang and Lindquist 6 , which included realistic back-flow conditions, did produce bubble dissolution. www.nature.com/scientificreports/ While we have examined the case where the gas phase is CO 2 , the analysis should hold for any other gas that preserves the assumptions used. Similarly, we have considered a liquid phase consisting solely of H 2 O, with dissolved CO 2 . This invites the question of how applicable our analysis is to a more complicated brine phase.
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 license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license 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 license, visit http://creat iveco mmons .org/licen ses/by/4.0/.