Connected Pathway Relative Permeability from Pore-Scale Imaging of Imbibition

Pore-scale images obtained from a synchrotron-based X-ray computed micro-tomography (CT) imbibition experiment in sandstone rock were used to conduct Navier–Stokes flow simulations on the connected pathways of water and oil phases. The resulting relative permeability was compared with steady-state Darcy-scale imbibition experiments on 5 cm large twin samples from the same outcrop sandstone material. While the relative permeability curves display a large degree of similarity, the endpoint saturations for the CT data are 10% in saturation units higher than the experimental data. However, the two datasets match well when normalizing to the mobile saturation range. The agreement is particularly good at low water saturations, where the oil is predominantly connected. Apart from different saturation endpoints, in this particular experiment where connected pathway flow dominates, the discrepancies between pore-scale connected pathway flow simulations and Darcy-scale steady-state data are minor overall and have very little impact on fractional flow. The results also indicate that if the pore-scale fluid distributions are available and the amount of disconnected non-wetting phase is low, quasi-static flow simulations may be sufficient to compute relative permeability. When pore-scale fluid distributions are not available, fluid distributions can be obtained from a morphological approach, which approximates capillary-dominated displacement. The relative permeability obtained from the morphological approach compare well to drainage steady state whereas major discrepancies to the imbibition steady-state experimental data are observed. The morphological approach does not represent the imbibition process very well and experimental data for the spatial arrangement of the phases are required. Presumably for modelling imbibition relative permeability an approach is needed that captures moving liquid-liquid interfaces, which requires viscous and capillary forces simultaneously.


A C C E P T E D M A N U S C R I P T 1 Introduction
For the description of multiphase flow in porous media on the Darcy scale, commonly the (phenomenological) extension of Darcy's law to multiphase is employed. The concept of relative permeability is introduced (Dullien, 1979;Corey, 1954) to account for the presence of other immiscible fluid phases (Niessner, 2011). Together with the standard parameters porosity , (absolute) permeability K, fluid viscosities  i , and the capillary pressure p c vs saturation S w function, the relative permeability k r is one of the key parameters to model the flow of multiple phases. For instance, in reservoir engineering, where on the field scale viscous forces dominate over capillary forces (Hilfer, 1996), the relative permeability vs saturation function together with absolute permeability has a large influence on the flux of fluid phases. Relative permeability can largely vary depending on rock type, wettability and other parameters, which can have a major impact on oil recovery by water flooding. At the same time, relative permeability is difficult and is experimentally time-consuming to obtain.
Since the multiphase extension of Darcy's law is phenomenological it does not within its own boundaries predict relative permeability. Relative permeability is traditionally determined experimentally using either the steady state (Dake, 1978;Kokkedee, 1996) or non-steady state methods (Johnson, 1959;Berg, 2013b), which are time-consuming and costly. Therefore within the last decade there has been increasing interest in predicting relative permeability by pore-scale simulation (Blunt, 2002), e.g. using pore network modelling or direct simulation techniques (Ferrari, 2013, Koroteev, 2014Landry, 2014;Gray, 2015). Generally these approaches can be categorized into pore network models (Lenormand, 1988Raoof, 2013Ghanbarzadeh, 2014) and direct approaches depending on whether the modelling is performed on a network model abstracted from the pore space or directly on the pore space obtained from X-ray computed micro-tomography (CT) (Prodanovic, 2006;Porter, 2009;Landry, 2011;Wildenschild, 2013;Landry, 2014;Manuel, 2014;Fusseis, 2014).
Pore-scale modelling approaches can be further categorized into quasi-static and dynamic approaches where capillary and viscous forces act at the same time (Leverett, 1941). Combinations of all categories exist, i.e. there are quasi-static (Blunt, 2002) and dynamic pore network models (Joekar-Niasar, 2012;Hammond, 2012;Prodanovic, 2014), but also quasi-static (Becker, 2008) and dynamic direct simulation methods (Coles, 1998;Jaqmin, 1999;Raeini, 2014a,b;Ferrari, 2013;Ramstad, 2009, Vogel, 2005Demianov, 2011;Boek, 2014;Gray, 2015). In practice, the most common approaches are quasi-static pore network models, which are computationally less costly than dynamic and direct approaches due to simplifications made. However, these simplifications introduce uncertainty. First, the abstracted pore networks are often subject to tuning on reference data (e.g. capillary-pressure data), which decreases a direct link between the actual pore space and the abstracted network (Raeesi, 2013). Elementary pore-scale displacement events (Lenormand, 1988) are captured only in a mechanistic way. Then the capillary-viscous flow is decomposed into capillary-dominated steps determining the pore-scale saturation distribution and viscous-dominated steps where the hydraulic conductivity of connected phases are determined. The main limitation of quasi-static pore network modelling is that only the connected pathway flow is captured. However, multiphase flow in porous media consists of connected flow and ganglion dynamics (Youssef, 2014; caused by the motion of disconnected fluid phases (Mason, 1983;Yadav, 1983;Datta, 2014). In reality there is mass exchange between connected and disconnected fluid phases (Hilfer, 1998;Hilfer, 2006a;Hilfer, 2006b), which not only adds to the net flux but also influences the time evolution of the pore-scale fluid distributions, including the configuration of connected phases.
Dynamic pore network models are more sophisticated and attempt to capture dynamics effects like capillary-viscous displacement including ganglion dynamics correctly (Joekar-Niasar, 2012;Hammond, 2012), but are more expensive from a computational perspective. Direct dynamic simulation approaches based on the lattice Boltzmann technique (Coles, 1998;Ramstad, 2009; A C C E P T E D M A N U S C R I P T 4 Chukwudozie, 2013; Boek, 2014;Gray, 2015) and Navier-Stokes approaches (Jacqmin, 1998;Demianov, 2011;Raeini, 2014a, Raeini, 2014bFerrari, 2013) are becoming more popular because no a priori choice between level of rigor and captured phenomena are made. Like other dynamic approaches, capillary and viscous forces act at the same time and, depending on the choice of flow parameters, both capillary and viscous dominated flow are captured (Raeini, 2014b). That allows for the description of a wide range of flow regimes and simulation of a wide range of pore-scale dynamics such as cooperative and/or non-local displacement processes (Holtzman, 2015). Such nonlocal processes have been observed during drainage in 2D micromodels (Armstrong, 2013; and during imbibition  in with fast CT (Berg, 2013a;Berg, 2013b;Youssef, 2014). However, from a computational perspective they are much more expensive than pore network models.
Along with the choice of the size of the computational domain and associated length scale, often implicitly a choice on the underlying flow regime is also made. Avraam and Payatakes (Avraam, 1995) demonstrated in 2D micromodels during fractional flow that a ganglion dynamics flow regime exists in combination with connected pathway flow. For ganglion dynamics the non-wetting phase is transported in the form of non-percolating clusters (Youssef, 2014). These move through the pore space in a sequence of breakup and coalescence processes , which requires direct numerical simulation to capture (Demianov, 2011;Gray, 2015). Clusters are ultimately mobilized directly or indirectly by viscous forces. Direct mobilization occurs e.g. during "capillary de-saturation" (Mason, 1983;Yadav, 1983;Lake, 1989) where the viscous pressure drop caused by the average wetting phase flow field over non-wetting phase clusters exceeds the capillary trapping forces. For field-relevant flow rates such non-wetting phase clusters can range from a few micrometers (Iglauer, 2012) to centimeters , which easily exceeds the simulation domain size of direct dynamic simulations. Therefore in a respective numerical experiment aiming to determine relative permeability, the contribution of large non-wetting phase clusters to the overall non-wetting phase flux is commonly dismissed.
Mobilization of non-wetting phase clusters can occur in a more subtle way, which has been discovered in a recent study . Cluster dynamics can also occur in capillary-dominated flow, i.e. where the viscous forces associated with the average flow field are not sufficiently large enough to overcome capillary forces and thus cause mobilization. However, the local flow field velocities associated with pore-scale displacement events can exceed the average flow field velocities by orders of magnitude. Haines jumps and snap-off processes are relatively fast events (Armstrong, 2013) reaching Reynolds numbers near unity . These events range over many pores, indicating that besides capillary and viscous forces also inertial forces become relevant (Ferrari, 2013) to describe the pore-scale displacement processes. In this way, snap-off processes can trigger coalescence in distant pores connected to the same non-wetting phase cluster , which adds to the viscous flow through connected pathways for the non-wetting phase flux. These processes occur on much shorter length scales than the one defining the balance of viscous and capillary forces . The associated processes are typically sufficiently captured by direct dynamic simulators using relatively small computational domain sizes (Koroteev, 2014;Raeini, 2014a,b).
Our hypothesis is that if all of the relevant length and time scales associated with the possible flow regimes are captured correctly, then a 1-to-1 match between experiment and simulation would become possible. However, for the computation of relative permeability the question arises of whether there is a significant difference in practice. Large clusters mobilized by viscous forces do contribute to the non-wetting phase flux, no doubt. But if they are infrequent (as percolation models suggest, see Iglauer, 2011 andGeorgiadis, 2013), then their contribution may become negligible. The same holds true for the cluster dynamics caused by non-local and cooperative displacement on smaller length scales. To answer this question, we compared Darcy scale steady-state relative permeability to the relative permeability computed  from pore-scale fluid A C C E P T E D M A N U S C R I P T 5 distributions obtained from a synchrotron-based CT imbibition experiment . This approach was already suggested by (Turner, 2004;Hussain, 2014), who performed quasistatic drainage CT experiments, where flow was stopped for imaging, and computed relative permeability from the connected pathway fraction of the resulting fluid distributions using a Stokes flow solver. We advanced this idea in several ways. First, we used pore-scale fluid distributions from a dynamic synchrotron-based CT imbibition experiment . Dynamic means that fluid distributions were imaged under flowing conditions, i.e. the visco-capillary balance was maintained. Second, while in previous studies the approach was used for drainage processes, we performed an imbibition study on sandstone rock. From many pore network modelling studies it is known that imbibition processes are less approachable by quasi-static modelling approaches than drainage processes . While quasi-static approaches consider only the connected pathway contribution, in reality snap-off (Raoof, 1970) and other mechanisms cause the non-wetting phase to be at least partially disconnected. Therefore, connected pathway flow and ganglion dynamics may co-exist . For validation of the pore-scale modelling, we then compared steady-state relative permeability measurements conducted on a larger-dimension twin sample from the same block to the pore-scale simulation data from the CT experiment. Lastly, to assess whether quasi-static modelling is also sufficient to describe the pore-scale fluid distribution, i.e. for situations where pore-scale fluid distributions are a priori not available, a morphological approach (Becker, 2008;Hilpert and Miller, 2001;Vogel, 2005;Silin, 2010) was used to simulate the capillary-dominated invasion process.

Porous Rock
Gildehauser sandstone (Dubelaar, 2015) was used, which has an average porosity of 20% and a permeability of 1.5 ± 0.3 D. It has low clay content and a narrow pore-size distribution (<35 m pore throat diameter), which makes the rock well suitable for CT imaging . Also, the relatively high permeability and low capillary entry pressure make the rock suitable for saturation with a lowpressure flow apparatus.

Pore-Scale Flow experiments and Pore-Scale Imaging
Cylindrical Gildehauser rock samples of 4 mm diameter (which is larger than the representative elementary volume -REV, see Figure A1 in Appendix A) and 10 mm length were embedded into a cylindrical polycarbonate flow cell by heat-shrinking. The assembly was mounted on a micro-pump assembly, which includes a micro-piston pump for fluid injection. More details about the setup are described in (Berg, 2013a). N-decane was used as non-wetting phase and brine doped with cesium chloride (10 wt-%) as the wetting phase. The interfacial tension between both fluids was = 35 mN m -1 .
The rock sample was first saturated with brine and then n-decane was injected at high rate (1000 µL/min) to reach a non-wetting phase saturation around 78%. Then the imbibition experiment was initiated by injecting brine at a volumetric injection rate of 0.1 µL/min, corresponding to a microscopic capillary number of Ca micro = 1.8·10 -8 which is defined as (Lake, 1989).
Here  w =1 mPa·s is the wetting phase viscosity, v the linear flow velocity of the injected wetting phase. The cluster-based or macroscopic capillary number micro p cl macro Ca r l Ca     /  was between 9.1·10 -7 and 1.3·10 -6 (l cl is the length of the non-wetting phase clusters determined from the CT images and r p is the average pore throat radius), indicating that the flow regime (based on average flow rates) is capillary-dominated. During the flow experiment, 60% of the non-wetting phase (in total 2.5·10 9 µm 3 ) was mobilized. The flow experiment is described in more A C C E P T E D M A N U S C R I P T 6 detail in , where saturation as a function of time and the (homogeneous) saturation distribution(s) are also shown.
During imbibition the pore-scale saturation distribution was imaged by synchrotron-based fast X-ray micro-tomography under dynamic flow conditions at the TOMCAT beamline at the Swiss Light Source, Paul Scherrer Institute, Villigen, Switzerland. The beam energy was 37 keV (monochromatic, with a multilayer filter). For flow experiments, 1500 projections were collected at 20 ms acquisition time using an edge camera (PCO AG, Kelheim, Germany). The spatial resolution was 2.2 µm/voxel and the temporal resolution for a 3D rendering was 45 s. Scans of the dry samples were obtained at 60 ms acquisition time to obtain a higher-quality image, as displayed in Figure 1A. It was used for more accurate segmentation of the pore space that is then later used as a mask for segmenting the fluid phases during the flow experiment. The respective greyscale image from the flow experiment is shown in Figure 1B, where the doped water phase appears in bright and the (undoped) oil in dark. Image processing (filtering and segmentation, see Schlüter, 2010 andSchlüter, 2014) was performed with the Avizo Fire 8.0 software package (Visualization Science Group), following a workflow previously described in (Berg, 2013a;. The reconstructed images were first filtered with a non-local means filter and then registered on the dry scan to apply the mask for the pore space of the rock. Then the fluid phases were segmented using a watershed-based segmentation method followed by a cleanup step using morphological operations. An example of the resulting segmented image is displayed in Figure 1C.  (bright) and oil (dark) in the pore space, and (C) segmented image from the flow experiment using the segmented pore space from the dry image to mask the pore space.

Pore-Scale Flow Simulations
The pore-scale imbibition experiment corresponds in principle to an unsteady-state experiment . However, the accessible length scale of the experiment, i.e. the length of the sample, is smaller than the capillary dispersion zone . Therefore, the Darcy-scale methodology like the analytical JBN inversion technique (Johnson, 1959) or history matching with a reservoir simulator (Berg, 2013a) cannot be used to determine the relative permeability from the experiment. The pore-scale fluid distribution was available from the CT imaging. Single-phase Navier-Stokes flow simulations were performed on the segmented fluid phases obtained from the CT imaging during the flow experiment using the software package GeoDict (version 2014; Fraunhofer ITWM, Math2Market GmbH, Kaiserslautern, Germany) directly on the voxel grid from the segmented CT images. Wetting and non-wetting phase flow were simulated separately for each time step i of the experiment (with respective wetting phase saturation ( pore space of the dry sample, yielding (in the connected pathway) relative permeabilities as for the wetting and non-wetting phases, respectively. The obtained relative permeabilities represent the connected pathway flow contribution, but are based on a saturation distribution obtained under dynamic conditions. They do not consider the flux contribution of non-wetting phase ganglion dynamics.
Capillary pressure vs saturation relationships were approximated with a morphological approach (Hilpert and Miller, 2001;Ahrenholz, 2008;Silin, 2010) implemented in GeoDict. Two approximations were applied, which were (i) the morphological approach approximating drainage or imbibition as a sequence of (Young-Laplace) equilibrium states, as opposed to dynamic approaches involving moving liquid-liquid interfaces, and (ii) the Young-Laplace equilibrium states found by a series of morphological erosion and dilation operations using a spherical structuring element. The drainage algorithm starts with a completely water-saturated porous medium and begins with a maximal pore radius R, which is then decreased step by step. In each step, a pore is drained if (i) the pore radius is larger than R, (ii) the entering non-wetting phase is connected by pores larger than R to the nonwetting phase reservoir, and (iii) the replaced wetting phase is connected to the wetting phase reservoir. Conditions (i) and (ii) describe the original algorithm as defined by (Hilpert and Miller, 2001) and used in (Becker, 2008;Vogel, 2005). Condition (iii) introduces a residual wetting phase and was added by (Ahrenholz, 2008). The imbibition algorithm in turn starts with a completely oilsaturated porous medium and minimal pore radius R, which is then increased step by step. In each step a pore is imbibed if (i) the pore radius is smaller than R, (ii) the entering wetting phase is connected to the wetting phase reservoir, and (iii) the replaced non-wetting phase is connected to the non-wetting phase reservoir (Ahrenholz, 2008). This approach is computationally less costly than the (capillary-dominated quasi-static) level-set method (Jettestuen, 2013) and also more efficient than simulating the Navier-Stokes equations for two-phase flow with mobile liquid-liquid interfaces. This approach was first used to obtain a simulated mercury-air (Hg-air) intrusion capillary pressure for comparison with experimentally obtained Hg-air intrusion data (often also referred to as mercury intrusion capillary porosimetry, MICP). In a second step the morphological method was also used to simulate drainage and imbibition. From the resulting pore-scale saturation distribution connected pathway flow relative permeabilities were also obtained in a similar way, as described above. In GeoDict this workflow is automated in the SatuDict module. As evident from approximations taken in the workflow the resulting relative permeabilities were entirely quasi-static.

Steady-State Relative Permeability Measurements
For validation of the pore-scale modelling we compare the pore-scale simulation data from the beamline experiment to steady-state relative permeability measurements conducted on a larger dimension sub-sample from the same block. Darcy-scale relative permeability was determined by a steady-state method (Dake, 1978) using a custom-built X-ray saturation measurement setup at Shell in the Netherlands (Kokkedee, 1996;Berg, 2008). Cylindrical rock samples of 5 cm height and 3.8 cm (1.5") diameter were first cleaned using Soxhlett extraction with toluene and, subsequently, an azeotrope mixture of chloroform (84%), methanol (14%), and water (2%). After drying, the porosity  was determined by helium porosimetry. The absolute permeability to brine K abs was determined in a flow experiment from the pressure vs flux relationship for multiple flow rates according to eq. 1: Then the sample was mounted inside a sleeve in an X-ray transparent Hassler-type core holder and transferred into the flow setup for steady-state relative permeability measurement by co-injection of saline brine (doped with potassium iodide for increasing the X-ray contrast) and n-decane. The measurements were conducted at constant flow rate, starting with primary drainage where the fractional flow f w (which is the ratio of water flux over total flux, see also eq. 3) was systematically changed from 100% brine (f w = 1) to 100% oil (f w = 0), followed by 1 st imbibition where the fractional flow f w was systematically changed from 100% oil (f w = 0) to 100% brine (f w = 1) in 10 saturation steps each. At each saturation step f w was kept constant and the saturation and phase pressures were recorded when steady-state was reached (constant phase pressures and saturations, i.e. in practice variation of the respective averages less than 1% per hour). Pressures were recorded with differential pressure transducers (Kokkedee, 1996). In addition to flow rates and phase pressures, the saturation profiles at 30 positions along the rock sample were also measured, by X-ray transmission measurements (which are calibrated to both 100% brine and oil saturation). The basic procedure and measurement cycles are sketched in Figure 2. For details of the equipment, operating procedures and accuracy we refer to (Kokkedee, 1996;Berg, 2008).
From the resulting phase fluxes q i , pressure drops i p  (using pressure data over length L of the central part of the sample only) and saturations, the relative permeability could in principle be directly determined using the two-phase extension of Darcy's law (equation 1; Dake, 1978). However, in order to account for potential capillary end-effects (Huang, 1998), all relative permeability data was obtained by numerical simulation and history matching using Shell's in-house reservoir simulator MoReS (Kokkedee, 1989;Regtien, 1995). More details on the experimental methodology are given in (Kokkedee, 1996;Berg, 2008).

Mercury-Air Intrusion Capillary Pressure Measurements
Mercury porosimetry was used to determine the mercury-air (Hg-air) capillary pressure p c (Stegemeier, 1977). Measurements were performed on smaller sub-samples from the same Gildehauser block using an automatic pore injectivity apparatus (Autopore 9520, Micromeritics). A packing or closure correction is applied during the mercury intrusion, which accounts for packing of mercury around the sample surface prior to entrance of into the pores. The porosity of the sample was determined from the volume of mercury intruded into the sample at 4100 bar, assuming that at this pressure the pore space is completely saturated with mercury. 3 Results and Discussion

Relative Permeability from CT Flow Experiment
The basis for the relative permeability calculation from the CT flow experiment are the (pore-scale) wetting and non-wetting phase saturation distributions obtained at each time step during the imbibition experiment. The respective connectivity of the wetting and non-wetting phases was determined using the Avizo software package. In our experiment both the wetting and non-wetting phases maintain connectivity between inlet and outlet, i.e. both are percolating. However, as the wetting phase saturation increases, the non-wetting phase becomes increasingly disconnected, which is visualized in Figure 3A. Initially the non-wetting phase is mostly connected (red). Then during imbibition the fraction of connected non-wetting phase (red) decreases and the fraction of disconnected non-wetting phase increases (yellow). The decrease in non-wetting phase saturation, loss of internal connectivity and its increasing level of fragmentation into disconnected clusters  causes the permeability of the oil phase to decrease as shown in Figure 3C. As the non-wetting phase saturation decreases, in turn the wetting phase saturation increases and consequently the wetting phase permeability increases, as displayed in Figure 3C. By normalizing the wetting and non-wetting phase connected pathway permeability to the absolute permeability K abs , the (connected pathway) relative permeabilities were computed and plotted against wetting phase saturation in Figure 4. Relative permeability-saturation functions often show power law-like behavior (Corey, 1954), with values being significant (experimental accuracy on flux measurement is about 10 -3 ) and relevant for practical field applications over three orders of magnitude. Therefore it is customary in the special core analysis community, which in addition to the hydrology community is also addressed in this paper, to always display relative permeability (and also compare different datasets, validate fits etc.) both on linear and logarithmic scales as is done in Figure 4 and also later figures.  Following the main objective for this study, i.e. to investigate to which extent quasi-static methods are sufficient to estimate relative permeability in imbibition, the results from the quasi-static simulations are validated against a reference dataset from steady-state imbibition experiments. In Figure 4 the relative permeabilities from the Darcy-scale steady-state imbibition experiments are superimposed for direct comparison. There is fair agreement between the two datasets in terms of magnitude of relative permeabilities, saturation range, general trends and shape of the curves. The only obvious difference is the saturation endpoint, which for the µCT experiments apparently lies at 10% higher water saturation than for the steady-state experiment.
However, in our experience, such differences in endpoint saturation lie well in the range of subsample heterogeneity, and can potentially also be an effect of sample length . To compensate for differences in endpoint saturation, we replotted in Figure 5 the relative permeability data from Figure 4 in terms of mobile wetting phase saturation representing the saturation interval between the irreducible wetting and non-wetting phase (endpoint) saturations, S w,c and S o,r , respectively. For the µCT experiment, S w,c = 0.21 and S o,r = 0.29, and for the steady-state experiment, S w,c = 21 and S o,r = 0.38 were used, respectively (see also Table  1). There is now a better agreement both on the linear and on the logarithmic scale within the experimental error (which for our equipment is approximately 3% for saturation and 10% for relative permeability, see Berg, 2008), except for the water endpoint. It is important to mention here that this comparison is very sensitive to the residual oil saturation S o,r in particular for the dataset obtained from the µCT experiment. S o,r is very difficult to determine in a reliable way in such dynamic flooding experiments. Centrifuge data would be required for a more reliable determination of S o,r , but which was not available.
Apart from the difference in endpoint saturation, i.e. after rescaling to the mobile saturation range, there is good agreement between the steady-state measurement and the quasi-static simulation of the beamline experiment. This means that for this particular case, where the non-wetting phase percolates over the whole saturation range, the contribution of ganglion dynamics (which has been observed in  to the overall non-wetting phase flux is minor.

Impact on Fractional Flow
The question is now whether the difference between the steady-state experimentally measured relative permeability and the quasi-static simulation from the beamline experiment have any significant implications for practical applications such as oil recovery. Without going into the details of fractional flow, a first good indicator for the amount of recoverable oil is the crossover point between k r,w and k r,o . The crossover point indicates the saturation at which the relative permeability of water and oil are approximately the same and it is a first indicator of the approximate economic limit for oil recovery. In Figure 4 the crossover point for the steady-state data is at S w = 0.54 whereas for the quasi-static simulation from the beamline experiment the crossover point is at S w = 0.61, which is significantly different. Using the normalized mobile saturation data in Figure 5, the crossover point for both datasets equals at S * = 0.80 (and the exact endpoint saturation S o,r can then be determined from a more applicable centrifuge measurement). For a more detailed analysis, we consider the fractional flow f w for the water phase, which is the ratio of the water flux q w over the total flux (q w + q o ). In the limit for vanishing capillary pressure (p c = 0) f w is expressed as (Buckley and Leverett, 1942): where the water viscosity  w = 1 cP, and the oil viscosity is  o = 0.92 cP. In the next step the fractional flow curves f w (S w ) were compared between the quasi-static simulation of the CT experiment and the experimental steady-state data.
In addition to f w directly computed from the individual data points from Figure 4, also smoother f w curves were computed from fits of the relative permeability data with Corey functions (Corey, 1954;Masalmeh, 2006): with the mobile saturation S * from eq. (2). The results of the fit are displayed in Figure 6, indicating a valid representation of the steady-state data by the Corey function. The respective parameters of the Corey fits, i.e. endpoint saturations and endpoint relative permeability  Table 1.
The Corey function parameterizes relative permeability with a smooth curve that follows a power law functionality and the last data point, i.e. the endpoint, should be consistent with the overall trend of the curve representing the other data points. Assuming that relative permeability is well described by Corey functions, the comparison presented in Figure 6 provides an additional level of consistency and reduces the uncertainty with respect to the endpoint saturation S o,r . For the simulation data of the µCT experiment, the oil relative permeabilities are represented well by Corey functions but the water relative permeabilities at low saturations clearly deviate from a Corey representation (deviation larger than the experimental error). In particular in the saturation range S * < 0.6 the agreement for k r,o is good but for S * > 0.6 the agreement is modest at best, in particular for the water endpoint. The origin of this particular discrepancy is not clear. However, the impact on the fractional flow curve in Figure 7 and the breakthrough saturations and recoveries are within the experimental uncertainty.  Figure 4. Here S w,c is the connate water saturation, S o,r the residual oil saturation, k r,o (S w,c ) the oil relative permeability endpoint (at connate water saturation), k r,w (1-S o,r ) the water relative permeability endpoint (i.e. at residual oil saturation), and n w and n o are the water and oil Corey exponents, respectively. S w @BT is the water saturation at breakthrough, f w @BT the respective fractional flow, and Recovery@BT the respective oil recovery.    Figure 7, the fractional flow curves of the quasi-static simulation of the µCT experiment were compared with the experimental steady-state data using both the absolute saturation S w and the mobile saturation S * scales. For 1D displacement in the absence of any capillary pressure, key parameters such as breakthrough water saturation (S w @BT), water cut at breakthrough (f w @BT) and recovery at breakthrough (Recovery@BT) of the shock front can be computed analytically from the fractional flow curve f w (S w ) using the Buckley-Leverett formalism (Buckley and Leverett, 1942;Dake, 1978). The graphical representation of these calculations is a tangent from S w,c to the fractional flow curve included in Figure 7. The saturation at which the tangent intersects with f w is the breakthrough saturation, and the respective value of f w is the water cut at breakthrough. The results of the Buckley-Leverett calculation for the steady-state experiment and simulation results on basis of the µCT experiment are listed in Table 1. Without normalizing to mobile saturation, there is a difference of 10% in saturation units for recovery at breakthrough between the steady-state experiment and the simulation results, which is roughly the difference in endpoint saturation. On the mobile saturation scale S * , the fractional flow curves and the respective tangents are very similar. Consequently, in the hypothetical case that the endpoints were identical, the respective breakthrough saturations and recovery at breakthrough would be identical within less than 5% for both the steady-state experiment and simulation results.

Capillary Pressure from the Dry Pore Space Using a Morphological Approach
A comparison between the experimentally measured (MICP experiment) and simulated (using the morphological approach in GeoDict, Becker, 2008;Hilpert and Miller, 2001;Vogel, 2005) Hg-air capillary pressure p c data is displayed in Figure 8. Note that both cases represent capillary-dominated techniques where the contribution of viscous forces is absent or at least negligible. There is a good agreement between the simulated and the MICP data in a wetting phase saturation range above 0.2, which is similar to observations presented by Hussain, 2014. The discrepancies at lower wetting phase saturations are related to limitations by the imaging resolution, where the pore space is not fully resolved . Figure 8: Simulated Hg-air intrusion capillary pressure with morphological approach using GeoDict and comparison with experimental MICP data for Gildehauser sandstone rock. In addition, the GeoDict oil-water "imbibition" capillary pressure simulation scaled to Hg-air by the interfacial tension ratio between oil-water and Hg-air is plotted for comparison.

Relative Permeability from Morphological Approach
Apart from end-point saturations, relative permeability obtained from quasi-static simulations based on experimentally measured pore-scale fluid distributions show good agreement with steady-state relative permeability reference measurements. One may argue that this is good news because for practical applications the computationally much cheaper quasi-static simulation could be used for predicting relative permeability. However, this approach is not practical because the pore-scale fluid distributions of wetting and non-wetting phases are required, which for the presented data requires synchrotron-based CT flow experiments. Therefore, as the next step a morphological approach (Becker, 2008;Hilpert and Miller, 2001;Vogel, 2005) was used to simulate the capillary-dominated invasion process. As already shown in Figure 8, this approach gives excellent agreement between simulated and measured mercury-air intrusion capillary pressure-saturation relationships. Similar observations have been made for Berea sandstone rock in a previous study . However, mercury-air intrusion is a drainage process. The important question is now how well this approach works for imbibition. Therefore in the next step we simulate pore-scale fluid distributions using the (capillary-controlled) morphological approach (Becker, 2008;Hilpert and Miller, 2001;Vogel, 2005) in "imbibition mode" as implemented in GeoDict. Then we compute the quasi-static hydraulic conductivities on the resulting pore-scale fluid distribution analogous to the approach taken with the CT data. The final results are compared to the steady-state relative permeability reference data.
The results show that the morphological approach using the "imbibition mode" setting provide relative permeability curves that are more representative of a drainage process than an imbibition process. As displayed in Figure 9, the resulting imbibition relative permeability shows very large discrepancies from the measured relative permeability and the simulation from the beamline data. However, the water relative permeability from the morphological approach matches well with the drainage steady-state relative permeability as shown in Figure 10. There are still discrepancies with the oil relative permeability, which indicates that the morphological approach has fundamental difficulties in "simulating" connected oil phases, in particular at low oil saturations. Also, the capillary pressure data from Figure 8, where the scaled oil-water imbibition capillary pressure superimposes with the Hg-air drainage capillary pressure, indicates that the morphological approach approximates more a drainage process with respect to the computed pore-scale fluid distribution. The reason why capillary pressure is predicted reasonably well by the quasi-static method is that the morphological method describes correctly in which pore the non-wetting phase is distributed in individual pores, also probing directly connected pore throats. Also, in a drainage process the conductivity is captured sufficiently well because it is an invasion process where, by design of the algorithm used, both phases are connected.
For imbibition, however, the long-range connectivity of the non-wetting phase is not described very well because processes like snap-off and corner flow are not captured. In reality, there is mass exchange between connected and disconnected fluid phases (Hilfer, 1998;Hilfer, 2006a;Hilfer, 2006b), which is not only adding to the net flux but also influences the time evolution of the porescale fluid distributions, including the configuration of the connected phases. The morphological approach used here employs local physics to determine the pore-scale saturation distribution (and hence is conceptually similar to quasi-static pore network modelling which uses local rules). Approaches based on local physics cannot capture the non-local and cooperative processes that have been recently observed in drainage (Armstrong, 2013; and imbibition Holtzman, 2015). However, our findings show that the local approximation has a more dramatic effect on imbibition than on drainage. Our work clearly shows that this has less to do with the quasi-static computation of the flow field based on a given saturation distribution but rather the approximation of the saturation distribution by a quasi-static morphological approach. A possible explanation can be offered by considering the different pore-scale displacement processes for drainage and imbibition. While for drainage frontal displacement and to some extent cooperative filling are the most relevant process, for imbibition the snap-off and corner flow processes are fundamental to what ultimately determines the connectivity of the non-wetting phase. Snap-off is a process that is influenced by both the local pore and fluid geometry in pore throats, but also the nonlocal, long-range connectivity of the non-wetting phase which is required to provide the volume of wetting fluid to allow the snap-off to occur (Mohanty, 1987). By applying only local rules, this longrange connectivity is not considered and only based on local rules snap-off processes are permitted even if the long-range connectivity may suppress them. Also, to capture the physics of corner the flow the tiniest crevices of the pore-structure must be imaged and simulated with sufficient resolution. Different realizations of the pore space corners would result in wetting film advancement along different paths and thus ultimately influence the sequence of pore filling events. Therefore the resulting pore-scale fluid distribution is different (and difficult to capture for imbibition), which ultimately leads to a significantly different connectivity and hence relative permeability of particularly the non-wetting phase.
With respect to the data from the morphological approach it is also important to note that the oil relative permeability is only accessible until saturation values of S*  0.4. This is because for larger water saturations, i.e. lower oil saturations, the oil phase is not percolating any more, as also illustrated in Figure 11, and thus the connected pathway hydraulic conductivity is zero. In the morphological approach the saturation value of S w  0.36-38 or S*  0.4 corresponds to a percolation threshold for the non-wetting phase whereas in the experiment the non-wetting phase percolates up to much higher saturations of S w  0.70. Figure 11: Comparison of the pore-scale fluid distributions from the beamline experiment (top) with the morphological approach from GeoDict in imbibition (bottom) as a function of the wetting phase saturation S w . Red indicates a connected and yellow a disconnected non-wetting phase. For the morphological approach the non-wetting phase disconnects i.e. stops percolating at much lower water saturations (around 36-38%) compared to the beamline data where the non-wetting phase percolates up to 70% in this particular experiment.
Currently, it is not fully clear whether it can be expected in all cases that connected pathway flow will dominate over ganglion dynamics and that the percolation threshold for the non-wetting phase ranges up to S w  0.70. Recent studies by Reynolds et al. (Reynolds, 2015) provide an example with lower percolation threshold and hence more contribution of ganglion dynamics. There the connected pathway flow approach would not succeed to compute relative permeability beyond the percolation threshold.

Summary and Conclusions
The ability to predict parameters for Darcy-scale flow from quasi-static pore-scale flow simulations conducted directly on the pore space imaged by CT was investigated. Previous studies have shown reasonably good agreement with Darcy-scale measurements for both drainage and imbibition when using a capillary-viscous flow simulation, which is computationally very expensive (Vogel, 2005;Ramstad, 2009;Ferrari, 2013;Koroteev, 2014). Therefore, the focus here is on quasi-static simulation, which comes at a substantially lower computational cost.
Previous studies (Hussain, 2014) have demonstrated that for drainage good agreement with Darcyscale measurements can be obtained using a morphological approach to estimate the pore-scale fluid distribution and then conducting quasi-static flow simulations on the connected fluid pathways. We find that for imbibition this approach works well if the pore-scale fluid distributions are available from a flow experiment with CT pore-scale imaging. Apart from a 10% difference in the endpoint saturation between twin samples, which can also have different origins than the applied methodology, the relative permeability computed on the connected pathway from a synchrotron beamline flow experiment matches closely with the Darcy-scale steady-state measurement. This means that for a given pore-scale fluid distribution where the non-wetting phase percolates over the whole saturation range, the contribution of ganglion dynamics  to the overall flux is very minor. However, it is not fully clear whether it can be expected in all cases that connected pathway flow will dominate over ganglion dynamics. Recent studies by Reynolds et al. (Reynolds, 2015) provide an example that seems to be more dominated by ganglion dynamics and there at least the connected pathway flow approach would not succeed to compute relative permeability, at least not at lower oil saturations.
In summary, we find that quasi-static flow simulations give a good approximation of relative permeability when the pore-scale fluid distributions at each saturation step are known and the oil is mainly connected. If the pore-scale fluid distributions are not available, the approximation with a morphological approach does provide a good basis for the computation of drainage capillary pressure (which matches with measured Hg-air intrusion data very well) but not for imbibition relative permeability. The reason is that quasi-static morphological approaches are not capable of predicting such pore-scale saturation distributions for imbibition. The relative permeability computed from the respective pore-scale fluid distribution using the connected pathway conductivity of the two fluid phases is much closer to a measured drainage relative permeability.

Acknowledgements
Ton Blok is gratefully acknowledged for the MICP measurements. Sebastiaan Pieterse and Niels Brussee are gratefully acknowledged for the steady-state relative permeability measurements. Fons Marcelis and Ab Coorn are acknowledged for sample preparation and Soxhlett cleaning. We acknowledge the Paul Scherrer Institut, Villigen, Switzerland for provision of synchrotron radiation beamtime at the TOMCAT beamline of the SLS and would like to thank Sally Irvine and Marco Stampanoni for assistance.

A Representative Elementary Volume
A representative elementary volume (REV) analysis with respect to porosity and permeability for the Gildehauser sandstone sample used for the synchrotron beamline flow experiments was performed using the CT image from Figure 1 following the methodology by (Al Ibrahim, 2012). For increasing window size, porosity and permeability values were determined and plotted as a function of window size as shown in Figure A1. In addition the mean + 1 standard deviation (standard deviation is obtained from multiple sub-samples at same window size) is plotted. According to (Al Ibrahim, 2012) the intersection between mean and mean + 1 standard deviation is an estimate for the REV. The results in Figure A1 show that the sample size of 4 mm diameter is larger than the (single phase) REV. Note that while for porosity the mean becomes independent of the window size for windows sizes larger than the REV, for permeability this is not the case because it is not a linear property. Also note that the 2-phase REV can be substantially larger than the single-phase REV (Georgiadis, 2013).