Pore-scale dynamics and the multiphase Darcy law

Synchrotron x-ray microtomography combined with sensitive pressure differential measurements were used to study ﬂow during steady-state injection of equal volume fractions of two immiscible ﬂuids of similar viscosity through a 57-mm-long porous sandstone sample for a wide range of ﬂow rates. We found three ﬂow regimes. (1) At low capillary numbers, Ca, representing the balance of viscous to capillary forces, the pressure gradient, ∇ P , across the sample was stable and proportional to the ﬂow rate (total Darcy ﬂux) q t (and hence capillary number), conﬁrming the traditional conceptual picture of ﬁxed multiphase ﬂow pathways in porous media. (2) Beyond Ca ∗ ≈ 10 − 6 , pressure ﬂuctuations were observed, while retaining a linear dependence between ﬂow rate and pressure gradient for the same fractional ﬂow. (3) Above a critical value Ca > Ca i ≈ 10 − 5 we observed a power-law dependence with ∇ P ∼ q at with a ≈ 0 . 6 associated with rapid ﬂuctuations of the pressure differential of a magnitude equal to the capillary pressure. At the pore scale a transient or intermittent occupancy of portions of the pore space was captured, where locally ﬂow paths were opened to increase the conductivity of the phases. We quantify the amount of this intermittent ﬂow and identify the onset of rapid pore-space rearrangements as the point when the Darcy law becomes nonlinear. We suggest an empirical form of the multiphase Darcy law applicable for all ﬂow rates, consistent with the experimental results.

fuel cells [5] and filtration, for instance. The saturation and flow-history dependence of the relative permeability and its impact on trapping, which governs the security of carbon dioxide storage and the efficiency of oil recovery, for instance, has been widely studied [6][7][8][9]. However, the link between pore-scale displacement and averaged flow properties is still not fully understood, particularly when, with increasing flow rate, viscous effects become significant.
In this paper we will consider one example of multiphase flow: the simultaneous flow of two fluid phases, oil, o, and water, w, of similar viscosity through a porous rock at steady state as the flow rate is increased. The oil will be the nonwetting phase, with a tendency to occupy the wider pores, while water will be wetting, and will occupy the smaller pores, as well as corners and roughness. We will observe that the average saturation is approximately uniform across the sample. In this case, the macroscopic capillary pressure P c (S w ) = P o − P w will be constant. Since there is flow from high to low pressure, there will be a pressure gradient in each phase. However, since the pressure difference between the phases is fixed, the pressure gradients in the oil and water phases will be the same (∇P o = ∇P w = ∇P); however, the flow rates of each phase and their saturations can be different. We assume that the sample size exceeds a representative elementary volume, REV, so that averaged saturation-dependent macroscopic flow properties can be defined [3]. We also consider fluctuations in pore-scale configurations that occur below the REV. In this case, and ignoring the effects of gravity, from Eq. (1) we can write the total Darcy flux, q t = q o + q w as where λ p = Kk rp /μ p is a mobility and the total mobility λ t = λ o + λ w . We will quantify the relationship between total flux, q t , and pressure gradient, ∇P. We will keep the fractional flow, the ratio of the volume of water injected to the total injected volume, f w = q w /q t , fixed. First, we theoretically examine the effect of flow rate, captured by the balance of viscous to capillary forces, using a capillary number, defined here as Ca = μq t /σ where μ is the average viscosity of the two fluids, σ is the interfacial tension between water and oil in what follows (σ is the interfacial energy per unit area of the fluid-fluid interface). Viscous forces dominate at the pore scale for Ca > Ca v ≈ K/lr, where l is a typical distance between pores and r is the mean radius of a restriction, throat, separating them [10,11]. The term K/lr is a geometric factor which means that the transition from capillary-controlled to viscous-controlled displacement does not occur when Ca ≈ 1, but for a much lower value. We will study Bentheimer sandstone, with K = 1.9 × 10 −12 m 2 , l = 79 μm, and r = 24 μm [ 4]. Hence, Ca v ≈ 10 −3 .F o rC a≪ Ca v capillary forces dominate at the pore scale and it is traditionally considered that viscous effects can be ignored. It is assumed that the fluids occupy fixed pathways, in that the pore-occupancy does not change in steady-state flow, and Eq. (1) is valid.
However, recent advances in pore-scale imaging using x-rays and confocal microscopy, combined with micromodel experiments and dynamic pore-scale simulation, have demonstrated that multiphase flow in porous media experiences a complex dynamics even for Ca < Ca v [12][13][14][15][16][17][18][19][20][21]. Specifically, it has been shown that there is an intermittent flow regime, where locally regions of the pore space are periodically occupied by one phase and then the other, even though averaged over many pores, the fluid saturations remain constant [14,17]. Furthermore, the degree of intermittency has been measured as a function of Ca and the mobility ratio between the fluids, with more intermittency observed for larger Ca and when the injected phase is less viscous [22,23]. However, the pore-scale behavior has not been linked to a quantitative assessment of macroscopic flow properties, and in particular, how the pressure gradient varies with Ca for a fixed fractional flow and how this is linked to the pore-scale dynamics.
Sinha and coworkers have suggested, based on experiments in bead packs [24] and associated simulations [24,25], that for Ca < Ca v there is a power-law scaling of flow rate with pressure drop, P: q t ∼ ( P − P y ) 1/a with an exponent a ≈ 0.5 and a yield pressure for either phase, P y : there is no flow for P P y . While a yield pressure may be needed to mobilize trapped ganglia for Ca > Ca v , Sinha et al. predict a finite threshold for the whole saturation range. This model also implies that at low flow rates the Darcy law, with a linear relationship between pressure gradient and flow rate, is not valid.
This paper aims to study the relationship between flow rate, pressure gradient and pore-scale fluid occupancy using synchrotron x-ray imaging with a high spatial and temporal resolution. We investigate the behavior when an equal volume of oil and water at steady-state is injected through a sandstone for a range of flow rates to determine flow regimes. While only a single rock type and fractional flow is studied, we seek to uncover the evolution of a complex pore-scale dynamics by quantifying time-dependent fluid rearrangements in the pore space. We propose a mathematical model for the rate dependence of the pressure gradient that is consistent with our measurements, bearing in mind that further work is required to test its general applicability.
The structure of the rest of the paper is as follows. We first describe the experimental methods in Sec. II. We then place the results in context in Sec. III A by providing a theoretical analysis for the threshold capillary numbers when viscous forces have an effect, based on a pore-scale energy balance. We then delineate our proposed flow regimes using the measured pressure gradients in Sec. III B. We follow this by a quantitative analysis of pore-scale fluid occupancy provided by our time-resolved imaging in Sec. III C.

II. IMAGING EXPERIMENTS
We imaged steady-state two-phase flow at the Diamond synchrotron in Oxfordshire, UK. The methods and apparatus are described in detail elsewhere [17,26]; using pink-beam radiation with a range of x-ray energies [14] images were acquired every 60 s. High-precision pumps injected both brine, the wetting phase (deionized water with 15 weight % potassium iodide, KI, dissolved to improve image contrast), and oil (n-decane), the nonwetting phase, at a fractional flow of 0.5 through a cylindrical sample of water-wet Bentheimer sandstone of diameter 6.01 mm and length 57.11 mm. The experiments were conducted at ambient temperature, 20 • C, and a fluid pressure of 2 MPa. The measured viscosity of the water was 0.83 ± 0.01 mPas, while n-decane had a viscosity of 0.838 mPas (PubChem, open chemistry database).
The flooding sequence was as follows: (1) A dry scan was taken with 2 MPa confining pressure, which compresses the Viton sleeve outside the core to avoid fluid bypass.
(2) CO 2 was injected into the sample for half an hour to displace air.
(3) The 15 wt% brine was injected to completely saturate the sample: any trapped CO 2 was dissolved. A differential pressure transducer was installed in the tubing (Keller PD-33X). The measurement accuracy was ±250 Pa up to 500 kPa. This pressure differential was continuously monitored with readings taken every 2 s.
(4) The brine-saturated sample was scanned. A back pressure of 2 MPa was set for the whole system. (5) Oil was injected at 2 ml/min for 30 min to reach the initial brine saturation. (6) Then brine and oil were injected for fractional flows f w of 0.15 and 0.3 by keeping the total volumetric flow rate fixed at 0.02 ml/min. Each fractional flow was injected for 90 min.
(7) Brine and oil were injected at an equal flow rate of 0.01 ml/min each: the fractional flow was 0.5. At the same time, the pressure drop across the whole sample was recorded during the experiment. Injection continued for at least two hours to allow the pressure to stabilize. Threedimensional images were taken every minute for the final 30 min of flow.
(8) The total flow rate was increased to 0.04 ml/min, 0.08 ml/min, 0.4 ml/min, 0.8 ml/min, 1.2 ml/min, and 4 ml/min step by step while the fractional flow was kept at 0.5. For each flow rate, we waited 2 h until steady state was reached. Scans were taken every minute for the final 30 min during injection at each flow rate. (9) After cleaning all the tubing and end fittings used in the experiment, steps 5 to 8 were applied to the whole system without the rock sample. The pressure drops of the apparatus itself were recorded at all flow rates. 013801-3 (10) The pressure differences were calculated by subtracting the pressure drops along the tubing measured at step 9 from the total pressure drops measured from steps 7 and 8. The pressure gradients were computed as the average of these pressure differences during the final 30 min at each steadystate divided by the length of the sample, 57.11 mm.
Images were acquired using 1 000 projections with an exposure time of 0.04 s. The total acquisition time for each image was 1 min: 40 s to take the projections and around 20 s to transfer the data from camera to computer and rotate the sample back to the starting position. The number of projections was chosen to optimize the resolution of the final images without significantly increasing the acquisition time. The lowest flow rate corresponded to Ca = 2.1 × 10 −7 , while the highest rate had Ca = 4.2 × 10 −5 .

A. Critical capillary numbers derived from energy balance
To provide a theoretical foundation to identify flow regimes, we first perform an analysis of porescale dynamics based on an energy balance and then compare with experiment. When a nonwetting phase is displaced by a wetting phase (oil by water in this paper), the nonwetting phase can become trapped, meaning that it is completely surrounded by wetting phase [4]. For flow rates below Ca v viscous forces are insufficient to mobilize these trapped ganglia, that is, to push them through the pore space [27]. Instead viscous effects allow interfaces to be formed and then destroyed to facilitate flow, which can be understood from an energy balance [14,28].
Instead of calculating the pressure necessary to force a fluid through a throat, we find the interfacial energy needed to create an interface at the pore scale and relate this to the energy provided by fluid injection. This interfacial energy is approximately σ r 2 (the energy per unit area times a typical area of an interface in a single pore), while the energy added to the system due to the flow PdV is, over a length l, equal to −∇Pφl 4 , where φ is the porosity (the pore volume dV = φl 3 while P =−l∇P). Therefore, the interfacial energy needed to create an interface is balanced by the work done by the flow when Then using ∇P ∼−μq t /K fromfromEq.(2) this occurs at a threshold capillary number The intermittent creation and destruction of interfaces can lead to two effects which we will describe in more detail later, based on an analysis of our pore-scale images. The first one is a rearrangement of fluids in the pore space, which may lead to a shift in saturation, but without changing the flow conductance: we define Ca * as the capillary number at which this occurs. The second effect is to cause a noticeable change in both saturation and conductance, facilitating the flow of both fluids for Ca i Ca * . We hypothesize that Ca i has a value less than or equal to the threshold given by Eq. (4) and hence For Bentheimer sandstone Ca i 10 −4 for φ = 0.2, at least one order of magnitude lower than Ca v due to different geometric scaling factors. Note the extreme sensitivity to l in Eq. (5): we provide a conservative limit on the threshold capillary numbers, but cannot be more quantitative without direct experimental measurements, which we present next.
This interpretation of pore-scale dynamics is related to the concept of a macroscopic capillary number [29]: viscous forces allow the periodic connection and disconnection of clusters of the nonwetting phase. In this work we study steady-state flow and quantify the relationship between flow rate and pressure gradient: A theoretical framework to describe a wider range of phenomena including the reduction in trapped saturation with flow rate and its relationship to displacement cycles is presented elsewhere [30,31].
This analysis has only provided an upper bound on the threshold capillary numbers. Next we present the pressure measurements which we use to identify flow regimes and to quantify the values of Ca * and Ca i . This analysis is then further supported through a quantification of fluid occupancy based on pore-space images.

B. Measurements of pressure gradient
At the lowest flow rates Ca is proportional to the pressure gradient, ∇P,Eq.(1),follo wedbya transition to the intermittent flow regime with a nonlinear relationship between pressure gradient and Ca, ∇P ∼ Ca a , as demonstrated in Fig. 1. In the Darcy regime we consider the five points with the lowest flow rates and find the best-fit straight line on linear axes with an intercept of zero, where R 2 = 0.995. We estimate Ca i = 10 −5 as the value when we first see a deviation from the linear Darcy law, and we fit a power-law to the three points with the highest flow rates to find a = 0.60 ± 0.01.
The measured Ca i is ten times lower than the upper bound given by Eq. (5). However, note the fourth-power scaling with l. As shown later, intermittency allows a significant rearrangement of fluid configurations over a few pore lengths -this regime starts when the viscous dissipation is equivalent to the interfacial energy over approximately 2 pore lengths.
In the capillary-controlled limit (1), the pressure differential shows only tiny deviations of around 0.5 kPa, the sensitivity of the transducer, Fig. 2. However, at higher flow rates we see the onset of pressure fluctuations [16,17,24,32]; their magnitude, around 10 kPa, is consistent with the independently measured capillary pressure [26]. These fluctuations start at around Ca = 10 −6 , below Ca i but initially do not affect the flow conductance. At the highest rates, the fluctuations are smaller: The fluid rearrangements may be too rapid for the pressure transducer to respond. The typical period of these oscillations was around 10 s; however, we cannot capture fluctuations that occur more rapidly than the 2 s sampling time. Since these rearrangements can occur anywhere in the pore space, their frequency is dependent both on the number of pores (and hence the size of the system) and on the period of the individual events themselves.
Based on these measurements our suggested flow regimes are shown in Fig. 1. We assume that transient effects have dissipated [12,32].
(1) For Ca Ca * we have capillary-dominated flow described by Eq. (1): ∇P ∼ Ca. This regime is identified by a constant pressure differential as a function of time, shown in Fig. 2.
(2) For Ca i Ca > Ca * we retain a linear Darcy law, Eq. (1), for fixed fractional flow, but with some rearrangement of fluids in the pore space which may change the saturation. However, there is no significant impact on conductance. We call this the onset of dynamics. In the experiments Ca * ≈ 10 −6 identified as when the first pressure fluctuations are observed, which indicate the presence of pore-scale displacements even at steady state, Fig. 2.
(3) For Ca v > Ca > Ca i we observe similar behavior to that described by Refs. [24,33], which we call intermittent flow, albeit with no yield pressure. ∇P ∼ Ca a with an exponent a = 0.60 ± 0.01 based on the data shown in Fig. 1.C a i is identified when there is first a deviation from a linear relationship between Ca and ∇P: in our experiments this occurs at Ca i ≈ 10 −5 as shown in Fig. 1.
In regime 3 for fixed fractional flow and ignoring gravity we propose that the scaling between total Darcy flux and pressure gradient is written as follows: where ∇P w = ∇P o ≡ ∇P. ∇P i is the pressure gradient at the beginning of the intermittent flow regime when Ca = Ca i : the term in the denominator ensures that Darcy flow, Eq. (2), is obtained smoothly for Ca Ca i . We assume that λ t = λ 0 t , its value in the slow-flow Darcy regimes (1) and (2) for a fixed value of f w . We encapsulate the rate dependence through a nonlinear scaling of the pressure gradient, see Fig. 1, rather than writing the total mobility as a function of flow rate.
For Ca Ca v there is the viscous limit when capillary forces are weak at the pore scale, which has been studied previously [10,24,25,34]. Here again the flow rate is proportional to the pressure gradient and Eq. (2) is obeyed, with q t =−λ v t ∇P where λ v t is the total mobility in the viscous limit. Capillary forces still play a role: for instance, the fluid patterns are sensitive to wettability [35].
However, the emphasis of this paper is to study the flow regimes for which viscous forces are small at the pore scale, Ca < Ca v , yet influence the fluid dynamics.
The exponent a in the scaling law, Eq. (6) is an empirical fit to the measurements. The exponent can be obtained by considering the behavior for the capillary and viscous limits. For Ca < Ca i the total mobility is λ 0 t while for Ca > Ca v the total mobility is λ v t . Then for continuity when Ca = Ca v , and the pressure gradient is ∇P v ,Eqs. (2) and (6)imply which, since in the viscous and capillary limits Ca v ∼ λ v t ∇P v and Ca i ∼ λ 0 t ∇P i respectively, Eq. (7) becomes All the terms in Eq. (8), except a, are fixed by the porous medium geometry, the viscous limit and the Darcy flow configurations. Hence, a is an empirical fit, extrapolating between well-established limits.
While the evidence for regime 3 is based on only a few data points in Fig. 1, it is possible to argue that it has to exist, even if a strict power-law dependence between flow rate and pressure gradient is not observed. This is because the total mobility in the capillary limit is lower than in the viscous limit: λ v t >λ 0 t [10,36]. Hence, the total mobility has to vary between these two values as flow rate is increased which requires a nonlinear scaling between Ca and ∇P.
The flow regimes displayed in Fig. 1 are supported not just by the pressure measurements, but by direct pore-space imaging of the displacement, shown next. However, these are hypothesized based on simple scaling arguments and energy balance: it is likely that if the full range of flow conditions, namely, different fractional flows, viscosity ratios, rocks types, and wettability, an even richer behavior may emerge. We cannot, for instance, in these experiments determine the saturation dependence of fractional flow and relative permeability.

C. Pore-scale imaging
Three-dimensional images with 10 9 voxels of size 5.2 μm were acquired every minute. The final 30 images at each flow rate were used for quantitative analysis. These tomograms were segmented to show the solid, brine, and oil phases using differential imaging; the details are provided elsewhere [17,37]. The saturation averaged in slices perpendicular to the flow direction was constant in all cases: we saw no capillary end effect.
At the lowest flow rates, the fluid configurations remained fixed in time, even when the flow rate was increased. As illustrated in Fig. 3, the oil may follow tortuous paths through the largest regions of the void space. In the intermittent regime, there is sufficient energy to allow the oil to short circuit some of these pathways to decrease the overall flow resistance. However, these connections are not fixed. Oil moves, not as disconnected ganglia [38] but along pathways that intermittently connect and reconnect at critical junctions, like cars controlled by traffic lights [14].
From the final 30 images at each flow rate the average saturation was recorded. The results are showninT ableI together with the observed flow regimes. At the lowest flow rates, under capillary control, regime (1), the saturation remains constant, within the measurement uncertainty, while the pressure trace, see Fig. 2, was constant. Then we observe a shift in both saturation and the onset of dynamic effects, regime (2), with pressure fluctuations and local rearrangements of the pore space, but without a deviation from Darcy flow. At the highest flow rates the saturation change is more significant, indicating the intermittent flow regime (3).
We measured the fraction of the pore space in which the fluids rearranged. We observed two types of fluctuation: type 1, is when we see a change in the phase occupying a voxel from one scan to the next, representing variations in configuration that occur over a minute or more. Type 2 is when the fluctuations are sufficiently rapid that they occur within 1 minute: we cannot explicitly see the changes; instead a voxel has an intermediate gray-scale value between water and oil. The fraction of the number of void-space voxels that are type 1 or type 2 are shown as a function of capillary number in Fig. 4(a). For type 1 we measure the fraction of the voxels that change occupancy from oil to water or vice versa between two consecutive scans, averaged over the final 30 scans at each flow rate. For type 2 we record the fraction of intermediate gray-scale voxels. The error bars indicate the standard deviations in these measurements from image to image.  4. (a) The fraction of void-space voxels for which the fluid occupancy changes from one scan to the next (type 1, left axis) and during a scan (type 2, right axis) as a function of capillary number Ca. The scans were taken every 60 s. The values are averaged over the final 30 images at each Ca. The error bars indicate the standard deviation. Regime 1 has no rapid oscillations (type 2) and very few type 1 changes. Regime 2 again has no type 2 intermittency, but more sustained low-frequency rearrangements, while in regime 3 we see significant fluctuations of both types, leading to a change in flow conductance. (b) The fraction of void-space voxels for which the fluid occupancy changes from one scan to the next (type 1) as a function of time during the final half hour of the scans. In regime 1 there are periods with no changes in fluid configuration, while in regime 2 there are always some low-frequency periodic rearrangements.
In regime 1, capillary controlled, we see no rapid type 2 fluctuations and only occasional, low-frequency alterations, as shown in Fig. 4(b) which shows the fluid occupancy as a function of time during the final 30 min of the scans; in fact, there are periods when there is no change in the pore-scale occupancy at all. The fluids effectively move through fixed pathways. Regime 2 is characterized by more frequent type 1 rearrangements that are seen between every scan, but there are no high-frequency fluctuations and there is no apparent impact on the linearity of the Darcy law. In regime 3, rapid, type 2, changes in occupancy occur. These are sufficiently extensive to alter the flow conductance.

IV. DISCUSSION AND CONCLUSIONS
We have studied the relationship between flow rate (total Darcy flux) and pressure gradient for the simultaneous flow of oil and water of similar viscosity through a small sandstone sample. We have also simultaneously used fast synchrotron imaging to study the fluid configurations at the pore scale.
We observed three flow regimes.
(1) The first regime is capillary-controlled flow where, for capillary numbers Ca < 10 −6 in these experiments, our results are consistent with the traditional multiphase Darcy law, showing a linear relationship between flow rate and pressure gradient where the phases follow fixed pathways. At the pore scale we see only rare changes in fluid configuration with significant time periods with no rearrangements at all.
(2) The second regime marks the onset of dynamics, where at a capillary number, Ca * ≈ 10 −6 , fluctuations in local pressure are observed with an overall shift to lower oil saturation at constant fractional flow, while the total mobility is approximately unchanged. This is accompanied, at the pore scale, by low-frequency changes in fluid configuration in a small fraction of the pore space.
(3) At a critical Ca i ≈ 10 −5 the third regime, intermittent flow, leads to an increased flow conductance and the onset of a nonlinear relationship between flow rate and pressure gradient. Here we see rapid, subminute, changes in fluid configuration in a significant fraction of the pore space. The results were fitted to a power-law relationship between flow rate and pressure gradient.
We suggest that the same flow regimes are seen for different porous systems, fluid pairs, and over the whole range of fractional flow, albeit with different threshold capillary numbers and exponents, a. We do not expect to see a threshold pressure gradient for flow at low flow rates if the fluids are connected.
Future work could study microporous carbonates [39], different fluid viscosity ratios [22,23], the full range of fractional flows and mixed-wet media [39][40][41][42], where intermittent flow may be more significant than in the water-wet rock studied here [43,44]. Further experimental work could also determine the saturation and rate dependence of relative permeability for the different flow regimes outlined here. We propose that accurate pressure measurements during steady-state flow, combined with pore-scale imaging, allows a detailed characterization of multiphase flow under the conditions encountered in the subsurface [26,45,46].