An Automatic In Situ Contact Angle Determination Based on Level Set Method

A novel level set‐based approach is presented to calculate in situ contact angle distribution, θ, from pore‐scale immiscible fluids and rock configuration directly imaged with micro‐computed tomography (micro‐CT) techniques. We first identify interfaces of the fluid‐fluid and fluid‐solid by their zero level set functions. This is accomplished by reinitializing the level set functions with a signed distance function. Then the three‐phase contact line is determined at the crossover of the two zero level set functions that represent the two interfaces. The normal vectors of both surfaces are calculated directly using the two level set functions, and the contact angle is found from the dot product of these vectors where they meet at the contact line. We first validated our newly proposed method for the semianalytically calculated fluid configurations in a 2‐D tube and then tested the algorithm on a synthetic spherical oil droplet residing on a tilted flat solid surface where the contact angle is analytically defined. It was then used to measure the in situ contact angle of droplets directly imaged by micro‐CT, and the results are compared with the manually and other available automatically measured results. Compared with other available automatic approaches, our approach is mathematically well defined, and it does not require any other complicated tuning procedures for surface smoothing. This proposed approach allow us to accurately characterize local in situ pore‐scale wettability, which is essential to model multiphase flow in porous media and eventually help us to design and assess optimal processes, such as hydrocarbon recovery and carbon dioxide storage.


Introduction
Wettability is a key property that defines the tendency of a fluid to spread on a solid surface in presence of two or more immiscible fluids, which plays a critical role in multiphase flow in porous media as it controls pore-scale fluid configurations and displacement scenarios (Frette & Helland, 2010;Ma et al., 1996;Zhou et al., 2014), which in turn influences the macroscopic flow function, such as capillary pressure and relative permeability curves (Anderson, 1987;Jettestuen et al., 2013;Zhou et al., 2012Zhou et al., , 2014. It is relevant to a number of processes, such as extraction of hydrocarbon from subsurface reservoir, CO 2 geological storage and sequestration, and fuel cell efficiency (Blunt, 2017;Liu & Wang, 2020). Wettability could be quantified by averaged wetting index at core scale, for example, Amott index and USBM index derived from capillary pressure curves (Amott, 1959;Donaldson et al., 1969), and by contact angle defined at pore scale with measurements performed on smooth pure mineral surface (Blunt, 2017;Morrow, 1990). Besides conventional lab-measured contact angle, quite a few numerical simulation approaches, for example, molecular dynamic simulations (Derksen, 2015;Stukan et al., 2010;Tian & Wang, 2017;Yong et al., 2019), disjoining pressure modeling using Frumkin-Derjaguin equation (Zhou et al., 2017) have also proposed to simulate the contact angle under different contexts. However, these measurements and numerical approaches fail to capture the pore surface roughness and mineral heterogeneities at the pore space.
Beside the ex situ measurements, advanced imaging technologies provide another promising approach to study the in situ contact angle based on imaged fluid configuration directly at pore space. The idea was initially proposed by Andrew et al. (2014), who manually measured the contact angle directly on a raw image of the plane perpendicular to the three-phase contact line. This image-based manual method has been adopted to study in situ contact angle under various conditions (Khishvand et al., 2016;Singh et al., 2016); however, the manual measurement is time-consuming and may also introduce human bias, and thus it is limited to relatively small samples (AlRatrout et al., 2017). In order to overcome the limitation of this manual in situ ©2020. The Authors. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. contact angle measurement method, Klise et al. (2016) and Scanziani et al. (2016) proposed automated approaches to determine in situ contact angle based on fitting planes/lines to voxelized images. As pointed by AlRatrout et al. (2017), these attempts have been proved to be unsatisfactory when measuring in situ contact angle for low-resolution micro-computed tomography (micro-CT) images and for complex rocks, so they developed another automated method to measure the in situ contact angle at pore scale. Unlike the planes/lines fitting approach, this recent approach proposed by AlRatrout et al. (2017) is based on discretizing the fluid-fluid and fluid-solid and define the vectors that have a direction perpendicular to these surface, and the contact angle defined at three-phase contact line is defined as the dot product of vectors describing the fluid-fluid interface and solid surface as given below: where n ! f f and n ! f s denotes the vector normal to the fluid-fluid and fluid-solid interface, respectively. This approach was implemented in the OpenFoam platform and requires tuning processes to generate smooth surfaces (AlRatrout et al., 2017). What's more, the algorithm is complicated and not well explained mathematically.
In this work, we developed a level set-based method to calculate the effective contact angles directly from imaged fluid distributions in pore space. The principle of our approach is to represent the fluid-fluid and pore-solid with two different level set functions. The normal vectors of the two interfaces (represented by the zero level set function) is calculated based on the level set functions, and then the contact angle defined at the three-phase contact points is calculated using the dot product of the two vectors that are normal to the interface. The paper is organized as follows: Section 2 describes the level set-based method devised to obtain contact angle distribution automatically. The detailed numerical implementation is presented in Section 3. In Section 4, we first validate the method on semianalytically calculated fluid configurations in a straight capillary tube and then in synthetic 3-D images of a spherical droplet of oil residing on a tilted flat solid surface. Section 5 demonstrates the application of our method on in situ wettability measurements in X-ray images of water-wet and mixed-wet porous media.

Method Description
The level set based automatic in situ contact angle estimation method includes three main steps: First, establish level set function for both solid phase and nonwetting phase; second, calculate normal vectors for both fluid-fluid and fluid-solid interface; and then estimate contact angle as the dot product of vectors describing the fluid-fluid interface and solid surface at the three-phase (wetting, nonwetting, and solid) triple line.

Level Set Function for Solid and Nonwetting Phase
The level set method is a numerical interface tracking method (Osher & Fedkiw, 2003;Sethian, 1999), and its zero level set of a function ϕ describes the interfaces implicitly, which is one dimension higher than the interface. The sign of level set function ϕ determines the interior and exterior regions of that separated by the interface, and the interface could move normally to itself with a velocity. The motion of the level set function is governed by the given equation where V n is the speed function that is the normal vector of the velocity. In order to maintain numerical stability, the ϕ is reinitialized occasionally during evolution of the level sets by equations 2. In order to reinitialize ϕ to a signed distance function, we reinitialize ϕ to a stationary state with the method proposed by Sussman et al. (1994) using the following equation: Here SðϕÞ ¼ ϕ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi ϕ 2 þ j∇ϕj 2 ðΔxÞ 2 q is a sign function suggested by Peng et al. (1999) to reinitialize the level set function, and Δx is the spatial step.

10.1029/2020WR027107
Calculating the in situ local contact angle for the segmented three-phase (includes solid, wetting, and nonwetting phases) image requires the detailed distribution of the solid and pore interface, wetting and nonwetting interface, and the vector normal to the interfaces defined at the three-phase contact line. In this work, we define a function ψ to represent the pore and solid space locations: ψ > 0 represents the pore space, ψ < 0 defines the solid matrix, and consequently ψ = 0 defines the interface between the solid and pore space. Similarly, the nonwetting/wetting fluid locations are defined by the function ϕ, in which ϕ < 0 represents nonwetting fluid and ϕ > 0 represents wetting fluid. The fluid interfaces are described by ϕ = 0. Figure 1 gives a schematic of the setup.

In Situ Contact Angle Defined at Three-Phase Contact Line
The zero level set function of the solid-pore, ψ, and wetting-nonwetting, ϕ, explicitly express the solid-pore interface and fluids interface, and the three-phase contact points could be determined by the crossover points of the two interfaces. Due to discretization and finite resolution in the segmented fluid-rock images, the smeared-out Heaviside step function and the smeared-out delta function (H(ϕ)) and δ(ϕ)) were used to extract the interfaces from the stationary level set function. These two functions are defined as (Osher & Fedkiw, 2003) : and Figure 1. The figure shows an illustration of the setup. The shaded regions represents the solid, θ is the contact angle, the red vectors, n ! φ and n ! ψ represents the surface normals of φ and ψ at the pore boundary, and β is the angle between them. : where ϵ is the parameter that determines the width of the interface. In this work, we choose ϵ = 1.5Δx. The three-phase contact points could be determined by δðϕÞδðψÞ > 0: The contact angle defined at the three-phase contact line could be estimated directly using the two level set functions as given below: The advantage of the level set function is that the normal vector for both interfaces could be calculated from the level set function automatically.

Numerical Implementation
Images of the segmented porous structure and fluid distribution are taken as input to the model. The corresponding level set functions for these solid-pore and wetting-nonwetting, ψ and ϕ, respectively, are computed by conventional reinitialization using equation 3. We used the mirror-reflected boundary conditions at the edges of the computational domain, that is, the stencils involving boundary voxels in the discretization schemes are mirrored across the boundary. In this procedure, the gradients of the level set functions are discretized using the third-order weighted essentially non-oscillatory (WENO) scheme and Godunov's method for upwinding (Osher & Fedkiw, 2003), and the detailed scheme is given as below:  10.1029/2020WR027107
The time is updated using the third-order Runge-Kutta method (Osher & Fedkiw, 2003), and it is started with Euler step to advance the solution to time t n + Δt, and then followed by a second Euler step to advance the solution to time t n + 2Δt, and then obtain the solution to time t n þ 1 2 Then the last Euler step to advance the solution to time t n þ 3 2 Δt is The final solution for next time step is then given as The two level set functions, ϕ and ψ, are reinitialized independently and are considered to be stationary when their error of the last reinitialization are sufficiently small. This convergence criterion is expressed as follows: Here Ω is the entire computational domain and tol is a tolerance value which is set equal to 0.002 in the simulations presented in this work.

Model Validation
We first validated our method for the fluid configurations simulated at a 2-D cross section using the semianalytical method (Frette & Helland, 2010); then we also tested the algorithm on a synthetic spherical oil droplet residing on a tilted flat solid surface where the contact angle is analytically defined.

Fluid Configuration Simulated in Capillary Tube
The pore space is represented by a capillary tube with constant cross section, and the geometry of the cross section is extracted directly from 2-D scanning electron microscopy (SEM) image of Bentheim sandstone. In such a tube, the fluid configuration at given capillary pressure and contact angle during drainage is simulated using the semianalytical method (Frette & Helland, 2010). For example, Figure 2 presents a cross section of the simulated fluid configuration in a capillary tube (the contact angle used for the semianalytical mode is θ = 20°), and its dimension is 2,056 × 1,288. As shown in this figure, there are four fluid-fluid interfaces (eight associated contact points), and the red and blue regions denote the nonwetting and wetting phases, respectively. In the test case, we constructed the capillary tube with 10 identical cross section along the z direction, so the size of the computational domain is 2,056 × 1,288 × 10 voxels.
As expected, we have identified 80 contact points using our level set-based automatic approach in this test case, and the contact angle is in a range of 15°to 25°. The measured contact angle deviates from the analytical one that has been used to generate the equilibrium fluid configuration in the capillary tube; the potential source of the uncertainty in contact angle in this capillary tube mainly originates from the process of image discretization and solving the equation numerically.

Static Spherical Oil Droplet Residing on a Tilted Flat Solid Surface
A benchmark validation with a static spherical oil droplet residing on a tilted flat rock surface, where the analytical solution is available, is performed to test the accuracy of our level set-based automatic in situ contact angle measurement method. In the test case, the solid-pore interface is represented by a tilted plane with a slope of y = x, and the spherical oil droplet is sliding on the tilted plane and surrounded by brine (see Figure 3). As shown in this figure, the voxel resolution is represented by the length scale of 1 voxel, which is presented as 1/R of the diameter of the spherical droplet (here R is the diameter of the spherical droplet, and fine resolution is represented by a larger value of R). As no gravity and other external force is considered in this validation case, and the contact angle defined at the contact line could then calculated by It has to be noted that for a given R, the different contact angles defined at the contact line could be achieved by moving the blue line using different c values in equation 15.  resolution 1/6, middle intermediate voxel resolution 1/15, and lower fine voxel resolution 1/25. The analytical contact angles for the three different voxels resolution case are 55°, 80°, and 86°. As shown in this figure, the contact angles are clearly defined at the three-phase contact points. Due to the stair-step voxel configuration of the fluid and rock phases, the contact angle estimated from our approach is a distribution. For the coarse resolution, 38 contact points has been identified, and the estimated contact angle is in a range of 56°to 70°and its mean contact angle is 66°. For the intermediate resolution, 72 contact points has been identified, and the estimated contact angle is in a range of 57°to 86°and its mean contact angle is 76°. For the fine resolution, 84 contact points has been identified, and the estimated contact angle is in a range of 60°to 90°a nd its mean contact angle is 80°.
The influence of voxel resolution on the accuracy of estimated contact angle is further investigated by assigning the length scale of 1 voxel, which is set to 1/12, 1/17, and 1/27 of the diameter of the spherical droplet. Their fluid-rock configuration are shown in Figures 5-7, respectively. The analytical contact angle in these tests varies from 5°-175°. As shown in Figure 5, the analytical contact angles defined at the solid-fluids contact points for the coarse voxel resolution case are θ = 64°, θ = 85°, and θ = 95°. The analytical and simulated contact angles for the coarse case are presented in Figure 8. As shown in this figure, the analytical solution is well captured by the automatic measurement, and it is located in the error bars, which show the difference between the maximum and minimum value in the contact angle distribution estimated from our level set based approach. A similar situation was also found for the test cases with intermediate and fine resolution, as shown in Figures 9 and 10 (please also refer the relevant fluid-rock configurations for intermediate and fine resolution in Figures 6 and 7, respectively). Note that our approach tends to slightly overestimate contact angle at intermediate wetting condition, θ changes from 75°to 105°, and underestimate the contact angle at both strongly water-wet or oil wet conditions. The difference between the analytical and our simulated contact angles becomes smaller when the resolution is improved. Similar test cases have been used to validate the most recent automatic approach (AlRatrout et al., 2017), which, however, requires to run a significant number of sensitivity tests to find the best tuning parameters to achieve a good in situ contact angle simulation. One strength of our approach is that it does not require to run larger number of presimulations to choose the best tuning parameter, and the only uncertainty in our approach is ϵ that defines the interface thickness. Previous studies have shown that ϵ = 1.5Δx is a good approximation to capture the interfaces.

Results and Discussion
We now apply our algorithm on two sets of high-resolution microtomography data sets of a ganglion in water-wet and altered wettability systems. The in situ contact angle distributions of the two different wetting systems are presented and compared with the most recent automatic approach proposed by AlRatrout et al. (2017).

In Situ Contact Angle of Oil Ganglia in Water-Wet Pore Space
The high-resolution (voxel resolution is 2 μm) fluid-rock images were obtained using the X-ray microtomography techniques for a water-wet Ketton carbonate rock sample after a drainage and imbibition cycle. The original data were prepared by Singh and Blunt (2018) and are available from the Digital Rock Portal. The two oil ganglia are presented in Figures 11 and 12. The dimension of the first microtomography images  is 365 × 255 × 225 voxels, and the oil ganglion is relatively smaller and regular (see Figure 11), while the second microtomography image is larger, 630 × 410 × 510 voxels, and a more complex oil ganglion pattern could be identified (see Figure 12).
A total of 2,324 and 16,239 three-phase contact points are identified using the criteria described in equation 6 for the high-resolution images, and the calculated in situ contact angle distributions are presented in Figures 13 and 14, respectively. Similar to Scanziani et al. (2016), we then interpret the estimated contact angle with an empirical probability density function through the truncated Gaussian model as follows: Here μ and σ are mean of the calculated contact angle and standard deviation, and a and b are set to be 0°and 180°, respectively. p(x) is the normal Gaussian distribution function, and P(x) is the Gaussian cumulative distribution function; they are given as The matching parameters of the modeled in situ contact angle probability function for the first case are μ = 51.6°, σ = 13.2°and for the second case are μ = 48.9°, σ = 17.4°. The two means are very similar, and this suggests that the two systems are very similar from a probabilistic point of view. The estimated mean contact angle for the first oil ganglion with the most recent available approach by Scanziani et al. (2016) is μ = 37.5 and by AlRatrout et al. (2017) is μ = 51.3. The standard deviation from two above-mentioned in situ contact angle measurements is σ = 13.7°and 12.9°. These three automatic measurements all suggest the system is weaker than slightly water-wet system, and these automatic measured contact angles are consistent with the ones measured manually using the approach proposed by Andrew et al. (2014). It was found that the contact angles were in the mid 40s (Scanziani et al., 2016).

In Situ Contact Angle of Water Ganglia in Mixed-Wet Pore Space
A sandstone sample was aged with crude oil over 2 weeks, and then it was flooded after drainage and imbibition cycle; the detailed fluid and rock configuration was then imaged with X-ray microtomography with voxel resolution of 2 μm. The dimension of the microtomography images is 132 × 122 × 22 voxels, as shown in Figure 15, and there is one larger water ganglion sitting in the middle of the pore space and amount of small oil ganglia attached to the rock surface.
A total of 2,712 three-phase contact points are identified using the criteria described in equation 6 in a mixed-wet rock sample, and the calculated in situ contact angles distributions are presented in Figure 16. The matching parameters of the modeled in situ contact angle probability function of this oil wet system are μ = 96.9°, σ = 39.6°. It has to be noticed that compared with water-wet case, the standard . The error bars show the difference between the maximum and minimum value in the contact angle distribution estimated from our level set-based approach. Figure 11. The first trapped oil ganglion isolated from the water-wet system.

10.1029/2020WR027107
Water Resources Research deviation of the contact angle distribution is much larger; this implies that contact angle spread in mixed-wet rock sample is much wider.
As discussed in the above cases, our level set-based approach has been proven to be effective in measuring in situ contact angle in both water-wet and mixed-wet system for brine-oil system. The method itself should also be applicable to other systems, such as air-brine, air-oil, and other immiscible multiphase systems.

Summary and Conclusions
In this work, we have introduced a level set-based automatic method to calculate the in situ contact angle from imaged three-dimensional pore-scale fluid configuration. The main idea is to represent the fluid-fluid and pore-solid interfaces with two level set functions and then determine the contact angle at the three-phase contact points by calculating the dot of the two vectors normal to the two interfaces.

Water Resources Research
The newly proposed approach was first validated using fluid configuration at straight capillary tube calculated using the semianalytical approach and single spherical oil droplet on a tilted plane. The algorithm was then applied to estimate in situ contact angle distribution on micro-CT images of fluids and rock configurations in water-wet and mixed-wet reservoir rock. The distribution of contact angles measured by our method is consistent with manual measurements (Andrew et al., 2014) and the most recent automatic approaches (AlRatrout et al., 2017;Scanziani et al., 2016). As demonstrated in the in situ contact angle distribution, both water-wet and mixed-wet contact angles could be expressed as normal distribution probability density function. The standard derivation of contact angle distribution function in mixed-wet sample is much larger than that of water-wet samples. Compared with other automatic measurement approach, our newly proposed level set-based approach is mathematically well defined and requires less tuning parameters to reduce the uncertainty in the measurements. However, one limitation of the current version of our level set-based method is relevant computational demand during reinitialization to compute level set function compared with the method proposed by AlRatrout et al. (2017). One natural extension of our level set-based approach could be extended to measure the local curvature to estimate the local capillary pressures for multiphase flow in porous media.
In situ contact angles measured at different conditions in porous media make it possible for researchers to define more realistic contact angle distribution for reservoir and fluid characterization. This will help us to reduce one of the main uncertainties for multiphase flow modeling at pore scale.