Quantification of Nonlinear Multiphase Flow in Porous Media

We measure the pressure difference during two‐phase flow across a sandstone sample for a range of injection rates and fractional flows of water, the wetting phase, during an imbibition experiment. We quantify the onset of a transition from a linear relationship between flow rate and pressure gradient to a nonlinear power‐law dependence. We show that the transition from linear (Darcy) to nonlinear flow and the exponent in the power‐law is a function of fractional flow. We use energy balance to accurately predict the onset of intermittency for a range of fractional flows, fluid viscosities, and different rock types.

where f w = q w /q t is the fractional flow. For a displacement with fixed fractional flow, but where the flow rate varies, Equation 2 implies a linear relationship: ∇P ∼ Ca.
It is well known that this linear law breaks down when, at the pore scale, viscous forces become comparable to capillary forces with a threshold value of Ca of around 10 −3 (Blunt, 2017;Lake, 1989): in this regime, the two-phases flow together through the pore space with a fractional flow which is proportional to saturation. Recently, however, a body of research has demonstrated that multiphase flow in porous media has a complex dynamics even at low Ca where displacement is still dominated by capillary forces at the pore scale (Armstrong & Berg, 2013;Berg et al., 2013;Datta et al., 2014aDatta et al., , 2014bReynolds et al., 2017;Rücker et al., 2015). It has been proposed that there is a transition from capillary-dominated to an intermittent flow regime with a power-law dependence between pressure gradient and flow rate   a P Ca (4) with 1 > a > 0 for Ca > Ca i . Tallakstad et al. (2009) conducted steady-state two-phase simultaneous flow experiments in a quasi-two-dimensional porous medium and found a ≈ 0.5 in Equation 4. Rassi et al. (2011) measured an exponent a between 0.3 and 0.45 depending on the fluid saturation from steady-state twophase flow in bead packs. Sinha et al. (2017) also conducted experiments and simulations and proposed a power-law flow regime with an exponent a = 0.5.
High-resolution X-ray tomography has allowed the fluid configurations to be imaged inside porous materials to interpret the physical origin of this nonlinear regime, which is due to the intermittent occupancy of regions of the pore space by the two phases that facilitates flow (Gao et al., 2017(Gao et al., , 2020Spurin et al., 2019aSpurin et al., , 2019b. From a theoretical perspective, Sinha and Hansen (2012) suggested that the exponent a was around 0.5 and confirmed this behavior using a dynamic pore-scale model. Roy et al. (2019) proposed a size dependence of the non-linear regime, such that for larger systems, the transition from intermittency to linear flow occurred for lower Ca in a high flow rate regime. In the limit of an infinite system, this suggests that there is no linear Darcy regime at all. Gao et al. (2020) quantified the threshold capillary number for the onset of intermittency Ca i as ∼10 −5 and found a = 0.6 from experiments of steady-state flow on water-wet Bentheimer sandstone for f w = 0.5; however, there were only eight data points. Overall, despite this body of theoretical, numerical, and experimental work, there is no accurate quantification of when the transition to nonlinear flow occurs and the relationship between pressure gradient and flow rate in this intermittent regime.
In this paper, we study steady-state immiscible two-phase flow through a water-wet Bentheimer sandstone sample with different flow rates (capillary number varies from ∼10 −7 to ∼10 −4 ) and fractional flows during an imbibition displacement (f w = 0.2, 0.4, 0.5, 0.6, 0.8, and 1). We quantify the threshold capillary number for the onset of the power-law non-Darcy regime, Equation 4, and measure the exponent a, which we find to be a function of fractional flow. We use energy balance to predict accurately the onset of non-Darcy flow as a function of capillary number and fractional flow, and for different rocks and viscosity ratios.

Materials and Methods
We performed experiments on a water-wet Bentheimer sandstone sample (5.97-mm diameter and 27.88mm length) mounted in a specially designed core flooding system: details of the experimental apparatus can be found in the literature (Gao et al., 2017(Gao et al., , 2020. A long sample was chosen to allow an accurate measurement of pressure drop. The diameter was chosen to be consistent with previous work (Gao et al., 2020) where pressure measurements were combined with pore-scale imaging. Further details about Bentheimer sandstone and associated pore-scale images can be found in Muljadi (2015). The wetting phase was 15 wt % KI (potassium iodide) brine, and the nonwetting phase was n-decane, both injected by high-precision ISCO pumps. To ensure uniform flow at the inlet, we used a dual injection port (Gao et al., 2019); this provides a more even flow than T-junction injection used in previous studies which may have over-stated the degree of intermittent flow as a consequence (Gao et al., 2020).
During the flow, the pressure gradient between inlet and outlet of the sample was measured by a high-precision pressure transducer Keller PD-33X. The viscosity of the brine μ w = 0.821 mPa s, the n-decane viscosity μ nw = 0.838 mPa s (PubChem, open chemistry database), while the interfacial tension was measured to be σ = 47 mN/m. The absolute permeability for the sample was K = 1.85(±0.02) × 10 −12 m 2 , measured from the relationship between flow rate and pressure gradient when the core was fully saturated with brine. We then injected n-decane at 3 ml/min for 30 min to reach the initial brine saturation.
We started the two-phase flow experiment by injecting the two phases at a fractional flow f w of 0.2 at 0.04 ml/min total flow rate (Ca = 4.2 × 10 −7 ); the pressure gradient was recorded after 10 h when it stabilized, and then we gradually increased the flow rate from low to high and recorded the pressure gradients at steady state. A stable pressure difference was used to determine that the system was at steady state. The time for the pressure gradient to become constant depended on the flow rates: it was up to 10 h for the low flow rates (<0.04 ml/min) but as little as 5 min for high flow rate flooding (greater than 3 ml/min). The highest flow rate was 4.25 ml/ min (Ca = 4.5 × 10 −5 ).
Once we had reached the highest flow rate for a given fractional flow, we injected n-decane at 3 ml/min for 30 min again to return to the initial saturation. We repeated the injection sequence at a series of increasing flow rates for other fractional flows: 0.4, 0.5, 0.6, 0.7, 0.8, and 1; we studied seven fractional flows with 20-30 measurements at different rates for each fractional flow, making a total of 178 measurements.

The Nonlinear Flow Regime
The results, Figure 1, show that Ca is proportional to the pressure gradient ∇P at low flow rates: ∇P ∼ Ca with a = 1 is valid for all the fractional flows studied for sufficiently small Ca. However, in all cases, this is followed by a nonlinear regime, Equation 4. The exponent a is a function of the fractional flow, see Table 1: the highest value a = 0.74 ± 0.02 occurs when f w = 0.2, which displayed the lowest degree of intermittency defined as the deviation from a linear Darcy law; f w = 0.6 had the lowest exponent a = 0.44 ± 0.02 indicating a strong deviation from linear flow.
These observations are consistent with previous experiments where the rock and fluids were imaged at micron resolution, which concluded that intermittency, or a fluctuating pore occupancy, only constituted a small percentage of the total pore volume for f w < 0.5 (Spurin et al., 2019a). Moreover, in our experiments, when f w < 0.5, the pressure gradient increased with fractional flow indicating that the total mobility decreases; for f w > 0.5, a more complex behavior emerged, where curves of ∇P as a function of Ca for constant f w crossed each other, indicating a significant increase in mobility with flow rate for the intermediate fractional flows, were both phases compete to occupy high-conductivity paths through the pore space (Gao et al., 2020).
The threshold capillary number Ca i between the linear Darcy regime and intermittent flow, Table 1  decreased from ∼10 −5.1 to ∼10 −6.4 as the fractional flow f w increased from 0.2 to 1. Thus, higher fractional flows have a lower critical flow rate i t q for the onset of intermittency. Note that the results for a fractional flow of 1 represent the capillary desaturation process (Lake, 1989); at steady state, there is only one mobile phase.

Quantification of the Transition From Linear to Intermittent Flow
We assume that we see the onset of intermittency when the energy, associated with the injection of fluids over a characteristic length l, related to the distance between pores, is first matched by the change in surface energy needed to create a fluid meniscus (Blunt, 2017;Cueto-Felgueroso & Juanes, 2016;Gao et al., 2020). The interfacial energy required to form an interface between the fluids inside a single pore is of order σr 2 where r is a typical pore radius. The mechanical (PdV) work occurs over a volume l 3 or a pore volume of dV = ϕl 3 where ϕ is the porosity. The change in pressure, P, across this length is −l∇P where ∇P is the pressure gradient. The mechanical work PdV is therefore −ϕl 4 ∇P. Hence, intermittency first occurs when the interfacial and mechanical energies match, or We now need to estimate ∇P over a typical pore length. We can assume a Darcy-like flow, although we need to note that we apply it at the pore scale. In imaging experiments, intermittency is caused by the nonwetting phase periodically finding more conductive pathways through the pore space (Gao et al., 2020;Spurin et al., 2019b). We assume that since both phases have to move, the pressure gradient necessary to allow intermittency is controlled by the total threshold flow rate i t q . The limiting mobility is assumed to be (1 − f w )/μ w , governed by the flow of the wetting phase into and out of pores filled with nonwetting phase with an effective relative permeability, at least in the viscous-flow limit, of 1 − f w , or the nonwetting phase fractional flow. Hence, from Equation 1, we estimate  and where the dimensionless number Y i is defined by Note that Y i incorporates information on both pore structure and viscosity ratio. In our experiments, for Bentheimer sandstone, the porosity ϕ is 0.2, the mean pore radius r is 24 μm (Blunt, 2017), while to capture the onset of intermittency l has a value of ∼150 μm (Gao et al., 2020): l is the mean pore-to-pore distance obtained from a pore-network analysis of the pore structure (Raeini et al., 2017). Then using the fluid viscosities in this experiment we calculate Y i ≈ 10 −5 .
In Figure 2, we delineate the threshold between Darcy-like and intermittent flow using Equations 7 and 8 which provides an accurate prediction compared to the values in Table 1. Note, however, that Equations 7 and 8 are only valid when we strictly have multiphase flow with 1 > f w > 0 and cannot be applied for f w = 1 or 0.
Our experiments were performed for one set of fluids with almost equal viscosities. To test the applicability of our theoretical result for different fluids, we apply Equation 9 to the data from Reynolds et al. (2017). As in this work, the rock studied was Bentheimer sandstone, so that the geometric terms in Equation 9 are the same. We show again that we make accurate predictions where we calculate Y i = 10 −6.5 , for nitrogen and brine injection: the threshold line (Figure 3) between Darcy and intermittent flow is consistent with the X-ray pore-scale images described in Reynolds et al. (2017).
Furthermore, we also apply our theory to a different rock sample, Estaillades carbonate with both low (nitrogen and brine, μ nw /μ w = 0.026) and high (hexadecane and brine, μ nw /μ w = 4.21) viscosity ratio fluids; here ϕ = 0.295, r = 7.5 μm, K = 1.49 × 10 −13 m 2 , and l has a value of ∼63 μm (here the pore length is assumed to be the inverse of the cube root of the number of pores per unit volume, while the radius is found from a pore-network analysis) (Blunt, 2017): Equation 9 is used to calculate Y i = 10 −7.3 and 10 −5.1 , respectively. The predicted thresholds are consistent with the results of Spurin et al. (2019b) where X-ray scanning results showed that the low viscosity ratio fluid system displayed strong intermittency, but no intermittency was observed at the pore scale for the high fluid viscosity ratio case although the total capillary numbers were similar. Estaillades has a more heterogeneous pore space with a wider range of pore size than Bentheimer sandstone (Blunt, 2017), and yet we still accurately capture the onset of intermittency for this rock.
We finally apply our theory to other data in the literature for Bentheimer sandstone, Estaillades limestone and glass bead packs (Datta et al., 2014a;Gao et al., 2017Gao et al., , 2019: details are provided in Figure 3 and Table 2. For a bead pack, we estimate l = Dϕ/(1 − ϕ) (Blunt, 2017;Porta et al., 2015;Whitaker, 2013), where D is the bead diameter; the radius, r, is assumed to be l/2. The permeability K can be estimated from the Kozeny-Carman equation K = ϕ 3 D 2 /180(1 − ϕ) 2 (Blunt, 2017). In all cases, we predict the onset of intermittency, as demonstrated in Figure 3, with only a few discrepancies for the bead pack data. Note, however, that in the experiments pressure gradient was not measured, and so the onset of intermittency was determined from an analysis of the pore-scale dynamics as described in the legend of the figure.

Conclusions
We have measured the pressure gradient as a function of flow rate for different fractional flows of oil and brine through a small sample of Bentheimer sandstone. We have observed the Darcy flow regime and the transition to non-Darcy or intermittent flow. At low flow rates, the flow follows a standard linear Darcy law, while for higher flow rates we see ∇P∼Ca a . We proposed a relationship for the threshold capillary number for the onset of intermittent flow, Equations 7 and 8, which accurately matched the experimental results and is applicable to the different viscosity ratio fluids and different rock types, reconciling a large body of experimental results in the literature.
Future work could include the study of flow at different viscosity ratios, system lengths, mixed-wet and oil-wet media, and different types of porous material. We could also study what controls intermittency and the ZHANG ET AL.
10.1029/2020GL090477 5 of 7   (Datta et al., 2014a;Gao et al, 2017Gao et al, , 2019Reynolds et al., 2017;Spurin et al., 2019b) where different viscosity fluids were injected into different rocks; full details can be found in Table 2. In these experiments, the pressure gradient was not measured; instead the onset of intermittent flow was estimated from the pore-scale dynamics as described in the legend. The dashed line is the predicted threshold from Equations 7 and 8. power-law exponent in Equation 4, including the balance of viscous, inertial and capillary forces (Ferrari & Lunati, 2013), and shear-stress at the fluid-fluid interface (Roman et al., 2020;Zarikos et al., 2018).