Bubble Dynamics in Stationary Two-phase Flow Through Disordered Porous Media

Two-phase flow through porous media leads to the formation of drops and fingers, which eventually break and merge or may be trapped behind obstacles. This complex dynamical behavior highly influences macroscopic properties such as the effective permeability and it also creates characteristic fluctuations in the velocity fields of the two phases, as well as in their relative permeability curves. In order to better understand how the microscopic behavior of the flow affects macroscopic properties of two phases, we simulate the velocity fields of two immiscible fluids flowing through a two-dimensional porous medium. By analyzing the fluctuations in the velocity fields of the two phases, we find that the system is ergodic for large volume fractions of the less viscous phase and high capillary numbers Ca. We also see that the distribution of drop sizes m follows a power-law scaling, P ( m ) ∝ m − ξ . The exponent ξ depends on the capillary number. Below a characteristic capillary number, namely Ca* ≈ 0.046, the drops are large and cohesive with a constant scaling exponent ξ ≈ 1.23 ± 0.03. Above the characteristic capillary number Ca*, the flow is dominated by many small droplets and few finger-like spanning clusters. In this regime the exponent ξ increases approaching 2.05 ± 0.03 in the limit of infinite capillary number. Our analysis also shows that the temporal mean velocity of the entire mixture can be described by a generalization of Darcy’s law of the form v ̄ ( m ) ∝ ( ∇ P ) β where the exponent β is sensitive to the surface tension between the two phases. In the limit of infinite capillary numbers the mobility term increases exponentially with the saturation of the less viscous phase. This result agrees with previous observations for effective permeabilities found in dissolved-gas-driven reservoirs.


INTRODUCTION
Due to its technological application, two-phase flows in porous media have been an active subject of research for decades [1][2][3][4]. It is well known, for instance, that the flow of a water/oil mixture through a porous rock strongly depends on the volume ratio of the two fluids, the surface tension between the two phases and their wetting interactions with the solid walls [5][6][7][8]. Phenomenologically, this behavior has been described by relative permeability curves, which need to be determined experimentally. Over the years, many models have been proposed to explain those experimental measurements, or to approximate them when experimental gauging curves are not available [9][10][11].
Most of these studies focused on the displacement of an interface moving in a transient regime, such as in gas/oil drainage or water/oil imbibition. Understanding this process is of great practical relevance in petroleum engineering where the main goal is to increase production by delaying the breakthrough of the injected fluid, i.e., to increase the amount of the displaced fluid by extending the time the injected fluid takes to reach the producing well [12].
Immiscible two-phase flows in porous media commonly exhibit intriguing phenomena that arise due to the competition between capillary and viscous effects. Generally speaking, the resulting interface formation is controlled by two dimensionless numbers, namely, the capillary number Ca, which describes the balance of viscous forces to capillary forces, and the ratio of the two viscosities. By carrying out multiple experiments, Lenormand classified the emerging patterns in two-phase flow into three major regimes, depending on these two dimensionless parameters. This classification is today known as the so-called Lenormand's phase-diagram [13], which was recently extended including effects of surface adhesion, namely, also considering contact angles [14]. While the general partitioning of the different flow regimes is characterized by Lenormand's phase diagram, the transitions at which one regime changes to another are not universal and depend on the particular pore geometry and the dimensionality of the system [7,15,16]. Low capillary numbers and high viscosity ratios are frequently associated with a "capillary fingering" regime, while high capillary numbers and high viscosity ratios to a "stable displacement." In addition, when a low viscous fluid displaces a high viscous one (low viscosity ratio), the resulting interface forms unstable patterns, known as "viscous fingering" [17][18][19], which can be obtained for any value of the capillary number. This flow regime can be described in terms of a Laplacian growth problem [20], which is also equivalent to diffusion-limited-aggregation (DLA) systems [21]. Experimental studies have also revealed that the capillary number influences the effective permeability of two immiscible phases in terms of a powerlaw k eff~C a α [22,23] and numerical studies suggest that the exponent α depends on the relative fraction of the two phases in the mixture [24].
The description of two-phase flow at the Darcy scale as an effective medium has a long history [25,26], but the dynamics of the interface between the different phases at the pore scale still represents a challenging and scientifically important problem. In this case, the geometric properties of the substrate create heterogeneous flow paths which can be invaded by one fluid or the other or both. In such mixing flows, drops of many sizes and shapes naturally emerge. These drops may split and re-merge while they are dragged by the flow through the porous medium or maybe trapped behind obstacles for some time. Numerical simulations [5] and micro-tomography [6,27] have been used to understand the interplay between two phases inside a porous medium on the scale of individual pores, but the connection between this micro-dynamic behavior and macroscopic quantities such as permeability or displacement efficiency is still poorly understood [28]. Over time, the dynamic of a twophase flow inside a heterogeneous medium eventually reaches a stationary regime, where certain variables fluctuate around a long-time averaged value [8]. In this regime the application of novel techniques from non-equilibrium statistical physics [29] proved to be very helpful in establishing a proper connection between micro and macro scale. Recent studies of fluctuations in mesoscopic properties, such as relative permeability and flow velocity, paved the way to a new perspective on the conceptual description for two-phase flow at low Reynolds numbers [30][31][32] with focus on its statistical properties. Here, we study through two-dimensional simulations the characteristic fluctuations in the velocity time series of two immiscible fluids flowing through an irregular porous medium.
The paper is organized as follows: In Section 2, we present the details of our pore geometry, the mathematical model and numerical technique used to calculate the velocity fields of the two phases. Section 3 covers the analysis of the numerical results. More specifically, we analyze temporal correlations in the spatially averaged time series of the velocity fields of the different phases and the mixture in the stationary regime and apply Onsager's reciprocal relations in order to investigate time reversibility. Moreover, we show that the drop size distribution follows a power-law scaling and propose a generalization of Darcy's law with a non-linear coupling between flow rate and pressure drop in Section 3.4. We close our analysis with discussions and conclusions in Section 4.

MATHEMATICAL MODEL
The pore geometry of our system consists of a two-dimensional "Swiss cheese" type of porous medium [33][34][35], which is made of circular obstacles that can overlap with each other. More specifically, the porous domain is composed of a 50 mm × 50 mm square, which is iteratively filled with randomly placed discs of 1 mm in diameter until a desired porosity ϱ is reached. For all simulations performed here, ϱ = 0.8. Periodic boundary conditions are applied in both directions in order to avoid finitesize effects. The pore space is filled with two immiscible Newtonian fluids with a surface tension, γ, acting at the interface between them. A global pressure gradient, ∇P, in the horizontal direction (x-direction) drives the flow. Figure 1A shows the initial condition of the system, where the filling fraction S 1 of the low viscous phase (blue) is set to 0.2. The two phases flowing through the pore space are solved numerically by a Volume-of-Fluid (VoF) formalism [5,36], which is an adaptation of the Navier-Stokes equations for multi-phase flow as implemented in the Ansys Fluent ™ software [37]. This numerical technique has been previously validated in several studies through numerous distinct applications [38][39][40][41][42][43][44]. In particular, the surface tension and contact angle in the model of Fluent's VoF scheme have been successfully applied to quantitatively describe the droplet pinch-off dynamics in a microfluidic step emulsification device [45].
In the Volume-of-Fluid formalism, the Navier-Stokes and continuity equations describe the conservation of linear momentum and mass of the entire mixture, where ρ, μ, p and u are density, viscosity, pressure, and the fluid velocity, respectively. The term u › u is an outer product, and s(x) is the fraction function of phase 1 in a given control volume. Due to mass conservation, the fraction of phase 2 is given by 1 − s(x). Because Eq. 1 describes the motion of both phases, ρ and μ are not constants, but scalar fields, which depend on the fraction function s(x). More specifically, a linear mixing rule μ(x) = μ 1 · s(x) + μ 2 · [1 − s(x)] is used to define a local effective viscosity of the mixture, where the constants μ 1 and μ 2 stand for the viscosities of phase 1 and 2, respectively. The same holds for the density in case of phases with different densities, ρ 1 and ρ 2 . The term γκ∇s in Eq. 1 describes the interfacial tension force and is proportional to the gradient of s normal to the interface. Here, the variable κ = ∇·n is the curvature of the interface. Because ∇s is nonzero only along the interface, interfacial tension forces vanish in the bulk of phase 1 and phase 2, where s is constant. As boundary conditions on the solid walls, we limit our study to neutral wettability (contact angle α = 90°) and no slip. The fraction function s(x) has a sharp interface at the border between the two phases, and is advected by the flow via In this analysis, we keep the density of the two fluids equal ρ 1 = ρ 2 = 1.0 kg/m 3 , to avoid disturbance by inertial effects. Moreover, we set μ 2 = 1.0 Pa·s and keep the viscosity ratio M = μ 2 /μ 1 at either 1 or 10. During a simulation, the integration time step dt is kept fixed and small enough to avoid high Courant numbers, generally dt ≈ 10 −4 s. If not mentioned otherwise, our fluid mixture consists of 20% of phase 1 (blue, lower viscous phase) and 80% of phase 2 (red, high viscous phase). Initially, the two phases are placed in vertical strips (see Figure 1A) in the fluid domain with their interface perpendicular to the pressure gradient. Different initial configurations have been tested in multiple simulations to ensure that the stationary regime does not depend on this initial condition. Using this setup, we study the behavior of the mixture as we change the two major interactions between the phases, namely, the viscous forces, which is controlled here by the global pressure gradient, ∇P, and interfacial forces, which depends on surface tension, γ. We also ran simulations with vanishing surface tension and equal viscosities in order to test whether our two-phase model recovers the properties of a singlephase flow. Simulation data are only analyzed after the stationary state is reached, in order to produce statistically meaningful time series. This highly increases the computational costs as the initial transient part of the calculation has to be discarded.

RESULTS
Shear and surface tension exert forces on the fluid while it flows through the porous medium. These forces eventually lead to the breakup of connected components of a phase and, thus, to the generation of drops. These drops may flow separately through the

Velocity Time Series
In the VoF method, the volumetric average velocity of a phase, at time t, may be written in terms of the fluid velocity, u, and the fraction function s(x, t). For the x-component of phase 1's velocity we find where Ω 1 = Ω s(x)dΩ is the domain occupied by phase 1, Ω is the total porous space occupied by the fluid mixture and u x (x, t) is the x-component of u, solved with Eqs 1, 2. Equivalently, the mean xvelocity of phase 2 is given by with Ω 2 Ω [1 − s(x)]dΩ. Since the two phases are incompressible and mass is conserved, Ω 1 and Ω 2 are constants. Finally the mean x-velocity of the entire mixture is given by The temporal mean values of these time series are computed through v (j) 1/T t 0 +T t0 v (j) dt, where j = 1, 2, m stands for the two phases (1,2) or the mixture (m). The variable T corresponds to the time window spanning the stationary regime. The capillary number describes the ratio between viscous forces and surface tension being defined here in terms of the mean mixture velocity as The snapshots in panels Figures 1C,F show the distribution of phase 1 (blue) and phase 2 (red) in the pore space for the two cases Ca → ∞ (equivalent to γ = 0) and Ca = 0.006. At high capillary number, the flow is characterized by a large number of small droplets with a highly active merging and splitting dynamics. In some of these cases, the formation of one or more finger-like cluster spanning from left to right along a major flow path can be observed. In the case with low capillary number, the phases stick together and create more cohesive drops, preventing split-merging events. In cases of very small capillary numbers, the pressure gradient may even not be sufficient to overcome the interfacial forces and thus the formation and motion of drops may be suppressed. In this situation, permanently trapped drops are observed in the porous medium. The velocity maps corresponding to the snapshots shown in (C) and (F) are plotted in panels (D) and (G). Darker colors represent regions with low velocities, while brighter colors indicate fast flowing regions.

Time Correlations in Two-phase Flow
In order to develop a description of the flow as a stochastic dynamical system, it is important to know whether the system is ergodic or not. Although originally applied to thermal fluctuations on the molecular scale, Onsager symmetries [46,47] have been successfully applied to describe the behavior of multi-phase flows in macroscopic porous media [48] and in pore-network models [30][31][32]. For this purpose, we apply Onsager's reciprocal relations in order to study ergodicity. Onsager's reciprocal relations are based on time correlations which can be calculated as follows Here ΔA C AB and C BA are equal [49]. This condition is satisfied when C AB (τ) collapses on C BA (τ) at all times, except for some very small timescales in which the contribution to the integral is negligible [32]. Autocorrelations, on the other hand, measure the randomness in the time series and the rate at which fluctuations change with respect to the time scale τ. It can also be interpreted as the memory of a stochastic process, determining whether successive observations are independent or not.
In the following, we analyze the correlations in the time series of the spatially averaged x-velocity of phase 1, phase 2, and the entire fluid as defined through Eqs 3-5. Here, C 12 (τ) and C 21 (τ) describe the cross-correlation between phase 1 (2) and phase 2 (1), where the evolution of phase 2 (1) is shifted forward in time relative to phase 1 (2). C mm (τ) is the autocorrelation of the mixture. Figure 2 shows C 12 (τ), C 21 (τ), and C mm (τ) for ∇P = 1.0 kPa/m and M = 10. Generally, for small delays, phase 1 and phase 2 tend to be anti-correlated and weakly correlated for increasing τ, i.e., if one phase speeds up (comparing to its average velocity) the other phase slows down. In the limit of very large delays, the cross correlations of C 12 (τ) and C 21 (τ) both vanish. Conversely, the autocorrelation functions show the opposite behavior with strong correlations for small delays, slight anticorrelations in a intermediate range, and vanishing correlations for τ → ∞. Figure 2 shows how the filling fraction of phase 1 and the capillary number affect the different correlations in the flow. Panels (A), (B), and (C) show the cross-correlation functions for decreasing Ca, keeping a constant filling fraction of S 1 = 0.2. For Ca → ∞ and S 1 = 0.2, C 12 (τ) and C 21 (τ) seem to collapse only for large values of τ. As Ca decreases, the cross-correlation functions approach zero more rapidly, but the fluctuations also increase and thus the collapse of C 12 with C 21 is not completely clear for small Ca. For capillary-dominant systems, Ca → 0, time-reversibility is less clear, similar to a recent experiments with two-phase flows in three-dimensional heterogeneous systems [30]. The figure also shows the effect of increasing filling fraction of the less viscous phase. Here S 1 is varied from 0.2 in panel (A) to 0.4 in (D), to 0.6 in (E), and to 0.9 in (F) while all other parameters are kept fixed. For filling fractions of 0.4, 0.6, and 0.9 the cross correlation functions almost perfectly collapse onto each other suggesting a more time reversal dynamics compared to the case with only 0.2 filling fraction, as shown in panel (A). In fact, by varying the saturation S 1 from 0.1 to 0.9, a good collapse of the two cross-correlation function is observed for S 1 ≥ 0.4 at all time scales. Regarding the autocorrelation functions (black lines in Figure 2), our results show that they slowly become decorrelated, which is similar to Brownian motion [50]. However, as the surface tension increases, C mm (τ) vanishes sooner, at smaller values of τ.

Drop Size Distribution
In order to study how the agglomerates of drops may influence the mesoscopic behavior of the two-phase flow, we determine the sizes of the different drops for every time step during the stationary regime and then calculate their distributions. Figure 3A shows the drop size distribution P(m) as a function of capillary number for ∇P = 2.0 kPa/m and M = 10, where m is the drop size normalized by the whole fluid volume. The distribution follows a power law, P(m) ∝ m −ξ , whose exponent ξ depends on the capillary number. For Ca ≲ 0.046, the exponent ξ is roughly constant with a value of ξ ≈ 1.23 ± 0.03. In this capillary range, the drops are mostly large and the splitmerging events (as well as the spanning channels) are rarely observed. For Ca > 0.046, however, the drops are shattered into many smaller droplets and few fingers (thin elongated clusters) along higher speed channels. In this regime, ξ increases from 1.73 ± 0.06 to 1.92 ± 0.06 for Ca = 0.078 and 0.150, respectively, and then converging to 2.05 ± 0.03 as Ca approaches infinity. The characteristic capillary number, Ca* ≈ 0.046, separates the two regimes, which are marked by light blue and cream background in Figure 3. The heavy-tail scaling of the drop sizes in the large capillary regime may be associated with the emergence of longrange correlations, similar to those found in anomalous diffusion [51].

Generalized Darcy Law
Next, we analyze the stationary state from a single-phase Darcy flow perspective, namely, how the averaged fluid velocity depends on the pressure gradient. Figure 4 shows the xvelocity of the mixture averaged over space and time, v (m) , as a function of the applied pressure drop ∇P. As depicted, the temporal average velocity follows a power-law scaling, v (m) ∝ (∇P) β , for different values of the surface tension γ. This equation can be interpreted as a non-linear form of Darcy's law for two-phase flows. A similar generalization between flux and pressure drop has been successfully applied to non-Newtonian flows in pipes and pore networks [52,53]. The exponent β varies as a function of surface tension. As In regimes of very high capillary numbers, the flow behavior is dominated by the presence of many small droplets. Figure 5 shows how the filling fraction of the less viscous phase (phase 1) impacts the average velocity for ∇P = 1 kPa/m and M = 10. The plot shows that v (m) increases exponentially with filling fraction v (m) ∝ e (δS 1 ) with an exponent δ = 2.34 ± 0.07, where e (δS 1 ) is a mobility coefficient for the mixture [54]. In a single phase flow, the mobility term is the ratio between the permeability and the viscosity. In a two-phase flow system, the mobility term is related to relative permeability curves [10]. The effective permeability in the form of an exponential increase with saturation fits particularly well with measurements of gas percolation in dissolved gas-driven reservoirs [55], where the oil phase is saturated with dispersed small bubbles, the so-called "foamy oil" [56,57], in order to enhance the recovery rates. Despite the complexity involved in the two-phase flow dynamics, our results suggest that the mean flow velocity of the mixture can be described by a simple function of the saturation and the gradient of pressure.

DISCUSSION AND CONCLUSION
We investigated the stationary flow regime of two immiscible and incompressible Newtonian fluids in porous media by solving Navier-Stokes equations in multi-phase flow in two dimensions. The behavior of the time series of the fluid's velocity is influenced by the complex dynamics of drops which form as the two phases interact with each other and the heterogeneous pore space. Despite the apparent disorder, in the stationary regime, the drop size distribution follows a well-defined power law, P(m) ∝ m −ξ , whose exponent ξ depends on the capillary number. This exponent is roughly constant ξ ≈ 1.23 for Ca ≲ 0.046, where the drops are mostly large and cohesive, and splitting and merging are less common. For Ca > 0.046, however, ξ increases systematically, reaching 2.05 ± 0.03 for Ca → ∞. This regime is characterized by the presence of a large number of small droplets and few finger-like clusters. Fluctuations in the time series are analyzed via cross-correlation functions between the x-velocities of the two phases showing that Onsager's reciprocal relations and time reversal symmetry are fulfilled for volume fraction above 0.4 and high capillary number. At lower capillary numbers Ca the time reversibility of the flow is less clear which is consistent with observations made in three dimensional two-phase experiments. Finally, we study the macroscopic scaling of the average velocity. Our results show that two-phase flows can be modeled by an effective Darcy type of description, namely, v (m) ∝ (∇P) β . In this generalization of Darcy's law, the exponent β depends on the surface tension between the phases (and, thus, on the capillary number). When the surface tension is neglected and the viscosity ratio is unity, the traditional Darcy relation is recovered, β = 1. For Ca → ∞, when the systems is dominated by the presence of many small droplets, the averaged fluid velocity increases exponentially with the filling fraction of the lower viscous phase, v (m) ∝ e (δS1) , where this term represents the mobility coefficient and δ is a constant. This behavior is similar to effective permeabilities found in dissolved-gasdriven reservoirs. We believe that our results have direct applications in the behavior of the mesoscopic flow (at the level of the pore) in several real situations, and can help in the description of the macroscopic propagation of the invasion front in oil reservoirs.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
HS, CO, and JA designed the research. JS performed the simulations with input from HS and CO. All authors contributed to the discussion and interpretation of the results. HS, CO, and JA wrote the paper with input from JS. All authors approved the submitted version.