ELSEVIER

Contents lists available at ScienceDirect

Journal of Colloid and Interface Science

journal homepage: www.elsevier.com/locate/jcis

Automatic method for estimation of in situ effective contact angle from X-ray micro tomography images of two-phase flow in porous media

Alessio ScanzianiaÄ*, Kamaljit Singh b, Martin J. Blunta,b, Alberto Guadagninia,c

a Politecnico di Milano, Dipartimento di ¡ngegneria Civile e Ambientale, Piazza Leonardo da Vinci 32, 20133 Milan, Italy b Imperial College London, Department of Earth Science and Engineering, London SW7 2AZ, UK c University of Arizona, Department of Hydrology and Atmospheric Science, Tucson, fZ 85721, USA

GRAPHICAL ABSTRACT

ARTICLE INFO

ABSTRACT

Article history: Received 2 December 2016 Revised 4 February 2017 Accepted 6 February 2017 Available online 8 February 2017

Keywords: Contact angle Multiphase flow Porous media Wettability Micro-CT imaging

Multiphase flow in porous media is strongly influenced by the wettability of the system, which affects the arrangement of the interfaces of different phases residing in the pores. We present a method for estimating the effective contact angle, which quantifies the wettability and controls the local capillary pressure within the complex pore space of natural rock samples, based on the physical constraint of constant curvature of the interface between two fluids. This algorithm is able to extract a large number of measurements from a single rock core, resulting in a characteristic distribution of effective in situ contact angle for the system, that is modelled as a truncated Gaussian probability density distribution. The method is first validated on synthetic images, where the exact angle is known analytically; then the results obtained from measurements within the pore space of rock samples imaged at a resolution of a few microns are compared to direct manual assessment. Finally the method is applied to X-ray micro computed tomography (micro-CT) scans of two Ketton cores after waterflooding, that display water-wet and mixed-wet behaviour. The resulting distribution of in situ contact angles is characterized in terms of a mixture of truncated Gaussian densities.

Crown Copyright © 2017 Published by Elsevier Inc. This is an open access article under the CC BY-NC-ND

license (http://creativecommons.org/licenses/by-nc-nd/4XI/).

Abbreviations: micro-CT, micro computed tomography; 2D, 3D, two-three dimensional; USBM, United States Bureau of Mines; SEM, ESEM, scanning electron microscopy, environmental SEM; ML, maximum likelihood; RMSE/RMSD, root mean square error/difference; ROI, region of interest.

* Corresponding author at: Politecnico di Milano, Dipartimento di Ingegneria Civile e Ambientale, Piazza Leonardo da Vinci 32, 20133 Milan, Italy.

1. Introduction

Multiphase flow in natural or engineered porous materials is ubiquitous in modern applications involving e.g., chemical transport, geochemical reactions, seawater intrusion, gas diffusion in fuel cells, conventional and unconventional oil recovery, and

http://dx.doi.org/10.1016/j.jcis.2017.02.005 0021-9797/Crown Copyright © 2017 Published by Elsevier Inc.

This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).

carbon dioxide capture and storage [1-8]. The pore scale arrangement of multiple fluid phases in such porous materials is controlled by the topological and geometrical characteristics of the pore space and the wettability of the system. Wettability is represented through the spatial distribution of contact angle at the three-phase contact between two residing fluids and the host solid matrix, and its in situ characterization in complex pore spaces is still an open challenge. The arrangement of fluids in a porous medium is governed by the Young-Laplace equation, which stems from energy balance considerations [1] and defines the equilibrium capillary pressure Pc, related to the interfacial tension y and the principal curvature radii (ri and r2) of the interface:

pc - y(n+r2) -Ky' (1)

where k is the total curvature of the interface. It can be noted that Pc, and hence the curvature, is constant when the system is at equilibrium, i.e., the fluids are at rest: this will be the key insight used in our approach to characterize the wettability.

The contact angle Q between two-phases (w and n) is related to the interfacial tension between the two fluids through the Young equation [9]:

Cs-n = Cs-w + Cn-w cos Q; (2)

where ys-i, (i = n, w) and yn-w respectively are the interfacial tension between the solid surface and the two fluid phases, and the interfacial tension between the two fluid phases themselves. For example, in a water-oil system, wettability is classified on the basis of the contact angle, as measured through the water phase. When considering the flow of water and oil in a porous medium, water- and oil-wet systems are characterized by values of Q which are respectively less than or greater than 90° [10].

A series of studies have documented that the spatial distribution of wettability - the pore-scale distribution of contact angle - strongly affects fluid displacement and recovery in rocks [10-19] as well as the performance of fuel cells [20,21]. For example, fluid snap-off is enhanced in water-wet systems and causes high residual saturation of the non-wetting phase [1,2,22-24]. As another example, the wettability influences the relationship between capillary pressure and liquid saturation which in turn governs gas diffusion in the porous electrodes of fuel cells [7,25].

Measurements of contact angles have been typically performed on flat surfaces at ambient conditions through sessile drop and/or captive bubble approaches [11,13,26-28], for different pressure and temperature conditions [29]. In this context, micro-CT imaging has been used to observe drops of fluid [30]. Bachmann et al. [31] modified the sessile drop method for the assessment of the contact angle of powered or granular material. While these methods yield valuable information, they do not yield indication about the in situ contact angle within a porous medium, as roughness of the rock surface and the irregular shape of pores are not taken into account. Another key drawback of these approaches is that they usually consider pure mineral surfaces. The ensuing contact angle estimates are then seldom transferable to settings typical of engineering applications, such as hydrocarbon reservoirs, where the porous media have a mix of mineralogy and small-scale pore texture and morphology.

In current petroleum engineering practice the contact angle, or wettability, is indirectly inferred from capillary pressure curves [14], through the Amott or USBM (United States Bureau of Mines) indices [32,33]. These are also used for the characterization of gas diffusion of fuel cells [7]. Other indirect approaches are documented in the literature, such as the Wilhelmy plate method [34], the water drop penetration test [35] and the capillary rise

method [36]. These have been applied on different granular soils to obtain estimates of receding and advancing contact angles, without the direct visualization of the microscopic interfaces [37]. All the indirect methods provide an appraisal of the average behaviour of macroscopic samples without providing insights on the distribution of local contact angles.

Scanning electron microscopy (SEM) methods can also be applied to study contact angles between fluids and a surface at high resolution [38]. Cryo-SEM methods establish a representative in situ distribution of fluids, and then freeze and cleave the sample followed by imaging to determine contact angle [38,39]. However, the freezing process might alter the contact line between the phases: the method does not therefore directly yield contact angle in situ at the conditions under which the fluids flow in reservoir settings. Environmental SEM (ESEM) methods can study contact angles under a range of temperature and pressure conditions [40], but again do not directly observe these angles during displacement [40,41]. Instead what is needed is a way to observe the contact line inside a rock at representative conditions during, or at the end of, a multiphase displacement.

Knowledge of effective in situ contact angle is needed for pore-scale models of multiphase fluid flow to quantify the fluid configurations and threshold capillary pressures. Averaged quantities, such as relative permeability and capillary pressure can than be predicted using pore-scale modelling - see, for instance [42-44]. However, at present - in the absence of direct measurements - a distribution of contact angle is simply guessed, possibly to match measured capillary pressure or residual saturation, which adds a significant uncertainty to the predictions of these pore-scale models.

The advent of pore-scale imaging has allowed the determination of contact angle in situ from imaged rocks with a micron-scale resolution. This information can constitute a critical input to advanced direct pore-scale numerical modelling of multiphase flow because it can account for the effect on rock wettability of solid matrix roughness and mineralogical composition. Andrew et al. [45] were the first to illustrate a procedure to perform such direct measurements of rock wettability for a supercritical CO2/ brine system. This approach has been applied to different fluid systems in porous media by several authors [42,46-49] and here it will be used for comparison with the new proposed algorithm.

The approach of Andrew et al. [45] is based on a manual measure of contact angles. As such, it might be prone to bias/subjectivity and is time-consuming, thus hampering the possibility of acquiring extensive data sets to capture spatial distributions of contact angles. Khishvand et al. [46] applied this method on segmented images to estimate advancing and receding contact angles under two- and three-phase flow conditions. Klise et al. [50] proposed an automatic method, which was tested only on simple settings: bead packs comprised of beads of one or two uniform wettabilities. Their method relies on the availability of accurate imaging of the three-phase (e.g., water-oil-solid) contact line, a feature which is very rarely available due to the difficulty of obtaining a clear and accurate discrimination of the three phases in realistic rocks where there is inevitable trade-off between resolution and system size.

To overcome the challenge posed by the need to accurately discriminate each phase at the three-phase contact between two fluids and the solid, we propose an approach which is grounded on the physical insight that the fluid-fluid interface has a constant curvature. This will enable us to infer the location of the interface as well as the contact angle with the solid without the need to ground the entire analysis on very precise observations of the fluid-fluid interface.

The estimated angle is effective in the sense that it is the angle that the fluid-fluid interface would have in the case of a locally

smooth surface: this is the angle we require for pore-scale calculations of fluid configuration and threshold capillary pressure. Indeed, while the "true" contact angle, defined at the atomic scale, provides valuable information on the balance of surface forces, it is not the value we require for pore-scale modelling and interpretation: rather we wish to find an effective angle that incorporates the effects of roughness, flow direction and local pore geometry.

The method is first tested on simple synthetic two-dimensional (2D) images where the contact angle is known. It is then compared against the manual method [45] and is finally applied to assess the spatial distribution of contact angles associated with oil ganglia trapped within water-wet and mixed-wet rock samples.

2. Materials and methods

2.1. Materials

The two experimental datasets used to establish and test our approach were obtained from X-ray micro-CT scans of two samples of Ketton carbonate that experienced a complete drainage-waterflooding cycle. The samples were 25 mm long with a cross-section diameter of 4 and 5 mm, respectively for a water-wet and an altered wettability medium. Decane was the oil phase to reproduce a water-wet system while a solution 0.01 M of oleic acid in decane (hereafter called doped-decane) was used in the second Ketton micro core to permit a wettability alteration after drainage due to the sorption of the oleic acid on the solid surface yielding a mixed-wet system. The aqueous phase (or brine) was a solution of 7 wt% potassium iodide salt in deionized water. This setting was conducive to a clear contrast between the two fluid phases and the solid matrix. The core-flooding apparatus and the experimental procedure follow the same protocols described in [15,45], to which the reader is referred for additional details.

1. The samples were first flushed with CO2 to displace air. Brine was then injected to displace CO2, thus ensuring that the sample was fully saturated with brine.

2. Oil phase was injected with a flow rate corresponding to a low capillary number (Nc ' 10-7), for a total of 15 pore volumes (drainage).

3. For the mixed-wet system, the oil phase - doped-decane - was left in the core overnight to allow the wettability change to occur.

4. A total of 20 pore volumes of brine were then injected so that residual oil saturation was attained in the sample.

5. The cores were scanned with a Zeiss Xradia 500 Versa 3D X-ray microscope, obtaining the three-dimensional (3D) tomograms with high resolution (voxel size: 2 im for water-wet and 2.5 im for the mixed-wet system) that will be used to estimate contact angles.

2.2. Imaging process

Raw three-dimensional images were filtered with a non-local means filter to remove noise and obtain improved demarcation between phases [15,45,51,52]. The non-local means filter adopts as the grey-scale value of a certain voxel the average of those of neighbouring voxels. In addition, a value proportional to an index that measures the similarity between the neighbouring and the voxel taken into consideration is used to weight the value taken from each neighbour. This value, called the similarity value, ranges between 0 and 1 and should be carefully chosen as it is a trade-off between acceptable noise removal (for which it needs to be quite high) and a clear definition of the borders between the phases, that is crucial to correctly determine the contact angle (in this case, we need a low value of the similarity value). In the cases studied, the similarity value is chosen as 0.6, to obtain filtered images where rock, water and oil phases can be easily distinguished. A range of values between 0.45 and 0.7 all gave acceptable results. The image is then segmented, i.e., each voxel is assigned to either the brine, oil or solid phases. In recent years segmentation methods have improved significantly [53,54] to yield a level of accuracy which enables one to crisply distinguish between phases even in small corners of the complex pore space. The approach we used is the seeded watershed algorithm [55], that in our case was able to distinguish voxels in the image belonging to the different phases with an acceptable degree of accuracy [15]. When applying the seeded watershed algorithm, the minimum and maximum values for grey-scale intensity and gradient magnitude of each phase need to be specified; these values are quite case-dependent and the procedure is iterated until acceptable results are obtained with visual inspection. Fig. 1 illustrates, as an example, the raw, filtered and segmented images of a slice of 3D data. Calculations are performed using Avizo (Visual Sciences Group, www.vsg3d.com) and Matlab (MathWorks, www.mathworks.com) software.

Fig. 1. Illustration of the stages of data processing: the raw reconstructed data obtained from the scan (A) are first filtered (B) and then segmented (C; rock, oil and brine phases are respectively depicted in grey, green, and blue). Phase interfaces are well preserved during the process, even in small corners. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

2.3. Manual method for the estimation of wettability

Andrew et al. [45] introduced a method for the identification of the local contact angle which is based on the series of steps depicted in Fig. S1 of the Supplementary material and listed in the following:

1. Raw data are filtered and segmented, as described in Section 2.2.

2. The contact points between the three (fluid-fluid-solid) phases are automatically selected from the segmented three-dimensional images. The set of all these points defines the three-phase contact line.

3. For a given contact point identified at step 2, a slice of the grey-scale 3D image is taken along the direction perpendicular to the three-phase contact line.

4. The contact angle is then manually measured on this slice.

2.4. New automatic algorithm

The new automatic algorithm we propose differs from the procedure of Andrew et al. [45] in the way steps 3 and 4 are performed, as described in the following.

direction obtained with the procedure described above is compared against the plane determined manually according to step 3 of the manual method illustrated in Section 2.3 [45]. Differences between the results obtained are quantified through the angle / between the two normal directions. We assess the consistency of our results by considering the moving-average smoothing procedure of Eq. (3) with N = 4 or N = 6, i.e. using the coordinates of 4 or 6 points around the one to be replaced. Table 1 lists values of the mean angular difference (mean(/)) between the automatically and manually computed normal directions at 23 three-phase points and calculated by considering the three moving average options N = 0 (no smoothing), N = 4 and N = 6. We also list the resulting mean contact angle (mean(h)) and its standard deviation (Std dev) when completing the algorithm with these parameters. It can be seen that the best results are obtained using N = 4, with an associated mean angular difference of / = 17.3°. Previous work has shown that if the plane is tilted from the true value by less than 40°, then the measured contact angle is not affected significantly [45]. This is confirmed by results of mean(h) in the table, that are close to each other. For this reason, we do not expect our selection strategy of a perpendicular plane, which is hereafter performed by relying on N = 4, to cause a critical impact on the final estimate of the local contact angle.

2.4.1. Identification of planes locally normal to the three-phase contact line

Perpendicular planes of the kind described in Section 2.3 are identified using the three-dimensional coordinates of each three-phase contact-voxel, defined as a voxel characterized by the presence of two different neighbouring phases; an example is a water voxel having both rock and oil voxels as neighbours. By definition, one neighbour has at least one of its six faces in common with the considered voxel. As such, each three-phase physical contact point is identified by more than one voxel in the imaged system, i.e., at least one voxel for each phase. This observation is consistent with the result that the curve connecting all of the contact points is associated with a complex spatial pattern (see the red curve in Fig. S2 of the Supplementary material). The resulting three-phase contact curve is then smoothed through a moving average procedure (see the green and blue curves in Fig. S2 of the Supplementary material). The latter is a technique that is commonly used for signal smoothing [56]; it transforms the original data replacing each value xi with an average yi which is calculated through the data comprised within a given interval centred at xi:

yi = N^

where N is the number of points used to compute the moving average. The normal direction of each plane perpendicular to the three-phase contact curve, where the contact angle has to be estimated, is assessed with vectors vnk. The components of these are the differences of the coordinates of two subsequent points along the smoothed curve, i.e.,

vn,k = (Xk+1 - xk)i +(yk+1 - yk)j +(zk+1 - zk)k.

In Eq. (4) xk,yk and zk are the three Cartesian coordinates of the kth contact point, and vn k is a three-component vector with the desired direction, i.e., aligned with the three-phase contact curve. The normal unit vector nk is then found dividing v„yk by its modulus | vn k |:

„*-\v nj

To check the accuracy of this approach, some points are randomly selected along the curve, and the plane associated with the normal

2.4.2. Contact angle

As described in Section 2.4.1, 2D images are obtained as slices of segmented 3D images whose normal is aligned with the three-phase contact line delineated with the procedure described above. Fig. 2 depicts an example of one of these extracted 2D slices. The

Table 1

The set of three-phase contact points (three-phase contact line) (in yellow in Fig. 2) is smoothed using a moving average: using 4 points to compute the moving average, the mean angular difference (mean(/)) between perpendicular planes automatically and manually found is the lowest. mean(0) and Std dev are respectively the mean contact angle and standard deviation obtained applying the algorithm described in the main text to the dataset SSa of Fig. 5 with the three values listed for Npoints employed to smooth the three-phase contact line.

Npoints mean(/) (°) mean(S) (°) Std dev (°)

0 34.7 38.8 19

4 17.3 37.7 14

6 28.0 36.6 13

Fig. 2. Visualization of a 2D image obtained as a slice of a segmented image taken along the plane locally normal to the three-phase contact line. The latter is depicted in yellow. The slices are selected within cubic subvolumes of side b around each three-phase contact point. Green, blue and grey in the 2D slice respectively correspond to oil, brine and rock phases. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

estimated value of the local contact angle is then obtained from such slices according to the following protocol:

1. Rock surface pixels (circled in black in Fig. 3A) are identified. These are defined as rock-phase pixels with a neighbour which is not rock.

2. Since we aim at characterizing an effective contact angle, we consider a locally smoothed rock surface. We do so through linear regression (identified by a solid black line in Fig. 3A) performed upon rock surface pixels around the three-phase contact point. Results of a sensitivity analysis on the size of the region within which the above mentioned regression is performed is presented in Section 3.3.

3. Fluid-fluid interface pixels (highlighted in white in Fig. 3A) are defined as pixels belonging to one fluid phase with a neighbouring pixel classified as the other fluid phase. For the estimation of contact angle, we are interested in interface pixels located in a region around a three-phase contact point, hence we partition the 3D data in cubic subvolumes of side length b around each three phase contact point (refer to Fig. 2). We will fit the fluid-fluid interface to a circle, whose curvature is related to the local capillary pressure, Eq. (1). The appropriate selection of the size b of the subvolume is critical. If it is too small, then the local curvature of the interface cannot be well approximated, as we do not have enough points to fit the circle. On the other hand, the selection of a wide search region can violate the constraint of constant curvature. As we move further from the contact point we may deviate from the principal plane, as we can approach another rock surface, that can have a different orientation and hence a different curvature of the two-fluids interface. In our study we select b — 40 pixels (80 im for a voxel size of 2 im). Results of a sensitivity analysis on the effect of this choice are illustrated in Section 3.3.

4. A constant curvature for the interface between two fluids under equilibrium conditions is defined through the Young-Laplace Eq. (1). This constraint enables us to construct a circle (in white in Fig. 3A) constrained by the selected fluid-fluid interface points. The coordinates of the centre of the circle and its radius are fitted using the least square method [57].

5. The contact angle is computed as the angle between two straight lines: the tangent to the circle at the three-phase contact point (in yellow in Fig. 3A) and the rock surface line (in black in Fig. 3A). Since the contact angle is conventionally defined as measured through the denser phase, one needs to select one of the two angles that are defined by these two secant straight lines. We term the acute angle as a and construct the line normal to the circle at the three-phase contact point. The result of this procedure is depicted in Fig. 3B and C, where this normal line is dashed. Upon examining the pixels crossed by this line, we can see the following:

• If these pixels do not belong to the denser phase, the contact angle to be selected is acute, i.e. a, as shown in Fig. 3B.

• If some of these pixels belong to the denser phase, the contact angle is obtuse, i.e. p - a, as we can see in Fig. 3C.

2.5. Method of analysis of the results

Application of the approach described in Section 2.4 yields a set of spatially distributed values of contact angles. These are then interpreted in a probabilistic framework. A probability distribution of the contact angle is inferred and interpreted according to a truncated Gaussian model or a mixture thereof, depending on the experimental scenario investigated. The parameters of the interpretive models are estimated via a maximum likelihood (ML) approach [58,59]. A comparison between the sample probability distribution of contact angle obtained for two different ganglia from the same water-wet core is performed upon relying on an established metric such as the Kullback-Leibler divergence (DKL) [60-62]. This quantifies the discrepancy between two probability distributions P and Q. For continuous variables it is defined as:

l(P||Q ) = £ p(x) log m dx.

2.5.1. Interpretive model for the probability distribution of contact angles

The empirical probability density function, f, of contact angles for the water-wet system analyzed is interpreted through a truncated Gaussian model. The latter is defined on the support [a — 0°, b — 180°] and is characterized by a location and a scale parameter, respectively denoted as i and r, i.e.,

f (x; I, r, a, b) =

where /(x) is the normal Gaussian distribution function, and U(x) the Gaussian cumulative distribution function, respectively defined

Fig. 3. A: Representation of the 2D slice of 3D segmented images where contact angle Q is computed. Oil is green, brine blue and rock grey. A circle (in white) is fitted on coordinates of the interface between two fluids. B is the angle between the yellow tangent to the circle in the contact point and the black solid line representing the rock surface. The correct angle is selected as the one measured through the denser phase (water, in blue). The dashed line is normal to the fluid-fluid interface at the contact point and is used for the identification of the correct angle in panels B and C. In this case B — a in B and B — p - a for case C. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

exp(- 2n2

U(x) =

p22n .

-t2=2 dt.

For the dataset with altered wettability, we will show that the sample density of contact angles is well approximated by a mixture (with weight w) of two truncated Gaussian distributions, i.e.,

g(x; l1, ff1, l2, 02, w, a, b) = wf (x; i1, r, a, b)

+ (1 - w)f (x; 12, r2, a, b). (10)

where li and ri (i = 1 , 2) respectively are the location and scale parameters of the two distributions constituting the mixture.

3. Testing of the approach

3.1. Analysis of synthetic 2D images

We validate the method we propose by relying on simple 2D images where the contact angle is known a priori. All of these are characterized by the presence of a circular droplet of oil phase located at the centre of the image. For the first set of images a horizontal solid surface scans the droplet from top to bottom, defining contact angles ranging from 0° to 180°. The second set of images is characterized by a solid phase which is tilted of 135° from the horizontal. The two sets of figures are produced with different voxel sizes relative to the diameter d of the droplet, and the method is applied to each of them.

The results shown in Fig. S3 of the Supplementary material as scatter plots of actual versus calculated angles demonstrate that increasing the resolution produces results closer to the true value. We used images with increasing resolution from 14 to 112 pixels/d and compute the root mean square error RMSE for each set of results as

RMSE =y1 Xn=1 (y> - y)2, (11)

where n is the number of measurements, y and yt respectively being the estimated and the true (analytical) value of the contact angle.

As for the 3D images, we need to partition a region around each three-phase contact point. In these 2D cases, the side length b of the cubic subvolume is represented by the side of a squared region of interest (ROI). In our tests on 2D synthetic images we use a default value of ROI size corresponding to 50% of d, and the results show an RMSE lower than 4° for resolutions with at least 28 pixels/ d. From Fig. S4 of the Supplementary material it is clear that increasing resolution and ROI dimension, RMSE lowers down to 0.5° (for resolution higher than 56 pixels/d and ROI size higher than 50% of d).

3.2. Comparison between manual and automatic methods on the same slices, from images of real systems

The tomograms obtained from the X-ray micro-CT scan of the Ketton-brine-decane system with the procedure described in Section 2.1 are post-processed with the Avizo software: the images are filtered and segmented with the techniques mentioned in Section 2.2 and illustrated in Fig. 1. One ganglion of oil is selected from the water-wet dataset and the manual method described in Section 2.3 is used to obtain contact angle estimations in a total of 44 of its three-phase contact points, on grey-scale images. The automatic algorithm is applied on slices of segmented images with the same orientation as those used for the manual method, on the very same points. A comparison between contact angle values obtained through the manual and automated methods is rendered by the scatterplot depicted in Fig. 4. Overall, the agreement between the results obtained with the two methods can be considered as mutually consistent.

3.3. Sensitivity analysis

Our automatic algorithm requires tuning two key parameters:

u 20 40 60 80

Manual estimate [°]

Fig. 4. Scatter plot comparing manual and automatic measurements of contact angles associated with the water-wet system considered.

1. The length ' of the line representing fluids-rock interface.

2. The dimension b of the side of the cubic subvolume considered within which the fluid-fluid interface is represented by a circle through a fitting procedure against the experimental data.

An analysis on the impact of b and ' on the quality of the results is performed by (a) considering the results depicted in Fig. 4 and associated with the manual method as reference values, (b) sampling b and ' on a regular grid in the space within which these parameters are defined, and (c) computing the resulting root mean square difference (RMSD) between manual and automatic measurements, with the same Eq. (11) used for RMSE, considering yt as the manual measurements. The results of this analysis depicted in Fig. S5 of the Supplementary material suggest that considering b ranging between 30 and 50 pixels and a value of ' ranging between 100 and 400 pixels yields the lowest values of RMSD. For our applications, we use b = 40 and ' = 200.

3.4. Limitations of the method

As already discussed in Section 3.1, the resolution of the imaged rock affects the accuracy of the results. The comparison with the manual method illustrated in the previous section demonstrates that, in our case, the results are acceptable with a 2 im voxel size. It is currently possible to produce micro-CT images on mm-sized samples with a 1.5 im voxel size [46], and according to the tests (refer to Fig. S4 of the Supporting material) this would yield in even more accurate results (lower RMSE). In addition, a lower radius of curvature may rise RMSE, increasing the ratio pixels/d. The curvature of the interface is influenced by several factors, including the size of the throats and the wettability of the system. Further tests of the method with different rock types and wettability grades would give indications about the robustness of the new algorithm.

The quality of the image segmentation also influences the results. The seeded watershed segmentation technique adopted in our work gives accurate results, as long as the interfaces are well preserved (Fig. 1). Due to the key role of segmentation in many fields, improvements are continuously being studied, so that future work will likely have access to higher quality segmented images [63-65].

4. Results

4.1. Water-wet system

We consider here the water-wet system illustrated in Sections 2.1 and 2.2. The two oil ganglia depicted in Fig. 5 are identified

Fig. 5. Trapped oil ganglia SSa (C) and SSb (B) isolated from the water-wet system with the procedure described in Sections 2.1 and 2.2 and employed for the demonstration of the new automated method. In panels B and C, the three-phase contact line is depicted in yellow. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

and isolated for the application of our automated methodology. These are characterized by different shape, dimension and topolog-ical structure. Ganglion SSa (Fig. 5C) is relatively regular, as opposed to the bigger ganglion SSb (Fig. 5B), which is characterized by a complex three-dimensional pattern.

Application of the automatic algorithm leads to the identification of the spatial distribution of the contact angle at a set of points located along the contact lines highlighted in yellow in Fig. 5B and C. A total of 1652 and 4025 points are selected along the contact lines of samples SSa and SSb, respectively. As an example of the results obtained, Fig. 6 depicts a three-dimensional view of the contact angle values obtained from the analysis of sample SSa. These values are then interpreted through the truncated Gaussian model, Eq. (7), yielding ML estimates of the distribution parameters as well as associated confidence intervals. The modelled probability densities associated with the two ganglia are compared through the Kullback-Leibler divergence, Eq. (6). The very small

Fig. 6. Three-dimensional view of the contact angle values evaluated by our automated algorithm at Npoints - 1652 points identified along the contact lines of sample SSa.

Table 2

Samples considered for the application, number of points (Npoi„ts) at which the contact angle is assessed, ML estimates of parameters (Par) 1 and r (see Eq. (7)) and w, 11, r1,12 and r2 (see Eq. (10)). The width of the 95% confidence interval associated with the estimate of each parameter is also listed.

Dataset Npoints Par Value Confidence interval

SSa 1652 1 37.5° [36.9° 38.2°]

r 1 3.7° [13.2° 14.2°]

SSb 4025 1 42.2° [41.7° 42.6°]

r 14.6° [14.3° 14.9°]

SSa U SSb 5668 1 40.8° [40.4° 41.2°]

r 14.5° [14.2° 14.8°]

Mixed-wet 2657 w 0.75 [0.73 0.77]

57.5° [56.4° 58.5°]

r 21.5° [20.8° 22.2°]

85.5° [85.3° 85.6°]

r 2 1.84° [1.71° 1.97°]

value of the associated Kullback-Leibler divergence (i.e., Dm — 0.06) suggests that the two samples are very similar from a probabilistic standpoint. The results of the analyses are collected in Table 2. The latter lists the samples considered for the application, the number of points (Npoints) at which the contact angle is assessed, as well as the ML estimates 1 and r (see Eq. (7)) and the associated 95% confidence intervals. Results obtained by interpreting the empirical probability density function obtained by jointly considering the two datasets are also listed. Fig. 7A complements our results by depicting the empirical and modelled density of the contact angle sample obtained by jointly considering measurements associated with samples SSa and SSb.

Fig. 7. Final result of contact angle distribution for the water-wet case (A), merging data from the two different ganglia SSa and SSb, and for the mixed-wet case (B). Dots represent empirical probability density distribution values, while solid lines are defined by function of Eq. (7) and (10) for the water- (A) and mixed-wet (B) system respectively.

4.2. Mixed-wet system

The automated method is also applied to the analysis of a ganglion extracted from the dataset obtained from the wettability-altered Ketton carbonate sample. The affected contact angle is here estimated at a set of 2657 contact points. Fig. 7B depicts the empirical and modelled density of the resulting contact angles. In this setting we observe the clear effect of ageing, as the mixed condition is represented by a bi-modal distribution. The latter is well interpreted by a mixture of two Gaussian distributions, whose estimated parameters are listed in Table 2 together with the associated 95% confidence intervals. Here we find a wider distribution of contact angles with values both less than and greater than 90°, meaning that parts of the pore space are water-wet while others are slightly oil-wet: this is what we mean by a mixed-wettability. There is not a single contact angle with some variation around a typical value, but two peaks for areas that have been affected to different degrees by the ageing procedure.

5. Conclusions and future work

The new algorithm for the estimation of contact angle presented in this work, based on the constraint of constant curvature defined by Young-Laplace Eq. (1), allows the automatic computation of thousands of reliable estimates of contact angle from a rock core. The method is validated with synthetic images and through comparison with results obtained using the manual method proposed by Andrew et al. [45]. Manual and automatic estimates are in good agreement.

Several methods available in the literature infer wettability indirectly from capillary pressure curves [14,32,33] or with capillary rise and Whielmy plate methods and water drop penetration test [7,34-37], without the direct visualization of contact angle. Direct measurements of contact angle are widely performed with sessile drop method on simple settings [11,13,26-31], while the only direct assessment of in situ contact angle within the complex pore space of natural rocks was a manual method proposed by Andrew et al. [45]. The new automatic algorithm adds a step forward with an automation of the process using the physical constraint of constant curvature of the fluid-fluid interface.

This method is applied on two different ganglia of a water-wet Ketton-brine-decane system, obtaining effective contact angle distributions modelled as truncated Gaussian probability density functions (Eq. (7)) with estimated location parameters l = 37.5° and l = 42.2° respectively. These are slightly lower than a mean value of 47° found for a similar system in [15]. Lastly, a ganglion from a mixed-wet Ketton-brine-doped-decane system is considered, obtaining the more spread distribution of Fig. 7, shifted towards higher values of effective contact angle [15,66,67]: this sample has contact angles that are both less than and greater than 90°, indicating the presence of both water-wet and oil-wet surfaces.

These results lead us to the following conclusions:

1. The automatic algorithm is a reliable method to estimate the distribution of in situ effective contact angle.

2. This new approach returns a number of estimates that is large enough for an interpretation of results in a probabilistic framework.

3. The probability density distribution of effective contact angles of a water-wet system can be modelled as a truncated Gaussian density function, characterized by a location parameter l and a scale parameter r.

4. The data from two ganglia of a Ketton-brine-decane water-wet system can be merged together with a resulting overall

distribution of effective contact angle with location and scale parameter l = 40.8° and r = 14.5° respectively. 5. The use of decane doped with an oleic acid can alter the wettability through ageing of the rock, producing a different distribution of contact angle, bimodal and shifted towards higher values, that can be modelled with the weighted mix of two truncated Gaussian density functions of Eq. (10).

Future work is planned to further test the robustness of the algorithm and extend this method to the study of different rock types, rocks with spatially variable wettability and different mineralogies, and to measure contact angles during a displacement.

Acknowledgements

The authors thank Dr. Qingyang Lin at Imperial College London for the useful help with TCL encoding.

Appendix A. Supplementary material

Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.jcis.2017.02.005.

References

[1] M.J. Blunt, Multiphase Flow in Permeable Media: A Pore-Scale Perspective, Cambridge University Press, 2017.

[2] G.F. Pinder, W.G. Gray, Essentials of Multiphase Flow and Transport in Porous Media, John Wiley & Sons, Inc., 2008.

[3] J. Bear, C. Tsang, G. de Marsily, Flow and Contaminant Transport in Fractured Rocks, Academic Press, 1993.

[4] W. Plug, J. Bruining, Capillary pressure for the sand-CO2-water system under various pressure conditions. application to CO2 sequestration, Adv. Water Resour. 30 (2007) 2339-2353.

[5] I. Gaus, Role and impact of CO2-rock interactions during CO2 storage in sedimentary rocks, Int. J. Greenhouse Gas Contr. 4 (2010) 73-89.

[6] L. Moghadasi, A. Guadagnini, F. Inzoli, D. Colapietro, M. Bartosek, D. Renna, Laboratory-scale investigation of two-phase relative permeability, Proc. Environ. Sci. 25 (2015) 166-174.

[7] J.T. Gostick, M.A. Ioannidis, M.W. Fowler, M.D. Pritzker, Wettability and capillary behavior of fibrous gas diffusion media for polymer electrolyte membrane fuel cells, J. Power Sources 194 (2009) 433-444.

[8] M.D. Simoni, J. Carrera, X. Sanchez-Vila, A. Guadagnini, A procedure for the solution of multicomponent reactive transport problems, Water Resour. Res. 41 (2005) W11410.

[9] T. Young, An essay on the cohesion of fluids, Philos. Trans. 95 (1805) 65-87.

[10] N.R. Morrow, Wettability and its effect on oil recovery, J. Petrol. Technol. 42 (1990) 1476-1484.

[11] S. Iglauer, C. Pentland, A. Busch, CO2 wettability of seal and reservoir rocks and the implications for carbon geo-sequestration, Water Resour. Res. 51 (2015) 729-774.

[12] S. Saraji, L. Goual, M. Piri, Adsorption of asphaltenes in porous media under flow conditions, Energy Fuels 24 (2010) 6009-6017.

[13] D.N. Espinoza, J. Santamarina, Water-CO2-mineral systems: interfacial tension, contact angle, and diffusion-implications to CO2 geological storage, Water Resour. Res. 46 (2010) 729-774.

[14] A.L. Herring, A. Sheppard, L. Andersson, D. Wildenschil, Impact of wettability alteration on 3D non-wetting phase trapping and transport, Int. J. Greenhouse Gas Contr. 46 (2016) 175-186.

[15] K. Singh, B. Bijeljic, M.J. Blunt, Imaging of oil layers, curvature and contact angle in a mixed-wet and a water-wet carbonate rock, Water Resour. Res. 52 (2016) 1716-1728.

[16] Y. Tanino, B. Akamairo, M. Christensen, S.A. Bowden, Impact of displacement rate on waterflood oil recovery under mixed-wet conditions, Soc. Core Anal. 11 (2015).

[17] H.S. Rabbani, V. Joekar-Niasar, N. Shokri, Effects of intermediate wettability on entry capillary pressure in angular pores, J. Colloid Interf. Sci. 473 (2016) 3443.

[18] H. Geistlinger, I. Ataei-Dadavi, Influence of the heterogeneous wettability on capillary trapping in glass-beads monolayers: comparison between experiments and the invasion percolation theory, J. Colloid Interf. Sci. 459 (2015) 230-240.

[19] S. Iglauer, M. Fern0, P. Shearing, M. Blunt, Comparison of residual oil cluster size distribution, morphology and saturation in oil-wet and water-wet sandstone, J. Colloid Interf. Sci. 375 (2012) 187-192.

[20] V. Gurau, M.J. Bluemle, E.S. De Castro, Y.-M. Tsou, J. Adin Mann Jr., T.A. Zawodzinski Jr., Characterization of transport properties in gas diffusion layers for proton exchange membrane fuel cells: 1. Wettability (internal contact

angle to water and surface energy of {GDL} fibers), J. Power Sources 160 (2006) 1156-1162.

[21] K. Jiao, B. Zhou, Effects of electrode wettabilities on liquid water behaviours in {PEM} fuel cell cathode, J. Power Sources 175 (2008) 106-119.

[22] M. Rucker, S. Berg, R.T. Armstrong,A. Georgiadis, H. Ott, A. Schwing, R. Neiteler, N. Brussee, A. Makurat, L. Leu, M. Wolf, F. Khan, F. Enzmann, M. Kersten, From connected pathway flow to ganglion dynamics, Water Resour. Res. 42 (2015) 3888-3894.

[23] M. Andrew, B. Bijeljic, M.J. Blunt, Pore-by-pore capillary pressure measurements using X-ray microtomography at reservoir conditions: curvature, snap-off, and remobilization of residual CO2, Water Resour. Res. 50 (2014) 8760-8774.

[24] N. Alyafei, M.J. Blunt, The effect of wettability on capillary trapping in carbonates, Adv. Water Resour. 90 (2016) 36-50.

[25] P.K. Das, X. Li, Z. Xie, Z.-S. Liu, Effects of catalyst layer structure and wettability on liquid water transport in polymer electrolyte membrane fuel cell, Int. J. Energy Res. 35 (2011) 1325-1339.

[26] N.M. Dingle, M.T. Harris, A robust algorithm for the simultaneous parameter estimation of interfacial tension and contact angle from sessile drop profiles, J. Colloid Interf. Sci. 286 (2005) 670-680.

[27] M. Schmitt, K. Groß, J. Grub, F. Heib, Detailed statistical contact angle analyses; slow moving drops on inclining silicon-oxide surfaces, J. Colloid Interf. Sci. 447

(2015) 229-239.

[28] A. Promraksa, L.-J. Chen, Modeling contact angle hysteresis of a liquid droplet sitting on a cosine wave-like pattern surface, J. Colloid Interf. Sci. 453849 (2012) 172-181.

[29] M. Arif, A.Z. Al-Yaseri, A. Barifcani, M. Lebedev, S. Iglauer, Impact of pressure and temperature on CO2-brine-mica contact angles and CO2-brine interfacial tension: implications for carbon geo-sequestration, J. Colloid Interf. Sci. 462

(2016) 208-215.

[30] M. Santini, M. Guilizzoni, S. Fest-Santini, X-ray computed microtomography for drop shape analysis and contact angle measurement, J. Colloid Interf. Sci. 409 (2013) 204-210.

[31] J. Bachmann, R. Horton, R Van Der Ploeg, S. Woche, Modified sessile drop method for assessing initial soil-water contact angle of sandy soil, Soil Sci. Soc. Am. J. 64 (2000) 564-567.

[32] E. Amott, Observations relating to the wettability of porous rock, Petrol. Trans. AIME 216 (1959) 156-162.

[33] E. Donaldson, R. Thomas, P. Lorenz, Wettability determination and its effect on recovery efficiency, SPE J. 9 (1969) 13-20.

[34] L. Wilhelmy, Ueber die abhängigkeit der capillaritäts-constanten des alkohols von substanz und gestalt des benetzten festen körpers, Ann. Phys. 195 (1863) 177-217.

[35] L. Dekker, C. Ritsema, Wetting patterns and moisture variability in water repellent dutch soils, J. Hydrol. 231 (2000) 148-164.

[36] A. Siebold, A. Walliser, M. Nardin, M. Oppliger, J. Schultz, Capillary rise for thermodynamic characterization of solid particle surface, J. Colloid Interf. Sci. 186 (1997) 60-70.

[37] J. Bachmann, S. Woche, M.-O. Goebel, M. Kirkham, R. Horton, Extended methodology for determining wetting properties of porous media, Water Resour. Res. 39 (2003).

[38] M. Robin, E. Rosenberg, O. Fassi-Fihri, Wettability studies at the pore level: a new approach by the use of cryo-scanning electron microscopy, SPE Formation Eval. 10 (1995) 11-20.

[39] J. Schmatz, J.L. Urai, S. Berg, H. Ott, Nanoscale imaging of pore-scale fluid-fluid-solid contacts in sandstone, Geophys. Res. Lett. 42 (2015) 2189-2195.

[40] R. Combes, M. Robin, G. Blavier, M. Aldan, F. Degreve, Visualization of imbibition in porous media by environmental scanning electron microscopy: application to reservoir rocks, J. Petrol. Sci. Eng. 20 (1998) 133-139.

[41] T.M. Okasha, J.J. Funk, H.N. Rashidi, Fifty years of wettability measurements in the arab-D carbonate reservoir, in: SPE 105114, Proceedings of the SPE Middle East Oil and Gas Show and Conference, Manama, Kingdom of Bahrain, 11-14 March, 2007.

[42] A.Q. Raeini, M.J. Blunt, B. Bijeljic, Direct simulations of two-phase flow on micro-ct images of porous media and upscaling of pore-scale forces, Adv. Water Resour. 74 (2014) 116-126.

[43] M.J. Blunt, B. Bijeljic, H. Dong, O. Gharbi, S. Iglauer, P. Mostaghimi, A. Paluszny, C. Pentland, Pore-scale imaging and modelling, Adv. Water Resour. 51 (2013) 197-216.

[44] A.Q. Raeini, B. Bijeljic, M.J. Blunt, Modelling capillary trapping using finite-volume simulation of two-phase flow directly on micro-ct images, Adv. Water Resour. 83 (2015) 102-110.

[45] M. Andrew, B. Bijeljic, M.J. Blunt, Pore-scale contact angle measurements at reservoir conditions using X-ray microtomography, Adv. Water Resour. 68 (2014) 24-31.

[46] M. Khishvand, A. Alizadeh, M. Piri, In-situ characterization of wettability and pore-scale displacements during two- and three-phase flow in natural porous media, Adv. Water Resour. 97 (2016) 279-298.

[47] X. Li, X. Fan, Effect of {CO2} phase on contact angle in oil-wet and water-wet pores, Int. J. Greenhouse Gas Contr. 36 (2015) 106-113.

[48] A. Aghaei, M. Piri, Direct pore-to-core up-scaling of displacement processes: dynamic pore network modeling and experimentation, J. Hydrol. 522 (2015) 488-509.

[49] C. Chen, J. Wan, W. Li, Y. Song, Water contact angles on quartz surfaces under supercritical {CO2} sequestration conditions: experimental and molecular dynamics simulation studies, Int. J. Greenhouse Gas Contr. 42 (2015) 655-665.

[50] K.A. Klise, D. Moriarty, H. Yoon, Z. Karpyn, Automated contact angle estimation for three-dimensional X-ray microtomography data, Adv. Water Resour. 95 (2016) 152-160.

[51] A. Buades, B. Coll, J. More, Nonlocal image and movie denoising, Int. J. Comput. Vis. 76 (2008) 123-139.

[52] M. Andrew, H. Menke, M.J. Blunt, B. Bijeljic, The imaging of dynamic multiphase fluid flow using synchrotron-based X-ray microtomography at reservoir conditions, Transp. Porous Media 110 (2015) 1-24.

[53] H. Ng, S. Ong, K. Foong, P. Goh, W. Nowinski, Medical image segmentation using k-means clustering and improved watershed algorithm, IEEE Xplore (2006) 61-65.

[54] G. Hamarneh, X. Li, Watershed segmentation using prior shape and appearance knowledge, Image Vis. Comput. 27 (2009) 59-68.

[55] A.C. Jones, C.H. Arns, A.P. Sheppard, D.W. Hutmacher, B. Milthorpe, M.A. Knackstedt, Assessment of bone ingrowth into porous biomaterials using MICRO-CT, Biomaterials 28 (2007) 2491-2504 (Imaging Techniques for Biomaterials Characterization).

[56] D.M.B. Vandeginste, S. Deming, Y. Michotte, L. Kaufman (Eds.), Chemometrics: A Textbook, Data Handling in Science and Technology, 1st repr. 1990 ed., vol. 2, Elsevier Science, 1988.

[57] V. Pratt, Direct least-squares fitting of algebraic surfaces, ACM SIGGRAPH Comput. Graph. 21 (1987).

[58] S.M. Ross, Introduction to Probability and Statistics for Engineers and Scientists, third ed., Academic Press, 2004.

[59] J.P.M. de Sa, Applied Statistics Using SPSS STATISTICA MATLAB and R, second ed., Springer, 2007.

[60] S. Kullback, Information Theory and Statistics, Dover Publications, 1987.

[61] I. Sraj, O.P. Le Maitre, O.M. Knioc, I. Hoteita, Coordinate transformation and polynomial chaos for the bayesian inference of a gaussian process with parametrized prior covariance function, Comput. Methods Appl. Mech. Eng. 298 (2016) 205-228.

[62] M. Panzeri, M. Riva, A. Guadagnini, S. Neuman, Theory and generation of conditional, scalable sub-gaussian random fields, Water Resour. Res. 52 (2016) 1746-1761.

[63] O.P. Verma, M. Hanmandlu, S. Susan, M. Kulkarni, P.K. Jain, A simple single seeded region growing algorithm for color image segmentation using adaptive thresholding, in: 2011 International Conference on Communication Systems and Network Technologies (CSNT), IEEE, pp. 500-503.

[64] J. He, C.-S. Kim, C.-C.J. Kuo, Interactive image segmentation techniques, in: Interactive Segmentation Techniques, Springer, 2014, pp. 17-62.

[65] J. He, C.-S. Kim, C.-C.J. Kuo, Conclusion and future work, in: Interactive Segmentation Techniques, Springer, 2014, pp. 75-76.

[66] W. Anderson, 1987, Wettability literature survey - Part 1 to part 6, J. Pet. Technol. 4 (1986) 1125-1144.

[67] J. Buckley, Y. Liu, S. Monsterleet, Mechanisms of wetting alteration by crude oils, SPE J. 3 (1998) 54-61.