Evaluation of methods using topology and integral geometry to assess wettability

a r t i c l e i n f o


Introduction
The recent development of advanced imaging and analysis tools has enabled the sophisticated description of flow in porous https://doi.org/10.1016/j.jcis.2020.04.118 0021-9797/Ó 2020 The Author(s). Published by Elsevier Inc. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). materials at the pore, or micron, scale [1,2]. One topic of active interest is the assessment of wettability, contact angle, within the pore space [3], often associated with estimates of interfacial curvature [4]. Recent work has developed a suite of analysis techniques to measure contact angle automatically [5][6][7][8][9][10] and, combined with the determination of curvature, has been used to quantify wettability and capillary pressure [4,[11][12][13].
Two approaches can be used to determine contact angle. The first is geometric, in that the angle that the fluid/fluid meniscus makes with the solid is determined directly from the image [3]. This method has the advantage of being straightforward whilst automatic methods can generate millions of point values from billionvoxel images [14]. There are two limitations, however. The first is that at the three-phase contact line, the identification and segmentation of the three phases (two fluids and the solid) is uncertain, which can lead to errors in the estimates, depending on image resolution, with a tendency for the values to approach 90 (see, for instance [15][16][17][18]). The second concern is that, even if accurately determined, these angles may represent hinging values on a rough or altered-wettability surface so are not necessarily the correct angles to use to characterize displacement or as input into porescale models [15].
The other approach is to use an energy balance [15,19]. This has the advantage of finding the contact angle associated with displacement, but also suffers from two possible drawbacks. Firstly, it is based on quantifying differences in fluid configurations to find an average value, rather than finding local variations in wettability. However, simulation studies have shown that it is possible to find the contact angle on a pore-by-pore basis when individual displacement events are confined to a single pore region [20]. Secondly, the energy balance does not consider viscous dissipation and is therefore inaccurate for drainage processes with Haines jumps [21], although this effect can be included in the analysis [20].
Recently, concepts in topology and integral geometry have suggested a powerful approach to the analysis of multiphase flow in porous media [22,23]. Sun et al. [17,24] have used the Gauss-Bonnet theorem to derive an effective macroscopic contact angle h macro to study the relationship between intrinsic, advancing and receding angles on rough surfaces. Sun et al. [17] showed that h macro was close to the average contact angle measured automatically [5] within a water-wet sandstone. This approach is, potentially, a major advance, allowing an accurate estimation of contact angle inside porous materials, grounded on fundamental topological principles.
In this paper we first explore the use of the Gauss-Bonnet theorem to determine contact angle, inspired by the work of Sun et al. [17,24]. It is not possible to uniquely associate an average angle with surface integrals of Gaussian curvature of fluid interfaces and the Euler characteristic of fluid clusters. This is because the orientation of the three-phase contact line relative to the fluid/fluid and fluid/solid interfaces is, in general, not uniform along the three-phase contact line. However, we do suggest an approximate relation, consistent with the work of Sun et al. [17], under certain assumptions. We then explore the use of this relationship using simulated distributions of fluids in the pore space, for which the configurations can be determined accurately and where the wettability is known.

External angles of a polygon and curvature
We will start with a pedagogic presentation of topology relevant to flow in porous media for readers who are unfamiliar with the subject (see, for instance, [25] for further reading). Fig. 1 shows a polygon: the sum of the external angles is 2p : where we have n external angles a i . Note that this relationship holds even for a concave polygon, where the angles in the concavity are negative (for example a 4 in Fig. 1). If we have a smooth shape, without corners, there is a similar relationship, if we introduce the concept of curvature, j. At any point L on a smooth curve, l, a circle can be fitted: the curvature is the inverse of its radius where u is an angular coordinate along l: this will decrease along l in concave regions where the curvature is negative. The final case is a shape which is a combination of smooth and angular regions. The general formula is

The Gauss-Bonnet equation in three dimensions
We now consider the extension of Eq. (1) to three dimensions. Here, we start with a smooth surface bounding a threedimensional object. The Gauss-Bonnet theorem states [26] For a twodimensional shape bounded by straight sides, sharp angles and smooth lines, the generalization of the relationship between curvature and angles is given by Eq. (1).
which provides a relationship between the topology of an object and the integral of the Gaussian curvature across the bounding surface. At any point on a surface we can define two principal curvatures in orthogonal directions, j a and j b with corresponding radii R a and R b respectively. If we consider a sphere of radius R, these two radii are equal: R a ¼ R b ¼ R; for a cylinder of radius R; R a ¼ R while R b ¼ 1. If the surface is saddle-shaped one of the radii of curvature is negative. The Gaussian curvature, j G ¼ j a j b , which appears in Eq. (2), is the product of the two principal curvatures: j G ¼ 1=R 2 for a sphere, j G ¼ 0 for a cylinder, while j G < 0 for a saddle. The other term in Eq. (2) is v, the Euler characteristic which is a measure of the topology of an object: 1 plus the number of holes in the structure minus the number of loops. Thus Eq. (2) relates the integral of the Gaussian curvature of an object to its topology.
Let's take a simple example. A solid sphere has an Euler characteristic of 1 (it has no holes or loops). The surface area of a sphere of radius R is 4pR 2 , so the integral of the Gaussian curvatue (which is constant) over this surface is 4p in agreement with Eq. (2). In general though, Eq. (2) is a remarkable result, as it states that however much we distort a sphere (or indeed any object), as long as we do not create holes or loops, the integral of the Gaussian curvature over the surface remains constant.
We will now apply the Gauss-Bonnet theorem to study fluids in a porous medium. First, we will describe the fluid arrangement that we will consider, illustrated in Fig. 2. We have two fluid phases, 1 and 2, separated by a meniscus, described by the surface (or interface) S 12 . The phases also contact the solid: the three-phase contact line, l, between phases 1, 2 and the solid forms a closed loop, and we define a contact angle h measured through phase 1. The solid surface in contact with phase 2 is S 2s . The object we will consider is a blob, cluster or ganglion of phase 2: there is a smooth surface between the two fluid phases, another surface between phase 2 and the solid, and a discontinuity in orientation of the surface at the three-phase contact line -this is encapsulated by the contact angle. We can consider applying Eq. (2) to the cluster of phase 2, where there is an additional contribution to the curvature from the contact line.
Sun et al. [17] derived a version of the Gauss-Bonnet theorem, Eq. (2), as follows: Here we integrate over both the fluid/fluid meniscus, S 12 , and the contact of phase 2 with the solid, S 2s . There is a final term, k d , called the deficit curvature which quantifies the additional contribution from the contact line.
Before continuing, we will explain the purpose of this analysis. The direct evaluation of k d involves an integral along the contact line that -as we show later -requires an estimation of the orientation of this line relative to the fluid/fluid and fluid/solid interfaces, which is similar to what is done to find contact angle directly [5], although even more involved. Hence, it suffers from the same uncertainties and errors in segmentation. What we will look for, as in [17], and discussed in detail later, is a general analytical relationship between k d and contact angle (or an average value), and then use Eq. (3) to estimate this contact angle from integral of Gaussian curvature over the fluid/fluid meniscus, which can be determined with reasonable accuracy from pore-space images [27,18].

Evaluation of the deficit curvature
Physically the deficit curvature k d in Eq. (3) accounts for the fact that there is a kink in the surface of the cluster of phase 2 at the three-phase contact line, or a discontinuity in the orientation of a tangent plane to the surface. This is represented geometrically by the contact angle. If the contact angle is 0 (complete wetting) the entire surface (the union of S 12 and S 2s ) remains smooth, k d ¼ 0 and Eq. (3) reduces to relating the surface integral of the Gaussian curvature of an unbounded shape to its Euler characteristic, Eq. (2). If the contact angle is not 0, there is a sharp change in orientation of the surface of phase 2 at the contact line, as it passes from contact with phase 1 to the solid.
To evaluate k d , which from Eq. (2) represents the integral of the Gaussian curvature in the neighbourhood of the contact line, we mathematically smooth the surface in an infinitesimal region as shown in Fig. 3. Locally the surface near the contact line is a patch of a torus: we call this surface S l . The radius of curvature of the contact line is R at point L, while r is the radius of curvature in the plane perpendicular to the contact line: R ) r and we can take the limit that r ! 0.
We first define the geometry local to a point on the three-phase contact line. For convenience we define the plane coincident with the contact line l as horizontal. The contact angle is h. A line in the vertical plane perpendicular to the meniscus, the surface S 12 , makes an angle a to the horizontal, while the perpendicular to the solid S 2s makes an angle Àb (as drawn in Fig. 3 it is oriented below the horizontal). Then h ¼ a þ b.
On the local torus we define a coordinate system u and v such that a distance increment around the torus in a plane perpendicular to the contact line is rdv, while it is Rdu parallel to the contact line. A surface element dS l ¼ rR dudv. Then the deficit curvature k d , which obeys Eq. (3), can be written as a surface integral of the Gaussian curvature: where j Gt is the Gaussian curvature of a torus. This will not be derived here, but it can be shown that [28]: We consider the limit r ( R such that j Gt ¼ cos v=Rr and hence Eq. (4) becomes: We cannot trivially perform the integral over u since a and b may be functions of u in general -the meniscus and solid do not necessarily maintain a constant orientation at the contact line.
Eq. (7) shows that there is no unique relationship between the deficit curvature and the contact angle; it also depends on b, the orientation of the contact line and the solid surface. This equation is similar to the expression in [17] (see their Eq. (S9) in the Supplementary Material).

Testing the Gauss-Bonnet equation on simple examples
To demonstrate the analysis, we will consider three examples illustrated in Fig. 4; in each one the contact angle and orientation b will be constant around the contact line. In this case, from Eqs. (6) and (7), we have: We have labelled the cases in Fig. 4 such that the angles a and b correspond to the same angles shown in Fig. 3.
The first, easiest, example 1 is a droplet on a flat surface: v ¼ 1. Here b ¼ p=2 and the Gaussian curvature of S 2s is zero. Hence we only consider S 12 in Eq. (3) to find: using a ¼ h À p=2 in this case, which corresponds to Eq. (8).
The second example is when both the meniscus and the solid form parts of spheres. Using simple geometry, the surface area of a sphere of radius r a to an angle p=2 À a to the horizontal, as shown in Fig. 4, is 2pr 2 a ð1 À sin aÞ. Similarly the area of the portion of the sphere representing the solid is 2pr 2 b ð1 À sin bÞ, where r b is the radius of the sphere. The Gaussian curvature of the surface, S 12 , is 1=r 2 a and the surface S 2s is 1=r 2 b . Hence we can write Eq. (3) as: The Euler characteristic of the cluster of phase 2 is 1 (the same as a sphere). Hence to agree with Eq. (3): which is correct -that is, consistent with Eq. (8) for constant a and b.
In the third example the solid bulges into phase 2: a < 0 and b > p=2. We can define a 0 ¼ Àa and b 0 ¼ p À b. Then the areas of S 12 and S 2s are 2pr 2 a ð1 þ sin a 0 Þ and 2pr 2 b ð1 À sin b 0 Þ respectively. As before we obtain from Eq. (3): which again is consistent with Eq. (8).

A simple relationship between Gaussian curvature and contact angle
These examples suggest an approximate yet general relationship between Gaussian curvature and contact angle through making three assumptions. (i) We assume that an arbitrary rough surface has, on average, a zero Gaussian curvature, which means that the integral over S 2s in Eq. (3) vanishes. (ii) Furthermore, we consider that, on average, the orientation of the contact line is such that b has a mean value of around p=2. Under these circumstances we have a situation which is a simple generalization of a drop on a flat surface, example 1, with k d given by Eq. (9). We may have more than one contact line loop, and each one will add to the deficit curvature [17]. This leads to the equation, combining Eq. (9) with Eq.
where n is the number of contact line loops where the single cluster of phase 2 touches the solid. However, as we show later, this does not lead to accurate results. We make one final assumption: (iii) we treat each cluster of phase 2 to be topologically equivalent to a sphere and so v ¼ 1. This appears to be a strong assertion for complex clusters with holes and loops. We will, in the next section, test these approximations and discuss their implications. If we make these three simplifications, we further simplify Eq. (13) to: We can use Eq. (14) to estimate an average contact angle for every discrete cluster of phase 2 in a pore-scale image. Since the orientation of the surface is not always flat, the right-hand side of Eq. (14) will tend to be lower than if a flat surface were present: this means that the estimated value of 1 À cos h may be smaller than its real value, yielding estimates of h erroneously close to p=2 (90 Þ, as encountered when the geometric contact angle is measured [18,14]; on the other hand, consideration of small loops on the surface may tend to underestimate the contact angle. Sun et al. [17], using a sphere on a flat surface (our example 1), identified a macroscopic contact angle h macro with pð1 À cos hÞ=2. They found good agreement with the analytical result when a sufficiently refined surface mesh was used.

Assessment and validation using direct numerical simulation
We will demonstrate the use of Eq. (3) and test the accuracy of Eq. (14) using lattice Boltzmann simulations of multiphase flow in porous media. The advantage of using a simulation method is that the results are unaffected by errors or uncertainty in image acquisition and segmentation. We will compute the terms on the righthand-side of Eq. (14) to estimate the contact angle and compare with the input to the simulations. The lattice Boltzmann method is used with a precise wetting boundary condition to specify local contact angles: details of the approach and validation are presented elsewhere [29,30].
We will consider four cases, Fig

Simulation conditions
The simulation conditions for the four cases are summarized in Table 1. We used densities and viscosities of q ¼ q 1 ¼ q 2 ¼ 1000 kg/m 3 and l ¼ l 1 ¼ l 2 ¼ 1 mPaÁs respectively, while the interfacial tension between the fluid phases was r ¼ 25 mN/m. For the first two cases, we used simple solid geometries where the Gaussian curvature of the interface is analytically determined. For case 1, in the z-direction, a wall boundary was placed at the top and bottom of the domain, while a periodic boundary condition was applied in the other directions. Initially, a semi-spherical droplet of phase 2 with a radius of 20 lattice nodes was placed at the bottom of the domain, while the rest of the domain was filled with phase 1. For this case we simulated a range of contact angles, similar to the example studied in [30].
For case 2, a spherical solid with a radius of 20 lattice nodes was placed at the lower part of the domain. A periodic boundary condition was applied to the boundaries in all the directions. Initially, a known volume of phase 2 was placed on the spherical solid, while the rest of the domain was filled with phase 1. The simulations were performed until they reached an equilibrium state. For this case we studied contact angles of 60°and 120°, similar to the example studied in [30].
For case 3, a cylindrical pore structure with an isosceles triangular cross-section was used of length 175.42 lm (49 lattice nodes), while the triangle side lengths were 71.6 lm, 78.6 lm and 78.6 lm. A periodic boundary condition was applied to the direction perpendicular to the cross-section. Initially, a known volume of phase 2 was placed in the centre of the pore structure, while the rest of the domain was filled with phase 1. The simulations were run until they reached an equilibrium state. For this case we imposed a contact angle of 45°, similar to the example studied in [31].
For case 4, we used two pore structures: a synthetic image of a beadpack and a micro-CT image of a Bentheimer sandstone. The porosity of these structures were 36% for the beadpack, and 21% for Bentheimer sandstone. Porous plates wetting to phases 1 and 2 were placed at þx and Àx faces of the porous domain respectively, while the boundaries in the other directions were closed with a solid boundary. We applied a contact angle of 45°. Drainage simulations were performed in the Àx-direction by gradually increasing the phase 2 pressure through a constant pressure boundary condition at the phase 2-wet porous plate region, while the pressure of phase 1 was maintained constant through a constant pressure boundary condition at the other porous plate. After the simulations reached an irreducible saturation of phase 1, phase 1 was injected in the þx direction by gradually decreasing the pressure of phase 2 while keeping the pressure of phase 1 constant. The simulations were performed until the system reached the residual saturation of phase 2. This simulation condition mimics porous plate capillary pressure measurements in the laboratory. Further details can be found in [20,32].

Estimation of contact angle
In the simulations, the fluid occupancy in each lattice node was obtained based on a colour function, q N , which takes a value of 1 in Table 1 Details of the simulation conditions for the four cases studied.  phase 2 and À1 in phase 1, while it takes the value of À1 < q N < 1 in the fluid/fluid interface region. Our lattice Boltzmann model is classified as a diffusive interface model which produces a slightly diffusive interface with a thickness of 2 to 3 lattice nodes [20,32]. The exact location of the fluid/fluid interface, S 12 , is given by a contour surface corresponding to q N ¼ 0, which has smoothness below the resolution of the simulation grid. The solid surface covered with oil was extracted from the simulation results as a triangulated surface. Then we applied 600 iterations of Laplacian smoothing to this surface. This provides sufficient smoothness while maintaining the original shape of a surface [27,31].

Ganglion configuration
The fluid/fluid and fluid/solid interfaces were modelled as triangulated surfaces on which the Gaussian curvature, j G , was computed using commercial image analysis software (Avizo). This was performed by fitting a quadratic form to the elemental triangulated surfaces, then the magnitude and direction of the principal curvatures were obtained from the eigenvalues and eigenvectors of the fitted quadratic form [4]. The integral of the Gaussian curvature was obtained by taking an area-weighted summation of j G for all the triangulated surface elements. Table 2 shows a contact angles obtained for cases 1 to 3 where relatively simple geometries were considered. There is good agreement between the input contact angle and the value determined using Eq. (14) when the approximations used in the analysis are valid ( R j G dS 12 ¼ 0 and b ¼ p=2), namely a flat solid surface oriented in the plane of the contact line, cases 1 and 3. The contact angle can be evaluated within 5°.
For case 2, when we have a curved solid boundary, representing a ganglion on a spherical grain, the use of Eq. (14) does not work, as the assumption of a flat solid surface oriented with the contact line is no longer valid. However, when the full equation is employed, Eq. (3), accounting for the solid curvature, R j G dS 2s , and its orientation. b, with Eq. (8) to find k d , we can obtain contact angles with similar accuracy to that found for cases 1 and 3.
For case 4, contact angle was determined for residual clusters of phase 2 in the beadpack and Bentheimer sandstone. There were 47 and 78 disconnected ganglia in the beadpack and Bentheimer sandstone, respectively. We chose clusters whose volume was greater than 10 4 voxels for the analysis, resulting in 7 clusters for the beadpack and 17 clusters for Bentheimer sandstone, which accounted for 99% and 97% of the total volume of residual phase 2, respectively. Fig. 6 illustrates the steps taken to perform the analysis for this fourth case, both on a single cluster and for all the clusters in the image. The integral of the Gaussian curvature was calculated on the simulated interface, Fig. 6b and e, while the number of closed contact loops were obtained from the extracted three-phase contact line (boundary voxels of different phases), Fig. 6c and f. As discussed above we assumed v ¼ 1 for all the clusters. Table 3 shows the estimated contact angles for the pore-space images on a cluster-by-cluster basis. It can be seen that for individual ganglia there can be significant errors in the measurement, since the assumptions made in Eq. (14) -a flat surface oriented with the contact line -may not be valid. There is no clear relationship between cluster size and degree of error. However, when averaged over all the ganglia, the contact angle is determined to within 3 , which indicates that this simple approach may indeed provide a valuable average assessment of pore-scale wettability.

Evaluation of the approximations made
In this section we will critically examine the approximations made in the use of Eq. (14) to estimate contact angle. The first and most obvious assumption is that the Euler characteristic, v, is 1. As evident from visual inspection of Fig. 6 some ganglia contain loops and hence do not have v ¼ 1 (the topological equivalent of a solid sphere). We analysed every ganglion and counted the number of loops and holes to calculate v, which is 1 -the number of loops + the number of holes. We then replaced 4p in Eq. (14) with 4pv, Eq. (13): the estimated contact angles are shown in Table 3. What we find is that for all three cases where the Euler characteristic is different from 1, the calculated contact angle is undetermined: this means that j cos h j> 1 and so a physical value of h cannot be found. The reason for this is that we cannot, in these cases, ignore the contribution of the fluid/solid curvature. We can see why this is the case if we consider an example of a cluster containing a hole -a hole that surrounds a solid grain. Then the Euler characteristic is 1 + 1 = 2, adding an additional 4p to the righthand-side of Eq. (3), which has to be exactly balanced by the integral of the Gaussian curvature around the grain which, from Eq. (2), is 4p. Similarly, for a cluster with a loop, the solid-fluid Gaussian curvature has to provide a negative contribution.
This analysis suggests that we should, instead, consider Eq. (3) in its entirety. This is presented in the final column of Table 3. However, we find that in all but two cases the angle is undetermined; this approach simply does not work on pore-space images. The problem is not the evaluation of the orientation angle b, but the computation of the solid/fluid Gaussian curvature: the undetermined values shown in Table 3 represent cases where regardless of h and b no physical values can be assigned that enable Eq. (3) to be obeyed.
In our simulations, and indeed from experimental pore-space images, the fluid/fluid interface is, on physical grounds, smooth. Hence it is possible to extract this meniscus and compute its total and Gaussian curvature with reasonable accuracy, as quantified in Table 2 The results of contact angle obtained for cases 1 to 3. In case 2, the orientation angle of the three phase contact line, b, and the integral of the Gaussian curvature, R jGdS2s, were analytically obtained based on the input contact angle.

Case
Input   previous work [27]. On the other hand, the solid/fluid interface, in both these simulations, and in imaging experiments, is segmented as a voxelized surface which can be rough at all length scales. Even the beadpack, for which the grains are smooth, has sharp contacts with divergent curvature. While the surface is smoothed to remove voxelization artefacts, this does not guarantee that the Gaussian curvature can be found accurately. Hence for both of the samples studied, standard image analysis methods were unable to compute the integral of the Gaussian curvature reliably, largely due to distortion near the three-phase contact line.

Conclusions
We have presented concepts in topology and integral geometry relevant to the study of multiphase flow in porous media. We have shown that the Gauss-Bonnet theorem does not provide a unique relationship between the integral of the Gaussian curvature of the surface of a fluid phase cluster, its topology and contact angle. We have presented, however, an approximate relationship under the assumption that the fluid-solid surface has a zero average Gaussian curvature (it is topologically flat) and is oriented, on average, in the plane of the contact line: this extends the work of [17,24] to provide a simple-to-apply relationship between topology and wettability.
We tested our theoretical analysis against lattice Boltzmann simulations of multiphase flow in porous media. We showed that when the assumptions of the analysis are valid, the contact angle can be determined accurately. When applied to simulations of residual non-wetting phase saturation in pore-space images, the average contact angle was accurate to within 3 for the cases studied. However, individual estimates for single clusters could be significantly in error because of the approximations used in the analysis.
We also showed that attempting to calculate an average contact angle from the full form of the Gauss-Bonnet theorem, including the Gaussian curvature of the solid/fluid interface was not successful, because of uncertainties in the evaluation of this curvature over a voxelized interface, even after smoothing.
Overall, the use of the Gauss-Bonnet theorem, as introduced by Sun et al. [17], does provide a useful methodology to determine an average contact angle. It is a complement to methods that evaluate curvature and contact angle geometrically from pore-space images [4,3,5]. It provides a topological counterpart to the thermodynamic contact angle that characterizes a displacement process [15].
Future work could be devoted to extending this analysis to a wider range of porous media, wettabilities, and to the analysis of experimental datasets. We could consider the analysis of wettability in mixed-wet media where the Gaussian curvature of the fluid menisci is negative [12] and for intermediate saturations where both phases are continuous throughout the image domain. It would also be of value to examine and reduce the errors in the calculation of the Gaussian curvature of solid surfaces from porespace images [31]. Finally, it would be instructive to relate the contact angle found using our proposed expression, Eq. (14), to macroscopic flow properties, such as relative permeability and capillary pressure [24,22].