Evaluation of Capillary Pressure Methods via Digital Rock Simulations

Capillary pressure measurements are an important part of the characterization of petroleum-bearing reservoirs. Three commonly used laboratory techniques, namely the porous plate (PP), centrifuge multi-speed experiment (CM), and mercury intrusion (MICP) methods, often provide nonidentical capillary pressure curves. We use high-resolution μ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mu }$$\end{document}-CT images of Fontainebleau and Bentheimer sandstones to derive saturation profiles numerically in 3D at the pore scale through morphological distance transforms to simulate the above experiments. In the invasion simulation, the capillary pressure is realized by using as structural element a ball—whose diameter is a function of a local pressure potential, which in turn is a function of radial distance in centrifuge experiment and constant throughout the sample in MICP and porous plate measurements. To assess the effect on heterogeneous rock samples, we compare the computed saturation profiles of the relatively homogeneous sandstones to a highly heterogeneous numerical model of rock generated by a mixture of a Gaussian random field approach for the large-scale features and two Poisson particle processes at the small scale. The comparison of the image-based pore-scale numerical interpretation to capillary drainage experiments reveals their match and demonstrates the influence of boundary conditions and heterogeneity on the resulting saturation profiles and capillary pressure curves. Simulated centrifuge experiments may assist in estimation of experimental equilibrium times and provide a useful tool in speed schedule design.


Introduction
In petroleum and groundwater engineering applications, capillary pressure relationships are among the most vital inputs required for reservoir description and modeling. Information I. Shikhov · C. H. Arns (B) School of Petroleum Engineering, The University of New South Wales, Sydney, NSW 2052, Australia e-mail: c.arns@unsw.edu.au about capillary pressure is used to estimate reservoir initial fluid saturations, fluid contacts, and transition zones; evaluate caprock sealing ability and displacement pressure; provides supplimentary data for relative permeability estimates and enables to predict dynamic behavior. One or few special core analysis laboratory (SCAL) techniques are normally used to determine capillary pressure curves: the porous plate (PP), centrifuge (CM), and mercury intrusion porosimetry (MICP). Due to the substantial difference in the underlying physics, these methods often provide sufficiently different pressure curves and consequently residuals, which are difficult to reconcile. In this paper, we consider primary capillary drainage only.
Capillary pressures across the rock sample arise between interfaces of two immiscible fluids. Usually, one phase is considered as a wetting phase and the other as a non-wetting phase. However, intermediate cases complicating the picture are common. The drainage case, i.e., a non-wetting phase displacing a wetting phase applies to hydrocarbon migrating into a previously water-saturated rock. Thus, the drainage data can usually be used to predict nonwetting fluid saturation at various points in a reservoir, and the imbibition data can be useful in assessing the relative contributions of capillary and viscous forces in dynamic systems.
Of the three common methods (PP, MICP, CM) to obtain capillary pressure curves for a rock core, the porous plate technique (Leverett 1941) is usually considered the most accurate since it avoids a saturation gradient and operates with native fluids. It can be combined with other measurements, e.g., with resistivity index, leads to a relatively homogeneous saturation profile, and allows further analysis with the core. It is a slow method-20 to 30 weeks are not uncommon to obtain an oil-water drainage curve (Wilson et al. 2001).
Mercury intrusion capillary pressure measurements (MICP) (Purcell 1949) are fast and also provide access to very small pores. Because of the strong non-wetting behavior of mercury, the mercury invasion is close to ideal drainage and therefore provides an excellent measure of the connectivity of the pore structure. MICP also has a number of drawbacks: contamination of the sample preventing its further use, non-representative fluids, limitations on sample size, and safety/environmental issues.
Evaluation of petroleum reserves requires analysis of a vast number of cores, making time and accuracy equally critical. The multi-speed centrifuge method (CM), first proposed by Hassler and Brunner (1945), is the most obvious option for routine core analysis enabling both speed and potentially, accuracy. In a centrifuge drainage experiment, a core saturated with wetting fluid is rotated while stepwise increasing rotational speeds. The core holder contains non-wetting fluid, which is allowed to drain into the core. The amount of displaced fluid provides the estimate of average core saturation at each rotational speed ω. The raw average saturation data should be converted to the capillary pressure P ci at a certain point, normally at a core plug inflow face (inlet). Correspondingly, an average saturationS needs to be converted to a saturation at that point S(P ci ). That transform (S, ω) → (S, P c ) is the major complication of the centrifuge technique as no analytical solution exists.
Numerous approaches have been proposed to transform average saturation of a centrifuged core to a saturation at the inlet: Hassler and Brunner (1945) approximation (HB), "type curves" by Bentsen and Anli (1977), "an exact equation" of van Domselaar (1984), "theoretically correct analytical solution" by Rajan (1986), and Ruth's (1988) iterative match, to mention just few. All those techniques have been found to be approximations. A good account of comparing 13 data reduction techniques can be found in Seth (2006).  who analyzed 19 various techniques pointed out that for typical centrifuging conditions, the HB approximation can provide a poor estimate, while other known interpretation techniques provide estimates with associated error of ±3 saturation units at best (in addition to experimental error).
One long-standing problem associated with centrifuge capillary pressure measurements is an inaccurate estimate of average saturation at equilibrium for each rotational speed (O'Meara et al. 1992). Furthermore, the lack of knowledge about an appropriate speed schedule may lead to a low-quality data set (a few useful saturation points at equilibrium) even for extremely long measurements (Ruth 2006;Ferno et al. 2007). An equilibrium state for a given rotational speed can be defined as the saturation where no additional fluid production is observed. The difficulty in estimating the time required to reach an equilibrium in centrifuge experiments has been discussed by Slobod et al. (1951), Hoffman (1963, Szabo (1972), Ward and Morrow (1987), Fleury et al. (2000), Seth (2006), Ferno et al. (2007). From a practical point of view, the criteria of equilibrium include either the absence of observable production over a certain period of 1-24 h, (Hoffman 1963;Seth 2006;Ferno et al. 2007) or a fixed time step (Szabo 1972). The time required to establish equilibrium state of two phases saturating the rock sample mainly depends on permeability and wettability conditions. While generally 24 h is accepted as an appropriate time step between the rotational speed increments, it obviously may not work equally well for all rocks, considering that permeability may vary by as much as eight orders of magnitude (between shale to benchmark rocks). Another possible explanation for the variable time to reach equilibrium, and why equilibrium is difficult to foresee in samples with comparable permeability, was proposed by Ward and Morrow (1987). They emphasized the influence of fluid connectivity provided through liquid wedges retained at the corners of pores and through surface roughness of natural porous rocks. The difficulty in estimating the degree of connectivity in natural rock may explain why it is yet not resolved how long a core plug needs to be spun to reach saturation equilibrium. The proposed numerical model does not answer directly the question about the production rate during centrifuging; instead, it inherently reproduces a static final equilibrium state. This may be practically employed to estimate production between centrifuge speed increments allowing to apply time required for eqilibration, thus enable a meaningful speed schedule design.
In recent years improvements in the interpretation of centrifuge experimental results have been achieved by the development of techniques enabling registration of saturation profiles, including CT-assisted imaging (Wunderlich 1985), nuclear tracers imaging (Graue et al. 2002), NMR transverse relaxation (Kleinberg 1996), and MRI (Green et al. 2008). For cores exhibiting higher permeability, this poses the problem that fluids tend to redistribute quickly and the actual saturation profiles are difficult to establish.
While the previously proposed approaches rely on experiments, we develop a pore-scale technique characterizing the saturation distribution within the centrifuged rock plug based on mathematical morphology and a realistic digitized representation of the sample. In the present paper, we aim to examine the three main capillary techniques using numerically simulated drainage on realistic rock morphologies derived from 3D digitized tomographic images. Focus is given to centrifuge technique as it has not been simulated before using a combination of μ-CT and morphological transforms. The organization of the paper is as follows: In the next section, the theoretical background to the centrifuge drainage capillary pressure technique is given and associated data reduction problems are discussed. In the following section, we describe in detail the numerical approach and describe the simulations on Fontainebleau and Bentheimer sandstones, and GRF model structure. Comparisons of the simulated capillary drainage curves to experiment are given thereafter, followed by conclusions. Fig. 1 Schematic diagram of the centrifuge drainage method. In rotating Cartesian reference frame (RCRF), the axis of rotation is aligned along the Z-axis and the long axis of a core plug is along the transverse axis X. In the text, a radial distance r is a distance along the X-axis in RCRF; r i denotes a discrete value of r

Theoretical Summary
At hydrostatic equilibrium, fluids saturating a porous medium subjected to centrifuging are distributed following centrifugal force field (and its history). Below we express equations relating to capillary drainage in terms of centrifugal force field potential. This enables generalization and makes accomodation of gravity and radial effects easier. The classical expression of radial centrifugal force acting on an element of mass m is following: where r notes the distance from the axis of rotation, a c (r ) centrifugal acceleration, ω the angular frequency, and U c (r ) the potential energy in the centrifugal force field (see Fig. 1). Correspondingly, integration of (1) with respect to r provides the expression of a potential (potential energy per unit mass) which is defined by the following: It is easy to see that at hydrostatic equilibrium, capillary pressure P c and potential Ψ c are connected through the simple relationship: where ρ is the difference in density. A general expression for the potential combining gravity, centrifugal force, and a background energy independent on position (e.g., stepwise applied pressure) is: where g is gravitational acceleration. Here, we neglect the gravity term, so that for CM experiment dΨ/dr = −ω 2 r . For PP and MICP, the potential is independent on radial position along the core; dΨ/dr = 0.
Assuming that a core is homogeneous, the grain matrix incompressible, and interfacial tension constant (isothermal conditions), the pressure potential along a sufficiently short core subjected to spinning can be described by the following simple integral (Hassler and Brunner 1945): where r 2 is the radial distance from axis to the outer face of the core plug and r the radial distance from the axis to an arbitrary point of the plug (see Fig. 1). By applying the limits, the change in potential energy for the particle moving from arbitrary position r to r 2 in the direction of the centrifugal acceleration is given by: For the Hassler-Brunner outflow boundary condition of Ψ c = 0 (for a discussion see e.g., O'Meara et al. 1992), negligible equilibration time, absence of end piece and radial effects, absence of bubble formation, and short core approximation: r 1 /r 2 ≈ 1, Eq. (6) may be reduced to: where r 1 is the distance from the rotational axis to the near end (inlet face) of core. We use the expression for potential Ψ to define a local capillary radius R to set an interface of fluids saturating a core (see next section). In the experimental part (see Sect. 6), we use a traditional interpretation which involves match of observed average saturation to inlet capillary pressure. The most relevant equations are summarized below. The general mathematical expression describing a centrifuge experiment, where observed average saturation (S) is connected to pressure potential (Ψ c ) through the various parameters, can be represented as an integral sum of the local saturations along the length of the core and given by the following:S Making the usual assumption that the projection (intersection) of a sample on equipotential surface is slender enough (Ayappa et al. 1989), the average saturation can be represented as the integral sum of the local saturations along the core: By manipulating (9) and neglecting radial and gravity effects, the average saturation for drainage in the centrifugal field in respect of the inlet is given by Hassler and Brunner (1945) in integral form:S where ξ = P c /P ci is a dimensionless variable and B = 1 − r 1 /r 2 . The Hassler-Brunner expression suggests a simple solution for B → 0, which practically implies short cores: This expression may be written in a differential form as (referred in the text as HB approximation): We use both the HB approximation and Forbes method  for experimental data interpretation.

Numerical Model
We use mathematical morphology and appropriate morphological operations to describe the spatial positions of labeled fields representing phases saturating a porous medium. The entire sample volume may be divided into three domains: (1) pore space a ∈ A; (2) solid phase c ∈ C; (3) subspace x ∈ X, X ⊂ A representing a result of a set of morphological operations reflecting propagation of a draining fluid. The properties of the latter are represented by structuring element-a covering sphere, a radius of which mimics local experimental conditions of a particular capillary technique (PP, MICP, CM). First, we briefly introduce the basic morphological operations and notations following (Najman and Talbot 2010).
More complete descriptions of the basic concepts and techniques can be found elsewhere (Serra 1982;Hilpert and Miller 2001;Thovert et al. 2001). The connectivity tests on the set of invading spheres are implemented as in .
In Euclidean 3D space, the translate of B by a vector R ∈ E is the set B R , Here, R defines a translational vector. The complement of a set A is the set of elements; it does not contain, Dilation and erosion are dual by complimentation: the dilation of a set A by B is the erosion of its complementary set A c using the symmetric structuring element of B. The morphological dilation (or Minkowski sum) of a pore phase domain A by a structural B and A are subsets of underlying space E and index R denotes vector R-the translate of B by R, ∩-denotes intersection and ∪-stands for union. The resulting dilation is the union of the B R such that R belongs to A: The erosion is the dual operation corresponding to dilation of the complement A c of A and corresponds to all points in A not covered by a sphere B R centered out of A: The erosion of A by B is the locus of the points R such that B R is entirely included in A.
The morphological opening O of a geometrical object A by a sphere b of radius R is the domain swept out by all the translates of R that belong to A where O R (A) is the set of points in A that can be covered by a sphere of radius R contained in A. Opening filters out small convexities and removes smaller isolated clusters (Fig. 2). The numerical simulation of MICP experiments on micro-CT images based on morphological transformations was successfully demonstrated by . A good account of the method is given by Hilpert and Miller (2001). Here, we consider for all following Fig. 2 A series of centrifuge saturation RPM-maps attained on a beeds pack model. a RPM-map attained in 4,000-6,000 rpm interval; b expanded speed interval 4,000-8,000 rpm; c a mapped rotational speed interval further expanded to 4,000-21,000 rpm. In case of MICP and PP, the fluid would propagate through all the throats (white circles) once the far-left one is percolated, while CM requires further increase in rotational speeds due to the pressure potential Ψ c drop toward the outlet (right end) numerical drainage experiments an idealized wettability model: Solid is presumed to be totally wet to a wetting fluid, i.e., a contact angle of θ w−s = 180 • , while perfectly non-wet to the invading phase. Fluids are immiscible and incompressible, residual saturation due to electrokinetic interaction of wetting fluid with solid surface is ignored as well as snap-off effect at the outlet boundary. All three drainage experiments (PP, MICP, CM) are simulated in 3D by considering the invasion of a digitized image by a continuous structure using as structuring element B R -a ball of a radius R inversely proportional to pressure potential Ψ c (r ). In a general case, the radius of a structuring element B R (ω, x), where x is a radial coordinate and ω-radial speed, may be expressed as following: with receding contact angle θ r , pressure potential Ψ c , surface tension σ , and density difference ρ. In PP and CM experiments, boundaries other than inlet and outlet faces are completely closed and the invasion occurs just from one side (inlet, Figs. 1, 3c-d, g-h, left), while for MICP, invasion occurs from all sides (Fig. 3e-f). The outlet boundaries for PP and CM are given by semi-permeable membranes which only the wetting fluid can traverse (without resistance), while for MICP, the defending phase is considered to be vacuum. The porous plate (PP) drainage experiment is characterized by a slow propagation of invading nonwetting fluid from the inlet side of a sample toward the opposite outlet side, displacing wetting fluid initially saturating the sample. Wetting fluid is allowed to go through a fluid plate or membrane attached to the outlet impermeable for the draining fluid, so the outer boundary is effectively closed. Following the assumptions used in the present paper at any given pressure step, the pressure potential Ψ c is constant all across the sample. The spatial position of draining and displaced phases is described by a maximum inscribed sphere B R of a radius R, which is placed on one side of the medium and allowed to progress through the system until it reaches a pore that cannot be traversed without overlapping the interface. The saturation map relative to the particular pressure step can then be produced. The radius is decremented and the traverse step is repeated until the sphere is becoming trapped. This process continues until the discretized sphere macroscopically traverses the medium. For PP and MICP, a series of capillary pressures is considered by increasing the pressure of the invading phase, leading to a particular saturation map for each step in capillary pressure (R = const) with a constant structuring element B R . These saturation maps are combined into a single field by storing for each voxel the largest invasion radius when the voxel is drained (Fig. 3c-f). A saturation map for a particular capillary pressure  (g, h). The region close to the outlet (g, h) desaturates only at high rotational speeds can then be extracted by thresholding this field with a desired radius cutoff corresponding to a target capillary pressure.
For the centrifuge capillary pressure simulation, the approach is similar, with the difference that for a given rotational speed, the size R of B R is now a function of radial distance x from the inlet and R = R(x) (see Fig. 1). The size of the structuring element B R can be calculated via equation (6), (19) by converting the capillary pressure at location x to an equivalent capillary pressure, leading to a structuring element B R (x). A saturation map is now defined for each chosen rotational speed, and the natural way of accumulating all saturation maps is by storing the rotational speed at which a particular voxel is drained (Fig. 3g, h). A saturation map for a particular rotational speed, representing a whole range of capillary pressures, can then be recovered by thresholding this RPM-map (rotations per minute).

Rock Sample Representation
The numerical model of capillary drainage techniques was first tested on two idealized samples in this study. The first represents a homogeneous sample and consists of a 15 mm long and 3.75 mm wide reconstructed Fontainebleau sandstone using a stochastic geomet-rical model (Latief et al. 2010). We consider the discretisations provided at resolutions of 1.83 µm (512 × 512 × 2, 048 voxel), 3.66 µm (1, 024 × 1, 024 × 4, 096 voxel), and 7.32 µm (2, 048×2, 048×8, 192 voxel). These reconstructions were shown to reproduce the Minkowski functionals of the original imaged rock samples as well as e.g., local percolation properties (Latief et al. 2010). We use the different discretisations to discuss discretisation effects for the CM method in the discussion section. The bulk porosity of the sample is φ = 13.6 %, and we calculated the permeability to be about k = 1.18 ± 0.04 D using a MPI parallel LBM method with Cartesian decomposition based on Arns (2004). The permeability calculations were carried out on 1, 024 3 blocks of the reconstructed Fontainebleau image at a resolution of 1.831 µm, e.g., on blocks of sidelength 2 mm, significantly above the size considered in Arns (2004). Compared to the experimental data in Fredrich et al. (1995), Arns (2004), this is slightly higher than that expected for the porosity of 13.6 %.
The second sample is constructed as a dual-scale medium by using a Gaussian Random Field (GRF) to spatially separate two independent Poisson particle placement processes of spheres (r = 12 voxel) and oblate disks (half-axes a = b = 24, c = 3 voxel), discretised at a resolution 4.688 µm (1, 600 × 1, 600 × 3, 200 voxel). We follow the technique of Cahn (1965), Roberts (1997),  to generate a GRF and use the field-field correlation function where r c = 0.4033, ξ = 0.4031, d = 7.7069 in voxel units (see Roberts 1997;). We take a symmetric 1-level cut through the resulting field to separate the space into two equal partitions. The particle density of the Poisson processes is chosen such that the porosity in both partitions is φ = 30 %. The length of the samples (15 mm) is selected to be sufficiently short in respect of standard Beckman centrifuge distance to inlet to test the Hassler-Brunner approximation (Hassler and Brunner 1945), which is supposed to work well if r 1 /r 2 > 0.7. Here, we have r 1 = 7.1 cm and r 2 = 8.6 cm and therefore r 1 /r 2 = 0.83. Cross sections of the considered samples are depicted in Fig. 3, first row. Comparison between simulation and experiment was performed on Bentheimer sandstone (see Sect. 6).
To aid the comparison of samples and provide a measure of heterogeneity, we derive the variograms of the digital representations of all samples; Fontainebleau sandstone at 3.66 µm resolution, GRF/Boolean (4.688 µm resolution), and Bentheimer (2.89 µm resolution). We use the approach taken in Marcotte (1996),  to calculate smoothed surrogate measures of the Minkowski functionals. The Minkowski measures itself are derived using the algorithm introduced in Arns et al. (2001). In the notation of , we define the We can then define the curvature measures C X,i of a given (and then fixed) random structure X for a variable set B in R 3 as C X,i (B) = V i (X ∩ B). We use as Ba sphere b of radius R = 5 voxel which provides the support of the measurements. The sphere is moved over the structure X to derive a field of curvature measures Z i (x) = C i (b(x, R)) for x ∈ R 3 , from which the variograms are derived; γ i (r ) = 1 2 (Z i (o)−Z i (r)) 2 . The variograms are depicted in Fig. 4 and normalized by their respective variances. Mean and variance for each measure with respect to the macro-solid (grain) fraction are given in Table 1. The GRF/Boolean model was constructed with constant porosity across both largescale partitions generated by the GRF model. Accordingly, we see only the short correlation length of the Boolean model reflected in γ 3 (r )/σ 2 3 (Fig. 4d). However, both the length scale of the GRF coarse-scale structure controlling the distribution of the Boolean sphere models and  Table 1) of the intrinsic volumes and curvature measures C X,i (B) of the reconstructed Fontainebleau sandstone, the GRF/Boolean model, and Bentheimer sandstone for R = 5 voxels following the notation of . a Euler characteristic χ V , b integral mean curvature M V , c surface area, S V , and d volume V V Table 1 Mean μ and variance σ of the curvature measures C X,i for reconstructed Fontainebleau sandstone, the GRF/Boolean model, and Bentheimer sandstone Sample and moment ( the short length scales associated with the sphere radii are reflected in γ 1 (r )/σ 2 1 and γ 2 (r )/σ 2 2 . A single-phase REV based on the two-point correlation function is typically reached at about 5-10 times the correlation length (Joshi (1974)). This would imply an REV of less than 1 cm 3 for all structures. For two fluid phases, it is known that the REV can be substantially larger (Arns et al. 2003). Given that the correlation lengths for some of the structural measures are an order of magnitude larger for the GRF/Boolean model than for both Fontainebleau and Bentheimer according to Fig. 4, we expect more heterogeneity effects for the GRF/Boolean model. This is in particular the case since the generating grain for the GRF/Boolean model is a sphere for both partitions; the S/V ratio of each partition roughly reflects the ratios in critical radius of the two domains, which has a significant impact on capillary drainage as it controls percolation.

Numerical Results
The simulated drainage experiments on Fontainebleau sandstone and the dual-scale GRF/Boolean model structure provide the saturation maps (Fig. 3) and profiles (Fig. 5) over a broad range of capillary pressures and rotational speeds. Figure 3c-h depicts the x-z cross sections of the drainage maps for the central y-slice. The invasion from all sides for MICP is clearly visible, and for the GRF, some pores in the center of the structure are still invaded at lower pressure. Furthermore, there is a clear ordering in the invasion of small and (b) (a) Fig. 6 Comparison of the simulated capillary pressure curves with the Hassler-Brunner approximation. The point sets depict the capillary pressure curves for a selection of equally spaced slices along the core for the centrifuge drainage simulation case, thus enabling a comparison of the average capillary saturation relationship with local measurements, which for the homogeneous case (top) coincide. a Fontainebleau sandstone. b GRF/Boolean structure large porosity for the GRF/Boolean system. Figure 3 also illustrates the relative homogeneity of the saturation profiles for PP and MICP away from the boundaries. Figure 5 shows the relative homogeneity of the saturation profiles for the PP method. In comparison, the saturation profiles for the centrifuge method show strong trends reflecting the spatial gradient in the potential function. The maps as well as the saturation profiles illustrate clearly the difference in propagation of a saturation front a) for the different capillary drainage techniques and b) between homogeneous Fontainebleau sandstone and the heterogeneous dual-scale GRF/Boolean model system. The dual-scale structure of the GRF model already illustrated using the variogram technique (Fig. 4) leads to larger-scale variability in the saturation profile (Fig. 5b, d, f). The resulting capillary pressure curves according to (7) are depicted in Fig. 6.

Validation by Experiment
We compare numerically simulated MICP and centrifuge drainage to experiment using Bentheimer sandstone samples and their digitized representation.
The CM and MICP experiments were performed on Bentheimer core plugs of (2.5 cm diameter, 2.5 cm length) and (0.5 cm diameter, 1.5 cm length), respectively. The brine porosity was measured as φ = 23.3 %. The MICP experiment was performed using a Quantachrome PoreMaster PM33-17 porosimeter. Acquired data were reduced assuming a surface tension value σ Hg−vac = 480 dynes/cm 2 and a contact angle σ = 140 • . CM experiments were performed using Beckman L60 M/P ultra-centrifuge (PIR-20 rotor) with 12 h between rotational speed increments (22 steps in total). Increments were 50 rpm at the low-speed region (below 1200 rpm) then gradually increased. The displaced volume of brine was visually registered with the aid of stroboscope and graduated pippets of the bucklets containig the centrifuged samples following the standard procedure. Raw data were reduced using Forbes method ) accounting for gravity and radial effects.
The numerical part was performed on a 2.89 µm resolution μ-CT image of 1, 800 × 1, 024 × 1, 024 voxels. The image was mirrored twice, resulting in dimensions of 5, 400 × 1, 024 × 1, 024 voxels corresponding to a real physical length 1.5 cm (see Fig. 7, top). The corresponding centrifuge drainage RPM-map is given in Fig. 7, bottom. Note that the outlet is on the left side. Furthermore, the highest RPM simulated is lower than in Fig. 3, thus Pressures are converted to air/water a significant part of the sample has not been invaded. The propagated draining fluid (air) as a function of rotational speed (rpm) is coded by different colors. Note a perceivable heterogeneity of Bentheimer expressed as fingering of propagating drained fluid at 750, 950, 1,100, and 1,500 rpm. Furthermore, the induced symmetry by mirroring the sample is reflected as saturation spikes in Fig. 8.b. This is likely due to a shift in connectivity across the mirrored boundaries, where some pore space can only be drained at relatively elevated RPMs. We compare in Fig. 9 the simulated capillary pressure curves to the experimentally derived ones. Both simulated and experimental MICP capillary pressure curves are converted to airwater fluid pairs (σ w−air = 72 dynes/cm 2 ) to enable comparison to PP/CM data. There is very good agreement between the MICP simulation and experiment. The match between experiment and simulation for CM is acceptable.

Discussion
Both simulated MICP and CM primary drainage capillary pressure curves attained on the Bentheimer μ-CT image being adjusted for air-water pair almost perfectly coincide, except for the very high saturation region where the observed difference is due to the specific boundary conditions (see Fig. 9). In the low-saturation region, simulated curves deviate substantially from experiment due to the image resolution limit: the unresolved clay region, which has a volume fraction of 1.5 %, contributing about 0.7 % porosity, would in principle account for the steep end of the curve together with fine-scale surface roughness. In the medium saturation region, the simulated MICP agrees well with experiment; however, the measured CM capillary pressure curve deviates from the other curves. Simulations were performed in assumption of perfect wetting conditions, i.e., receding angle θ r = 0 • and uniform and known surface tension values of σ w−air and σ Hg−vac . While 0 • is often used in numerical models (Øren and Bakke 2003;Zhao et al. 2010), in a real experiment the contact angle of natural rock (without plasma cleaning or firing) may depart quite significantly from that ideal value. In the future, we plan to extend this work to include nonzero contact angles, which can, e.g., be introduced morphologically using the level-set method (Sethian 1999;Prodanović and Bryant 2006).
Other factors include sample shape and size, which result in different entry effects. Simulations were performed on samples represented by a rectangular slender prism, 5 mm long (mirrored twice in case of CM for a total length of 1.5 cm), while cylindrical samples of 25 mm length were used in the experiments.
We finally consider discretization effects in the numerical approach using the reconstructed Fontainebleau sandstone sample at different resolutions for the case of the CM method. Figure 10a, b shows the saturation profiles of Fontainebleau at higher ( x = 1.831 µm) and lower ( x = 7.324 µm) resolution than in Fig. 5e, as well as the respective differences to x = 3.662 µm. A positive difference indicates higher saturation for higher resolution. The largest absolute errors occur closer to the saturation front, where only the higher curvatures occur and where the saturation profile is steepest. This is associated with the discretization of a sphere on the lattice becoming increasingly inaccurate for small sphere diameters affecting both percolation through constrictions and into crevices. The error is in general well below 0.1 saturation units for the steepest part of the curve or below about 10 % relative error. In terms of the resulting centrifuge capillary pressure  (Fig. 10e), the difference between the x = 1.831 µm and x = 3.662 µm curves is marginal and the difference between x = 3.662 µm and x = 7.324 µm is well below 10 % except for the high-pressure region where one enters the resolution limit.

Conclusions
We present a numerical approach enabling comparisons of the three main capillary drainage experimental techniques. The impact of heterogeneity on capillary drainage estimates is demonstrated by simulation of three systems with different degree of heterogeneity. The numerical results have been experimentally validated. It was demonstrated that for the sufficiently homogeneous system of appropriate size, the CM drainage agrees well with porous plate and MICP results. Furthermore, a simulated centrifuge experiment is free from main uncertainties specific to CM-unknown equilibrium time between speed increments and a transform between average core saturation to saturation at the inlet face provides a tool for CM experimental design. In addition to comparisons between different capillary pressure measurement techniques, this morphological method allows to design more targeted speed schedules for complex workflows involving centrifuge experiments. It may also assist in the typing of microporosity by comparing CM simulation with actual measurements for rocks exhibiting significant porosity below imaging resolution.
The current work does not consider the simulation of more complex wettability conditions including more advanced level-set morphological techniques or the extension to relative permeability prediction. This will be the subject of a future article.