# Robust detection and characterization of a bifurcated bi-phasic interface propagating through a granular medium: physically and morphologically

Rana Al Nemer, Julien Réthoré, Giulio Sciarra

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

## To cite this version:

Rana Al Nemer, Julien Réthoré, Giulio Sciarra. Robust detection and characterization of a bifurcated bi-phasic interface propagating through a granular medium: physically and morphologically. Experiments in Fluids, In press. hal-04778816

HAL Id: hal-04778816

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

Submitted on 12 Nov 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Robust detection & characterization of a
bifurcated bi-phasic interface propagating through
a granular medium: physically and morphologically

Rana Al Nemer, Julien Réthoré† and Giulio Sciarra†

Nantes Université, Ecole Centrale Nantes, CNRS, GeM, UMR 6183,
F-44000, Nantes, France.

*Corresponding author(s). E-mail(s): rana.al-nemer@ec-nantes.fr; rana.alnemer96@gmail.com;

Contributing authors: julien.rethore@ec-nantes.fr; giulio.sciarra@ec-nantes.fr;

†These authors contributed equally to this work.

Abstract

The detection of an interface separating two zones of low contrast and evolving in a deformable medium is very challenging. The low-contrasted images describing a partially saturated granular medium are replaced by high contrasted images generated from the correlation error maps issued from a standard Digital Image Correlation (DIC) calculation. A robust algorithm has been developed to capture automatically the sought interface evolving in 2D space through time. Capable of detecting all morphological reliefs, this approach is completed by a sequence of steps to remove aberrant contours and track the main advancing interface. Consequently, adequate algorithms have been developed to determine physical and morphological properties that will give insights on the factors that monitor the interface propagation.

**Keywords:** interface detection, bi-phasic flow, imaging techniques, deformable granular medium
# 1 Introduction

Hydraulic instabilities establishing as a finger-like pattern have been determined to occur in a drainage process depending on the governing forces: capillary or viscous [1]. When the fluid is injected with high flow rate, the flow regime is viscous-dominated and the established displacement pattern is described as viscous fingering; in reverse, when the flow rate is slow, the flow is capillary-controlled and the displacement pattern is described as capillary fingering [2, 3]. Overall, the topology of the developed pattern is characterized by a very ramified and heterogeneous structure with various local reliefs especially when this interface is travelling a granular porous medium which may complexify its detection.

Extensive studies have investigated these hydraulic instabilities and focused on the detection of the interface separating two immiscible fluids mainly in a 2D configuration using the so-called Hele-Shaw cell [4–6] but also in 3D configuration using X-ray tomography [7] or the planar laser-fluorescence imaging [8, 9]. Focusing on the bi-phasic flow in 2D space, for instance in [10], when air is injected into glycerol-water saturated glass beads, the interface was detected using image segmentation: a suitable threshold is fixed in order to separate the two main clusters of gray level. Another technique used in [11–13] merely consists in subtracting the gray level of reference image from that characterizing the current state, neglecting possible grain reorganization caused by the infiltration process. These conventional techniques have turned out to be ineffective and inaccurate under conditions of low contrast and grains rearrangement. In addition, regardless the approach used for the interface detection, the analysis of the detected contour has been limited to macroscopic description without focusing on local characteristics such as local velocity and local characteristic size of the interface roughness which are essential to figure out the factors that could monitor the interface evolution.
The aim of this work is to develop robust algorithms capable of (i) detecting automatically the interface evolving within a low-contrast 2D granular medium taking into account the interactions between the flow and the grains and (ii) extracting information regarding the motion of the interface and its topological properties. This characterization is of paramount importance in order to get a full analysis of the flow within a granular medium and to correlate the interface propagation with the response of the skeleton.

The paper is organized as follows. First in 2, the experimental setup generating the hydraulic instabilities within a granular sand medium along with the acquisition system are briefly described. The observations reported in the paper correspond to the Taylor-Saffman fingering, the characteristics of the two fluids being consistent with an intermediate behavior between capillary and viscous fingering according to Lenormand's classification [1]. Then in 3, the technique used for a robust detection of unstable interface is revealed, followed by the physical and morphological properties that characterize the unstable flow in 4. Finally, the conclusions are addressed in 5.

## 2 Experimental & acquisition systems

The aim of this work is to characterize the unstable bi-phasic interface evolving in a 2D granular medium under a given hydro-mechanical loading. A concise description of the experimental and acquisition systems is presented here to highlight the way of nucleation of such unstable percolation and the recording process.

The bi-phasic interface in this work consists on a moving air/water front drying an initially water-saturated Fontainebleau sand medium (mean diameter of sand grain = 210 µm). The surface tension between air and water is $7.28 \times 10^{-2}$ N.m$^{-1}$ and the viscosity of water is $1.008 \times 10^{-3}$ Pa.s while the viscosity of air is $1.788 \times 10^{-5}$ Pa.s. This flow occurs in a prismatic sample of 5cm height, 4cm width and 1.1cm thickness, confined by two transparent windows (biaxial conditions) providing the access for a
<figure><img src="image_3.png" /><figcaption>**Fig. 1** Schematic of the sample cell showing (i) in the close-up, the frontal window in semi-transparent gray characterized by a circular sapphire glass ensuring the visual access to the surface of the specimen, (ii) in the middle distance, the different parts of the assembly where the purple pieces represent the bottom and top porous stone ensuring a uniform injection through and an extraction of the fluids from the sample, while the yellowish pieces represent the membranes applying the lateral confinement to the specimen and (iii) in the background another window similar to the frontal one.</figcaption></figure>

high-resolution optical system to film the infiltration process at 30 Hz. The thickness of the sample has been chosen to be as small as possible to approach the conditions of two dimensional fluid flow typical of the Hele-Shaw cell as well as to provide a minimal surface capable to sustain the mechanical loading applied to the sample. A schematic of the sample cell is presented in Figure 1 illustrating its different parts. Three pressures are controlled: the mechanical confining pressure and two hydraulic pressures: the pressure of the fluid in-place (e.g. water) and the pressure of the injected fluid (e.g. air). Consequently, the air infiltration is driven by the hydraulic loading, called as capillary pressure $P_c$ and defined as the difference between air and water pressure but it also happens under a mechanical loading acting on the grains, known as initial effective stress $\sigma'_0$ and calculated as the difference between the applied confining pressure and the water pressure. An example of a sequence of images acquired through the drainage process occuring under $P_c = 30$ kPa and $\sigma'_0 = 120$ kPa, is illustrated in
Figure 2 (see also video (Online Resource 1) in the supplementary material). More details regarding the set-up and all the preparation steps are provided in [14].

<figure><img src="image_4.png" /><figcaption>**Fig. 2** Raw images (1950×6600 pixels) captured during the drainage process: bright color represents the region invaded by air and dark color represents the region that remains saturated by water. The time-step between these images is 0.5s. A video (Online Resource 1) showing the sequence of the acquired raw images is provided in the supplementary material.</figcaption></figure>

## 3 Technique for the detection of unstable interface

The detection of the interface between the defending wetting and the invading non-wetting phase during the considered drainage process is of paramount relevance in order to have a characterization of the gas infiltration process as much detailed as possible. In particular measuring, all along the time evolution process, the interface length, the interface velocity, the interface corrugation and the interfacial length density, quantitatively characterizes the gas fingering structure and enlightens the characteristic scales of the infiltration of gas through the considered granular medium.
Taking into account the specific characteristics of the drainage process of an initially saturated granular medium, say: (i) the existence of a diffuse transition zone between the fully wet and the fully dry material, as thick as fine grains, (ii) the ramification of the transition zone and the possible pinch off of the branches, (iii) the coalescence of branches and (iv) the trapping of residual spots of the wetting fluid behind the non-wetting one, a new approach to the detection of the interface between the invading non-wetting and the defending wetting phase is developed. Based on a similar functional as the well-known Mumford-Shah one [15] for image segmentation, the proposed method allows to robustly account for the abovementioned topological transition of the flow. However, image segmentation works well when a significant contrast of gray levels between the different regions of the image is available and it only provides a description of the front in the deformed configuration thus ignoring the deformation of the grains. In our case, the raw images acquired during the infiltration process of air, lack in contrast between the two phases as reported in Figure 2. Further, if a significant remodeling of the granular medium occurs during the infiltration, the detected interface, defined in the deformed configuration, has to be mapped back to the reference configuration. As an alternative, the raw images are therefore replaced by pixelwise correlation error maps that are reconstructed from the full-field measures of displacement: first, a global Digital Image Correlation (DIC) approach is applied with *f* illustrated in Figure 3a as the reference image. In *g*, illustrated in Figure 3b as the deformed image, the gray level variation occurs not only due to the motion of the grains but also to the infiltration of gas. As DIC assumes that the gray level variation is only due to particle motion, it is run in the following, between *f* and *g* that is an element-wise corrected deformed image based on matching the brightness and contrast parameters with the reference image, in a similar way as reported in [16], more details are provided in [17, 18]; consequently, the correct displacement *ũ* is measured (Figure 3c shows the map of the correlation error defined as |*f*(*x*) − *g*(*x* + *ũ*)| which is
<figure><img src="image_5.png" /><figcaption>Fig. 3 Insight on the process towards generating amplified correlation error maps: (a) reference image f where sand medium is fully water-saturated, (b) raw deformed image g where sand medium is partially saturated with poor contrast between air/water, (c) correlation error map obtained from the difference between the raw gray level of the reference image f and the modified gray level $\tilde{g}$ of the deformed image taking into account the output displacement $\tilde{u}$ from the DIC calculation and (d) correlation error map obtained from the difference in raw gray levels between the reference image f and the deformed image g but always using the displacement $\tilde{u}$ to highlight the change of contrast due solely to air infiltration and not to grain rearrangements; this map is used to generate the amplified correlation error map.</figcaption></figure>

negligible almost everywhere, indicating the convergence of the DIC approach). Then, the correlation error $|f(x)-g(x+\tilde{u})|$, illustrated in Figure 3d, is taken as an alternative to the raw image (Figure 3b) for further analyses in particular for the interface detection. The $\tilde{g}(x + \tilde{u})$ is in particular defined at each iteration of the process conduction at the evaluation of the displacement $\tilde{u}$, via the following equation:

$$ \tilde{g}(x + \tilde{u}) = C(g(x + \tilde{u}) - \mu_d) + B \quad (1) $$
with $C = \sigma_{ref}/\sigma_{def}$ being the ratio of the standard deviations of the gray level within the corresponding elements and $B = \mu_{ref}$ the mean value of gray level in the reference element, similarly $\mu_d$ is the mean value of gray level in the deformed element.

The reliance on correlation error maps for interface detection is inspired from a similar approach adopted in the identification of fracture paths starting from the field of correlation errors, see [19]. To further increase the contrast between the saturated and unsaturated zone, k-means clustering algorithm [20], with two clusters, is first applied to the map shown in Figure 3d, to retrieve the intensity of the invaded zone, then the correlation error is normalised with respect to this intensity and suitably amplified (by 100 in the specific case). The final correlation error images, defined in the reference image frame, are illustrated in Figure 4 (see also video (Online Resource 2) in the supplementary material).

<figure><img src="image_6.png" /><figcaption>**Fig. 4** Correlation error images (1950 × 6600 pixels) with amplified contrast: bright color represents the region invaded by air and dark color represents the region saturated by water. The time-step between these images is 0.5s. A video (Online Resource 2) showing the sequence of the correlation error images is provided in the supplementary material.</figcaption></figure>
As previously mentioned, a similar approach as the one proposed in the semi-
nal paper by Mumford-Shah, and later developed by Ambrosio and Tortorelli [21],
is adopted introducing the order parameter $d$ ($0 < d \le 1$) in order to track and
quantitatively analyse the interface between the saturated and the dry zone. The
order parameter $d$ regularizes the sharp transition from the aqueous to the air phase
through a narrow strip of porous space. The $d$ spatial distribution is identified as that
minimizing the following energy functional:

$$
E(d) = \frac{1}{2} \int_{\Omega} (1-d)^2 E_r \, d\Omega + E_r^c \int_{\Omega} \frac{d^2}{2} + 2\ell^2 \nabla d \cdot \nabla d \, d\Omega \quad (2)
$$

where $E_r$ is the pixel map of the absolute value of the amplified correlation error, $\ell$ the
length over which the transition from 0 to 1 of $d$ is smoothed, $E_r^c$ the threshold value of
$E_r$ above which $d$ is expected to be higher than 0 and $\Omega \subset \mathbb{R}^2$ the domain that encloses
transient gray level change. This method was originally thought by Mumford & Shah
[15] to detect discontinuities due to edges and shadows in a 2D image. This is obtained
minimizing a suitable functional with respect to $u$, representing the approximation of
the image $g$ within a set of disjoint connected open subsets $R_i$ of a planar domain
$R$, each one with a piece-wise smooth boundary, and with respect to $\Gamma$ which is the
union of the part of the boundaries of the $R_i$ inside $R$, that represents the line through
which $u$ is allowed to be discontinuous. The minimization with respect to $\Gamma$ was later
regularized by Ambrosio & Tortorelli [21] introducing a scalar variable, the above
mentioned order parameter $d$, so that the minimizers of the functional with respect to
$d$ is different from zero only in a small neighborhood of $\Gamma$ when the characteristic length
$\ell$ tends to zero. The formulation proposed in [21] was devoted to capture sharp (but
smoothed over an internal length parameter) transition between two material orders.
The energy that promotes this transition can be related to strain energy or any other
quantity driving the change of order. Herein, we used the correlation error map with
an enhanced contrast indicating that the material point has been invaded by gas to promote the change of order (d changing from 0 to 1). This is exactly what Equation 2 indicates considering that the discontinuity of the Mumford-Shah/Ambrosio-Tortorelli problem is in our context the interface between the saturated and the dry zone.

This minimization is solved in its discrete format by using a finite element discretization inherited from the DIC measurements (i.e. a regular grid of square-shaped quadrilateral first order elements). The term involving $E_r$ is integrated using a sum over the pixels inside the elements because $E_r$ is defined pixel-wise. Conversely, a standard Gaussian quadrature using 4 Gauss points per element can be used for the second term in Equation 2.

<figure><img src="image_7.png" /><figcaption>Fig. 5 Output from the interface detection algorithm for unstable interface: a map of the scalar $d$, issued from $ℓ = 40$ pixels.</figcaption></figure>

After writing the stationarity of $E$ with respect to $d$, $d$ is obtained solving a linear system. It provides a regularized representation of the saturated/unsaturated interface which will be further exploited to analyse quantitatively the properties (length, corrugation,...) of this interface and the evolution of these last properties. This analysis
will be conducted using the correlation error maps for all the images in the sequence
during the unstable gas infiltration. Figure 5 represents the output maps of d, issued
from ℓ = 40 pixels for the considered sequence of images; the length of the transition
zone ℓ over which the analysis is performed is slightly larger than the characteristic
size of a FB sand grain (35 pixels) and is spread over 8 square finite elements of 5 pixels
each. Incidentally it is worth to notice that performing this analysis using the corre-
lation errors allows for detecting the sought interface in the reference configuration of
the analyzed sample while it would have been obtained in the deformed configuration
if performed using the raw images. This avoids solving an inverse mapping problem
for retreiving the interface position in the reference configuration. It is worth to point
out that this process implicitely reveals the fully desaturated region as in the case of
the bottom of the finger. The gray level in the desaturated zone is homogeneous and
brighter than in any other region associated to d smaller than 1.

Using this interface detection algorithm, for each time-step, a *d* map is calculated as shown in Figure 5, then a set of contours is identified by selecting a value for *d*. In the presented cases, a *d* = 0.25 is chosen as a compromise for having an interface describing the reliefs of the drained zone and located within the partially saturated transition zone. Figure 6 illustrates the set of contours detected at each timestep, for *d* = 0.25. As can be seen, the raw output from this selection requires some post-treatment in order to distinguish the advancing front (in blue) from non-simply connected or even non-connected parts (in black). The observed fronts can be characterized as fractal. Their fractal dimension is estimated as 1.29 by the box-counting method [22].

## 3.1 Post-treatment process for interface detection

As in segmentation process, where choosing a threshold value on the gray level implies
that a binarized image corrupted by noise shows up, similarly, choosing a value of d
does not result in a single contour separating two distinct zones, but into non-simply
<figure><img src="image_8.png" /><figcaption>**Fig. 6** The contours from the interface detection procedure, obtained for a $d = 0.25$. The blue lines represent the evolving main interface through time, called as the root interface while the black ones representing the closed contours are going to be discarded.</figcaption></figure>

connected and/or possibly non-connected domains. In order to overcome the difficulties intrinsic at the treatment of multi-connected surfaces, an algorithm is developed in 2D to distinguish automatically the unconnected open/closed contours from the main root. The output of this algorithm is illustrated in Figure 4 where the black lines indicate the fragments that will be discarded and the blue one denotes the main root. To explain the way the algorithm treats artefacted images and to test its capability to identify internal boundaries or boundaries of non-connected parts of the domain, a complex scenario issued from the case illustrated in Figure 7a and which has led to a set of detected contours, as that depicted in Figure 7b, is addressed.

The main outlines of the algorithm are reported in Algorithm 1. Being deduced from a discretization grid, the contours map associated to a given value of the scalar $d$, can be converted into an explicit geometrical representation of the interface. The first step is to select the open contour that has the maximum $y$ position. This condition can be further restrained to avoid tricky open-fragments or to collect the two “river
<figure><img src="image_9.png" /><figcaption>**Fig. 7** Example illustrating a correlation error image of complex, ramified interface in (a) and main steps to sculpt the instability, at a given time *t*, with all the details: (b): contours from a selected *d* value, (c): detection of the top contours in blue, (d): binarized image initiated from the top contour, (e): further separation of contours into sub-categories and (f): binarized image for (white) drained / (black) undrained zones.</figcaption></figure>

bank" of the finger, if this last is crossed by the top border of the ROI (see Figure 7c). Then, a closed contour is constructed by connecting the blue contour with the base of the ROI, serving for the generation of a binarized image, see e.g. Figure 7d. A sorting is then conducted on the rest of the red contours of Figure 7c to divide them into two categories: the first one completes the fingering morphology while the second one is located outside the target and so is neglected (colored in magenta in Figure 7e). The first category contributes either in spotting a (red) hole inside the finger to account for trapped drops of the wetting fluid against the non-wetting fluid invasion or in adding (light blue) fragments to the main interface to be counted later on in the total length, as illustrated in Figure 7e.
**Algorithm 1** Resume of interface detection

Select a *d* value for contours detection

Pick the contour(s) having maximum *y* position

Construct a binary image from the previous top contour(s) selected

Sort the rest of the contours according to their location with respect to the white
zone of the binary image:

**if** Inside **then**

**if** On the border of the ROI **then**

Contribute to the formation of the finger border

**else**

Contribute to the localization of undrained zones

**end**

**else**

Do not contribute

**end**

# 4 Physical & Morphological properties

The main finger being detected, further operations can be performed to retrieve con-
crete and quantitative information about the moving interface for all the aquired
images such as the interface length, the invaded area, the interface velocity, the inter-
facial length density, as defined by [23] and the interface corrugation, as defined by
[24].

## 4.1 Physical properties

A propagating interface is characterised by a set of physical properties that are mainly
driven by the applied hydro-mechanical loading. In the following, the interface length,
the underlying area and the growth property are addressed.
4.1.1 Interface length & invaded area

Once the interface of the finger is identified via the previous algorithm, its discretiza-
tion is assumed based on a set of points of known cartesian coordinates. In order to
calculate the length of the unfolded interface, a sequence of local curvilinear coordi-
nates is then defined based on the length of segments connecting each two consecutive
points, all along the interface. A reference point (s = 0) is set for one extremity of the
blue interface in Figure 6. In the case where the root interface is separated in discon-
nected branches by the image border as in Figure 7e, a reference point (s = 0) is defined
for each blue part. Finally, the total length of the finger is the sum of the lengths of the
different segments constituing the interface or merely the final curvilinear coordinate.

Concerning the invaded area, as already mentioned, the propagating air does not
succeed in draining all the saturated zone; some spots remain saturated; using the
previous algorithm the distinction between the invaded zones and the resisting ones
is ensured, as shown in Figure 7f where the black spots inside the large white domain
hide the undrained zones. The invaded area is computed by counting the number of
white pixels in the binarized image.

4.1.2 Interface velocity

The interface velocity depends both on the micro-structure of the medium and on
the applied hydro-mechanical loading: the hydraulic loading monitored by either the
imposed capillary pressure or the gas inflow and the mechanical loading controlled
by the confining pressure and indirectly by the effective stress. In the literature, the
interface tip velocity, as proposed in [4], is a macroscopic quantity describing only
the advancing velocity of the tip of the front, in the vertical direction. Recording the
maximum y position of the interface of each acquired image allows to deduce the jump
of the tip (Δy) between two consecutive images taken at gap of 1/ν [s], ν being the
acquisition frequency. Thus the interface tip velocity is ν × Δy. However, the velocity
of a curve moving in the 2D space or of a surface moving in the 3D space is correctly
defined as a normal velocity, say a vector field oriented as the outward unit normal
to the considered object of co-dimension 1, as reported e.g. in [25]. In addition, the
interface does not evolve in a unique direction: the front tip can cease and an expansion
of the interface can be triggered simultaneously at the most vulnerable spot of the
front. This can be clearly seen in Figure 8 that summarizes the history of the invading
front by assigning a color for the growth areas at each time step; the non-invaded zone
is represented by the white color while the invaded zone is painted by the gradually
colored time-steps. A measure of the local velocity is quantified to highlight the local

<figure><img src="image_10.png" /><figcaption>**Fig. 8** Element-scale map (1 element = 5 pixels) showing the temporal evolution of the finger formation: the gradually coloured palette indicates the time-step of the growth areas.</figcaption></figure>

jumps of the interface via a customized algorithm. Two consecutive interfaces, obtained
at $t_{i-1}$ and $t_i$ are scanned. It is important to note that the two interfaces, are extremely
heterogeneous, in terms of the intrinsic fluctuations of the profile, and significantly different one from the other; however the level of image resolution does not allow to capture the infiltration of the gas through the pores of the granular medium and the retraction of liquid meniscii. So this measurement is neither at the representative elementary volume (REV) scale nor at the pore scale. It solely measures the normal displacements between two consecutive interfaces at the considered meso-scale.

The interface velocity is calculated based on the following criterion: any point of the interface at time $t_i$ is regarded as the transported along the normal to the interface at time $t_{i-1}$ of a point staying on the latter. However, the local properties such as the curvature and the normal direction influence this searching process. Therefore, in order to further enrich the collected data, a similar process can be done but in the reverse direction by following the normal direction drawn on the (i) interface. These two

<figure><img src="image_11.png" /><figcaption>Fig. 9 Tracking process for the quantification of a velocity profile, where in (a) the normal is plotted on the (i-1) interface while in (b) the normal is plotted on the (i) interface; the time-step between the two interfaces is 1/30 [s]. In (c), the velocity magnitude data acquired from both analysis is plotted as function of the curvilinear abscissa of the (i-1) interface.</figcaption></figure>

ways complement each other and contribute to capturing the spontaneous appearance of the irregularities especially when sharp changes of curvature occurs between $t_{i-1}$ and $t_i$. An example issued from the direct application of the aforementioned process
is illustrated in Figure 9 where the tracking process in both directions provides a
complete description of the burst: in Figure 9a, the tracking process is based on the
field of normals (black lines) to the (i-1) interface while in Figure 9b, the normals are
those to the (i) interface. The magnitude of the interface velocity computed using both
the abovementioned approaches and expressed as function of the curvilinear abscissa
of the (i-1) interface, $s_{i-1}$, are plotted in Figure 9c: this last indicates negligible
velocity everywhere except at the bulge formation which occured at a speed of 60
mm/s. The calculated displacement behind this local velocity is nothing else than
the depth of this bulge, which is around 300 pixels; a similar result will be obtained
following a completely different approach in Section 4.2.1.

The algorithm developed to fit the normal velocity consists on finding the intersec-
tion point between each normal plotted on the (i − 1) interface with the (i) interface
when considering the combination (i − 1) → (i), and conversely when considering the
combination (i) → (i − 1). In the following, Algorithm 2 describes an iterative search-
ing process for the intersection points, based on the combination (i − 1) → (i). First, a
window of a certain size is defined and centered on each discretised point on the (i − 1)
curve to enclose the moving part of the interface. Within this window, a set of points
P of the discretised (i) interface is collected. The equation of the normal line is then
computed at each point of the (i − 1) interface; however to capture its intersection with
the fragment of the (i) interface enclosed within the window, this last is not defined by
a continuous equation but by the set of discretising points P . Therefore, the equation
of the line connecting any two consecutive points of P is defined and its intersection
with the normal line is computed. If the intersection point falls outside the connecting
segment formed by the considered couple of points of P , the iterative process con-
tinues with the next couple of points and so on. If the intersection point is located
within the connecting segment, its curvilinear abscissa is interpolated and recorded to
be used as the starting point in the searching process within the next window.
**Algorithm 2 Resume of interface velocity**

Consider two consecutive detected interfaces $(i-1)$ & $(i)$

Select the size of the sliding window centered over the $(i-1)$ interface of $\mathcal{N}$ discretised points $\mathbb{D} = \{\mathcal{D}_1, \dots, \mathcal{D}_N\}$

for each $\mathcal{D}_i \in \mathbb{D}$ do

Find the set of points of the $(i)$ interface within each window $\mathcal{P} = \{P_{i1}, \dots, P_{ik}\}$

Define the equation of the normal line generated at $\mathcal{D}_i$

$j=1$

while $j \le k-1$ do

Define the equation of the line connecting $P_{ij}$ and $P_{i(j+1)}$

Compute the intersection point $\mathcal{P} \cap$ between the two lines

if $\mathcal{P} \cap [P_{ij}P_{i(j+1)}]$ then

Compute the velocity $V = \nu \times [\mathcal{D}_i \mathcal{P} \cap]$

break (solution is found)

else

$j=j+1$

end

end

end

At the end of the process, the two normal displacements fields obtained, from both directions when needed, are assigned to the $(i-1)$ interface and therefore the velocity profile of the $(i-1)$ interface is computed as can be seen in Figure 9c. In parallel, a second example illustrating the outcomes from the developed interface velocity algorithm is provided in Figure 10 treating another two consecutive interfaces within the same experiment. After the points allocation process between $(i-1)$ and $(i)$ interfaces in Figure 10a, the velocity vector distribution is plotted in Figure 10b along with the velocity profile in Figure 10c. One can conclude that the velocity of the bursts
<figure><img src="image_12.png" /><figcaption>**Fig. 10** Another example for interface velocity quantification, in (a) the points allocation from (i-1) interface to (i) interface, in (b) the velocity vector distribution and in (c) the velocity profile.</figcaption></figure>

is variable all along the interface and from one time-step to another which confirms the relation between the interface velocity and the micro-structure of the medium under the constant imposed capillary pressure. Moreover, the case reported in Figure 9
occurs after around one minute from that in Figure 10 which underlines the fact that
the interface can cease and destabilize far from the front tip.

## 4.2 Morphological properties

Morphological properties describe the reliefs of an interface by highlighting the depth of the local bifurcations via the interfacial length density analysis, the macroscopic corrugation amplitude and the ondulation frequency. These properties are mainly influenced by the heterogeneity of the granular medium and are detailed in the sequel.

### 4.2.1 Interfacial length density

In 3D images of a partially saturated porous medium, the interfacial area is defined as the area of the meniscus between the non-wetting and wetting phase per unit volume [23, 26]. This quantity is an indicator for the compactness of the ganglion formations [27] and improves in the description of partially saturated medium besides the saturation degree parameter [28].

Imaging the 2D surface of the sample and tracking the moving interface on low-contrasted patterns does not allow to capture the separation between the solid and fluid phases. Therefore, the separation line identified as in Section 3.1 is considered as the interface between the two fluid phases and an analysis aiming at characterizing the interfacial length distribution along the air-water interface is developed. The goal is to bring out the characteristics of the structure of the finger front, and in particular to identify (i) the placement of bulges, indicating possible ramification of the finger, as the case in [29] reporting the sequence of bursts during a drainage process, together with (ii) their dendritic depth. To this aim a square window, of size $W$, is moved translating its centroid along the detected ondulated interface. The fragments of the interface residing within each cell are identified and their cumulative length $L_{cum}$ per cell is calculated and normalised with respect to cell's diagonal $D_W = W\sqrt{2}$. The interfacial length $I_L$ per cell is defined herein as $I_L = L_{cum}/D_W$. It is to be noted that
due to the vicinity of the two “river banks” of the finger, in particular at the finger tip,
this study is conducted, separately, on the two branches which are labelled later-on by
*left* and *right* parts, in order to avoid swallowing non-concerned fragments within the
sliding cell and to distinguish the contribution of the bulges which have nucleated on
the *left* and *right* borders. The *left* part is colored in magenta in Figure 11 and the *right*
part in cyan and they are split based on the maximum *y* position on the whole curve.
Different window sizes have been tested; apparently depending on the heterogeneity
of the depth of the bulges, with respect to the finger body, windows of different sizes
are needed to classify the heterogeneities. The aim of this study is to carry out this

<figure><img src="image_13.png" /><figcaption>**Fig. 11** Ondulated interface obtained from *l* = 40 pixels and *d* = 0.25, divided into two parts: *left* and *right* parts. A particular attention is drawn to the *left* part where local analysis is carried out in the sequel on two different sets of bulges: the first set indicated by the dashed black rectangles, is characterized by small size of bulges with respect to the second set which is defined by points {a,b,c} from one side and {d,e,f} from other side.</figcaption></figure>

classification. Mean and standard deviation of the interfacial length distribution are
therefore calculated assuming the window characteristic size as a parametrization.
Their evolution as functions of the window size is reported on Figure 12 for the two parts of the finger. $\mu_{IL}$ exceeding 1 indicates that local patterns, longer than $D_W$, are dominant. This interpretation is complemented by the evolution of the standard

<figure><img src="image_14.png" /><figcaption>**Fig. 12** Mean and standard deviation of the interfacial length function of window size for the *left* and *right* parts of the finger, in (a) and (b) respectively.</figcaption></figure>

deviation. When it increases along with the window size, this reflects a significant fluctuation in the interface morphology, say a transition from smooth fragments into sharply bended fragments, whereas when it decreases along with the window size, this means that the characteristic size of the bulges has been exceeded and as long as the cell size increases, the cell grabs almost the same texture while travelling on the interface. A qualitative description of the interfacial length as a function of the window size all along a defined bulge shape is summarized in Figure 13. When a small window size is employed, as in panel 1, a constant interfacial length is going to be found, which in the limit corresponds to zero standard deviation. If the window size increases, the texture of the bulge starts to be more valued, the interfacial length will exceed 1 and the standard deviation grows, as in panel 2. However the excessive increase of the window size will drive again to a constant interfacial length and a drop in the standard deviation, see panel 3.
<figure><img src="image_15.png" /><figcaption>**Fig. 13** Qualitative characterization of the interfacial length as a function of the window size, 1: small, 2: medium and 3: large. *s* indicates the abscissa along the interface.</figcaption></figure>

The result obtained for the *left* finger border in Figure 12a shows three different slopes of the standard deviation within the intervals of the window characteristic size *W* (100, 300), (300, 550) and (550, 750) pixels. The window sizes *W*={100, 300, 550, 750} pixels are considered as critical since they are the turning points for the change of slope in the standard deviation evolution. On the other side, the result obtained for the *right* border in Figure 12b, describes a different behavior. A decrease in the standard deviation is observed after the peak at *W* = 300 pixels. So the largest window needed to capture the heterogeneities on the *right* side is 300 pixels.

Focusing on the *left* border where the variation of the standard deviation as a function of *W* indicates the presence of heterogeneities of different characteristic sizes, the complete evolution of the interfacial length along the curvilinear abscissa is drawn in Figure 14a only for the critical window sizes. For the first window size tested, *W* = 100 pixels, the profile indicates the showing up of mild fluctuations that turn out to sharp fluctuations with the increase of *W*. One can notice that at the positions *s* pointed by the dashed black rectangles, the interfacial length *Ic* exceeds 1 for a certain *W* and then drops when larger cells are used. The identification of these positions on the finger structure presented in Figure 11, provides a justification of this observation
considering that at these locations the size of the heterogeneities is relatively small compared to other irregularities and so it is tightly enclosed within a cell of $W = 300$ pixels whereas the use of larger cells provides a measure which is not affected by these small details. However, the measurement of the interfacial length at positions {$d, e, f$}

<figure><img src="image_16.png" /><figcaption>**Fig. 14** Interfacial length function of curvilinear abscissa for several window sizes applied on the *left* border in (a) and on the *right* border in (b). Some specific bulges on the *left border*, indicated in Figure 11 by the two black rectangles and the letters {a,b,c,d,e,f} are also shown in panel (a), to emphasize the change in interfacial length as function of window size at these particular locations.</figcaption></figure>

provides a behavior contrary to the previous one: the bigger the cell is, the larger the amount of included bended fragment and the higher the interfacial length will be. This interpretation is valid up to a certain limit of $W$ which once exceeded, implies a backward trend. The morphology of the corresponding bulge, identified in Figure 11 explains the reason behind the “efficiency” of large cells. An in-depth investigation on the interfacial length as a function of the window sizes at specific positions where bold bumps are formed and denoted by {$a, b, c$} and {$d, e, f$} in Figure 11, is presented in
Figure 15. For the first bulge, a plateau is reached at *W* ranging from 300 to 550 pixels

<figure><img src="image_17.png" /><figcaption>**Fig. 15** Interfacial length function of window size at critical locations defining the bulges on the *left* part of the finger, which are depicted in Figure 11 by {a,b,c,d,e,f}.</figcaption></figure>

while for the second bulge, a plateau starts at *W* = 400 pixels followed by a peak at
*W* = 850 pixels. One can conclude that the first bulge is the main contributor for the
first change of slope of the standard deviation evolution reported in Figure 12a, while
the second bulge is the principal cause of the second change of slope. Therefore, the *left*
border of this finger can be characterised by three scales of the bulge characteristic
size: small size detected at *W* = 300 pixels corresponding to a bulge amplitude of
150 pixels (4 grains of FB), regular size detected at *W* = 550 pixels corresponding
to an amplitude around 300 pixels (8 grains of FB) and finally big size detected at
*W* = 750 – 800 pixels corresponding to an amplitude around 400 pixels (12 grains
of FB). Using a window size *W* = 750 px allows consequently to observe the main
characteristics of the micro-structure of the left part of the interface. The amplitude of
the bulges corresponds to the half of the window size, as the window is drawn centered
on the discretised interface. Contrary to what has been seen for the *left* part, the use
of large cells is not necessary for the *right* part as reported in Figure 14b, which is
consistent with the fact that the largest characteristic size of the bulges which form on
the right hand side of the finger is smaller than that relative to its left part. Therefore,
for the right border, one can consider a cell of characteristic size $W$ of 100–150 pixels corresponding to 3–4 grains of FB.

### 4.2.2 Corrugation amplitude

The granular medium through which the bi-phasic flow percolates is responsible of the dendritic structure of the finger sides. The definition of a corrugation amplitude $A$, describing the perturbation of the current ondulated interface with respect to a reference flat 1D form, has been proposed in [24, 30]; this last being defined as the standard deviation of the normal coordinate $y_ξ$, of each point $ξ$ of the ondulated surface measured with respect to the reference curve. Inspired by these studies as well as by the methods proposed in [31] for 2D interface, an estimation of this quantity is proposed herein to evaluate the ondulation of the air-water interface, evolving in a 2D space. The idea is to generate a “smooth” 2D contour $\bar{C}$, that preserves the shape of the corrugated one $\tilde{C}$ but with less relief-details, as an alternative to the referred flat interface in the 1D problem. This is achieved by manipulating the parameter $ℓ$ that controls the length of the regularization in the computation of $d$ via the minimization of the functional in Equation 2 and thus the smoothness of the detected contour: the largest the $ℓ$ is, the less details are detected. However, generating a smooth interface by increasing the regularization parameter $ℓ$ tends to shift outward the level-sets of $d$. To avoid the contact of the selected smooth interface with the borders of the ROI and subsequently its undesired fragmentation, the algorithm of interface detection is supplied by an enlarged input domain of correlation error image generated by virtually replicating the column of border pixels on the left and right hand sides of the finger as well as the row border pixels at the top, as illustrated in Figure 16a. So, in addition to the first front detection that delivers the corrugated interface $\tilde{C}$ for $ℓ_{\tilde{C}}=40$ pixels, a second one is launched with the same parameters except a larger $ℓ$ and on virtually larger input correlation error images in order to extract a smoother interface. The
<figure><img src="image_18.png" /><figcaption>**Fig. 16** Steps of the construction of smooth interface and its signed distance map, (a): enlarged correlation error image generated from the initial ROI pointed by the dashed rectangle, (b): level-set $d$ map calculated for $l = 200$ pixels and the smooth interface $\bar{C}$ in red is selected for $d = 0.05$ and (c): signed distance $\Psi$ map generated with respect to the smooth interface in red and the corrugated interface $\tilde{C}$ in blue resides within the signed distance area. The first image is in pixels scale while the last two maps are illustrated in element scale: (732 x 1434) elements of 5 pixels each.</figcaption></figure>

coarse regularization length is $l_{\bar{C}} = 200$ pixels. The $d$ value to be selected should provide a smooth interface that incorporates the whole ondulated interface with a particular attention for the tip of the finger. Herein, $d_{\bar{C}}$ is 0.05 (smaller than $d_{\tilde{C}}=0.25$) and the corresponding smooth interface is colored in red on the level-set map of $d$ in Figure 16b. Once the smooth interface is provided, the goal is to determine the distance between this latter and the corrugated one. A signed distance $\Psi$ is defined over $2 \times \eta$ elements in the normal direction to the curvilinear abscissa along $\bar{C}$, using a fast marching method on the regular grid defined by the DIC finite element mesh [32]. In this way, the first $\eta$ elements outside the domain delimited by $\bar{C}$ are assigned with negative values according to their distance from the smooth interface and similarly, the last $\eta$ elements inside the domain delimited by $\bar{C}$ with positive values; $\eta$ being defined as the minimum offset to $\bar{C}$ delivering a zone of signed distance covering all the reliefs of $\tilde{C}$, see Figure 16c. Having a defined signed distance map with respect to $\bar{C}$, an
interpolation is applied on the contour $\tilde{C}$ to deduce its signed distance $\Psi_{\tilde{C}}$ to $\bar{C}$. Indeed not all the points of $\tilde{C}$ contribute to the calculation of the corrugation amplitude. A sorting is applied: for each point on the smooth interface, a corresponding point on the corrugated interface is chosen if the three following conditions are satisfied: (i) the searched point has to be on a straight line normal to $\bar{C}$ with a given tolerance, (ii) in the case where more than one point belongs to the same normal direction, the one with the smallest distance is picked; this occurs when a concavity inversion of the interface profile exists on the normal direction and (iii) the index of the point to be picked has to be larger than that of the predecessors. Similarly to what already done in Section 4.2.1, a distinction between the left and right borders is taken into account for both interfaces. The result of this sorting is shown in Figure 17a where the black lines link the considered couples of points: each point of $\bar{C}$ with its counterpart on $\tilde{C}$. Figure 17a also provides a visualization of the appropriate selection of the value of $\ell_{\bar{C}}$ to

<figure><img src="image_19.png" /><figcaption>**Fig. 17** Estimation of the corrugation amplitude: (a) the output of the points selection process where the black lines link the couple of points: the one of $\bar{C}$ with its corresponding on $\tilde{C}$, this plot is in element scale, (b) the signal showing the evolution of the signed distance with respect to the curvilinear abscissa defined on the smooth interface. The magenta and cyan scars delimit the intervals considered for the calculation of the corrugation amplitude.</figcaption></figure>
be used so that the sought interface is not only smooth but it is capable to accurately
describe the reliefs of the corrugated curve. In addition, this approach allows us to
transform the geometrical interface to exploitable data. The evolution of the signed
distance with respect to the curvilinear abscissa defined on the smooth interface is
represented on Figure 17b for the left and right parts. Based on the definition of the
corrugation amplitude, the standard deviation of the signed distances is calculated
over the intervals of points delimited by the magenta and cyan scars on Figure 17b to
eliminate the uncertain measurements of the borders. The estimated corrugation for
the left part $A_{left}$ is around 90 pixels while for the right part $A_{right}$ is around 78 pixels.
It is worth to mention that the calculated corrugation amplitude underestimates the
actual ondulation, being defined from a reference non-rectilinear interface that holds
the signature of the ondulated interface. This amplitude reflects solely a macroscopic
value of the ondulation magnitude especially in the case of 2D heterogeneous ondulated
interface, where arbitrary major singularities break the somehow “periodic” topology
dominated by the minor fluctuations. Therefore, this can justify the fact that the
corrugation of the left border is slightly higher than the right one since the left border
is marked by the appearance of more significant irregularities.

This approach based on the definition of the level-set with respect to a “smooth”
interface has allowed to extract the distance between two interfaces, the smooth one
and the corrugated one. Conversely, calculating the distance between two consecu-
tive corrugated interfaces labelled (*i* − 1) and (*i*) as per Section 4.1.2, using the same
approach, is not an option because the definition of the signed distance will be lim-
ited over a small η due to the significant curvature of the corrugated interface. This
procedure is therefore effective for characterizing the topology of the interface but not
for quantifying the interface velocity.
### 4.2.3 Ondulation frequency

The signed distance profile reported in Figure 17b can be analysed as a signal via the Fast Fourier Transform algorithm (FFT), in order to retrieve the most redundant periods and deduce the approximate distances separating the main bulges formations. The graph showing the amplitudes of each period (T) of the Fourier domain is represented in Figure 18a for the left and right borders in logarithmic scale. A very dense set of small periods can be observed indicating the noisiness of the signed distance signal and indirectly that of the detected interface. A truncation is applied on the delivered period domains by referring to the previous analysis performed on the interfacial length density of Section 4.2.1. The window size for which the maximum standard deviation of the interfacial length is calculated (see Figure 12) is considered as the lower boundary of this truncation. This window size values the largest texture of the interface which could possibly be at least equivalent to one period. To delimit the upper limit, periods close to the total length of each smooth border, are as well discarded. The retained period intervals are shown in Figure 18a with solid and dashed lines for the left and the right border respectively while the removed sets of periods are pointed by square and round symbols. The interval of periods taken into account for the left border ranges between 750 and 3000 pixels (equivalent to 21 to 85 grains of FB) while for the right border between 300 and 2800 pixels (equivalent to 8 to 80 grains of FB). By comparing these values to the size of the reliefs on the ondulated interface, one can notice that the irregularities on the left border are further apart but more pronounced than those on the right border which are closer together. However, for both cases, the wavelength of irregularities is on the order of a few dozen grains at most, highlighting the fact that the interface morphology is driven by the granular structure.

The signals corresponding to the retained or removed domains are generated via the Inverse Fast Fourier Transform algorithm (IFFT) and are referred as filtered or
noisy signals respectively. In Figure 18c, the filtered signals in magenta and cyan for the left and right borders respectively and the noisy signal in black are plotted together with the original signal in blue. Using the filtered signal and the normal vectors to the smooth interface, one can reconstruct the filtered interface reported in Figure 18b. Apparently it describes the main reliefs of each border that are well detected by the considered interval of periods.

# 5 Conclusions

Being interested in the tracking of a bi-phasic flow percolating a granular medium, an effective algorithm capable of dealing with highly perturbated interfaces undergoing bifurcations/coalescences through the infiltration process, has been established. This algorithm coupled with a post-treatment operation has led to the automatic detection of the history of the main interface. Concerning the input to this algorithm, the acquired images through the air-water infiltration has led to low contrasted images which are critical for any segmentation-like approaches. For this reason, a turn around has been found which consists on taking advantage of the result of a DIC calculation, the correlation error maps, to feed the interface detection algorithm. Consequently, the detected interface is characterized by a geometrical description determined in the reference configuration. Then, an in-depth investigation on the detected interface is performed to extract the maximum information regarding the physical and morphological properties that will serve in future works for the analysis of the coupling between different hydro-mechanical loadings and the interface behavior. Particularly, considering a given hydro-mechanical loading, it was obvious that the local velocity of the interface characterizing the nucleation of local bursts depends solely on the microstructural factors such as the heterogeneity of the medium. In addition, referring to the morphological analysis performed on the left and right borders of the developed finger,
<figure><img src="image_20.png" /><figcaption>Fig. 18 Post-treatment of the signed distance signal with the FFT: (a) log-scale graph showing the output of FFT with periods sorting in two parts for both borders: the main part contains the main periods that will contribute to the reconstruction of a filtered interface and the removed part contains the set of small and large periods that will not be involved in the contour reconstruction, (b) the main parts represent the filtered signal of signed distance while the removed parts account for the noisy signal that has been removed; in blue is the original signal and (c) in magenta and cyan are reported the reconstructed borders of the contour based on the set of main periods selected; this figure is in element scale where each element is 5 pixels.</figcaption></figure>
the characteristic sizes of the bulges have also shown to be driven by the microstruc-
ture of the medium. And finally, as mentioned previously in Section 4.1.2 & 4.2.1, the
interface detection analysis is impacted by the presence of the solid grains so it will be
worthwhile to conduct experiments where the refractive index of the dyed-fluids [33]
matches that of the granular material in order to reconstruct fluid percolation paths
and grain reorganization through the thickness of the sample. Further developments
will be devoted to exploiting this technique.

**Supplementary information.** The Supplementary Material for this article can be
found online at: (link)

**Acknowledgements.** The authors would like to acknowledge the support of the French National Research Agency (ANR), project STOWENG (Project-ANR-18-CE05-0033) and the Carnot MERS institute, project TCIN (investigation de la sTabilité des sols dunaires sous chargement environnementaux par Corrélation d'Image Numérique).

**Declarations**

**Authors contribution**

R.A.N., J.R. and G.S. contributed in the analysis of the results. J.R. developed the raw
interface detection algorithm and R.A.N. developed all the post-processing algorithms
and prepared the figures. The first draft of the manuscript was written by R.A.N. and
all authors commented on previous versions of the manuscript. All authors read and
approved the final manuscript.

**Conflict of interest**

The authors declare no competing interests.
## Ethical Approval

Not applicable.

## Funding

This study is supported by the French National Research Agency (ANR), project STOWENG (Project-ANR-18-CE05-0033) and the Carnot MERS institute, project TCIN (investigation de la sTabilité des sols dunaires sous chargement environnementaux par Corrélation d'Image Numérique).

## Availability of data and materials

Data can be made available by the authors upon reasonable request.

## References

[1] R. Lenormand, Liquids in porous media. Journal of Physics: Condensed Matter 2(S), SA79 (1990)

[2] P.K. Sinha, C.Y. Wang, Pore-network modeling of liquid water transport in gas diffusion layer of a polymer electrolyte fuel cell. Electrochimica Acta 52(28), 7936-7945 (2007)

[3] F. Guo, S.A. Aryana, An experimental investigation of flow regimes in imbibition and drainage using a microfluidic platform. Energies 12(7), 1390 (2019)

[4] G. Løvoll, Y. Méheust, R. Toussaint, J. Schmittbuhl, K.J. Måløy, Growth activity during fingering in a porous hele-shaw cell. Physical Review E 70(2), 026301 (2004)

[5] R. Holtzman, M.L. Szulczewski, R. Juanes, Capillary fracturing in granular media. Physical review letters 108(26), 264504 (2012)
[6] F.K. Eriksen, R. Toussaint, K.J. Måløy, E.G. Flekkøy, Invasion patterns during two-phase flow in deformable porous media. *Frontiers in Physics* **3**, 81 (2015)

[7] K. Culligan, D. Wildenschild, B.S.B. Christensen, W. Gray, M. Rivers, Pore-scale characteristics of multiphase flow in porous media: A comparison of air-water and oil-water experiments. *Advances in Water Resources* **29**(2), 227–238 (2006)

[8] P. Sharma, P. Aswathi, A. Sane, S. Ghosh, S. Bhattacharya, Three-dimensional real-time imaging of bi-phasic flow through porous media. *Review of Scientific Instruments* **82**(11), 113704 (2011)

[9] M.J. Dalbe, R. Juanes, Morphodynamics of fluid-fluid displacement in three-dimensional deformable granular media. *Physical Review Applied* **9**(2), 024028 (2018)

[10] Y. Méheust, G. Løvoll, K.J. Måløy, J. Schmittbuhl, Interface scaling in a two-dimensional porous medium under combined viscous, gravity, and capillary effects. *Physical Review E* **66**(5), 051603 (2002)

[11] A.T. Krummel, S.S. Datta, S. Münster, D.A. Weitz, Visualizing multiphase flow and trapped fluid configurations in a model three-dimensional porous medium. *AIChE Journal* **59**(3), 1022–1029 (2013)

[12] J.M. Campbell, D. Ozturk, B. Sandnes, Gas-driven fracturing of saturated granular media. *Physical Review Applied* **8**(6), 064029 (2017)

[13] D. Ozturk, B. Sandnes, A study of three-phase fracturing in granular media using high-speed imaging. *Journal of Porous Media* **22**(8) (2019)

[14] R. Al Nemer, G. Sciarra, J. Réthoré, Drainage instabilities in granular materials: A new biaxial apparatus for fluid fingering and solid remodeling detection.
[15] D.B. Mumford, J. Shah, Optimal approximations by piecewise smooth func-
tions and associated variational problems. Communications on pure and applied
mathematics **XLII**, 577–685 (1989)

[16] M.A. Sutton, J.J. Orteu, H. Schreier, Image correlation for shape, motion and deformation measurements: basic concepts, theory and applications (Springer Science & Business Media, 2009)

[17] R. Al Nemer, J. Réthoré, G. Sciarra, Quantification of localised strains induced within a granular medium by a biphasic flow via digital image correlation. Strain p. e12480 (2024). https://doi.org/10.1111/str.12480

[18] R. Al Nemer, Effect of two-phase fluid percolation on remodeling of geo-
materials. Theses, École centrale de Nantes (2023). URL https://theses.hal.
science/tel-04353684

[19] S. Roux, F. Hild, H. Leclerc, Mechanical assistance to dic. Procedia IUTAM 4, 159-168 (2012)

[20] J.A. Hartigan, M.A. Wong, Algorithm as 136: A k-means clustering algorithm.
Journal of the royal statistical society. series c (applied statistics) **28**(1), 100–108
(1979)

[21] L. Ambrosio, V.M. Tortorelli, On the approximation of free discontinuity prob-
lems. Bollettino dell'Unione Matematica Italiana 7 pp. 105-123 (1992)

[22] K. Foroutan-pour, P. Dutilleul, D.L. Smith, Advances in the implementation of the box-counting method of fractal dimension estimation. Applied mathematics and computation **105**(2-3), 195–210 (1999)
[23] E. Dalla, M. Hilpert, C.T. Miller, Computation of the interfacial area for two-fluid porous medium systems. Journal of contaminant hydrology **56**(1-2), 25-48 (2002)

[24] A. Vreme, F. Nadal, B. Pouligny, P. Jeandet, G. Liger-Belair, P. Meunier, Gravitational instability due to the dissolution of carbon dioxide in a hele-shaw cell. Physical Review Fluids **1**(6), 064301 (2016)

[25] W. Mulder, S. Osher, J.A. Sethian, Computing interface motion in compressible gas dynamics. Journal of computational physics **100**(2), 209-228 (1992)

[26] K.A. Culligan, D. Wildenschild, B.S. Christensen, W.G. Gray, M.L. Rivers, A.F. Tompson, Interfacial area measurements for unsaturated flow through a porous medium. Water Resources Research **40**(12) (2004)

[27] K. Singh, B. Bijeljic, M.J. Blunt, Imaging of oil layers, curvature and contact angle in a mixed-wet and a water-wet carbonate rock. Water Resources Research **52**(3), 1716-1728 (2016)

[28] V. Joekar-Niasar, S. Hassanizadeh, A. Leijnse, Insights into the relationships among capillary pressure, saturation, interfacial area and relative permeability using pore-network modeling. Transport in porous media **74**(2), 201-219 (2008)

[29] A. Islam, S. Chevalier, I.B. Salem, Y. Bernabe, R. Juanes, M. Sassi, Characterization of the crossover from capillary invasion to viscous fingering to fracturing during drainage in a vertical 2d porous medium. International Journal of Multiphase Flow **58**, 279-291 (2014)

[30] C. Brouzet, Y. Méheust, P. Meunier, Co 2 convective dissolution in a three-dimensional granular porous medium: An experimental study. Physical Review Fluids **7**(3), 033802 (2022)
[31] M. Alava, M. Dubé, M. Rost, Imbibition in disordered media. Advances in Physics 53(2), 83–175 (2004)

[32] J.A. Sethian, Level set methods and fast marching methods: evolving interfaces in computational geometry, fluid mechanics, computer vision, and materials science, vol. 3 (Cambridge university press, 1999)

[33] A.N. Lembach, H.B. Tan, I.V. Roisman, T. Gambaryan-Roisman, Y. Zhang, C. Tropea, A.L. Yarin, Drop impact, spreading, splashing, and penetration into electrospun nanofiber mats. Langmuir 26(12), 9516–9523 (2010)