Contact Angles in Two-Phase Flow Images

In this work, we calculate contact angles in X-ray tomography images of two-phase flow in order to investigate the wettability. Triangulated surfaces, generated using the images, are smoothed to calculate the contact angles. As expected, the angles have a spread rather than being a constant value. We attempt to shed light on sources of the spread by addressing the overlooked mesh corrections prior to smoothing, poorly resolved image features, cluster-based analysis, and local variations of contact angles. We verify the smoothing algorithm by analytical examples with known contact angle and curvature. According to the analytical cases, point-wise and average contact angles, average mean curvature and surface area converge to the analytical values with increased voxel grid resolution. Analytical examples show that these parameters can reliably be calculated for fluid–fluid surfaces composed of roughly 3000 vertices or more equivalent to 1000 pixel2. In an experimental image, by looking into individual interfaces and clusters, we show that contact angles are underestimated for wetting fluid clusters where the fluid–fluid surface is resolved with less than roughly 500 vertices. However, for the fluid–fluid surfaces with at least a few thousand vertices, the mean and standard deviation of angles converge to similar values. Further investigation of local variations of angles along three-phase lines for large clusters revealed that a source of angle variations is anomalies in the solid surface. However, in the places least influenced by such noise, we observed that angles tend to be larger when the line is convex and smaller when the line is concave. We believe this pattern may indicate the significance of line energy in the free energy of the two-phase flow systems.

will depend on the interfacial tensions between the two fluids ( 21 ) and fluids with solid ( 2s , 1s ) through Young's equation (1), (Young (1805) and Laplace (1806)) The applicability of Eq. (1) is restricted to ideal smooth surfaces of a homogeneous solid material, and the corresponding contact angle is sometimes referred to as an intrinsic contact angle (Sun et al. 2020a). The effect of rough surfaces can be included in the surface tensions 1s and 2s , with a resulting effective contact angle . Even though the surfaces of glass beads in synthetic media are smoother than surfaces in natural porous media, the angles measured by means of micro-CT in the micron scale should be categorized as the effective angles. However, we expect the effective interfacial tensions to be fairly constant throughout the medium. We refer to Appendix 1 of reference (Berg et al. 2020), for further reading about sub-resolution roughness.
Accordingly, the contact angle is expected to be close to constant in a single material porous medium, such as sintered glass beads. Similarly, the mean curvature of fluid-fluid interfaces, H,f , depends on the pressure difference, P c , across the interface and the interfacial tension through Laplace's equation, and hence, the fluid-fluid interfaces are expected to have constant curvature.
The contact angle, as an important means of measuring wettability Anderson (1986), influences fluid distribution in two-phase flow. It is also an important input parameter in pore scale simulations (Ramstad et al. 2019). With recent advances in micro-computed X-ray tomography, configurations of two fluids in complex three-dimensional porous media have been imaged (Berg et al. 2013;Schlüter et al. 2017). Therefore, one wishes to use a contact angle measured in situ from images in the simulations. Andrew et. al. (2014) measured contact angles manually, whereas Klise et al. (2016), Scanziani et. al. (2017) and  performed automated measurements with different methodologies. However, the measured contact angles at multiple (thousands or tens of thousands) contact points in the X-ray images tend to show a wide distribution rather than a single value (AlRatrout et al. , 2018aAlhammadi et al. 2017). Mixed wetting, surface roughness, and advancing/receding contact angles are thought to be sources of the spread (AlRatrout et al. , (2018a. However, a constant contact angle and surface curvature is also not expected if the free energy contains additional contributions related to the length of three-phase contact lines or surface curvatures (Khanamiri et al. 2018), predicted by Hadwiger's characterization theorem (Hadwiger 1957).
User-biased image segmentation, the limited image resolution, and artifacts imposed by surface smoothing will contribute to a spread in the measured angles and local interface curvature. In order to answer questions related to the relative importance of surface roughness, advancing/receding angles, and additional free energy contributions for the contact angle and curvature variations, it is important to minimize and quantify these errors. The effect of advancing/receding angles was addressed by Mascini et al. (2020). They proposed an event-based approach where they investigated only receding contact angles in pore-filling events in the drainage process and observed a narrower distribution of the angles. In order to solve some of the mentioned issues, a number of averaging schemes have been proposed recently. For instance, Blunt et al. (2019) proposed a thermodynamic averaging scheme, which was updated later by including the dissipation Berg et al. 2020). In this work, we do not investigate the thermodynamic approach where the efficiency of the displacement process is a major unknown.
In an number of recent works, integral geometry is used as an alternative to the direct angle measurement. Sun et al. (2020b) introduced this approach by linking local curvatures and topology of a cluster using the Gauss-Bonnet theorem. Blunt et al. (2020) and Sun et al. (2020a, c) elaborated on the methodology. These methods often require calculation of local curvatures which are typically products of surface smoothing, a common part shared between the direct measurement and the integral geometry methods. Addressing the sources of angle distribution in the direct method, in this work, can potentially be used in improving accuracy of the methods based on integral geometry. Furthermore, our approach is entirely a mathematical solution to surface smoothing, without any assumptions based on the current understanding of a physical equilibrium, neither with any assumptions on the possible variables influencing the contact angles. We choose this approach in order to avoid smoothing artifacts possibly induced by such assumptions. The counterparts of our approach are for instance the work by Raeini et al. (2012) or smoothing by the Surface Evolver package. Raeini et al. (2012) ran two-phase flow simulation on the experimental images in order to correct the imaging artifacts where the simulated velocities were inconsistent with the physical model.
In this work, we present a method for extracting contact angles and curvatures from tomography images by smoothing triangulated surface representations of the solid-fluid and fluid-fluid interfaces (Khanamiri 2019). High resolution experimental data (Schlüter et al. 2016b) of sintered glass beads and simple fluids were investigated to minimize possibility of contact angle spread as a result of mixed wetting or other surface heterogeneities.

Triangular Mesh and Differential Geometry
In this section, a number of properties for a triangular mesh are introduced from differential geometry. The wire-frame in Fig. 1 shows 1-ring neighborhood N 1 (i) of vertex x i . N 1 (i) is the triangular domain confined by the immediate neighbor vertices of x i , for instance x j . Unit normal vector of each triangular face, placed at centroid of the triangle, is calculated by cross product of two edges of the triangle. The cross products should be consistently clockwise or anti-clockwise to ensure all normal vectors pointing from one phase-fluid or solid-to another in the entire mesh. Unit normal vector for x i is calculated by Eq. (3), where T k is a triangular face in N 1 (i) , with unit normal vector n(T k ).
The shaded region in Fig. 1 is Voronoi area of x i , which is a fraction of the area in the N 1 (i) neighborhood and defined by different methods. One simple method is to connect centroids of the triangles to midpoint of the edges. This results in a Voronoi area equals to one-third of the entire 1-ring neighborhood of the vertex, Desbrun et al. (1999), and can be calculated by, where angles ij and ij are the two angles opposite to the edge ij, as shown in Fig. 1. The definition by Meyer et al. (2003) is the same as Eq. (4) only when ij , ij and the angles at x i are nonobtuse. In case of obtuse angles, the portion of A voronoi (x i ) from the edge ij in triangle T is area(T )∕4 if T is obtuse at x i , and area(T )∕2 if T is obtuse in other angles. A visual inspection showed that the definition proposed by Meyer et al. (2003) ensures a robust smoothing in complex structures, whereas the method by Desbrun et al. (1999) can be unstable in neighborhoods with rapidly changing curvatures such as where two glass spheres are sintered.
Calculation of curvatures and smoothing a mesh requires finding Laplace-Beltrami operator, K for vertices. According to the derivation by Meyer et al. (2003), Equations (4) and (5) are known as cotangent discretization. For vertex x i , K(x i ) and the mean curvature H (x i ) are related by K(x i ) = 2 H (x i )n(x i ) . Therefore, using dot product of two vectors, we have For a mesh, say a fluid-fluid interface, the average mean curvature H can be calculated using where A is area of the interface and V is the number of vertices on the interface.
The contact angle of solid-fluid and fluid-fluid meshes on a vertex on the three-phase contact line x i,3 can be calculated by where n s and n f are unit normal vectors of x i,3 in the solid-fluid and fluid-fluid meshes, respectively.
Note that in the definition by Meyer et al. (2003), the sum of Voronoi areas for all vertices does not always add up to the total mesh area due to presence of the obtuse angles; and in order to obtain an accurate H , we should consider using correct mesh area in denominator of (7). Our calculations in geometries with known curvatures showed that this approach ensures correct result in calculating H using Voronoi region definitions by both Meyer et al. (2003) and Desbrun et al. (1999). Since there is no fundamental and unified definition of Voronoi regions, the point-wise mean curvatures tend to depend on the choice of Voronoi region and are unreliable. However, H obtained by integration in (7) is reliable as average mean curvature for a mesh. Furthermore, choice of Voronoi region does not have a large effect on point-wise contact angle calculation, as (x i ) in (8) is calculated using unit normal vectors which only weakly depend on the Voronoi region definition.

Mesh Generation
Segmented tomography images of fluids in porous media are volume data comprising voxels on a regular grid, and interfaces are defined on the voxel boundaries. A triangular mesh for the interfaces is created by defining triangles on each of the corresponding voxel boundaries. Every square face is symmetrically divided into four triangle faces as shown in Fig. 2. Using four triangles in the squares instead of two gives sufficient valence to all vertices for smoothing and curvature calculations. The resulting mesh will have a stair-case configuration with sharp edges and corners and must be smoothed in order to find contact angles and other related parameters such as curvatures.
Note that when building a mesh for an interface (say fluid-fluid interface) from a segmented image, it is straightforward to identify the triangular faces and also the orientation of these faces, i.e., the normal vectors of the faces. The edges and vertices defining these faces are, however, not always to be shared between faces even when they have the same coordinates. This occurs when an interface, due to insufficient resolution, touches itself or another interface at a vertex or along an edge as shown in Fig. 3. These vertices and edges must be identified and duplicated into separate entities for each part of the interface. The vertices that must be duplicated are called nonorientable since, for a self-touching interface, the surface orientation is undefined at these points if the vertices are not duplicated. Nonorientable vertices tend to be abundant near the three-phase contact lines and in narrow parts of the pore space, and incorrect treatment of these vertices will lead to errors in contact angle and curvature calculations. These errors are especially severe for properties of the wetting fluid trapped in narrow regions forming pendular ring structures.
Vertices in the mesh grid correspond to either a center or a corner of voxel boundaries in the voxel grid (see Fig. 2). The center vertices are always orientable, whereas the corner vertices can be nonorientable. Therefore, we use the correct orientation of center vertices in the 1-ring neighborhood of corner vertices to identify and correct the nonorientable ones. Figure 4 shows examples of nonorientable vertices from experimental images-both far from and in the vicinity of three-phase contact line-and how they affect the mesh smoothing with and without corrections. Figure 4d-f shows the nonorientable vertices close to the three-phase contact line that, if not corrected, can be falsely classified as contact vertices where the contact angles are measured. As an example, in an experimental image with 232882 vertices on the three-phase contact line, 15181 vertices were nonorientable, and 8813 of nonorientable vertices were either on the three-phase line or incorrectly considered as a part of it before corrections.

Mesh Smoothing
The generated initial mesh has sharp edges and corners and must be smoothed in order to find contact angles and other related parameters such as curvatures. The vertices on the threephase contact lines belong to both the solid-fluid and the fluid-fluid meshes which must be smoothed consistently through smoothing iterations. Different smoothing techniques have been proposed; for instance, Akai et al.  2017), in the sense that vertices are displaced to minimize the curvatures. However, weights of displacements and the way we use neighborhood of vertices are different in this work. Our methodology is mainly based on isotropic smoothing flow by Meyer et al. (2003). We displace the vertices along their unit normal vectors with mean curvature as the weight of disparagement. The constraint on the displacements Here f is the function expressing the surface. For the Laplace's equation, smoothing is also equivalent to Euler-Lagrange equation, H = 0 used for surface area minimization (Meyer et al. 2003). In absence of a function, the surface is discretized and the Laplace's equation is solved by a numerical method, e.g., finite difference. A 2-manifold-a mesh-is represented by vertices and discrete triangles rather than by an ordinary continuous surface. Generalization of Laplace's equation from surfaces to manifolds is known as Laplace-Beltrami method which can be discretized using cotangent scheme in order to calculate Voronoi area (see Fig. 1), Laplace-Beltrami operator and mean curvature in Eqs. (4)-(6).
Smoothing is an iterative procedure where in each iteration, all vertices are displaced simultaneously after calculating the displacements. A vertex x i is displaced a distance x i along the surface unit normal vector n(x i ) with H (x i ) as the weight of displacement. Smoothing the solid-fluid and the fluid-fluid meshes are performed in a similar manner, except for the vertices lying on the three-phase contact line.
The average mean curvature for individual fluid-fluid surfaces, H,f is updated for each iteration using Eq. (7). The coefficients 0 < t s < 1 and 0 < t f < 1 in (10) are tuning parameters used mainly to avoid too much displacement which can cause smoothing to diverge in meshes imitating complex porous structures. We start the smoothing with t s = 0.15 and t f = 0.3 and change the values in the iterations following a reward/penalty strategy, meaning that if for instance, solid-fluid mesh improves in one step t s is multiplied by 1.05, and if it worsens, t s is multiplied by 0.8 for the next step. Coefficient t f is tuned with the same rule. This procedure was fixed after several test runs, and there is no theoretical reason for the choice of values for the coefficients. One can perform the iterations with both different initial values and multipliers for t s and t f and obtain smooth meshes as long as the smoothing does not diverge.
The displacement of all vertices is subject to the following two constraints.
First constraint in (11) ensures that x i does not displace more than 0.2 voxel length at each iteration. The displacement is reduced to 0.2 voxel length, if it is larger than that. Second constraint ensures that x i does not move farther than √ 3 voxel length from its initial position. This value is the longest distance inside a voxel in the voxel grid. Maximum displacement of √ 3 is associated with the uncertainty in the experimental images such that x i could in reality be anywhere inside one voxel. In the images investigated, only a minor fraction of (9) the vertices are moved a distance close to √ 3 in voxel length. For instance, after smoothing a tomography image, mean displacements of vertices in solid-fluid surfaces, fluid-fluid surfaces, and three-phase contact lines were approximately 0.30, 0.30, 0.52 voxel length, respectively; vertices displaced more than 0.95 √ 3 were only 0.20, 0.03, 0.02 percent of all vertices in those entities.
As mentioned earlier, for a smooth solid-fluid surface, mean curvatures are being minimized H,s (x i ) → 0 . It should be noted that H,s (x i ) values are not becoming exactly zero because of the constraint in Eq. (11). The geometric implication is that for a solid-fluid surface throughout the iteration steps, the unit normal vectors of neighbor vertices align more parallel to each other once the surface is being smoothed, and their dot products approach unity. Therefore, we use s , as the target of smoothing for the solid-fluid mesh, where V s and E s are the number of vertices and edges in the solid-fluid mesh.
For a smooth fluid-fluid surface, the mean curvatures approach an unknown constant value H,f (x i ) → H,f . The standard deviation of H,f (x i ) should be locally minimized for each interface. Therefore, the target for the entire domain can be minimizing average of the standard deviations, std, for all fluid-fluid surfaces, and Eq. (13) serves as the target for smoothing the fluid-fluid mesh.
In Eq. (13), M is the number of fluid-fluid interfaces in the mesh grid.
We use and t variables to stop the iteration. Smoothing iterations continue as long as both s and f are reducing. If their minimum values are at the latest step, it means that they have not reached a minimum and smoothing is continued. The solid-fluid mesh often improves faster, and s can reach a minimum before f . At this point, f has not reached a minimum and the iterations must be continued. In these conditions, the quality of solid-fluid mesh may degrade as s will start to increase. Therefore, the coefficient t s in Eq. (10) reduces automatically to keep the changes in the solid-fluid mesh to a minimum. The reason we do not immediately set t s = 0 is that the solid-fluid mesh in every iteration needs to be readjusted in the vicinity of the three-phase contact line when fluid-fluid mesh improves. With t s = 0 and a changing fluid-fluid mesh, the line will be pulled out of the solid-fluid mesh. Once f has also passed a minimum and starts to increase, t f in Eq. (10) will automatically reduce. At last, when both s and f have passed minima, and one of the tuning coefficients is t < 10 −3 , the smoothing stops after five more iterations, and the results of iteration with least f is picked up as the final result. At the end, the number of iterations is approximately in the range of 50-100. Figure 5 shows variations of contact angle, s and f , as smoothing progresses in an analytical model of a sphere-plane geometry explained in the next section.
In a smooth mesh grid, s value is approximately in the order of 10 −3 with variations between two steps in the order of 10 −5 ; f is approximately in the order of 10 −2 with step variations in the order of 10 −3 . Smoothing several experimental images showed that although cotangent discretization smoothes surfaces, it is not convergent in the sense that despite the point-wise values of H become more uniform, they are not equal (Meyer et al. 2003;Xu 2013). Improvement slows down as smoothing progresses. The point-wise convergence may depend on the voxel grid resolution. In the next section, using analytical examples with known curvature and contact angles, we show that although point-wise H is not convergent, average mean curvatures H and point-wise contact angles are convergent when grid resolution is sufficient.

Analytical Examples
We attempt to verify the smoothing algorithm and to investigate the limitations by examining two examples with known contact angle and curvature. We create synthetic voxel images of two different configurations simply by intersecting sphere-plane equations where the sphere radius and the distance between the plane and the sphere center are known. Figure 6 shows the two examples where the solid-fluid interfaces are simply represented by planes, and the fluid-fluid interfaces are part of a sphere. For these two configurations, analytical H and are constant. Images with different resolution were created by increasing number of the voxels representing the geometry. Figure 7 shows the calculated error with respect to inverse of number of vertices in the fluid-fluid surface 1∕V f . The grid resolution is higher at smaller 1∕V f . Figure 7 shows that point-wise contact angles (Eq. (8)), average mean curvature H,f (Eq. (7)), and interfacial area are converging as 1∕V f → 0 , whereas point-wise mean curvatures H,f (Eq. (6)) are not converging.
V f = 148 ( 1∕V f = 6.7 × 10 −3 ) and V f = 700 ( 1∕V f = 1.4 × 10 −3 ) for the spherical cap, V f = 540 ( 1∕V f = 1.8 × 10 −3 ) and V f = 2226 ( 1∕V f = 4.4 × 10 −4 ) for the cluster between two parallel planes have not yet captured the contact angles with an acceptable accuracy. Furthermore, A f and H,f , which is the inverse of radius for the complete sphere in these examples, can also be used as further discriminants to determine if the resolution was adequate in order to find and H,f . Overall, a fluid-fluid interface composed of roughly 3000 Furthermore, for interfaces with sufficient resolution standard deviations of is roughly 1 − 2 • ; standard deviation of H,f is in the order of 10 −3 pixel −1 .

Experimental Examples
In this section, we investigate dependency of contact angles on the size of fluid clusters in a 3D experimental X-ray tomography image published by Schlüter et al. (2016b). The twophase flow experiments were performed with low flow rates at quasi-static equilibrium.  The solid was composed of sintered glass beads with two different diameters; the fluids were brine and n-dodecane (Schlüter et al. 2016a). In this paper, we use an image from the early main imbibition with water saturation of 15.1 percent. Comparison of mean contact angles for individual fluid clusters revealed that the calculated contact angles for small clusters are smaller than those for large clusters. Figure 8 illustrates six examples. The large clusters (a)-(d) with several thousand vertices in the brine-dodecane interface have mean contact angles m of 68 • , 69 • , 71 • and 74 • ; whereas smaller clusters with only a few hundred vertices have m of 49 • and 31 • . The cluster (f) in Fig. 8, has m of 31 • which is even smaller than 49 • for cluster (e). Cluster (f) is a droplet lying on the body of glass bead, while cluster (e) is one attempting to encircle the rim between two fused glass beads. We found that the small clusters of type (f) have often the least contact angles among the poorly resolved clusters. We believe this is a smoothing artifact. Smoothing without a constraint on the vertices displacement will eventually flattens a mesh. Displacement constraint of √ 3 , Eq. (11), is a large degree of change for a small surface which stretches considerably to be close to a flat surface. This results in underestimation of both the contact angle and the curvature. Figure 8 shows also number of three-phase contact lines L 3 , and Euler characteristic , of 2-manifold of the cluster surfaces. is calculated by Eq. (14), where V, E, and F are number of vertices, edges, and triangle faces, respectively, on the entire surface of fluid cluster regardless of how fluid-fluid and solid-fluid interfaces cover the surface. This should not be mistaken with Euler characteristic of the cluster as a volume which is half the value of . For clusters homeomorphic to a sphere and a torus is 2 and 0, respectively (see Fig. 8d, e). More complicated clusters will have values as negative even integers, as shown in Fig. 8a-c. On the size-biased contact angles for poorly resolved small clusters, results also showed that very small clusters with a volume of only one or a few voxels have a contact angle of near or equal to zero, i.e., smoothing collapses both fluid-fluid and solid-fluid interfaces on a plane. Figure 9a shows m reduces when size of the fluid-fluid, V f interface reduces. Although the number of three-phase contact vertices on the smaller interfaces is not large, prevalence of such interfaces will change the mean contact angle for the whole image. Figure 9d shows the histogram of angles for the whole image rises relatively sharply for angles close to zero. In the same figure, the blue histogram is drawn after eliminating contributions from interfaces with V f < 500 . The changes between the two curves for angles roughly below 45 • is noticeable. In Fig. 9d, the blue curve (unfiltered) is result of measurement from 376,827 three-phase contact vertices V 3 , with mean angle and standard deviation of 65 • and 17 • , respectively. The red curve (filtered) includes 354,362 vertices, 94 percent of all vertices on the three-phase line, with mean angle and standard deviation of 68 • and 11 • , respectively. Standard deviation of angles std( ) for the individual interfaces with respect to V f and V 3 are illustrated in Fig. 9b, c. std( ) of large clusters is roughly 10 • , whereas it is nearly any value from 0 to 25 • for smaller ones. Lower std( ) for a group of smaller clusters does not imply precise measurements. It only means the angles are closer to an average angle at vertices the experimental measurement is highly uncertain. As a rule of thumb and in order to be cautious, we recommend finding the mean angle from fluid It should be noted that because of high frequency of very small clusters ( V f < 100 ) and for the sake of better visualization, those clusters are not presented in Fig. 9a-c. The trends remain the same with/without such clusters, most of which are a few voxels or even one in volume. Furthermore, the largest fluid-fluid interface with V f = 1,602,085 and V 3 = 168,510 has m = 67 • and std( ) = 9 • which are similar to m and std( ) of large clusters and filtered results of the whole image. These high V f and V 3 are not presented in Fig. 9a-c, again for visualization purposes. The largest interface, which is the one separating the connected displacing and the connected displaced fluids, has m value close to the mean angles of the well-resolved trapped clusters in Fig. 8.

Contact Angle Variations Along three-Phase Line
We showed that contact angle is best to be measured or averaged for the larger fluid-fluid interfaces where there is less measurement uncertainty compared with the smaller ones.
In an example, the calculated for large clusters had similar values of m ≈ 68 • and std( ) ≈ 10 • . In this section, we investigate the local variations of contact angles along the three-phase lines. We show that 10 degrees standard deviation is partly random caused by imaging noise and partly a pattern that we believe is possibly related to the overlooked line energies.
We sort the values of the consecutive vertices along a three-phase line. In this way, angles resemble a time series, for which the auto-correlation function ACF can be calculated. Figure 10 shows the ACF, with 100 as number of the lag periods, for multiple lines on a single cluster. According to Fig. 10, ACF values are high for smaller values of lag, meaning that the variation of along the lines is smooth. We examine the hypothesis that whether or not the variations are related to the line curvature. Therefore, we calculate the angle created by the two edges neighboring each vertex on the three-phase lines. values are a measure of how much the line turns at a vertex and hence proportional to the line curvature. Figure 10 shows that lower values tend to occur when the three-phase line curvature is higher. The − plots have a high degree of noise and may be affected by the imaging noise and the smoothing artifacts in poorly resolved features. We further investigate the variations by visualizing the original and smoothed solid surfaces and the color-coded values in Fig. 11. As shown in Fig. 11, is lower at and around kinks (high ) in the line. The kinks occur at surfaces with anomalies in the form of a small whole or a small lump of voxels on surface of the glass spheres. These features are certainly not resolved, and the resulted and values are not reliable.
However, further investigation of variations along the line, with least noisy image features, revealed values tend to be larger at the convex parts of the three-phase line and smaller in the concave parts (see Fig. 12). This can also be explained by the generalized Young's equation Pethica (1977), Boruvka and Neumann (1977), where is the line tension and is the signed line curvature. Convex parts of the line with > 0 will have larger , whereas concave parts with < 0 will have smaller . This pattern in the experimental images suggests that the likelihood of line energy playing a role in twophase flow should not be ruled out. The variable contact angle correlating with line curvature hints that the free energy of the system should have an extra term for the line energy, in addition to contributions from the volumes and the surface areas.
In the following, we will estimate the line tension based on the changes in contact angle along the three-phase contact line. Given the rotation angle , in Fig. 10, at a vertex on a contact line, the line curvature at the vertex is Fig. 11 Color-coded-and size-codedvalues along three-phase contact line viewed from inside the solid; the noise in the solid surface represented by the kinks in the line results in smaller Fig. 12 -color-coded and size-coded-tends to be larger at convex parts of the three-phase line and smaller at the concave parts. This pattern was observed when the line neighborhood in the image was not noisy where s is the edge length between the two consecutive vertices. s is in the order of pixel length. Assuming effective average values for the tension terms in Eq. (15), one can construct the equation for two vertices on the convex and concave parts of a three-phase line.
Subtracting the two equations in (17), we have One can calculate using a pair of vertices with known and by mean of Eq. (18). Using five pairs of vertices, all on a single three-phase line, and an interfacial tension of 12 = 36 × 10 −3 N/m , measured by Schlüter et al. (2017), the line tension was estimated from Eq. (18) as ≈ 8.7 × 10 −7 N . Attempting the procedure on another line results in similar value of ≈ 9.3 × 10 −7 N . These two values were calculated from the two lines shown in Fig. 12.
The calculated is close to higher line tension values reported in the literature. Zhang et al. (2018) reported that the line tension of sessile droplets obtained from theories vary from 10 −12 N to 10 −10 N whereas experimental measurements vary from 10 −12 N to over 10 −5 N , depending on the studied system (Zhang et al. 2018). Further, it has been debated that while the line tension at atomic scale is on the order of nN, gravitational effects at mm scales can add significant contributions to the line tension (Bruce et al. 2017). In our microscale measurements, it is possible that is influenced by the sub-resolution surface roughness and wetting variations. Now we can estimate the interfacial energy F A 12 and the line energy F L for the fluid cluster for which the statistics is presented in Fig. 10.
The estimated line energy is nearly equal to the interfacial energy. The line energy will therefore influence the free energy of the two-fluid system in micro scale. The studied experimental system is composed of relatively large and smooth glass beads; thus, we can expect stronger effect in natural porous media such as sandstone rocks.

Summary
We calculated the contact angles and curvatures in X-ray tomography images of two-phase flow in porous material. Two analytical examples with known contact angles and curvatures were used to verify the mesh smoothing method. Analytical cases showed that average and point-wise contact angles, average mean curvature, and surface area converged with increased grid resolution. According to analytical examples, point-wise contact angles and average mean curvature for fluid-fluid surfaces resolved by 3000 vertices or more, (18) = 12 cos concave − cos convex − concave + convex .
(19) F A 12 = 12 A 12 = 1.8 × 10 −8 J F L = L ≈ 1.6 × 10 −8 J equivalent to an area of 1000 pixel 2 , could be calculated accurately. We addressed also the impact of poorly resolved small fluid clusters-which can often have high frequency in the experimental images-on the contact angle distribution and mean contact angle. We found that contact angles were underestimated for wetting fluid clusters where the fluid-fluid surface mesh was resolved with roughly less than 500 vertices. In addition, the cluster-based investigation showed that both mean and standard deviation of angles where the fluid-fluid surfaces had at least a few thousand vertices converge to similar values. We observed a standard deviation of approximately 10 degrees for large fluid clusters. Our investigation of the local variations of contact angles along the three-phase contact lines showed that a part of the variations is caused by measurement uncertainty. However, a part of variations represented a pattern in which contact angles tend to be larger at convex parts of the line and smaller at concave parts. This suggests that the free energy may have contributions from three-phase contact lines, in addition to the volumes and surfaces.