Event-based contact angle measurements inside porous media using time-resolved micro-computed tomography Journal of Colloid and Interface Science

Hypothesis: Capillary-dominated multiphase ﬂow in porous materials is strongly affected by the pore walls’ wettability. Recent micro-computed tomography (mCT) studies found unexpectedly wide contact angle distributions measured on static ﬂuid distributions inside the pores. We hypothesize that analysis on time-resolved mCT data of ﬂuid invasion events may be more directly relevant to the ﬂuid dynamics. Experiment: We approximated receding contact angles locally in time and space on time-resolved mCT datasets of drainage in a glass bead pack and a limestone. Whenever a meniscus suddenly entered one or more pores, geometric and thermodynamically consistent contact angles in the surrounding pores were measured in the time step just prior to the displacement event. We introduced a new force-based contact angle, deﬁned to recover the measured capillary pressure in the invaded pore throat prior to interface movement. Findings: Unlike the classical method, the new geometric and force-based contact angles followed plau- sible, narrower distributions and were mutually consistent. We were unable to obtain credible results with the thermodynamically consistent method, likely because of sensitivity to common imaging arti- facts and neglecting dissipation. Time-resolved mCT analysis can yield a more appropriate wettability characterization for pore scale models, despite the need

wettability is a function of the intermolecular forces between the fluids and the solid surface (including any coatings or impurities on it). In addition, most natural materials contain surface roughness from the nanometer scale upwards, which influences the effective contact angle observed at larger scales [8]. Due to local variations in mineralogy, surface roughness and coating, contact angles in porous media are often hysteretic, scale-dependent and variable throughout the pore space [5,[9][10][11][12][13][14][15][16][17][18]. One of the main open standing questions is therefore how to define and measure local wettability characteristics throughout the pore space with relevance to multiphase fluid dynamics. This is particularly important to inform pore-scale computational models [19][20][21].
Recent work has shown that contact angles can be measured by geometrical analysis on a 3D image of fluids in the pore space where the interfaces between the fluids remain static. Such images are typically acquired using micro-computed X-ray tomography (mCT), which has been established as an important tool to investigate multiphase flow at the pore scale in recent years [7,[22][23][24][25]. Local measurements of geometrical contact angles in the pore space can be made based on visual observation [22,26], automated algorithms [27,28] or methods based on the deficit curvature of the solid and fluid interfaces [29]. However, these measurements were shown to result in unexpectedly wide distributions of contact angles which are difficult to interpret and to use in pore scale models [30,31]. This is not fully explained by the significant uncertainty on determining the three-phase contact line and the normal to the rough solid surface, caused by partial volume effects and other imaging artefacts common in mCT [32].
None of the methods discussed so far take into account that contact angles measured on static fluid distributions may not be directly relevant to dynamic fluid displacements during drainage and imbibition. Due to unresolved roughness on the solid surface, contact angles can be hysteretic: at a pinned contact line in any pore, the contact angle can vary between an advancing and a receding value, at which the contact line finally starts to move. Furthermore, the location of the contact line in combination with the pore shape at the time of observation, as well as the equilibration time before imaging, likely also influence the observed contact angles [33,34]. Such effects cannot be discerned from an image of a fluid distribution at one specific time. Blunt et al. [31] addressed these concerns by indirectly estimating a thermodynamically consistent contact angle based on energy conservation. They compared mCT images of fluids in a rock sample at the start and end of imbibition, and then equated the pressure-volume work, estimated by measuring the curvature of the fluid-fluid interface and the saturation change, to the interfacial energy stored in the system. The latter can be expressed in function of interfacial areas and a ''thermodynamically consistent" contact angle. The sensitivity to contact line and solid surface normal estimation is therefore reduced or eliminated. Furthermore, this definition aims to take the effects of unresolved solid surface features into account, as it should yield an effective value related to the fluid displacement. Yet, two important issues remain. First, the method assumes that the invasion process can be approximated as being reversible, while this is unlikely to hold for general displacements. Seth and Morrow [35] found that up to 84% of the pressure-volume work during drainage of a limestone was dissipated. Second, the method provides a single contact angle value for the whole sample, rather than a localized measurement.
In this work, we propose to estimate pore-scale (receding) contact angles that are relevant to fluid displacement by analyzing time-resolved mCT datasets. First, we identified the time and location of pores in which fluid displacements took place during drainage. Then, geometric and thermodynamically consistent contact angles were computed on a pore-by-pore basis for each single displacement event at the time just before displacement. We introduce a new force-based receding contact angle definition derived from the measured curvature of a fluid meniscus which triggers it to move through a pore throat (i.e. a Haines jump). The method was tested on two publicly available drainage datasets [36,37].
The methodology and the experimental data are described in Section 2. In Section 3.1, the results from the detection of pore filling events are discussed, followed by a validation of the interfacial curvature analysis by experimental capillary pressure data in Section 3.2. In Section 3.3, the novel time-resolved contact angle measurements are compared to the prior static approach and to each other. This is followed by conclusions and discussion of the current limitations of the method in Section 4.

Materials and methods
In the following section, the conceptual framework for this study is first introduced, followed by a brief overview of the experimental data and a detailed description of the image analysis workflow used to identify filling events and calculate force-based, geometric and thermodynamically consistent contact angles.

Contact angles and displacement events
At low capillary numbers, drainage takes place as a sequence of unstable fluid redistribution events (Haines jumps) every few (tens of) seconds, interspersed by smooth reversible displacement [38][39][40][41]. During the smooth reversible displacement, the contact line can remain (nearly) static close to a local constriction ('pore throat') while the curvature of the interface increases due to the increasing capillary pressure (Fig. 1), as described by the Young-Laplace equation: where r is the interfacial tension, R 1 and R 2 are the principle radii of curvature of the fluid-fluid interface and j is its mean curvature.
When the pressure difference between the two fluids exceeds a certain capillary pressure threshold, an irreversible displacement takes place and the interface abruptly enters one or more neighboring pores (local dilations in the pore space). This threshold is determined by the geometry of the pore throat through which the interface has to pass -the narrower the throat, the higher the threshold -as well as by the local wettability. For example, the associated threshold curvature j thr in a cylindrical pore throat is given by the well-known Young-Laplace pressure: where h is the receding contact angle measured through the densest phase, and r is the throat radius. Therefore, the local wettability in a pore throat can be characterized by defining a force-based receding contact angle h f relevant to the displacement: Note that the force-based contact angle only yields the same value as the receding contact angle for perfectly cylindrical pore throats. However, the advantage is that it provides a direct link to the threshold capillary pressure without the need for measurements at or near the rough solid surface, which are highly dependent on the image resolution and the scale of surface roughness features. Furthermore, Eq. (3) can in principle be extended to arbitrarily complex geometrical pore throat models, e.g. pore throats with triangular [42] or hyperbolic polygonal cross-sections [43]. In this paper, we used the cylindrical model to obtain a first order estimate of the force-based receding contact angle for several hundreds of Haines jumps in a bead pack and a limestone. The following sections describe how h f can be estimated for each individual displacement event in time-resolved mCT datasets, as well as how geometrical and thermodynamically consistent contact angles were determined for comparison.

Experimental data
The proposed method places restrictions on the type of image data suitable for analysis. The flowrate used in the experiments needs to be sufficiently slow to be dominated by capillary forces (i.e. a capillary number of less than 10 À5 ). The temporal resolution needs to be fast enough to capture individual filling events and the spatial resolution needs to be sufficient to clearly resolve interfaces. In addition, pressure measurements made with an external pressure transmitter are useful to validate the obtained results. Based on this, two well-documented and publicly available primary drainage data sets were selected for analysis: one measured by Schlüter et al. on a sintered glass bead pack [36] and one mea-sured by Singh et al. on a Ketton limestone [37]. Both experiments are unsteady-state drainage experiments in which a non-wetting phase is injected into a cylindrical sample at a low capillary number to ensure capillary dominated flow. A summary of the experimental parameters is given in Table 1. The time resolved mCT experiments resulted in a series of segmented 3D images, each representing the fluid distribution in the pore space during a discrete, short time interval in the drainage. In our analysis, we use the segmented data by [36] and [37]. The former used a modified form of Markov random field (MRF) segmentation for this, the latter a seeded watershed algorithm.

Identifying fluid filling events in time-resolved mCT datasets
As a first step in the analysis, the pore space was divided into individual pores separated by throat surfaces using a watershed algorithm implemented in the open source algorithm ''pnextract" [44] (Fig. 2C). This algorithm determines the largest inscribed spheres in each pore and in each throat. Pore and throat radii were found as the radii of these inscribed spheres. To assess the fluid   ) is located in a pore throat with an interface with a certain radius of curvature (R*) and associated contact angle h*. Just before displacement (in red), the interface reaches a radius of curvature (R receding ) and an associated contact angle h receding (measured through the wetting phase and hence receding) upon which it will start to move into the next pore. R receding determines the threshold capillary pressure for the event. occupancy of the center of each pore and throat, the inscribed spheres can be overlain on the drainage dataset, which is registered to a high quality dry-state image from which the PNM is extracted [44]. The occupancy is then determined by determined if the majority of the voxels in the associated inscribed sphere were filled with wetting or non-wetting phase. This was done for each pore in each time step, and ''fluid filling" was consecutively found as the case where a pore changes its fluid occupancy (i.e. the fluid in the majority of its inscribed sphere's voxels) in consecutive time steps. Connected fluid fillings in the same time step were regarded to belong to the same filling event (Fig. 1). This was done by performing a graph-based connectivity clustering of the pores filled in each time step using MATLAB. The adjacency matrix was used to cluster pores which share a connection and change fluid occupancy in the same time step into one event [45,46]. In addition, the neighboring pores are recorded for each event. Finally, the result of this analysis was a list of filling events, detailing the time of filling, the pores that were filled (including location and characterization of the pores), and the pores and throats that neighbor each filling event. The source throat of each event was found by selecting the throat neighbor with the largest throat radius.

Force-based receding contact angles
As defined in Section 2.1, an estimate for the force-based contact angle can be obtained by linking the local capillary pressure which triggered the fluid redistribution to the radius of the source throat of each event. On 3D images of fluid distributions, the local capillary pressure can be estimated from the curvature of the fluid-fluid interface [45,46]. As the volume in space associated to each pore is known, it is straightforward to map point measurements of interfacial curvature to pores. At time step T-1, before a filling event, the fluid-fluid interface is located (partly) in the pores directly neighboring the event. Therefore, we determined the average fluid-fluid interfacial curvature both in pores that were invaded and in their neighbors during the time-step just before the associated filling event (Fig. 2D).
To obtain curvature values, the fluid-fluid interface was extracted from the segmented images by extracting a triangulated surface using the marching cube algorithm [47] in Avizo (Thermo-Fisher Scientific). Constrained surface smoothing using a Gaussian filter with an extent of 3 voxels was performed, analogous to the method described in [48]. Mean surface curvature was calculated using the eigenvalues and eigenvectors of a quadratic form fitted locally to the surface. The accuracy of this calculation is limited by the finite resolution of the images, especially in regions close to the three-phase contact line, where partial volume effects add to the uncertainty, and in arc-meniscus sections of the fluid-fluid interface [19,46,48,49]. The curvature data points were filtered based on several criteria to account for this. First, all data points with a curvature corresponding to a radius of curvature smaller than twice the reconstructed voxel size were omitted as these most likely represent noise. The non-wetting phase bulges into the wetting phase during drainage and terminal menisci thus have a positive curvature with respect to the non-wetting phase [46]. As this study is limited to drainage, all data points with a curvature equal to or smaller than zero, typically associated with arcmeniscus sections of the interface, were subsequently filtered out. Lastly, we followed the approach of [48] on the glass bead data set, to filter and weigh the curvature data points based [27,31,50] on their geodesic distance (20% of the maximum geodesic distance found in the image) to the edge of the surface. For the Ketton data set, we followed the approach of [49] by filtering data points with a Euclidian distance of 3 voxels from the pore wall.

Geometric (receding) contact angles
To reduce uncertainty related to the state of the interface (e.g. pinning), receding contact angles have to be measured at (or at least very near) the moment the fluid-fluid interface moves. Therefore, we propose to measure geometric contact angles at the appropriate time step in a small region directly surrounding the displacement event in time-resolved mCT data. This is in contrast to previous measurements of geometric contact angles that included all points on the three phase contact lines in a single mCT image of a static fluid distribution, typically acquired at the end of drainage or imbibition [27,31,50].
We compared both approaches of contact angle measurements using the fully automatic algorithm developed by [27]. The algorithm generated a smoothed mesh on which the three phase contact line is identified. Subsequently, it calculated geometric contact angles in each mesh point on the contact line based on the dot product of the vectors normal to the solid surface and the fluid-fluid interface. The method was applied to each time step image using the smoothing settings proposed by [27]. Event-based geometric contact angles were then determined by retaining and averaging only the geometric contact angle points measured in the event pores and their neighbors at the appropriate event time step (Fig. 2C). We compared this to the conventional static method by calculating the distribution of all the contact angle points in the last time step of drainage.

Thermodynamically-consistent receding contact angles
Here, we extend the work of Morrow [40] and Blunt et al. [31] to derive a thermodynamically-consistent contact angle, h t , for filling events during primary drainage: where j thr is the threshold mean curvature at invasion, DV nwp the change in volume of the non-wetting phase (positive for primary drainage), DA fluidÀfluid the change in fluid-fluid interfacial area, and DA nwpÀsolid the change in non-wetting-solid surface area (positive for primary drainage). Surface areas were calculated on triangulated surfaces extracted using the method described in Section 2.3.2. By applying this relationship, drainage is assumed to be an isothermal, reversible process, which has been shown not to be valid in general [35,38,40]. The calculated thermodynamically-consistent contact angle is therefore likely an underestimation, as part of the work done will be lost as viscous dissipation [31]. Similar to the calculation of the geometric receding contact angle, we compare the eventby-event based results to a static calculation using surface areas, curvatures and volumes of consecutive time step images.

Filling events
By detecting changes in pore occupancy during the experiments and clustering these into pore filling events (Section 2), 231 (glass bead pack) and 425 (Ketton limestone) filling events were detected and analyzed (Table 2). Up to 17 (glass beads) and 8 (Ketton limestone) different events were detected in a single time step (Fig. 3). The number of pore-filling events per time step increases over time (Fig. 3) while the volumes of these events decreases (Fig. 4). Well connected, larger pores were filled first, and smaller pores with fewer connections were filled later in the sequence, as expected from invasion percolation.

Validation of curvature measurements
The force-based contact angle introduced in this paper depends crucially on the measurement of interfacial curvature from the time-resolved mCT scans. Using Eq. (1), these measurements also yield the associated capillary pressure, which can be plotted against the wetting saturation determined from the image Number of pores in the largest event 18 46 Fig. 3. Bar charts showing the number of (multi-pore) events detected per time step. (Fig. 5). The resulting threshold capillary pressures can thus be compared and validated to capillary pressure curves measured with external pressure transducers [46]. Experimental pressure measurements were available for the glass bead pack dataset, but not for the Ketton limestone dataset. We validated the latter using mercury intrusion capillary pressure (MICP) data, measured on a different Ketton limestone sample [51]. We rescaled the saturation axis of the MICP data to include only the pore space larger than the voxel size of the images (3.28 lm). The Pc axis was rescaled using the interfacial tension measured by [37] and a range of contact angle values (20°-60°) expected for Ketton limestone. As can be seen in Fig. 5, the event-based calculations scatter around the externally measured pressures. The scatter is to some extent expected for drainage, since non-local effects can induce local pressure differences and the time scale of the experiments is likely too short to establish full capillary equilibrium [38,41].
Imaging and image analysis form additional sources of uncertainty. Limited spatial resolution causes uncertainty in the calculated mean curvature of the interface as the curvature is not fully resolved. Akai et al. [19] estimated that local capillary pressures can be estimated within 30% if the average radius of curvature is more than 6 times the image resolution. Applying these rules of thumb, a radius of curvature of 6 times the image resolution would be equivalent to capillary pressures of~0.7 kPa and~2.6 kPa for the glass beads and Ketton limestone respectively. As these values are close to the measured pressures, the spatial resolution is likely to influence the results, especially at lower wetting phase saturations.
The limited temporal resolution of the scans causes motion artefacts and decreases the signal to noise ratio, providing additional sources of uncertainty. Furthermore, this also means that the curvature is determined a few tens of seconds before the Haines jump. To estimate the uncertainty induced by this, we determined the average rate of change of the external Pc measurement in the glass bead pack. This was~43 Pa per time step, equating to an uncertainty of approximately 10% on the average measured event pressures. However, the uncertainty is less significant in the beginning of drainage, and more significant at the very end (when the Pc is changing more rapidly). The temporal uncertainty on the Ketton dataset can be assumed to be lower, as both the flow rate and the time per image are smaller (Table 1), yet there are no external pressure measurements available to estimate this quantitatively. Newer synchrotron beamline setups allow to further increase the temporal resolution by 2 orders of magnitude (1 s/scan) without compromising the image quality [52,53].

Force-based, geometric and thermodynamically consistent contact angles
The distribution of the force-based contact angles is compared to those of the static and event-based geometric contact angles and of the event-based thermodynamically consistent contact angles in Fig. 6. The average force-based contact angle was 63°a nd 56°in respectively the glass bead pack and the Ketton limestone datasets, compared to an average event-based geometric  contact angle of respectively 64°and 57°. The distribution of the event-based geometric contact angles matched the force-based angle closely. The event-based geometric contact angle has a notably narrower distribution than its ''static" counterpart (the distribution of all geometric contact angles measured in one time step image). This confirms that the unexpectedly wide contact angle distributions reported in literature are to a significant extent related to the dynamics of the interface motion, e.g. contact angle hysteresis, pinning, the location of the contact line at the time of fluid redistribution, or interface relaxation. The event-based thermodynamically consistent contact angles are more broadly distributed than the geometric and force-based methods. They had a lower mean value than the other methods for the glass beads (48°) and a higher value for Ketton limestone (66°). Similar values, 44°and 68°for the glass bead pack and Ketton limestone respectively, were found when using the approach of [31], using the full fluid distributions on consecutive images for the calculations. Fig. 7 shows how the determined contact angles varied over the course of the experiments. Event-based geometric and force-based contact angles remain almost constant over time for both data sets The calculated thermodynamically consistent contact angles show considerably more spread compared to the other two methods, which was found to increase over the duration of the experiment. This is likely related to increasing errors in determining the fluidsolid surface areas as smaller pores are invaded. The average thermodynamically consistent contact angle for the glass bead pack data set is negatively correlated with time.
The contact angle values for each event are cross-plotted in Fig. 8. The fairly good match between the force-based and the geo-metric methods suggests that geometric contact angles can be used to provide a reasonable prediction of the invasion capillary pressure in pore scale drainage models, yet only when measured on an event-by-event basis in a time-resolved dataset. This is crucial for e.g. pore network modelling studies, which predict the drainage filling sequence by sequentially invading pores in order of increasing threshold capillary pressure (invasion percolation). However, it is also clear that there was still a significant amount of scatter in the data, related to the spatial and temporal resolutions on the one hand and the irregular shape of the throats in realistic porous materials on the other hand. The latter could be improved by refining the definition of the force-based contact angle (Eq. (2)). Both sources of uncertainty are consistent with higher scatter in Ketton than in the simpler, wider pore space of the glass bead pack. The calculated thermodynamically consistent contact angles show significantly more scattering than the force-based contact angles (Fig. 8). In addition, the equation yielded imaginary values for the contact angle for 17 events in the glass bead pack data set and 125 in the Ketton limestone data set. The increased scatter could be attributed to both the added uncertainty in the calculations of fluid-fluid and fluid-solid surface areas and the assumption of no viscous dissipation.
Intuitively, receding contact angles in smooth, water-wet media are expected to be lower than the reported values. Manual observations using the method described by [54] confirmed contact angles around 50°À60°(e.g. Fig. 2) in both the glass beads and Ketton limestone data set, showing that these were likely not induced by the automatic image analysis method. As shown in Fig. 8, the calculated values were consistent with a description of the fluid  displacement at the scale of observation, and can thus be used directly for pore scale modelling. The discrepancy may be due to the converging geometry of pore throats, which was shown to increase the ''effective" contact angle (i.e. the one linked to the fluid-fluid interfacial curvature) [33]. Finally, it should also be noted that the limited temporal resolution of the mCT data tended to yield an underestimation of the fluid-fluid curvature, as the capillary pressure rises continuously during image acquisition, resulting in a slight overestimation of the contact angles.
As shown in Fig. 9, the values of event-based contact angles can be mapped back to the original pore space. This spatial data could be used to improve numerical simulations of multiphase flow by incorporating local information on the wettability of the sample, which could be especially valuable for samples with a mixedwettability.

Conclusion and outlook
This work aims to improve our understanding of wettability by calculating receding contact angles for individual pore-filling events, rather than for a static fluid distribution as a whole. The proposed method was applied to two unsteady state drainage data sets: a sintered glass bead pack and a Ketton limestone. Eventbased geometric contact angles show a distinctively narrower distribution than when these are calculated on the entire static fluid distributions [27,31,50], suggesting that wide contact angle distributions are likely caused by the unaccounted for dynamics of the interface. We introduce a force-based contact angle, which shows that event-based geometric contact angles produce plausible threshold capillary pressures for associated pore filling events. This suggest that event-based geometric contact angles may provide valid effective contact angles for the displacement process and are more appropriate for use in pore scale modeling. Despite these promising first results, we note the need for enhanced image quality and image processing methodologies to reduce the uncertainty of the proposed methods. Due to these uncertainties, we were not able to draw conclusions on the appropriateness of the thermodynamic contact angle as a concept.
Future work should point out if the described method can be used to quantify wettability for mixed-wettability systems and   9. Rendering of the 3D distribution of the event-based geometric contact angle for each pore. during imbibition. The success hereof would benefit of the increased temporal and spatial resolution available at synchrotron facilities optimized for fast imaging [52] and advances in iterative reconstruction techniques developed for low signal to noise ratios [55,56]. This enhanced image quality is crucial to distinguish different displacement mechanism during imbibition and to improve the accuracy in calculating interfacial curvature and area.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.