The development of intermittent multiphase fluid flow pathways through a porous rock

Intermittent fluid flow has recently been identified as an important transport mode for subsurface multiphase flow systems such as CO2 storage and natural gas production. However, due to experimental limitations, it has not been possible to identify why intermittency occurs at subsurface conditions and what the implications are for upscaled flow properties such as relative permeability. We address these questions with observations of nitrogen and brine flowing at steady-state through a carbonate rock. We overcome previous imaging limitations with high-speed (1s resolution), synchrotron-based X-ray micro-computed tomography combined with pressure measurements recorded while controlling fluid flux. We observe that intermittent fluid transport allows the non-wetting phase to flow through a more ramified network of pores, which would not be possible with connected pathway flow alone for the same flow rate. The volume of fluid intermittently fluctuating increases with capillary number, with the corresponding expansion of the flow network minimising the role of inertial forces in controlling flow even as the flow rate increases. Intermittent pathway flow sits energetically between transport through connected pathway flow and turbulent flow. While a more ramified flow network favours lowered relative permeability, intermittency is more dissipative than connected pathway flow, and the relative permeability remains unchanged for low capillary numbers where the pore geometry controls the location of intermittency. However, as the capillary number increases further, the role of pore structure in controlling intermittency decreases which corresponds to an increase in relative permeability. These observations can serve as the basis of a model for the causal links between intermittent fluid flow, fluid distribution throughout the pore space, and the upscaled manifestation in relative permeability. Preprint submitted to Elsevier November 18, 2020


Introduction
Accurate prediction of the propagation and trapping of multiple fluid phases in porous rocks of the Earth's subsurface is important to a number of engineered and geological applications; fluid flow is simulated to quantify the volume of CO 2 that can be stored safely underground (1; 2; 3) and to estimate the volume of hydrocarbon fluids in an oil field (4). These systems are modelled using a continuum framework, including Darcy's law extended to multiphase flow, to describe energy dissipation as fluid moves through the pore network of subsurface rocks. Darcy's law is predicated on the transport of the individual fluid phases through distinct pore networks, with static fluid-fluid interfaces. However, recent observations have demonstrated that multiphase flow through the pore networks of rocks is more complex, with modes of transport through both connected paths as well as the frequent rearrangement, or intermittency, of these flow networks (5; 6; 7; 8; 9).
Distinct modes of transport occur depending on the capillary number, the ratio of capillary to viscous forces, and the ratio of fluid viscosities, M = µ nw /µ w (6). At capillary numbers prevailing in the subsurface and with fluid pairs including medium crude oil-brine or CO 2 -brine that have less than unity viscosity ratios (M ≤ 1) there are two dominant modes of transport. Either fluid phases flow through connected pathways, with fixed pore networks for each fluid, or they flow through intermittent pathways, with the competition between the fluids resulting in pathways that are periodically broken and reformed in critical locations of the pore space. Intermittent pathway flow has been compared to cars at traffic lights; with the pore network being equivalent to the road network. Flow through the pore network, like traffic along a road, is not continuous (7). These intermittent flow pathways persist even under steady-state conditions whereby the macro-scale flow properties such as saturation and pressure are invariant when averaged over time (10).
The various features of the pore structure which control the particular locations of intermittency have been described in previous work (5). The intermittency occurs in intermediate sized pores, where the capillary pressure of the fluids is close to the entry pressure of the pore. Intermittency tends to occur in places where either the non-wetting phase or the pore space is poorly connected. Small variations in the initial invasion of a fluid phase can lead to variations in the location of intermittent pathways, but the characteristics of the intermittent pores remain the same. Observations also show that when a system is perturbed from steady-state, the oscillations in the pressure signal caused by intermittent pathway flow continue as the system transitions to a new steady-state (11).
However, due to experimental limitations, it has not been possible to identify what drives intermittency at the low capillary numbers characteristic of subsurface conditions, and how this impacts larger scale flow properties such as the relative permeability. In this paper we address these questions with observations of steady-state flow of nitrogen and brine through a porous carbonate rock. We observed intermittency in steady-state flow in unprecedented detail with high temporal resolution (1 s) synchrotron X-ray tomography for nine capillary numbers, with total flow rate and fractional flow varied between experiments. We revealed a previously unidentified scaling between capillary number and fluid flux through intermittent flow paths. We identified that as flow through certain areas of the pore space approached critical Reynolds numbers for the transition at which inertial forces are no longer negligible, more intermittent pathways over a more disperse area occurred, reducing the velocity of flow in a given region of the pore space. These observations can provide the basis for a model for the causal links between intermittent fluid flow, fluid distribution throughout the pore space, and the upscaled manifestation in relative permeability.

Methods
We observe the fluid distribution of nitrogen and brine with time as both are simultaneously injected into a carbonate rock. A total of 9 observations were collected for 3 different total flow rates. For each total flow rate, the fractional flow (which is the fraction of the total flow rate constituted by the brine flow rate) was changed 3 times. The list of experiments is given in Table 1. We use a capillary number which comprises properties of both fluid phases: where q is the total Darcy flux, σ the interfacial tension between the nitrogen and brine, and γ the total mobility (γ = fw µw + 1−fw µnw where µ is the viscosity and f w is the brine flow rate as a fraction of the total flow rate) ( (11). The processing is described in more detail in (11), but is repeated for completeness below, with Experiments 4-9 conducted with the same methodology. The observations were made on the same cylindrical Estaillades sample, 5 mm diameter and 20 mm length. The sample was initially saturated with brine (deionised water doped with 15% wt. KI to improve the X-ray contrast). The system was pressurised to 8 MPa to minimise compressibility effects, with an additional 2 MPa of confining pressure. Then both nitrogen and brine were injected simultaneously, with the flow rates of the nitrogen and the brine listed in Table 1. A differential pressure transducer (Keller, 300 kPa transducer with 1.5 kPa accuracy) was connected across the sample, to measure the pressure drop. Steady-state was identified as the point after which the differential pressure had plateaued. Between observations the sample was not resaturated with brine.
The X-ray imaging was conducted at the TOMCAT beamline at the Swiss Light Source, Paul Scherrer Institut, Villigen, Switzerland. The sample was exposed to filtered polychromatic X-ray radiation with a peak energy of about 26 keV originating from a 2.9 T bending magnet source. The filter was 2.3 mm-thick Silicon and 5 mm thick pyrolytic carbon. An in-house developed GigaFRoST camera (12) and a high numerical aperture white-beam microscope (Optique Peter) with 4x magnification (13) were used, yielding an effective pixel size of 2.75 µm. Each tomogram contained 1000 projections over 180 • rotation. Each scan lasted 1 s, with a further 1 s required for the sample to rotate back to its initial position before the next scan started.
Each image analysed was 4224 × 4298 × 4263 µm in size. The images were reconstructed from the X-ray projections using the propagation-based phase contrast method (14) and the gridrec algorithm (15). The images were then filtered with a non-local means filter to suppress noise whilst maintaining the information at phase boundaries (16). The first image was taken with just deionised water in the pore space. This image is used to segment the pore space from the rock grains using a watershed segmentation algorithm (17). Then the sample was saturated with the brine, and another image was taken; this is the brine saturated imaged. All subsequent images with the non-wetting phase (nitrogen) and brine present were subtracted from the brine saturated image; this results in a differential image whereby only the location of non-wetting phase remains. From this a simple greyscale value threshold was used to segment out the nitrogen. The pore space was overlain on this segmentation to locate the pore space occupied with brine.
The pore space morphology was extracted using a maximal inscribed spheres (maximum ball, MB) network extraction technique (18; 19). The fluid occupancy for each pore MB was assigned for every time step by overlaying the MBs on the segmented fluid distribution, and assigning the MB the fluid with majority occupancy.

The stability of intermittent fluid flow pathways
3.1. The fluid distribution and dynamics change to accommodate increasing fluid flux A total of 9 steady-state experiments were conducted across a capillary number range spanning two orders of magnitude (see Table 1). The average gas saturation for these 9 experiments was between 0.26 and 0.33, a difference of 27%. The average differential pressure across the sample increased from 48 kPa to 206 kPa, an increase of 429%.
The gas saturation remains relatively constant throughout the experiments, suggesting a similar capillary pressure for the observations across the range of capillary numbers. In contrast, the fluid phase distribution and the mode of transport for the gas phase varies considerably ( Figure 1). For experiment 1 (top left panel of Figure 1) the gas is highly channeled and is connected by a single pore in the center of the rock sample. This is also the case for experiments 2 and 3. However, for the group with the next highest total flow rate (experiments 4-6, central column in Figure 1) the gas distribution is more uniform throughout the pore space. The gas phase is more evenly distributed at the highest total flow rate (experiments 7-9, right hand column in Figure 1). As the total flux of the fluids increased, transport across the porous medium also changed. As we will show in subsequent sections, this corresponds to an increase in intermittency as a fluid transport mode. The brine and rock space are transparent in these images. Each distinct gas ganglion has been assigned a different colour. These images are for a single scan and although the ganglia do connect and disconnect with time, these images are representative for the connectivity.

Intermittent flow pathways develop in response to pathways with high
Re flow The highly channeled flow in the top left hand panel of Figure 1, whereby gas flow is forced through a single pore provides the opportunity for the Reynolds number to be calculated for that constriction.
The Reynolds number can be defined by: where ρ is fluid density, u is fluid velocity, D is the hydraulic diameter and µ is the fluid viscosity. The hydraulic diameter is calculated by: where A is the cross-sectional area of the flow path and W p is the wetted perimeter (the perimeter of the cross-sectional area of the gas). The fluid velocity u (m/s), is calculated from the injected volumetric flow rate v (m 3 /s) by Thus the Reynolds number can be rewritten as: For the experiments conducted at the lowest total flow rate (the left hand panel in Figure 1), the gas flow is confined to a single constriction at one point in the pore space. Here, the Reynolds number is calculated to be 27 ≤ Re ≤ 30 ( Table 2). The other experiments were conducted at higher flow rates, but there were more pathways available for gas to flow across the pore space (the middle and right hand panels in Figure 1). The fluid velocity cannot be as precisely constrained with many avenues for flow. To provide a comparison we calculate the wetted perimeter for the same slice and assume the gas flow in uniform across all pathways. Despite the higher total fluid flux for these experiments, the Reynolds numbers are 1-2 orders of magnitude smaller due to the altered fluid distribution, Re ≤ 1 ( Table 2).
This variation in Reynolds number is notable due to its proximity with the critical Reynolds number for flow in porous media. The critical Reynolds   The fluid redistribution in response to the increased total flux of experiments 4-9 reduces the influence of inertial forces, even as flow rate increases. As will be reported below, this also coincides with an increased flux through intermittent flow pathways. This suggests that while intermittent fluid flow may not be energetically favoured over laminar connected pathway flow, it provides more avenues for fluid flow which reduce the role of inertial forces in a given region of the pore space.

Intermittency scales with capillary number
The amount of intermittent fluid flow in the pore space has typically been quantified by the volume of pore space identified as intermittently occupied during imaging (5; 6; 9; 24; 8). A positive correlation exists between the percentage of the pore space intermittently occupied and the capillary number (shown in Figure 2a) for a given fractional flow. However, the trend does not scale across fractional flows. Thus, the percentage volume of the pore space intermittently occupied provides some insight into the dynamics, but it does not capture the full complexity of the dynamics.
With high temporal resolution synchrotron imaging, we were able to capture the movement of fluid-fluid interfaces. As a proxy for fluid flux through intermittent flow, we quantify the volume of gas fluctuating as the volume of an intermittently occupied pore multiplied by the number of times a change in occupancy occurs during an imaging window, averaged to the number of scans during the imaging window. The capillary number plotted against the average volume of gas fluctuating per scan is shown in Figure 2b. The inclusion of the number of fluctuations scales the points across all fractional flows towards a single curve in Figure 2b. The volume flux of fluid transport through intermittent flow pathways scales directly with capillary number.

The evolution of intermittent fluid flow pathways
4.1. Intermittent flow path evolution with fractional flow and capillary number By comparing observations between experimental stages, we tracked the evolution of flow pathways in pores, e.g, from intermittent to stable. At a given total flow rate, the majority of pore space that was intermittently occupied either remains intermittent or becomes a stable gas pathway as the fractional flow of gas increases, Figure 3a. Increasing the relative flow rate of the nitrogen induces an increase in the capillary pressure, allowing some regions of intermittent pathway flow to develop into stable flow pathways. The increased capillary pressure is also reflected in the trends towards the invasion of smaller pores with increased fractional flow, experimental number groups 1-3, 4-6, 7-9 in Figure 3b. In all cases, the pores with intermittent fluid occupancy are on average smaller than pores with stable non-wetting flow pathways.

The role of pore geometry in the development of intermittent pathways
Intermittent pathways occur in the intermediate sized pores, where the capillary pressure is close to the entry pressure required for non wetting phase invasion into the pore. The two phases compete for occupancy leading to intermittent occupation. It has been previously shown that intermittency is also more likely to occur in regions of the pore space that are poorly connected (5).
We observe that the amount of time gas spends in an intermittent pore is positively correlated with the pore size, with larger pores occupied with gas for longer (Figure 4a). The trend weakens as the total flow rate of the experiment increases. This shows a reduction of the role of pore size in controlling the intermittent events as the capillary number increases; pore geometry is less important at higher capillary numbers.
The declining role of pore structure with increasing capillary number is also reflected in the frequency of fluctuation of fluid phases within the pores. The role of pore radius in the average number of fluctuations per scan is shown in Figure 4b. For the lowest total flow rate (Experiments 1-3), the smaller pores fluctuate more. This trend weakens for the higher total flow rates, again showing that the pore structure is less of a control on intermittency at higher flow rates.

The impact of intermittency on the relative permeability
According to the two-phase Darcy equations, the relative permeability of the non-wetting phase is defined by: whereby q nw is the gas flow rate per unit area (m/s) calculated by dividing the non-wetting phase volumetric flow rate by the cross sectional area of the cylindrical sample, µ nw is the gas viscosity (Pa· s), L is the length of the sample (m), κ abs is the absolute permeability calculated from single phase brine flow through the sample (m 2 ) and ∆P is the pressure difference across the sample (Pa).
The gas relative permeability for all experiments is low due to the gas saturation being close to the percolation threshold (25) and the constricted nature of the gas pathway/s (see Figure 1). There is a correlation between the gas saturation and the gas relative permeability (shown in Figure 5) but the saturation covered in these experiments is small, while the range in relative permeability measured is large. Figure 5: The saturation of gas versus gas relative permeability. Note the "error bars" are in fact the maximum and minimum relative permeability or saturation observed during a 10 minute imaging window.
To illustrate the impact of increasing capillary number, we fit two Brooks-Corey models to the relative permeability observations using: While the individual parameters can correspond to physical properties of the rock material, we simply use the model empirically and have not made any attempt, e.g., to measure the irreducible water saturation (26; 27; 28). We hold all parameters constant except for the use of two different values of λ to model our data points. The relative permeability relationships appears to shift, increasing with increasing capillary number, although not in an obviously linear fashion (it remains relatively constant for the first two groups of experiments, and finally increases at the highest capillary number). It is notable that the shift can be captured with a variation in the parameter λ. This parameter reflects the controls of the pore structure on relative permeability. A decreasing value of λ with increasing capillary number suggests that intermittency leads to a weakening of the controls of the pore structure, i.e., the pore space is becoming effectively more uniform from the perspective of the flow process. This is consistent with the observations that intermittency permits fluid phases to access more of the pore space. It also suggests potential for a scaling of relative permeability with capillary number on the basis of the role of intermittency in controlling fluid distribution.

Conclusions
We have identified the links between fluid flow potential, transport through intermittent flow pathways, and fluid phase distribution within the pore space. We have also provided preliminary observations of the upscaled impact of these varying modes of transport on the continuum property relative permeability.
During transport through connected fluid pathways, the fluid phase distribution throughout the pore space is controlled by capillary pressure; nonwetting fluid cannot enter regions of the pore space without sufficiently high capillary pressure to overcome barriers to entry.
Fluid transport through intermittent pathway flow provides access to regions of the pore space previously inaccessible to the connected non-wetting phase due to entry pressure barriers. In this work we identify the development of intermittent pathways as we increase the total flow rate at fixed capillary pressure (saturation and fractional flow). As the capillary number is increased intermittent flow develops throughout a distributed network in the pore space. This happens as the flow velocity in the network of connected paths reaches the critical Reynolds number (where inertial forces are no longer negligible) in specific locations. Intermittent transport is evidently more energetically favourable than turbulent flow, and it provides a mechanisms to reduce the maximum flow velocity in the pore space through the generation of a more disperse flow network.
The volume of fluid fluctuating through intermittent flow pathways scales with capillary number irrespective of the fractional flow of the system. The particular scaling relationship is reflective of the characteristics of the pore structure. It shows an increasing trend with Capillary number towards transport through intermittent flow paths whereas increasing transport through connected flow paths could result in flow paths with high local Reynolds numbers and consequently increased viscous dissipation The fluid distribution throughout the pore space is heavily influenced by the transport of fluid through intermittent versus connected flow paths. It is notable that the relative permeability can be robust to changes in this underlying fluid distribution; two groups of relative permeability data with distinct underlying fluid distributions fall on the same curve. This probably reflects the energetic expense of intermittency: while intermittency permits a more disperse network of flow pathways, favouring an increase in relative permeability, energy is consumed in the creation and destruction of interfaces during this transport mechanism. Thus the relative permeability can remain almost the same in the initial transitions from pure connected pathway flow to intermittent flow.
At the highest capillary numbers, and greatest flux through intermittent pathway flow, we do see an increase in the relative permeability saturation relationship. This coincides with a weakening of the role of the pore structure in controlling locations where intermittency occurs. Whereas the location of intermittent transport is strongly correlated with pore size at low capillary number, intermittent transport occurs across the entire range of pore sizes at high capillary number. The corresponding shift in the relative permeability curve is consistent with a weakening of the role of pore structure.
Combined, these observations provide the basis for a model for the causal links between intermittent fluid flow, fluid distribution throughout the pore space, and the upscaled manifestation in relative permeability. It is possible that the locations of intermittent flow may be predicted a priori from knowledge of the pore structure, the capillary pressure, and the average fluid flow velocities; where turbulent flow is approached in a given pore, intermittency may develop in locations where the entry pressure is above the average capillary pressure in the system. Energy dissipation associated with distinct modes of transport, i.e., intermittent and connected pathway flow, may be summed to estimate the impact on upscaled relative permeability.

Acknowledgments
We acknowledge the Paul Scherrer Institut, Villigen, Switzerland for provision of synchrotron radiation beamtime at the TOMCAT beamline X02DA of the SLS. Catherine Spurin is grateful for her funding through the President's PhD Scholarship, Imperial College London. Tom Bultreys acknowledges the Research Foundation-Flanders (FWO) for his post-doctoral fellowship 12X0919N. We are grateful to all our colleagues within the Shell Digital Rocks Programme for their useful discussions and support. Vladimir Novak acknowledges funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sk lodowska-Curie Grant Agreement No 701647. We acknowledge John Spurin and Jane Spurin for transporting experimental equipment and Alessio Scanziani for help during the experiments.