A thermodynamically consistent characterization of wettability in porous media using high-resolution imaging

a r t i c l e i n f o

provide useful information on wettability on a pore-by-pore basis, there are some uncertainties in their interpretation. Firstly, the contact angle is defined at the three-phase contact line, which is where the uncertainties in image quality and segmentation are greatest [4,6]. Secondly, the angle is measured when the fluids are at rest and hence does not necessary represent the value for a displacement because of contact angle hysteresis [11]. The time-scale of pore-scale displacement is of the order of micro-to milli-seconds and so even the fastest synchrotron imaging-with a time resolution of a few seconds-cannot, currently, capture the dynamics of fluid movement in the pore space [12,13]. Finally, a distribution of angles is measured; however, in many places the interface is pinned and will not move during a displacement. It is not possible to determine a priori which three-phase contact points instead are about to move. As a result h g does not necessarily indicate the contact angle that controls displacement, or that should be used in a pore-scale modelling study to reproduce the behaviour. Indeed it has been shown that to match experiments where water displaced oil in a mixed-wet rock, the use of spatially-distributed values of h g at equilibrium after waterflooding was insufficient to model the behaviour accurately: instead, the contact angle in more oil-wet regions had to be increased, accounting for contact angle hysteresis, to predict the movement and connectivity of the oil phases accurately [14].
An alternative approach is to consider conservation of energy [11,15] while finding the contact angle from a generalized Young equation [16][17][18][19][20][21]. We consider an arbitrary finite control volume and calculate changes in the Helmholtz free energy. We do not consider the Gibbs free energy which is used to compute states of thermodynamic equilibrium, since for a displacement involving two immiscible fluids, we will consider the work done in a system of fixed volume and temperature [11]. Then the thermodynamically-consistent contact angle h t measured through phase 1 is defined as: where Áe is the change in surface free energy per unit area when phase 1 displaces phase 2. r is the interfacial tension between these two phases while r 1s and r 2s represent the interfacial tensions between the surface, and phases 1 and 2 respectively. The contact angle is an effective value, since the interfacial tensions represent energies per unit area, and this area-as we show later-is measured at a certain resolution over a rough surface. h t is simply the contact angle that provides the observed change in interfacial energy and does not represent the angle measured at the atomic scale.
The interfacial energy of the system-a porous medium containing two fluid phases-or any sub-set of the system, is given by: where A is an area, while the subscripts 1; 2 and s refer to phase 1, phase 2 and the solid respectively. The total solid area A s ¼ A 1s þ A 2s is fixed. Hence Eq. (2) becomes: Then from Eq. (1) we may write Eq. (3) as: Now consider the difference in energy for the transition from one local state of equilibrium to another-this can be an infinitesimal change, or a finite difference, but where the beginning and end points represent equilibria with the fluid interfaces at rest: For an isothermal process the Helmholtz free energy difference, DF, between these states should be zero [11,22,23]: where there is a pressure difference between phase 2 and phase 1, P c ¼ P 2 À P 1 , and an increase in volume DV 1 of phase 1. DV ¼ /VDS 1 where / is the porosity, S is the saturation and V is the volume of the region of the porous medium under consideration. We can find the capillary pressure from the interfacial curvature j using the Young-Laplace equation [24]: where r 1 and r 2 are the principal radii of curvature. Thus Eq. (6), substituting Eq. (5), becomes: We can define specific surface areas a ¼ A=V and rearrange Eq. (8) to obtain: Fig. 1 illustrates our approach for a bundle of capillary tubes. On rough or altered-wettability surfaces, the meniscus between the two phases can be pinned: the apparent geometric angle is not the angle necessary for displacement. In contrast, the energy balance approach only considers interfaces where displacement has occurred and ignores the pinned values, which may be larger or smaller than h t .
This treatment follows the work of Morrow [11], but is distinct conceptually from previous discussions of energy balance [22,23]. We assume that the free energy change between two equilibrium states at the same capillary pressure is zero. Hence the change in capillary pressure, or curvature, must be small during the displacement. h t is therefore the contact angle to use in, for instance, many pore-scale modelling studies that track a series of local equilibrium states where transient effects and viscous dissipation are ignored; hence the energy balance enshrined in Eq. (9) is assumed [13,25,26]. Specifically, it is the angle to use in Mayer-Stowe-Princen calculations of threshold capillary pressure [27][28][29][30]: this is the angle at the onset of displacement and is different from the angle for fluids at rest. The derivation is strictly correct for infinitesimal, reversible, changes in fluid configuration. In this paper, however, we will consider transitions between two states with a finite change in saturation. To do this we make two assumptions. Firstly, that the change in capillary pressure is small: Phase 2 Phase 1 Displacement θ t θ g = 66°P hase 1 Phase 2 θ g Fig. 1. An illustration of the thermodynamic approach to calculating contact angle for a bundle of capillary tubes. The bold lines indicate rough or altered-wettability surfaces where the contact angle for displacement is larger than the current, pinned or geometric value, hg whose values are indicated: these interfaces cannot move until the pressure in phase 1 is increased. The thermodynamic angle, ht , (indicated in bold) calculated using Eq. (9), only considers those regions of the pore space where a displacement has occurred and ignores the pinned values. In contrast, measuring contact angle directly on an image may over or under-estimate the contact angle, dependent on the nature of the pinned interfaces. The particular values of contact angle in this example correspond to the water-wet example shown later in the paper. the value of j used in Eq. (9) will be the average of j measured at the two states (before and after the displacement). Secondly, we ignore viscous dissipation during irreversible changes in fluid configuration. In experiments in bead packs it has been shown that up to 90% of the injection energy is converted to interfacial area [11]. For more complex porous media in drainage the conversion efficiency can be lower [22], but this was under the assumption that the contact angle was 0. Accounting for viscous dissipation would lower the calculated value of cos h t , making the system apparently less wetting to phase 1 (larger h t ). Despite this, for mixed-wet systems studied later, we demonstrate that the calculated value of h t is larger than the geometric angle, indicating that viscous dissipation, in this case, is less significant than the errors introduced through measuring angles at pinned interfaces.
The key insight of this work is that we can now measure interfacial areas [31][32][33], curvatures [34,35] and saturations from three-dimensional images to determine h t directly. Furthermore, unlike h g which is a point property, h t is the average contact angle that provides the correct energy balance for a displacement in some finite region of the pore space: it can be defined in any region where we see different fluid configurations containing identifiable menisci from one image to the next. Finally, note that the interfacial areas are measured at the resolution of the image and hence will ignore small-scale surface roughness, or menisci in this roughness [36]. In this respect, we compute an effective contact angle, valid to describe displacement at the scale at which we image, and not, as mentioned above, a chemically-consistent value at the atomic level.

Imaging experiments
We illustrate our approach by measuring h g and h t for two datasets. In both cases we consider the displacement of oil (phase 2) by brine (phase 1) through two cylindrical samples of Bentheimer sandstone arranged in series. These cores had a diameter of 6.1 mm and a combined length of 51.6 mm. Full details of the experiments are presented elsewhere [37][38][39]. The brine contained 3.5% by weight of KI-the addition of iodine aids the image contrast. The oil was decalin. The interfacial tension between brine and decalin for both experiments was measured to be 51:5 AE 1:6 mN/m by the pendant drop method [39]. Initially the samples were fully saturated with brine. In the water-wet case, oil was injected at a high flow rate to establish an initial brine saturation of approximately 0.14. In the mixed-wet experiment, crude oil from a producing reservoir was forced into the samples using a centrifuge leaving a brine saturation of approximately 0.11. The cores were then left at a temperature of 80 C and a pressure of 3 MPa for four weeks to allow direct contact between crude oil and the solid surface. This alters the wetting state of the rock in a process known as ageing [40,41]. After ageing, crude oil was replaced by decalin. For both experiments oil and water were injected simultaneously through the rock at a low flow rate (the capillary number Ca ¼ lq=r % 10 À7 where l is the brine viscosity and q the total Darcy flow rate). A fixed water fractional flow, f w (the volumetric flow rate of brine divided by the total flow rate) was imposed. f w was increased in a series of steps from 0 to 1 while keeping the total flow rate of brine plus oil constant. The flow was maintained for several hours until the pressure differential across the cores stabilized. This is the procedure used in a steady-state relative permeability measurement [42,37]. The experiments were conducted at ambient temperature, 20 C, and a fluid pressure of 3 MPa. During the final hour of each fractional flow step, the samples were imaged in three dimensions using a Zeiss Versa 510 X-ray microscope: a central region of size 1000 Â 1000 Â 3000 cubic voxels of size 3.58 lm was analysed in this study.
We will quantify the wetting state of the sample that had been in contact with crude oil, which we call the mixed-wet system in what follows, since-as we show below-we measured values of h g both above and below 90 [43]. The other sample had not been in contact with crude oil and remained water-wet: h g < 90 .
Every image was first segmented into brine, oil and solid phases using standard methods: in this case a seeded watershed algorithm was used [44]. The brine saturation was found from the number of brine voxels divided by the total number of oil and brine voxels. The interfaces defining A 1s ; A 2s and A 12 were identified and voxelization artefacts were removed using a boundary-preserving surface smoothing algorithm with 150 iterations [45]. The interfacial curvature j and its average value j were computed using commercial image analysis software, Avizo. Any fluid-fluid surfaces within three voxels of the solid were ignored in the estimation of curvature to avoid uncertainties in the segmentation close to the three-phase contact line. We also only considered interfaces where the oil phase was connected across the imaged volume; we ignored trapped phases in the computation of curvature, but not for the calculation of area or saturation. The reason for this is that a macroscopic capillary pressure is only defined for continuous phases: once trapped, the fluids can rearrange in the pore space in a new position of local capillary equilibrium with a pressure very different from the prevailing capillary pressure in the continuous fluids [48,13]. Finally, the geometric contact angle h g was estimated at every point on the three-phase contact line using an automated algorithm [3]: the oil/brine interface was further adjusted in this step to impose a locally constant curvature. Further details of the image analysis approach and a quantification of measurement uncertainties can be found in [38]. We have performed a detailed analysis of uncertainty in the estimation of curvature against accurate computations from direct pore-scale modelling: in this case the uncertainty was less than 20% [46]. Then h t was computed using Eq. (9) from two successive images and using the mean curvatures. Further details on the uncertainty of the measurements are provided in the next section.

Water-wet dataset
We first perform the analysis for the water-wet dataset. In this case we consider changes from the first non-zero fractional flow, where both connected oil and water can be imaged, f w ¼ 0:15, to the penultimate step with f w ¼ 0:85. The results are presented in Table 1. We do not use the final fractional flow, f w ¼ 1, since here all the oil was trapped and hence we cannot determine a continuous phase capillary pressure. For this experiment we have confirmed the validity of the curvature measurements, since they provide an estimate of capillary pressure (using Eq. (7) and the interfacial tension) which is very close to that measured independently on a larger sample of the same sandstone using the porous plate method [47,38].
h t was found from the difference between the measurements at the fractional flow shown and the preceding value: the j used was the average of j at these two fractional flows. Note that the j values are large in comparison to the interfacial areas. The expression for h t in Eq. (9) is therefore dominated by the first term on the right-hand-side; since the curvature is positive we predict water-wet conditions. We find h t ¼ 47:8 AE 5 as opposed to h g ¼ 66:4 , the average of the two values shown in Table 1. See the end of this section for a discussion of the uncertainty in h t . The geometric angle appears to be rather high, implying weakly water-wet conditions, even though the experiment was performed on a clean sandstone with a refined oil. However, as mentioned in the introduction, the estimated value could be affected by segmentation errors and, in any event, on a rough surface, some interfaces may be pinned with a large effective angle, while displacement occurs at other interfaces where the contact angles are lower, see Fig. 1. This is evident in Fig. 2 which shows the distribution of measured contact angles at equilibrium. A wide range of contact angle is estimated locally and it is not obvious what angle in this range determines displacement. In contrast, the thermodynamic angle is considerably lower and appears more plausible.
To make the comparison of thermodynamic and geometric angles more quantitative, we would like to know which angle better represents displacement in a pore-scale model: that is which angle should be input to predict the macroscopic properties of the rock correctly? In a recent pore-scale network modelling study, a comparison between predicted and measured relative permeabil-ities for water-wet Bentheimer sandstone, the system studied here, was presented. The best match to the data was obtained using an advancing contact angle of 46°, very close to the value we calculate for h t [26]. On the other hand, a contact angle of h g % 66 would provide an incorrect quantification of pore filling, specifically underestimating the degreee of snap-off and trapping, leading to poor predictions [26].

Mixed-wet dataset
In mixed-wet media the assignment of contact angle is even less certain. The measured values of h g shown in Fig. 2 again exhibit a wide range of values with an average of 80:8 . This is surprising, showing a rather modest change from water-wet conditions and indicating a system that is still weakly water-wet. This contrasts with the measured curvatures, see Table 2, which are negative, meaning that the capillary pressure is negative, Eq. (7), consistent with contact angles above 90°.
Before discussing the contact angle measurements themselves for the mixed-wet case, Table 2, we first place the work in context by showing portions of the three-dimensional images to illustrate key differences between the behaviour in the water-wet and mixed-wet systems. Fig. 3 compares the fluid configurations for the two cases. In the water-wet case, oil bulges into the brine with two approximately equal positive curvatures. This is as expected: near the end of the displacement the oil is trapped in quasi-spherical ganglia in the larger pores [38], which is similar to what has been observed in previous experiments [48,34]. In contrast, for mixed-wet rock, brine now tends to invade the larger pore spaces, with oil retained in connected layers and in the smaller regions of the void space [39]. While this has been observed previously in imaging experiments [49,9], the remarkable finding concerns the interfacial curvature. The mean curvature is approximately zero, see Table 2-the magnitudes are approximately 10 times smaller than in the water-wet case-but this does not mean that the interfaces are flat; instead the typical curvature in any direction is controlled by the pore geometry and broadly similar to that seen in the water-wet case [39]. However, one radius of curvature is positive, while the other is negative. Interfaces of both positive and negative curvatures have also been predicted in numerical simulations of pore-scale displacement, but for separate menisci, and not at the same location as here [50].
The oil/brine menisci are approximately minimal surfaces-interfaces of zero mean curvature which minimize their area subject to a specified boundary [51]. This implies that the majority of the contacts between the two fluids and the solid are pinned, which means that they are held in place during most of the displacement. Hence, to minimize energy, the interfaces find positions of lowest area subject to being stuck at the pinned contacts; this is the definition of a minimal surface. This observation reveals a strength of this approach, as the thermodynamic angle is not affected by these pinned or hinging contact angles that can vary continually but without movement, Fig. 1. Displacement proceeds as a series of lurches from one location of local equilibrium to another with very little change in overall curvature (and hence capillary pressure). Minimal surfaces also provide good connectivity of the phases Table 1 Data from the water-wet experiment. The porosity, /, was 0:207. The imposed fractional flow f w , the water (brine) saturations Sw, specific surface areas awo and aws, mean curvatures j, mean contact angles hg and the computed thermodynamic contact angle ht from Eq. (9) using the differences between the two results shown. Phase 1 is water (brine, w) and phase 2 is oil (o).   Tables 1 and 2 for the water-wet (top) and mixed-wet (bottom) samples. Shown as the vertical lines are the thermodynamic angles, ht, computed using Eq. (9). We estimate that the uncertainty in the calculations is AE5 . Note that the thermodynamic angle better discriminates between the two cases, consistent with very different pore-scale fluid distributions, Fig. 3. and provide a topological explanation for the favourable oil recoveries observed in these experiments [52][53][54][55]39]. Since, in the mixed-wet case, the curvature is almost zero, h t will be controlled by changes in the interfacial areas a wo and a ws .
Indeed in the limit that j ! 0 we can make the following approximation from Eq. (9): From Tables 1 and 2 we can see that the oil/brine interfacial area in the water-wet case, a wo , is largest at low saturation and then decreases, consistent with previous studies [31,32]. At low water (brine) saturation, there is a large area in arc menisci in the corners of the pore space; as water advances, oil becomes trapped in larger pores and the interfacial area is dominated by terminal menisci across throats. In contrast, a wo is generally larger in the mixed-wet case and reaches a maximum at an intermediate saturation when oil is connected in layers throughout the pore space, see Fig. 3, although the change in area with saturation is rather modest. In both cases, the area between the brine and the solid, a ws , increases approximately linearly with brine saturation: the values are lower for the mixed-wet case, as the brine tends to avoid those surfaces that have previously been in contact with oil. The solid/brine areas are always much larger than those between oil and brine, indicating that we have rough surfaces contacted by brine whose area is much larger than that of fluid-fluid menisci which only exist in a subset of the pore space.
We can now understand how the different terms in Eq. (9) contribute to h t . In Table 2 we only quote values where both brine and oil are connected across the system at the resolution of the images and only consider connected phases in the computation of interfacial curvature [39]. In the water-wet case, the value of h t in Eq. (9) is controlled by the large positive value of j, as mentioned previously. In the mixed-wet case, j is much smaller, and the principal determinant of h is the ratio of the changes in the oil/brine area, Da wo to the brine/solid area, Da ws , Eq. (10). The brine/solid area increases approximately linearly with saturation, as mentioned above. At the lower fractional flow values, when the interfacial area increases with saturation, this gives a contact angle slightly below 90 (cos h t > 0); at higher fractional flows when a wo decreases, we see h t > 90 . Again this makes physical sense: we would expect brine to invade locally more water-wet regions of the pore space initially followed by displacement in more oil-wet regions.
From Fig. 2, as discussed previously, h t is lower than the average h g for the water-wet case.
Here, it appears that the geometric angle measures some high values close to 90 which may represent contact lines pinned on a rough surface where a relatively large angle is required for the brine to advance: the displacement, in fact, occurs elsewhere at places where the contact angle is lower. h t is, however, larger than h g for the mixed-wet case. This means that the discrimination between water-wet and mixed-wet systems, which is a small difference in h g , is somewhat greater in h t , which is qualitatively consistent with the very different fluid configurations and curvatures evident in Fig. 3. In this case, we are likely to observe significant contact angle hysteresis: the brine can only advance over an altered wettability surface when the contact angle is large. h g captures smaller values at pinned interfaces where the brine cannot move [13].

Quantification of uncertainty
In the experiments, oil and brine were continually injected. We assume that when the images were acquired the system was at steady-state and in local equilibrium, in that the fluid interfaces remained fixed. This was confirmed by the stability of differential pressure measurements across the sample and from the images themselves [38]. This means that all the energy introduced from pumping the fluids was lost as heat by viscous dissipation. When the fractional flow was changed, there was displacement, and we have assumed that all the additional energy used to displace oil was converted into interfacial area.
The main uncertainty associated with the measurement of the thermodynamic contact angle comes from the segmentation of the raw X-ray images into brine, oil and rock phases, and in the estimation of interfacial curvature [46]. To compute the differences in saturation and area we used the same segmentation method for each image. We considered different plausible thresholds for the segmentation and recomputed the interfacial areas, saturations and curvatures. We calculated an estimated uncertainty in contact angle of approximately 5 : specifically different segmentation still provided a contact angle less than 50 for the water-wet dataset and a contact angle above 90 at the end of the displacement for the mixed-wet case. Table 2 Data from the mixed-wet experiment. The porosity, /, was 0.231. The imposed fractional flow f w , water (brine) saturations Sw, specific surface areas awo and aws, mean curvatures j, mean contact angles hg and the computed thermodynamic contact angle ht from Eq. (9) using the differences between two successive sets of images are shown. Phase 1 is water (brine, w) and phase 2 is oil (o).  The mixed-wet case when f w ¼ 0:8. Note that in the water-wet case, the two curvatures tend to be equal and positive where oil bulges into the brine. In contrast, for the mixed-wet rock, while the mean curvature is close to zero; in one direction the interface seems to bulge into the brine, whereas in a perpendicular direction, brine bulges into oil.

Conclusions and implications
We have defined a thermodynamically-consistent contact angle h t based on consideration of energy balance at states of local capillary equilibrium. It is the angle that correctly captures displacement in quasi-static pore-scale models. This angle was measured directly from high-resolution pore-space images and can be interpreted in terms of the mean curvature of the fluid-fluid interface and changes in interfacial area. We have shown that this angle provides a more sensitive characterization of wettability than the geometric angle, estimated at the three-phase contact line between two phases and the solid, and is consistent with the results of previously-published modelling studies. We suggest that this angle provides a robust characterization of in situ wettability and should be used in modelling studies to predict displacement. In this work, we studied a homogeneous sandstone and computed a single value contact angle; in systems that are more physically and mineralogically heterogeneous, the same approach could be extended to compute contact angle locally, or on a pore-by-pore basis.