Due to its technological application, two-phase flows in porous media have been an active subject of research for decades [1–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–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–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–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 power-law k _eff ∼Ca ^α [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 two-phase 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–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.

The pore geometry of our system consists of a two-dimensional “Swiss cheese” type of porous medium [33–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 finite-size 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 ₁ 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–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, ρ ∂ u ∂ t + ρ ∇ · u ⊗ u = − ∇ p + ∇ · μ ∇ u + ∇ u T + γ κ ∇ s , ∇ · u = 0 , (1) 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) = μ ₁ ⋅ s(x) + μ ₂ ⋅ [1 − s(x)] is used to define a local effective viscosity of the mixture, where the constants μ ₁ and μ ₂ stand for the viscosities of phase 1 and 2, respectively. The same holds for the density in case of phases with different densities, ρ ₁ and ρ ₂. 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 ∂ s ∂ t + u · ∇ s = 0 . (2)

In this analysis, we keep the density of the two fluids equal ρ ₁ = ρ ₂ = 1.0 kg/m³, to avoid disturbance by inertial effects. Moreover, we set μ ₂ = 1.0 Pa·s and keep the viscosity ratio M = μ ₂/μ ₁ 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⁻⁴ 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 single-phase 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.

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 porous channels or may be trapped for some time by the porous matrix. Moreover, colliding drops may also merge with each other forming larger clusters. This complex dynamics of merging and disconnecting drops dominates the flow’s microscopic behavior in the stationary state. As a result, the velocity fields of the two phases show typical transient fluctuations.

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 ( δ S 1 ) , where this term represents the mobility coefficient and δ is a constant. This behavior is similar to effective permeabilities found in dissolved-gas-driven 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.