# An open-source toolbox for curvature capture in porous media

Jacques Franc, Romain Guibert, Christophe Blondeau, Patrice Creux

► To cite this version:

Jacques Franc, Romain Guibert, Christophe Blondeau, Patrice Creux. An open-source toolbox for curvature capture in porous media. SoftwareX, 2025, 32, pp.102376. 10.1016/j.softx.2025.102376. hal-05482367

HAL Id: hal-05482367

https://hal.science/hal-05482367v1

Submitted on 28 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Original software publication

An open-source toolbox for curvature capture in porous media

J. Franc<sup>a ID*</sup>, R. Guibert<sup>b</sup>, C. Blondeau<sup>c</sup>, P. Creux<sup>a</sup>

ARTICLE INFO

Keywords:

Porous media

Interpolated curvature

Wettability

ABSTRACT

Tomography is a central and essential starting point for Digital Rock Physics (DRP) and more generally, in the study of multiphase flow in porous media. The flow of coexisting fluids in the porous media is highly sensitive to contact angles between phases. Consequently, it is then important to design a tool that is simple to use to enable the extraction of global and localized angles between any pairs of labeled phases.

Code metadata

Current code version

Permanent link to code/repository used for this code version

Permanent link to Reproducible Capsule

Legal Code License

Code versioning system used

Software code languages, tools, and services used

Compilation requirements, operating environments & dependencies

If available Link to developer documentation/manual

Support email for questions

v0.1

https://github.com/ElsevierSoftwareX/SOFTX-D-25-00519

https://codeocean.com/capsule/3828478/tree

GNU General Public License

git

C++

CGAL, cxxopts, CIMG, loguru

https://jafranc.github.io/sx-x-ccpm/

jacques.franc@univ-pau.fr

1. Motivation and significance

Medical diagnosis has made considerable progress with the develop-
ment of new imaging setups that provide high-quality, high-resolution
3D numerical datasets from biological samples. These datasets, which
range from the sub-cellular scale to organs and entire bodies, can be
obtained using a wide array of techniques, including various forms
of optical and electron microscopy, MRI, medical CT scanning, and
micro-XRM CT.

Among the many morphological features, curvature is one of the
most critical. Today, curvature analysis can contribute to identifying
at-risk patients, improving our understanding of disease mechanisms,
and ultimately aiding in the anticipation of future pathologies. A recent
review [1] examines various mechanisms that, at both the cellular and
organ scales, are associated with pathological developments. In the
field of porous media physics, curvature serves not only to characterize
the geometry of the porous structure but also to quantify capillary
forces. These forces, governed by the contact angle at the triple line

between the solid and the two fluid phases, reveal the affinity of the
less-wetting phase for the solid surface. This information is crucial for
predicting whether capillary forces will be significant, beneficial, or
detrimental in specific contexts, and it helps in understanding retention
phenomena in which portions of the non-wetting phase remain trapped
in the porous network.

In the early stages of X-ray micro-tomography, scientists measured contact angles manually using 3D geometric tools. However, it is important to understand that this angle is not unique and may itself reflect complex underlying mechanisms. These include variations in surface composition, surface roughness — which influences interfacial forces — and constraints related to the presence of a second fluid. All of these factors can significantly impact the local wetting behavior and the resulting capillary forces.

Among the current challenges in image exploitation, the ability to
access and analyze curvature through advanced imaging techniques
represents a major advantage—not only for medical diagnostics and

<footnote><sup>a</sup> Université de Pau et des Pays de l'Adour, E2S UPPA, CNRS, LFCR, Pau, France</footnote>

<footnote><sup>b</sup> Institut de Mécanique des Fluides de Toulouse (IMFT) - Université de Toulouse, CNRS-INPT-UPS, Toulouse, France</footnote>

<footnote><sup>c</sup> TotalEnergies, 64000 Pau, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail addresses: jacques.franc@univ-pau.fr, jacquesfrancdev@gmail.com (J. Franc).</footnote>
<table><caption>**Table 1**
Key differences between this work and the existing tool [6] in terms of functionality.</caption><thead><tr><th></th><th>ccpm</th><th>(Blunt et al.)[6]</th></tr></thead><tbody><tr><td>Surface Approximation</td><td>Surface Sampling</td><td>Castellated</td></tr><tr><td>Surface smoothing</td><td>(opt.) Shape regularization</td><td>Laplacian smoothing</td></tr><tr><td>Triple line capture</td><td>Projection Interpolation<br>Thresholding</td><td>Boundary voxels<br>flagging</td></tr><tr><td>Curvature computation</td><td>Interpolated corrected<br>curvature [4]</td><td>Quadratic polynomial<br>fitting (Avizo)<br>and eigen decomposition</td></tr></tbody></table>

biology but also for the field of geosciences in the context of understanding retention phenomena where portions of the non-wetting phase may remain trapped in the porous network.

The proposed C++ software offers a fast assessment of local curvatures and contact angles in images in which the phases are already labeled. Robust phase labeling is already a commonly adopted part of the post-treatment workflow of tomography-X imaging. This would be a valuable addition for rock process characterization.

The software exposes a simple two-stage interface. The first stage consists in taking a labeled 3D TIFF image and splitting it into a series of binary images for each phase and every connected component. The second stage internally produces a 3D surface representation of the binarized object, evaluates curvature at these discretization points and interpolates voxel-derived triple line onto these discrete points to allow for the derivation of contact angles there.

The image processing is deferred to the CIMG [2] library. It uses a classical algorithm [3] to label connected components. The computational geometry is handled through the multi-authored library CGAL [4] which uses surface sampling techniques and 3D Delaunay triangulation to build the surface mesh representation [5]. It employs local interpolation and polynomial smoothing to assess local curvature characteristics [4]. While pursuing the same goal as the works of [6–10], it offers a public implementation of contact angle assessment tool. To the best of the author’s understanding and knowledge, Table 1 reports differences in structures and tools used in the two approaches.

## 2. Software description

In fluid physics, interfacial curvatures play a major role as they result from the combined influence of various interfacial forces, often leading to complex deformations. A particularly critical case is that of capillary physics in porous media, where the solid matrix exhibits highly complex morphology and topology. As a result, under multiphase conditions, droplets or bubbles can adopt intricate shapes that are difficult to characterize, leading to poorly understood effects on fluid flow. The computational code developed and presented here represents a generalized approach in comparison to existing methods. Illustrative examples involving supercritical CO₂ bubbles in rock matrices in the presence of brine are provided to demonstrate the capabilities of this approach.

### 2.1. Software architecture

The C++ software, curvature capture in porous media (ccpm), is interfacing to two well-established image and computational geometry, respectively, CIMG [2] and CGAL [4,5]. It uses CIMG ability to read multiple formats, perform standard image processing transformation, and write to INR format. It uses CGAL’s ability to generate a surface mesh from a grayscale image and its ability to measure local curvatures. Both these tools are products from the well-known INRIA lab in Sophia-Antipolis and examples of the valorization of academic research. They are used in multiple open-source context such as MorphoGraphX, patternize, OpenMEEG and even OpenFOAM. However, alternative software exists that can provide similar tools. The present software could then have been integrated with those alternatives too.

It consists of two main templatized interface that utilizes CIMG [2] specifically Itf_CIMG.h, and CGAL [4,5] specifically Itf_CGAL.h, as its main drivers. Stages are triggered based on command-line arguments:

- `--image` specified the input TIFF image, either the labeled image (for stage 1) or the binarized connected component (for stage 2).

- `--isoval`, `-u` is followed by a list of values that need to be extracted and whose connected component will be used to generate a list of binary images.

- `--output`, `-o` specified by the output path, where the generated artifacts are produced (used in both stages)

- `--ncoordinates`, `-n` refers to the mappings of points and label values that are generated in stage 1 and re-used in stage 2 to represent the triple line.

- `--cutoff`, `-x` indicates the value from which to threshold to isolate the triple line in stage 2.

While the combination of I/O options (that are common to all stages) and isovalue flag trigger stage 1, whereas the I/O options, coordinates and cutoff, with the input pointing to a binarized connected component image, will trigger stage 2. Values of the angle will then be logged in the results of stage 2. Optionally, stage 1 offers an `--ncomponent`, `-c` option that is used to set the maximum number of connected components per phase to process. The default is set to 4. Additionally, the optional `--ard` in stage 2 sets the angle, radius and Haussdorff distance parameters that are used to generate the 3D surface. The default are set to 30° and $r = 2$.

For convenience, the `--image` parameters in stage 2 also accept surface description OFF format and thus bypass its own 3D surface generation process.

### 2.2. Software functionalities

The software will read in 3D segmented and phase-labeled TIFF-formatted images and will produce outputs in STereoLithography (STL) format as well as a series of Comma Separated Value (CSV) files that store probing points, their two main curvature values and derived angles. Local surface curvature is defined by the Frenet–Serret [11] formula as the 2D solution of in a local plane,

$$ \partial_3 \mathbf{T} = k_1 \mathbf{N}, \qquad (1) $$

$$ \partial_3 \mathbf{N} = -k_1 \mathbf{T} + k_2 \mathbf{B}, \qquad (2) $$

$$ \partial_3 \mathbf{B} = -k_2 \mathbf{N}, \qquad (3) $$

introducing $(\mathbf{T}, \mathbf{N}, \mathbf{B})$ the Frenet–Serret basis and $s$ the local arc-length. In differential geometry, Frenet–Serret [11] formula (see Appendix C) allows to study the local behavior of a curve provided that the curve is piecewise $C^1$. Then, in 3D space, the local normal and tangential directions and their variations form a differential system, whose coefficients are respectively the curvatures in both orthogonal directions $k_1$ and $k_2$.

To do so, the software utilizes the CIMG [2] library for its interface to libtiff and the ability to write in Inrimage (INR) format. This enables CGAL Surface Mesh [5] to generate a triangulated surface, constrained by three regularization parameters, respectively, angle, radius and
distance (ARD, see Appendix B) that can be specified as an input argument. They constrain the aspect ratio of triangle elements as well as their uniformity in size and their global curvature. It should be noted that the distance criterion is based on Hausdorff distance and does not play any role if it is larger than the square of the radius criterion.

The interpolated Gaussian curvature is then extracted using Polygon Mesh Processing CGAL library [4]. This allows to compute curvature averaged on an inclusive ball. This will be of interest when seeking to express the Gauss-Bonnet curvature balance on objects of interest.

The second mathematical basis of this work is the Gauss-Bonnet [12, 13] theorem, which states that for homotopical objects, the integral of the Gaussian curvature $κ_G = k_1k_2$ and the mean curvature $κ_m = \frac{1}{2}(k_1 + k_2)$ sums to an invariant constant $2πχ$, in which $χ$ is the topologically invariant Euler characteristic [14].

$$ \int_{\Omega} \kappa_G d\Omega + \oint_{\partial\Omega} \kappa_m d\Gamma = 2\pi\chi, \qquad (4) $$

using $Ω$ to denote the surface and then $∂Ω$ the boundary of this surface. A corollary of this theorem, established in [6], is that despite the approximated nature of surface mesh construction, the contact angle $θ$ can be retrieved as the angle defect from the real contact. On a topologically equivalent sphere object, it leads to:

$$ 2\pi(1 - \cos \theta) = 2\pi\chi - \int_{T_i} \kappa_G d\Omega, \qquad (5) $$

in which $T_i$ stands for the triple line, the locus of contact of all phases mentioned above.

# 3. Illustrative examples

## 3.1. Droplet with single and multiple triple lines

The first test mimics a droplet on a plane solid surrounded by another phase. In the generation of the image, the contact angle is controlled by the distance of the contact chords (see Fig. 1(a)). It introduces an $α \in [-1, 1]$ parameters that scale as $cos(θ)$. Indeed, using A to denote center of the clipped sphere and r its radius, it can be proven that the red angle sine is $α$. The sum of green and red angles is $π/2$.

This functionality is available in the C++ software through the options `-gen-drop`, `-test-angle` and `-test-rot` which correspond to the following parameters,

`--gen-drop`, the nominal pixel size $n$ of the full image,

`--test-angle` the comma-separated list of angles $α$,

`--test-rot` the comma-separated list of rotation angles.

Note that the rotation list must have one fewer element than the angle list. Indeed, the image for a droplet is generated, then rotated and stacked onto the previous ones. The first droplet image is considered unrotated.

The convergence of the total Gauss-Bonnet formula is reported, along with the convergence of the observed contact angle in table Fig. 1(d). There are nonetheless factors that can alter convergence behavior, e.g. the regularity of the sharpest edges and how the triple contact lines are interpolated onto the surface mesh.

The main concern is addressed by using 30 iterations of area and angle smoothing in order to obtain a more regular mesh. This last operation is information-destructive, as the construction of the surface mesh is altered; however, it is beneficial for a robust integration of the Gauss-Bonnet formula.

More generally, there is an optimum trade-off in the process of resolving these contact loops. Upon refinement, the local Gaussian curvature estimate is more precise, but integrated over fewer elements, and thus more sensitive to noise.

The convergence results are gathered in Fig. 1(d). Note that all the errors remain an absolute error threshold of 10°. The overall process is approximately $O(r)$, where $r$ denotes the element size

One of the limitations of the current algorithm is illustrated in Fig. 2. Indeed, as contacts are made on two different solids with different properties resulting in a single merged triple line, the retrieved contact angle is some average of the angle of both triple lines and the angle of intersection of angle loops. To the best of the authors' knowledge, this is a shared limitation with the alternative algorithm. It is though a very restricted case, that would not be impactful if angles are to be determined on a full core plug and a unique value retrieved.

Though appearing simplistic, those academic case covers the different configurations that would be encountered in real medium (see below) while preserving the ability to derive the contact angle analytically.

As we demonstrated the application on an academic case, let us see a use in the context of a micro-CT scan high-resolution image.

## 3.2. Real image case

In this last application section, a real image is input into the global workflow. It consists of a 800³ voxel image (with voxels measuring 4.4³ µm³ each) of brine and CO₂ co-existing in a sandstone rock [15].

### 3.2.1. Launching the case

In order to run such a case, it consists as mentioned in the online documentation (See code metadata table) in two major steps. From a clone of the repository:

```
git clone https://gihub.com/jafranc/sx-x-ccpm
cd sx-x-ccpm
git submodule update --init --recursive
git apply cgal.patch
```

The last line patches CGAL to output binary STL as wanted. Code compilation uses using CMake generator to handle the build stage.

```
mkdir build
cd build
cmake ..
make [-j nprocs]
```

Then, from a pre-existing labeled core-image, the software is first use to extract each label's connected components.

```
/path/to/bin/ccpm --image /path/to/full.tiff
-u 1,2,3
-o /path/to/output
```

with full.tiff being the input image and 1,2,3 the labels in it. This step will provide a series of images `isoVal_1_cc_1.tiff`, `isoVal_1_cc_2.tiff`, `isoVal_2_cc_1.tiff`, `isoVal_3_cc_1.tiff`, `isoVal_3_cc_2.tiff`... if the label 1 has 2 connected components, the label 2 only a single connected component and the label 3 more than 2. It will also produce a file `isoVal_mapping.csv` that maps every voxel with its label values and coordinates. It will be used as a source for projection onto the surface mesh.

The following step should then be integrated in a loop that browses all the connected components of one of the labels (e.g. the third).

```
/path/to/bin/ccpm --image /path/to/output/
isoVal_3_cc_1.tiff
-n /path/to/output/isoVal_mapping.csv -x 5.1
-o /path/to/output
```

The above single line will project values of the aforementioned mapping onto the generated surface using linear interpolation (-n). It will then isolate surface elements onto which the projected values exceed the values input by -x (e.g. 5.1). This will isolate the triple line for curvature extraction and integration. The log will provide as many surfaces and angles as were singled out.
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Geometric contact angle measurement Construction of the droplet and associated control parameters (a), and (b) illustrates of a surface mesh extracted from one of the generated images colored by phase boundary indicators; (c) shows a figure with a top and bottom triple lines (generated via a 180° rotation). In (d), a convergence study is presented with both single and multiple triple lines on a droplet. The reference linear slope is displayed as a thick black line for comparison. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Effects of intersecting triple lines. Here, the contacts with rocks of different properties result in a union of two triple lines. In our algorithm, it is then an average of both angle that is returned. Here, a 36° and a 63° triple line contact are generated, separated by two 90° angles. The resulting angles are converging towards 72°.</figcaption></figure>

### 3.2.2. Analysis of the case

Once each but one label has been isolated in *n* connected components, a surface mesh is built for every component. Then the three-phase triple lines are projected onto this 3D surface mesh. Finally, Gauss-Bonnet application Eq. (4) is used to derive angles.

On Fig. 3(f), the histogram of angles is displayed along with the surface area they represent. In this way, we can observe that the largest peak between 26° and 40° also corresponds for the largest area. This confirms that the statistics do not overestimate an angle due to its repetition on multiple small triple lines.

## 4. Impact

The generalization of the method to a wide range of configurations, combined with its potential implementation in an HPC (High-Performance Computing) environment, will enable the monitoring of capillary pressures — that is, the pressure differences between fluid phases within porous media — in 4D. This will provide access to the dynamic evolution of interfacial mechanisms over time. Another major objective is to extract quantitative insights from a single type of input data (e.g., 3D or 4D imaging), enabling for the reconstruction of multiphase interactions and the prediction of flow behavior without the need for additional experimental measurements. This capability is
<figure><img src="image_4.png" /><figcaption>Fig. 3. Workflow illustration on a full core scan. The labeled image (a) is loaded with two distinct phases (CO2: red and brine: white) and a solid phase (blue). Then (b) the largest component of the gas phase is extracted. Then (c), the three-phase line loops are calculated. Finally (d), the Gaussian curvature is produced and the angle determined using Eq. (4) on each triple line. On the last two lines, (e) the histogram of angle values and (f) their accumulated surface area are shown. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

especially valuable in complex porous media systems where in situ measurements are limited or inaccessible.

This modeling tool enables the quantification of capillary effects in complex porous media, facilitating a deeper understanding of how mineralogical heterogeneity influences two-phase flow dynamics.

The generalization of the method to all topological configurations commonly encountered in two-phase flow allows for the automation of tasks that were previously handled manually. This represents a major step forward in terms of efficiency and reproducibility. By providing

open access to the code, the project offers the broader scientific community a powerful tool to explore and quantify capillary phenomena in complex porous media. This is expected to foster new research directions and support advances in the modeling of multiphase flow, particularly in energy, environmental, and geoscience applications.

## 5. Conclusions

This C++ software, *ccpm*, is an open-source enabler for fast post-treatment of micro-CT porous media images to obtain the distribution
of contact angles [6–9]. It runs in approximately 2.5 h for a 800³ voxel image on a single 11th-gen i7-11700 core. It spends roughly half of the time loading the large labeled image and the other half splitting the connected component of the three phases in stage 1. In stage 2, it spends a third of the time generating a fine enough surface mesh while the remainder is devoted to interpolation of Gaussian curvature, its integration and I/O operations. Bridging the gap towards parallel implementation would be an interesting path to process large volumes of data, such as those produced by scanning campaign that aims at capturing dynamic effects. Integration over the successions of surface mesh generated in such a dynamic context can also provide insight on the capillary pressure evolution of the ganglia. The *ccpm* software is expected to significantly improve our ability to characterize mesoscopic interfacial forces that control fluid behavior in porous media. These forces become dominant in low-velocity flow regimes — such as those found far from wells — and are therefore crucial for understanding flow across most geological formations. By providing a more detailed and physically consistent representation of these phenomena, *ccpm* will contribute to more accurate predictions of subsurface flow, particularly in contexts involving multiphase systems and reactive transport. It can also be used as a direct solver for related inverse problems (e.g., nano-roughness characterization [16,17]). In order to tackle those new challenges, GPU implementation of the surface mesh construction can be triggered in the next release, following the CGAL guidance and existing structures. Another point for improvement in the development will be to introduce a retro-active re-meshing of the surface to increase element density in the area flagged as triple line. Eventually, a web-GUI front end and a dockerized back-end would be an enabler to provide users a better experience.

### CRediT authorship contribution statement

**J. Franc:** Writing – review & editing, Writing – original draft, Visualization, Validation, Software, Resources, Investigation, Formal analysis, Data curation, Conceptualization. **R. Guibert:** Writing – review & editing, Supervision, Methodology, Conceptualization. **C. Blondeau:** Writing – review & editing, Supervision, Project administration, Funding acquisition. **P. Creux:** Writing – review & editing, Writing – original draft, Supervision, Project administration, Methodology, Investigation, Conceptualization.

### Declaration of competing interest

The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: J. Franc reports financial support was provided by TotalEnergies SE. If there are other authors, they declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Appendix A. Gauss-Bonnet derivation

Gauss-Bonnet [12,13] formula is a key tool to link surface mean $κ_m$ and Gaussian $κ_G$ to its global topological invariant, Euler Characteristic $χ$ [14].

It states that in a compact Riemannian manifold $Ω$ with boundary $∂Ω$ integral of the Gaussian on the compact (from Gauss) and the mean curvature on the boundary (from Bonnet) sum up to the global topological invariant $2πχ$. That is a smooth manifold that locally resembles $R^n$ and possesses a Riemannian metric (e.g. an inner product).

$$ \iiint_{\Omega} \kappa_G d\Omega + \oint_{\partial\Omega} \kappa_m d\Gamma = 2\pi\chi. \qquad (A.1) $$

Now let us introduce the fact that at contact geometry might be smoothed out by design as mesh optimization tends to converge towards a moderate angle. On Fig. A.4, the water is light blue, the solid is

<figure><img src="image_5.png" /><figcaption>Fig. A.4. Geometric representation.</figcaption></figure>

black and the gas remains uncolored. The space partition resulting from surface derivation is split by colors while the contact should follow the $θ$ line (in the Riemannian hypothesis).

By introducing $α$ the $π/2$ complementary angle to the bisectrix from the real contact point to the center of curvature with respect to the gas-water interface and in the same way, $β$ the $π/2$ complementary angle to the same bisectrix with respect to the solid gas interface, it can be derived that $θ = α + β$.

Assuming this for the whole contact line will define a local torus in contact with the triple point, whose Gaussian curvature is given by:

$$ \kappa_G = \frac{\cos\nu}{r(R + r\cos\nu)} \qquad (A.2) $$

in which $ν$ and $u$ the torus' angular parametrization, $r$ the inner torus radius and $R$ the outer torus radius. In the limit of thin torus (i.e. $r \ll R$,

$$ \kappa_T = \int_{0}^{2\pi} \int_{-\beta}^{\alpha} \cos\nu r R d\nu du, \qquad (A.3) $$

which after some trigonometry can be rewritten as,

$$ \kappa_T = 2\pi(\sin\theta\cos\beta + (1 - \cos\theta)\sin\beta). \qquad (A.4) $$

On a plane surface, $β = π/2$, which simplify the above formula to:

$$ \kappa_T = 2\pi(1 - \cos\theta). \qquad (A.5) $$

Now as stated above by the Gauss-Bonnet theorem Eq. (A.1), the full curvature balance can be written as,

$$ \iiint_{\kappa_G} d\Omega + 2\pi(1 - \cos\theta) = 2\pi\chi \qquad (A.6) $$

if considering a borderless compact. This latter hypothesis is enforced by the design of the triple line's extraction.

## Appendix B. Sensitivity to the angle and radius meshing constraint

When generating the surface mesh [5], a triplet of parameters can be provided to control the resulting discretization. They noticeably introduce angle and radius constraints, defining respectively the maximum angle in the triangulation triggering splitting and the maximum inscribed radius.

In Table B.2, the sensitivity to angle and radius is explored on a droplet of 26° contact. The sensitivity to radius is the simplest to analyze as it directly controls the number of elements. It is then the trade-off parameter between memory/computation efficiency and accuracy. As illustrated in Fig. B.5, the variation in angles can result in less regular triangulation that can defeat the purpose of accurately
<figure><img src="image_6.png" /><figcaption>Fig. B.5. Illustration of (a) different radii constraint (respectively 4 and 2) (b) different angle constraint (respectively 20 and 30).</figcaption></figure>

<table><caption>**Table B.2**
Relative error for a 26° contact.</caption><thead><tr><th>Angle/radius</th><th>1.5</th><th>2</th><th>4</th></tr></thead><tbody><tr><td>30</td><td>0.0105</td><td>0.0245</td><td>0.0282</td></tr><tr><td>20</td><td></td><td>0.0834</td><td></td></tr><tr><td>15</td><td></td><td>0.0054</td><td></td></tr></tbody></table>

capturing the contact angles, while globally providing a more dense mesh. According to the author's experience a 30° is a good trade-off, while preserving the regularity.

## Appendix C. Frenet–Serret frame

In the oriented $\mathbb{R}^3$ Euclidian plane, an orthonormal plane is usually defined as $(O, x, y, z)$. Let us then define a curve $C$ that has $C^1$ regularity. Then we can define a parametrization of this curve $s(x, y, z)$ such that $C$ is exactly described by its value of $s$, for every point $M \in C$, there is a unique $s$.

Then, according to the cinematic definition of the curvilinear abscise $s$, the direction tangential to $C$ local to the point $M$ is the $T = \frac{dOM}{ds}$.

In the case of a parametric surface $S$, the normal direction is usually defined and noted $N$. Once discretized, this direction is linked to the elements edges. Eventually another tangential direction $B$ can be defined as the last direction forming a direct orthonormal frame. Those 3 directions defines local plane approximation of the parametric surface.

Likewise in a curve $C$, a single curvature can be described $k_1$ that is scaling as the inverse of the radius of the local circle tangential to $M$. It links the local variation of the tangential direction to the local normal and conversely.

$$ \frac{\partial T}{\partial s} = k_1 N \qquad \frac{\partial N}{\partial s} = -k_1 T \qquad (C.1) $$

On a parametric surface, the reasoning in the second tangential direction $B$ is duplicated. Hence,

$$ \partial_s T = k_1 N, \qquad (C.2) $$

$$ \partial_s N = -k_1 T + k_2 B, \qquad (C.3) $$

$$ \partial_s B = -k_2 N, \qquad (C.4) $$

using $k_2$ as the second curvature, scaling as the grand axis of the local ellipsoid. Indeed, the circle on a 1D curve is change for and ellipsoid on a 2D surface to account for non-isotropic curvature.

## References

[1] Schamberger B, Ziege R, Anselme K, Ben Amar M, Bykowski M, Castro AP, Cipitria A, Coles RA, Dimova R, Eder M, et al. Curvature in biological systems: its quantification, emergence, and implications across the scales. Advanced Materials 2023;35(13):2206110.

[2] Tschumperlé D. The cimg library. In: IPOL 2012 meeting on image processing libraries. 2012, p. 4-pp.

[3] Hesselink WH, Meijster A, Bron C. Concurrent determination of connected components. Science Computer Programming 2001;41(2):173-94.

[4] Lachaud J-O, Romon P, Thibert B, Coeurjolly D. Interpolated corrected curvature measures for polygonal surfaces. In: Computer graphics Forum, vol. 39, (5):Wiley Online Library; 2020, p. 41-54.

[5] Sieger D, Botsch M. Design, implementation, and evaluation of the surface mesh data structure. In: Proceedings of the 20TH international meshing roundtable. Springer; 2012, p. 533-50.

[6] Blunt MJ, Akai T, Bijeljic B. Evaluation of methods using topology and integral geometry to assess wettability. J Colloid Interface Sci 2020;576:99-108.

[7] Sun C, McClure JE, Mostaghimi P, Herring AL, Meisenheimer DE, Wildenschild D, Berg S, Armstrong RT. Characterization of wetting using topological principles. J Colloid Interface Sci 2020;578:106-15.

[8] Da Wang Y, Kearney LM, Blunt MJ, Sun C, Tang K, Mostaghimi P, Armstrong RT. In situ characterization of heterogeneous surface wetting in porous materials. Advances Colloid Interface Sci 2024;326:103122.

[9] AlRatrout A, Raeini AQ, Bijeljic B, Blunt MJ. Automatic measurement of contact angle in pore-space images. Advances Water Resources 2017;109:158-69.

[10] Janssen T, Llamas CG, Rieder D, Kuipers J, Baltussen M. A coupled local front reconstruction and immersed boundary method for simulating 3D multiphase flows with contact line dynamics in complex geometries. J Computational Phy 2024;519:113449.

[11] Serret J-A. Sur quelques formules relatives à la théorie des courbes à double courbure. J Mathématiques Pures Appliquées 1851;16:193-207.

[12] Gauss CF. In: Disquisitiones generales circa superficies curvas, vol. 1, DIETERICH; 1828.

[13] Bonnet O. Mémoire sur la théorie des surfaces applicables sur une surface donnée. deuxième partie. GAUTHIER-VILLARS; 1867.

[14] Euler L. Elementa doctrinae solidorum. Novi Commentar Academ Scient Petropolitanae 1758;109-40.

[15] Guibert R, Nazarova M, Voltolini M, Beretta T, Debenest G, Creux P. Influence of standard image processing of 3D X-ray microscopy on morphology, topology and effective properties. Energies 2022;15(20):7796.

[16] Ryan BJ, Poduska KM. Roughness effects on contact angle measurements. American J Phys 2008;76(11):1074-7.

[17] Yan L, Müller JC, van Noorden TL, Weigand B, Raoof A. Wettability-driven pore-filling instabilities: Microfluidic and numerical insights. J Colloid Interface Sci 2025;696:137884.