# Hysteresis phenomenon within unsaturated granular assemblies: Capillary forces and matric suction

Nabil Younes, Antoine Wautier, Richard Wan, Olivier Millet, François Nicot

► To cite this version:

Nabil Younes, Antoine Wautier, Richard Wan, Olivier Millet, François Nicot. Hysteresis phenomenon within unsaturated granular assemblies: Capillary forces and matric suction. Physical Review E, 2026, 113 (3), pp.035401. 10.1103/1qdn-hnns. hal-05534016

HAL Id: hal-05534016

https://hal.inrae.fr/hal-05534016v1

Submitted on 3 Mar 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
# Hysteresis phenomenon within unsaturated granular assemblies: capillary forces and matric suction

Nabil Younes,<fnref n="1" />,<fnref n="2" />,<fnref n="3" />,<fnref n="4" />,* Antoine Wautier,<fnref n="3" /> Richard Wan,<fnref n="4" /> Olivier Millet,<fnref n="2" /> and François Nicot<fnref n="5" />

In this paper, we present a numerical approach to simulate the condensation and evaporation processes of capillary bridges within granular materials. The formation and dynamics of capillary bridges are captured using the phase-field-based Lattice Boltzmann model (LBM) to solve Allen-Cahn and Navier-Stokes equations, while evaporation and condensation processes accounted by local changes in position of the capillary interfaces. Our model captures the emergence of hysteresis as a result of irreversible geometric transitions—such as bridge coalescence and snap-off—without prescribing any constitutive relation between suction and saturation degrees. This change in capillary regimes arises naturally from the interface dynamics and not solved by the LBM. In particular, we can capture and analyze the discontinuities in capillary forces, when capillary bridges merge or split within small elementary assemblies of three or four spherical particles. Having validated our numerical results for the above elementary assemblies, a poly-dispersed granular assembly composed of 1,000 spherical grains is next addressed. Our simulations capture well-known condensation and evaporation hysteresis phenomenon while offering the possibility to inspect the underlying topology of air/water cluster for the same water saturation along different hydraulic paths.

## I. INTRODUCTION

The presence of capillary bridges is ubiquitous in our daily life, especially within granular materials [1]. A capillary bridge takes place due to the interactions between liquid-gas-solid phases creating an apparent cohesión within a granular material. This apparent cohesión occurs because capillary bridges give rise to attractive forces between grains [2–6]. This is exactly the case for a beach sandcastle that gains strength in partially saturated conditions [7–10]. It is, however, well-known that the distribution of capillary bridges in unsaturated granular materials depends strongly on the hydraulic loading path—saturation/desaturation or condensation/evaporation. This phenomenon is known as hysteresis and arises from physical mechanisms at the pore scale, including the merging and splitting of liquid bridges, the nucleation of new bridges, and the rearrangement of grains [11–17]. It can be observed at any scale, even within small assemblies of grains where only merging/splitting conditions matter.

Vast numerical, theoretical, and experimental studies have been carried out to investigate the capillary forces as well as capillary bridges in small granular assemblies, i.e., two, three, and four spherical grains or between plane and plane, plane and sphere or poly-dispersed spheres [18–25]. From a numerical perspective, it is not easy to handle the merging of capillary bridges between three or more particles. For instance, Miot *et al.* (2021) [26] proposed to account for merging/splitting of capillary

bridges using an energy minimization principle and geometry constraints. Their findings were in accordance with experimental data reported in the literature. Nevertheless, the work is limited to a specific grain configuration and the coalescence of several capillary bridges could not be implicitly handled. Recently, Benseghier *et al.* (2022) [27] have instead used the Shan-Chen Lattice Boltzmann Method (LBM) model to calculate capillary forces between doublets and triplets of particles. The results reported in this work agreed with experimental and theoretical data, especially since they were able to model the merging of capillary bridges for the triplet. However, the use of the Shan-Chen model for capillary bridge problems requires simulations to be carried out in extreme temperature and pressure conditions, i.e., $T = 0.75T_c$, where $T_c = 373.946°C$ is defined as the critical temperature below which a liquid and its vapor coexist. Still in the same family of models, the phase-field-based LB models have proved to be powerful in dealing with multi-phase fluids [28–32]. The phase-field-based LB model was used by Younes *et al.* (2022) [33] to study capillary bridges for doublets and triplets of particles. Not only was good accordance reported between their findings and experimental results, but also the coalescence was implicitly handled without using any geometrical or energetic criterion.

It has been experimentally proven by Gras (2011) [34] that the merging and splitting of capillary bridges have a huge influence on the capillary forces. In the latter study, an experimental investigation in which a granular system composed of three particles and one merged capillary bridge is subjected to evaporation has been carried out. As the capillary bridge shrinks due to evaporation, it eventually breaks up into three small isolated capillary

<footnote n="1">¹LMGC, Univ Montpellier, CNRS, Montpellier, France</footnote>

<footnote n="2">²LaSIE, UMR CNRS 7356, University of La Rochelle, La Rochelle, France</footnote>

<footnote n="3">³INRAE, Aix-Marseille University, UMR RECOVER, Aix-en-Provence, France</footnote>

<footnote n="4">⁴Department of Civil Engineering, Schulich School of Engineering, University of Calgary, Calgary, Canada</footnote>

<footnote n="5">⁵University Savoie Mont-Blanc, ISTerre, Chambéry, France</footnote>

<footnote>* nabil.younes@umontpellier.fr</footnote>
bridges between each pair of particles. The transition from the funicular regime to the pendular regime results in a discontinuity in capillary forces applied on the top grain, with a reported decrease of 30%. Several years later, Gagneux and Millet (2016) [35] have theoretically demonstrated that once isolated capillary bridges of a triplet of mono-dispersed particles merge all together, a sharp increase of 30% is noted which is in accordance with the work of Gras (2011) [34]. In addition, Wang et al. (2022) [36] conducted an experimental investigation with three spherical particles to study the effects of sphere size, relative position, liquid volume, viscosity, and separation rate on the peak of capillary forces.

In addition to the geometric complexity of capillary bridge interactions, modeling the underlying phase-change processes—evaporation and condensation—poses its own set of challenges. These processes have been extensively studied using thermodynamically rich models, which often involve complex couplings and high computational costs. These approaches—such as those based on the modified Cahn-Hilliard equation [37, 38], pseudo-potential phase-field models [39, 40], and conservative Allen-Cahn formulations [41]—have successfully captured phase-change dynamics in various configurations. However, they often involve coupling with temperature fields, solving additional probability density functions, or managing spurious velocities at interfaces, which significantly increase computational cost and complexity. Motivated by the need for a computationally efficient yet physically meaningful approach, we propose a simplified model for simulating condensation and evaporation. In this work, our approach is based on the Allen-Cahn phase-field model coupled with the Navier-Stokes equations, without introducing thermodynamic variables or temperature coupling. This simplification allows us to efficiently simulate condensation and evaporation processes while capturing the key geometric transitions—such as bridge coalescence and rupture. To implement condensation and evaporation, we introduce a phenomenological modification of the equilibrium phase-field profile, whereby the interface is shifted by a lattice unit—Eqs. (9b) and (9a)—allowing us to simulate wetting and drying processes without thermodynamic coupling.

Unlike macroscopic models that attempt to eliminate hysteresis by enriching the set of state variables—such as interfacial area, mean curvature, or Euler characteristic [42–44]—our approach operates at the mesoscopic scale, where hysteresis emerges naturally from the physics of interface evolution. The phase-field formulation does not impose any constitutive relation between capillary pressure and saturation. Instead, it resolves fluid–fluid interfaces and their interactions with solid boundaries, allowing irreversible geometric transitions—such as snap-off, coalescence, and bridge rupture—to occur spontaneously. These transitions lead to path-dependent fluid configurations and discontinuities in capillary forces, which can be compared with experimental observations and inspected from a micro-mechanical point of view. Throughout this

paper, the term *hysteresis* is used in the conventional sense employed in soil and granular mechanics under unsaturated conditions, namely the path dependence that appears when quantities such as capillary force or suction are represented in reduced variables, e.g., force-volume or suction-saturation, rather than in the full Minkowski space description of the capillary state [42, 43]. As it will be shown in the paper, this hysteresis is attributed to the irreversible geometric and topological transitions of the liquid phase coalescence and bridge rupture resolved directly by the mesoscopic interface dynamics.

The primary objective of this paper is to investigate the sharp jumps and drops in capillary forces that occur when capillary bridges merge or rupture due to condensation and evaporation, respectively. We demonstrate that these transitions do not occur at the same saturation degree upon reversal of paths, leading to the emergence of the *hysteresis* phenomenon. This behavior is also observed in a larger fixed granular assembly, confirming its relevance beyond small-scale configurations.

The paper is organized as follows. In Section II we review the LBM model used for simulating capillary forces [33, 45]. We then address the new approach for condensation and evaporation processes in Section III. Next, the merging and splitting of capillary bridges are addressed for triplets of particles and tetrahedrons, in Section IV, showing that this model can capture the discontinuity in capillary forces. After validating the approach on small grain clusters, we apply it to a complex granular assembly composed of 1,000 spherical grains in Section V. Finally, we conclude with a summary and provide possible avenues for future research.

## II. NUMERICAL MODEL

In this paper, the Allen-Cahn phase-field-based LBM model is adopted given its performance in handling of liquid-gas interfaces [28–31]. The main idea of this model is to solve the modified Navier-Stokes equation whereby pressures and velocities are determined, coupled with the Allen-Cahn equation that tracks the capillary interface [28, 29, 33, 45]. The model was developed by the authors and implemented on GPU architecture, achieving approximately 50 times speedup on an NVIDIA Quadro RTX 5000 compared to an Apple M3 CPU for identical cases. As a reference, a full simulation involving the complete condensation or evaporation of a sample containing 1,000 grains—as described in Section V—takes approximately 24 hours on GPU.

### A. Governing macroscopic equations

The conservative Allen-Cahn equation reads

$$ \frac{\partial \phi}{\partial t} + \nabla \cdot (\phi \mathbf{u}) = \nabla \cdot [M(\nabla \phi - \lambda \mathbf{n})] \quad (1) $$
where **u** is the advection fluid velocity, *M* is the so-called mobility, **n** is the outward unit-normal vector as indicated in Figure (1b), and $\phi(\mathbf{r}, t)$ is the phase-field parameter that allows us to identify the regions occupied by the fluids. The fluid densities $\rho(\mathbf{r}, t)$ is linearly linked to the phase-field parameter following this relation

$$ \rho(\mathbf{r}, t) = \phi(\rho_{\ell} - \rho_{g}) + \rho_{g} \quad (2) $$

where $\rho_{\ell}$ and $\rho_{g}$ are the liquid and the gas densities, respectively. Therefore, $\phi(\mathbf{r}, t) \in \{0; 1\}$ corresponds to the bulk fluids (0 for gas and 1 for liquid). And, $\phi \in ]0; 1[$ corresponds to the interface. It is worth noting that at the mechanical equilibrium the distribution of $\phi(\xi)$ is a smooth Heaviside function, as depicted in Figure (1). The mathematical form of the phase-field parameter at the equilibrium is $\phi^{\text{eq}}(\xi)$

$$ \phi^{\text{eq}}(\xi) = \frac{1}{2} + \frac{1}{2} \tanh \left( \frac{2\xi}{W} \right) \quad (3) $$

where $\xi$ is the orthogonal local variable crossing through the interface—see Figure (1). In Eq. (1) the parameter $\lambda$ is defined as $\lambda = d\phi^{\text{eq}}/d\xi$.

The modified multi-phase Navier-Stokes equations for incompressible fluids read

$$ \rho \left( \frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} \right) = - \nabla p + \nabla \cdot \tau + \mathbf{F}_s + \mathbf{G} \quad (4) \\ \nabla \cdot \mathbf{u} = 0 $$

where $\tau = \rho\nu(\nabla\mathbf{u} + \nabla\mathbf{u}^T)$ is the viscous stress tensor ($\nu$ is the kinematic viscosity), $p$ is the pressure, $\mathbf{G}$ is the body force, and $\mathbf{F}_s$ is the surface tension term defined as follows

$$ \mathbf{F}_s = \mu_{\phi} \nabla \phi \quad (5) $$

where $\mu_{\phi}$ is the chemical potential $\mu_{\phi} = 4\beta\phi(\phi - 1)(\phi - \frac{1}{2}) - k\nabla^2\phi$ with

$$ k = \frac{3}{2}\gamma W; \quad \beta = \frac{12\gamma}{W} \quad (6) $$

and $\gamma$ is the surface tension and $W$ is the interface width. The detailed LBM schemes along with the boundary conditions (bounce-back and wetting condition) can be found in [28, 29, 33, 45].

## B. Capillary forces

For any shape of the capillary bridge, the associated capillary force can be calculated as follows [33]

$$ \hat{\mathbf{F}}_f = \underbrace{\int_{\Omega} \sigma_f n_{\Omega} dS}_{\mathbf{F}_p} + \gamma \underbrace{\int_{\Gamma} m d\ell}_{\mathbf{F}_{\gamma}} \quad (7) $$

<figure><img src="image_2.png" /><figcaption>FIG. 1. (a) The smooth evolution of $\phi$ through the interface whose thickness is $W$ (smooth Heaviside function). $\phi_{\ell}$ and $\phi_g$ correspond to liquid and gas phase parameters, respectively. (b) Scheme of the smooth interface. $\phi \in \{0; 1\}$ corresponds to the bulk fluid, gas, and liquid, respectively. And $\phi \in ]0; 1[$ corresponds to the thick interface of width $W$.</figcaption></figure>

where $\Omega$ is the integrated domain, $\sigma_f = -pI + \tau$ is the stress fluid with $I$ is the identity tensor in $\mathbb{R}^3$, $n_{\Omega}$ is the outward unit normal vector of the domain $\Omega$, $\Gamma$ is the intersection of $\Omega$ with the capillary interface, $\gamma$ is the surface tension, and $m$ is the tangent vector of the gas-liquid interface and perpendicular to the contour $\Gamma$. Readers can refer to Younes et al. (2022) (2023) [33, 45] for more details about the formulation and the discretization of capillary force calculations.
### III. EVAPORATION AND CONDENSATION PROCESSES

Evaporation and condensation processes occur at the interface separating water from air, leading to a change in the interface's position. Hence, in a phase-field model, where the interface corresponds to a threshold value for $\phi$, eliminating or adding matter may be achieved by increasing or decreasing the value of $\phi$ around the equilibrium position of the interface. Formally, this is done by applying an interface-shifting method of the $\phi^{eq}(\xi)$ to become $\phi^e(\xi)$ or $\phi^c(\xi)$ to simulate evaporation or condensation, respectively; see Figure (2). At mechanical equilibrium, the evolution of the phase-field parameter $\phi$ through the interface has to fulfill the equilibrium condition: $\phi(\xi) = \phi^{eq}(\xi)$. Based on Eq. (3), the differentiation of $\phi^{eq}(\xi)$ yields

$$d\phi = \frac{4}{W}\phi(1 - \phi)d\xi \quad (8)$$

where $d\phi$ is the infinitesimal variation of $\phi$ at equilibrium and $d\xi$ is the infinitesimal variation of $\xi$ within the capillary bridge interface. As such, Eq. (8) provides the relationship between the change in $\phi$ and the change in the position of the interface.

<figure><img src="image_3.png" /><figcaption>FIG. 2. Profiles of condensed and evaporated liquid-gas capillary interfaces $\phi^c(\xi)$ and $\phi^e(\xi)$, respectively for $\Delta\xi = 1$.</figcaption></figure>

The proposed numerical technique is such that evaporation or the condensation of fluid can be imposed as follows

$$\phi^e(\xi) = \phi^{eq}(\xi) - \Delta\phi(\Delta\xi) \quad (9a)$$

$$\phi^c(\xi) = \phi^{eq}(\xi) + \Delta\phi(\Delta\xi) \quad (9b)$$

where $\phi^e(\xi)$ and $\phi^c(\xi)$ are the profiles of the phase-field parameters after applying evaporation and condensation processes, respectively.

Figure (2) shows the equilibrium $\phi^{eq}(\xi)$ profile (black curve) and the new curves after condensation (red) and evaporation (blue) for an imposed displacement of the interface by 1 lattice unit $\Delta\xi = \pm 1$ lu. As can be seen from Figure (2), the variations of $\phi$ are significant only in the liquid-gas interface, whereas, in the bulk fluids the variation is virtually insignificant and still satisfies $\phi \in \{0; 1\}$. The modification of the local equilibrium value by an order parameter is a phenomenological strategy designed to mimic the growth and retreat of capillary bridges during condensation and evaporation, respectively. It is worth noting that this approach is not derived from thermodynamic first principles, but is physically motivated by the observation that increasing or decreasing in water volume occurs at the fluid interface that moves outward and inward. However, this approach cannot account for the nucleation of new capillary bridges. Our goal is not to model the full complexity of phase transitions, but to reproduce realistic interface dynamics under wetting and drying conditions so that hysteresis phenomenon can emerge.

### IV. JUMPS IN CAPILLARY FORCES AND HYSTERESIS EFFECTS ON SMALL GRAIN ASSEMBLIES

In this section, the condensation/evaporation technique is implemented to capture the hysteresis effects in small clusters composed of three and four spherical grains. For all the triplets and quadruplets applications, only monodispersed particles are considered with a radius of $R = 4$ mm, with a size mesh $\Delta x = 7.5 \times 10^{-5}$ m. For consistency in comparing results, the following dimensionless forms are used for force, volume, and pressure:

$$F^* = \frac{F}{\gamma R}; V^* = \frac{V}{R^3}; P^* = \frac{PR}{\gamma}; D^* = \frac{D}{R} \quad (10)$$

#### A. Equidistant triplet of grains

The first case study deals with a triplet of particles forming an equilateral triangle with a separation distance $D^* = 0.025$, as illustrated in Figure (3). This separation distance is the minimum distance we can simulate, which corresponds to 1 voxel. The latter distance (equivalent to $\Delta x = 7.5 \times 10^{-5}$ m) is relatively small compared to the size of the particles and hence they can be assumed to be almost in contact. We also show in the same figure how capillary bridges are initialized at the mid-distance of each pair of spherical grains. This configuration approximately matches the one used in Gagneux and Millet (2016) [35] and Gras (2011) [34], in which case the particles were all touching. After initializing capillary bridges, the proposed condensation/evaporation process is numerically enforced and the resulting capillary force acting on the top particle is calculated via Eq. (7).
<figure><img src="image_4.png" /><figcaption>FIG. 3. Three capillary bridges located at the mid-distance of each pair of particles.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIG. 4. Capillary force evolution in terms of the water volume as the system is being condensed and evaporated for $\theta = 50^{\circ}$ and for a separation distance of $D^* = 0.025$. The capillary force is calculated on the top grain in the y-direction, see Figure (3). It is important to emphasize that $\Delta\xi = 1$ was chosen, as demonstrated in Appendix A, since any value of $\Delta\xi$ less than or equal to 1 has no impact on the results, while any value of $\Delta\xi > 1$ is not acceptable. The red and blue dotted arrows correspond to changes in the capillary force during the condensation and evaporation processes, respectively.</figcaption></figure>

Figure (4) depicts the capillary forces evolution when the granular system is subjected to condensation followed by evaporation. As the volume of water increases, when the triple (three-phase contact) lines of the isolated capillary bridges touch each other, coalescence takes place, resulting in the transition of the system from the pendular to the funicular regime. The merging occurs for a volume of $V^* \approx 0.14$ as seen in Figure (4) from state (a) to (b). It is important to note that at the coalescence, a discontinuity in capillary forces is observed with a relatively sharp increase of 29.9% in capillary forces, which is completely in agreement with the theoretical calculation of 30% re-

ported by Gagneux and Millet (2016) [35] for the same contact angle $\theta = 50^{\circ}$.

To better understand the physics behind the sharp increase in capillary forces at coalescence, the evolution of Laplace $F_p^*$ and surface tension $F_\gamma^*$ forces are plotted in Figure (5). In the pendular regime and as the water volume increases, it can be observed that $F_p^*$ decreases continuously. By contrast, the surface tension force $F_\gamma^*$ term increases because the wetting surface (solid-water interface) expands gradually around the particle. At coalescence, it is observed that the Laplace force suffers a significant jump whereby it increases, whereas the surface tension force shows a slight decrease also through a jump. It is the two forces when combined that give an overall increase of $\approx 29.9\%$ in capillary force, which coincides with the value reported by [35]. It is worth noting that the discontinuities in $F_p^*$ and $F_\gamma^*$ have two origins. First, based on Eq. (7), the surface tension force $F_\gamma^*$ will increase in direct proportion with the size of the wet contour $\Gamma$. Figure (6) depicts how the contours evolve before and after coalescence. It can be observed that the contour before coalescence (blue contour) slides into the middle to fill the air between capillary bridges at coalescence, causing a slight decrease in the wet contour $\Gamma$ (red contour), hence a decrease in the surface tension force. Second, $F_p^*$ increases drastically despite the volume $V$ changes slightly at coalescence. This drastic increase in $F_p^*$ is due to the change in the suction $\Delta P^* = P_a^* - P_w^*$ when the isolated capillary bridges become one big cluster, which is illustrated in Figure (7), where $P_a^*$ and $P_w^*$ correspond to the air and water pressures, respectively. It is highlighted here that the change in the suction $\Delta P^*$ is caused by the change in curvatures upon the air-water interface after and before coalescence—according to the Young-Laplace equation [20, 46].

When the water volume exceeds $V^* \approx 0.46875$, we choose to begin the evaporation process by using Eq. (9a). The merged capillary bridge shrinks as evaporation proceeds until the water volume reaches $V^* \approx 0.08125$ at step (c) in Figure (4). At this point, the capillary bridge spontaneously splits up into three smaller capillary bridges of equal size as seen in Figure (4) at step (d). This rupture event causes a sharp decrease of 30.88% in the capillary force as also obtained by [34].

It is important to note that in Figure (5), the force acting on the top grain with one single bridge or three distinct bridges is the same during either the condensation or evaporation path, which is in agreement with the fact that the wetting angle remains constant during the simulation. Only the volume at which bridges merge or split is different. The same observation is true for the variation of suction as shown in Figure (7).

To understand the underlying physics of this hysteresis, we show the snapshots of the merged capillary bridge as it is being evaporated until the rupture phase—from (a) to (c) of Figure (8). As the merged capillary bridge is being evaporated, the thickness $z^*$ of the capillary bridge in the $yz$-plane is being narrowed, as illustrated in Fig-
ures (8(c.1), (c.2), and (9). Once $z^*$ reaches zero, the single bridge splits ultimately into 3 separate bridges as precisely illustrated in Figure (9). As previously mentioned, coalescence takes place when the triple lines of distinct capillary bridges touch each other. Therefore, during condensation and evaporation processes, the geometrical conditions for the switch from 1 to 3 bridges and *vice-versa* are forcibly different, which explains the hysteresis. In short, the rupture of a capillary bridge involves a thinning mechanism as shown in Figure (10), whereas the merging of capillary bridges involves the meeting of triple lines as shown in Figure (11).

<figure><img src="image_6.png" /><figcaption>FIG. 5. The evolution of pressure force $F_p^*$ and tension force $F_γ^*$ of Eq. (7) in terms of water volumes while subjected to condensation and evaporation processes. The red and blue dotted arrows correspond to changes in forces during the condensation and evaporation processes, respectively.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIG. 6. Air-water interface contour before and after coalescence for $\theta = 50°$ and $D^* = 0.025$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIG. 7. Pressure difference $\Delta P^* = P_a^* - P_w^*$ evolution in terms of water volumes $V^*$ during condensation and evaporation processes, where $P_a^*$ and $P_w^*$ are dimensionless air and water pressures, respectively. The red and blue dotted arrows correspond to changes in the suction $\Delta P^*$ during the condensation and evaporation processes, respectively.</figcaption></figure>

## B. Influence of separation distances and wetting angle

Figure (12) shows the sharp relative capillary force $|\Delta F|$ [%] in terms of the separation distances during the condensation and evaporation processes, respectively. It can be observed that the further apart the particles are, the lower the capillary force jumps/drops. These findings are consistent with the fact that when the separation distance increases, the suction $\Delta P^*$ decreases as already reported in Mielniczuk et al. (2018) [20]. Nevertheless, the *hysteresis* phenomenon is still observed even when increasing $D^*$, for the same reasons explained previously.

The effect of the wetting angle $\theta$ is next investigated. Multiple simulations are carried out with $\theta \in [10; 70]°$ and for $D^* \in \{0.025; 0.05; 0.075\}$. Figures (13)a and (b) show the absolute difference $|\Delta F|$ [mN] at coalescence and rupture, respectively, for different separation distances $D^*$ and wetting angles $\theta$. These findings show that increasing the wetting angle decreases the absolute jump/drop of capillary forces either at coalescence or rupture. It is interesting to note that during evaporation for $D^* = 0.075$ and $\theta = 70°$, the merged capillary bridge splits up into two capillary bridges instead of three. Therefore, no discontinuity in capillary force, applied on the top particle, has been observed for this particular point. The reason behind this revolves around the fact that just before the rupture, the bottom triple line interface is above the line joining the centers of the bottom particles as shown in Figure (23). It is believed that this is due to one of these two reasons: (1) the configuration is theoretically symmetrical. However, the horizontal meshing is different than the oblique one, and (2) the
<figure><img src="image_9.png" /><figcaption>FIG. 8. Evaporation process steps from 1 to 3 for a contact angle $\theta = 50^{\circ}$. (a) 3D view, (b) 2D section view in $xy$-plane, (c) 2D section view in $yz$-plane. The bottom particles in (c) are in transparent mode to better visualize the merged capillary bridge.</figcaption></figure>

configuration (shape of the merged capillary bridge) in Figure (23b) might not be stable. Indeed, this is a hypothetical assumption because there is no criterion indicating whether this configuration is stable or not. However, (c) is stable while ($z^*$) did not reach 0 yet. Therefore, the merged capillary bridge undergoes a spontaneous transformation into two stable capillary bridges. We note that another hysteresis effect has been reported for $\theta = 70^{\circ}$ between evaporation and condensation processes.

### C. Quadruplet of particles - Regular tetrahedron

From the planar configuration of the previous subsections, we now move on to 3D structures of grains. The question is whether the sharp jump and drop in capil-

<figure><img src="image_10.png" /><figcaption>FIG. 9. Evolution of the capillary bridge thickness in $yz$-plane $z^*/R$ in terms of water volumes during the evaporation phase.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>FIG. 10. Thinning of the merged capillary bridge to become three isolated capillary bridges during the evaporation process in (a) the $xy$-plane and (b) the $yz$-plane.</figcaption></figure>

lary forces remain ≈ 30% even for 4 particles (non-planar systems) for $\theta = 50^{\circ}$? To answer this question, in this section, we examine a quadruplet of particles forming a regular tetrahedron with different separation distances, e.g., $D^* \in \{0.025; 0.05; 0.075\}$, along with $\theta = 50^{\circ}$. Figure (15) illustrates the initialization phase of the capillary bridges created at the mid-distance between particles. Capillary forces are hereafter calculated on the top
<figure><img src="image_12.png" /><figcaption>FIG. 11. Growth of the isolated capillary bridges to become one capillary bridge cluster at the meeting of triple lines during the condensation process in (a) the xy-plane and (b) in the yz-plane.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>FIG. 12. The evolution of the relative difference in capillary forces $|\Delta F|$ [%] in terms of separation distance $D^*$ for both condensation and evaporation processes for $\theta = 50^{\circ}$.</figcaption></figure>

particle along the z-direction as shown in Figure (15).

The first configuration to be considered is similar to the one discussed in Section IVA, with $D^* = 0.025$ and $\theta = 50^{\circ}$. In the same way as before, the condensation and evaporation processes are imposed successively. Figure (16) depicts the capillary force evolution in terms of water volumes. It can be observed that the capillary force, in this case, follows roughly the same tendency as the case of a triplet shown earlier in Figure (4). Although it has roughly the same tendency, here, capillary force is larger than in the triplet case. This is due to the fact that in the four-particle case, 3 capillary bridges are connected to the top spherical particle, with 1 capillary bridge more than the top grain of the three-particle case. The jump and drop in capillary forces are found to be 58% and 52% for condensation and evaporation, respectively (compared with 30% in Section IVA). The reason

<figure><img src="image_14.png" /><figcaption>FIG. 13. The evolution of the absolute difference in capillary forces $\Delta F^*$ in terms of the wetting angle $\theta$ for different separation distances $D^* \in \{0.025; 0.05; 0.075\}$, during (a) condensation and (b) evaporation processes.</figcaption></figure>

why the jump and drop of capillary forces are considerably higher than in the triplet case is that, in this case, three capillary bridges are attached to the bottom grain instead of one as was the case of the triplet. This forms a larger volume of water that contributes to the capillary force applied on the top grain after coalescence. Once again, as expected, the hysteresis phenomenon is observed—coalescence and rupture do not occur for the same volume.

Figure (17) depicts the evolution of the relative difference in capillary forces as a function of separation distances during condensation and evaporation processes. These findings are in line with the triplet case in which the relative difference in capillary forces $|\Delta F|$ [%] decreases when the separation distance increases. However, contrary to the triplet case, $|\Delta F|$ [%] is larger than 30% ($\ge$ 40%) for all three separation distances $D^* \in \{0.025; 0.05; 0.075\}$. Two effects can be pointed out: (1) more “dead water” not contributing to the
<figure><img src="image_15.png" /><figcaption>FIG. 14. Profiles of capillary bridges. (a) Merged capillary bridge. (b) Merged capillary bridge just before the rupture. (c) After the rupture of the capillary bridge. In (b) the bottom air-water interface is at the same level as the horizontal line joining the particles' centers.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>FIG. 15. Initialization of capillary bridges at the mid-distance of spherical particles for a quadruplet of particles.</figcaption></figure>

capillary force on the top grain before coalescence, i.e., the horizontal capillary bridges in the triplet case, (2) rupture later because of "3D" effects (the planar bridge gets thinner quicker).

## V. HYSTERESIS EFFECTS FOR A LARGE GRANULAR ASSEMBLY

In the previous section, we have examined how the hysteresis phenomenon affects the capillary force jumps and drops within small granular assemblies during condensation and evaporation, respectively. Our devised approach can next be extended to a large system composed of 1,000

<figure><img src="image_17.png" /><figcaption>FIG. 16. Capillary force evolution in terms of water volumes for condensation and evaporation processes for a wetting angle of $\theta = 50^{\circ}$ and for a separation distance of $D^* = 0.025$. The capillary force is computed on the top grain in the z-direction. The red and blue dotted arrows correspond to changes in the capillary force during the condensation and evaporation processes, respectively.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>FIG. 17. Relative difference in capillary forces $|\Delta F|$ [%] in terms of separation distances $D^*$ for both coalescence and rupture for a wetting angle of $\theta = 50^{\circ}$ for the quadruplet case.</figcaption></figure>

particles, all of which are kept fixed, allowing us to separate the hysteresis phenomenon arising from capillary bridge topology from the particle rearrangement induced by capillary effects. In this section, the global suction $\Delta P^*$ of the granular assembly is defined as the average of all local suctions $\Delta P^i$ at capillary menisci weighted by their respective volumes as in [45], i.e.

$$ \Delta P = \frac{1}{V_w} \sum_{i=1}^{n} \overbrace{(P_a - P_w^i)}^{\Delta P^i} V_w^i \quad (11) $$
with $V_w$ is the total water volume in the sample, n is the number of capillary bridges, $P_w^i$ is the pore water pressure within the $i^{th}$ capillary bridge, and $V_w^i$ is the volume of the $i^{th}$ capillary bridge. Here, the air pressure $P_a$ is negligible when compared to the negative water pressures $P_w^i$ within capillary bridges. The suction $\Delta P^*$ takes then positive values. The $P_w^i$ is the average water pressure of each individual voxel that makes up the $i^{th}$ capillary bridge. Note that thermodynamic equilibrium is not necessarily achieved because vapor fluxes are not modeled. As a result, suction may vary between individual bridges rather than being uniform over the REV. This is justified when vapor diffusion occurs on a much slower timescale than liquid distribution during wetting/drying. We define the dimensionless form of $\Delta P^*$ as follows

$$ \Delta P^* = \frac{\Delta P d}{\gamma} \qquad (12) $$

where $d$ [m] is a reference length scale (here taken $D_{50}$) and $\gamma$ [N.m$^{-1}$] is the water-air surface tension. The Particle Size Distribution (PSD) is shown in Figure (18). Granular assembly characteristics are summarized in Table I.

<figure><img src="image_19.png" /><figcaption>FIG. 18. (a) An illustration of the granular assembly composed of 1,000 spherical grains. (b) The Particle Size Distribution (PSD) of (a). The largest and the smallest grain diameters are 64.48 µm and 30.45 µm, respectively.</figcaption></figure>

<table><caption>TABLE I. Granular system statistical characteristics: $N_{grains}$ is the number of grains, $D_{max}/D_{min}$ is the packing polydispersity, $D_{50}$ is the grain diameter at which 50% of cumulative percent is reached, $C_u = D_{60}/D_{10}$ is the uniformity coefficient, and $\phi^p$ is the granular media porosity.</caption><thead><tr><th>N<sub>grains</sub> [-]</th><th>D<sub>max</sub>/D<sub>min</sub> [-]</th><th>D<sub>50</sub> [µm]</th><th>C<sub>u</sub> [-]</th><th>φ<sup>p</sup> [-]</th></tr></thead><tbody><tr><td>1,000</td><td>2.11</td><td>50.35</td><td>1.22</td><td>0.4624</td></tr></tbody></table>

The first step consists of initializing equal-sized isolated capillary bridges at the mid-distance of each pair of grains. Only those that are separated by 1.5 µm or 1 voxel in the LBM discretization, capillary bridges are

initialized as shown in Figure (19). Consequently, 628 isolated capillary bridges are created for this case resulting in an initial degree of saturation of $S_r = 0.0414$ which is defined as the ratio of the water volume to the void volume. From this point, the condensation process begins, leading to the growth of capillary bridges that enable the water to transit from the *pendular* regime to the *funicular*, and eventually to the *capillary* regime.

Figure (21) shows the suction evolution during condensation and evaporation processes in terms of degrees of saturation $S_r$. A sharp decrease is observed from $S_r \in [0.04; 0.17]$. Afterward, a change of slope is then observed until $S_r \approx 0.6$. Finally, the suction tends to zero as the system enters the capillary regime. The evaporation

<figure><img src="image_20.png" /><figcaption>FIG. 19. (a) Initialization of the capillary bridges at the mid-distance of spherical grains. (b) A zoomed snapshot of (a). It is important to note that due to the large distance between its neighboring spheres, there is no capillary bridge initialized around the red-dotted-circle particle in (b).</figcaption></figure>

process begins once the granular system is saturated—blue curve of Figure (21). Approximately the same trend is observed between both curves, but the evaporation curve is shifted upwards with respect to the condensation one in both capillary and funicular regimes. In other words, the *hysteresis* phenomenon is also observed within the suction $s$ of the big granular assembly.

A thorough analysis of the air bubbles' topology can be carried out using the Flood-Fill algorithm to reveal the physics behind the difference in suctions. Figure (20) represents screenshots of the air bubbles distribution within the sample for $S_r \approx 0.4$. It is very important to note that in Figure (20) the big air cluster has been omitted for better visualization, and only the trapped air bubbles are displayed. It can be easily seen that in (a) the number of air bubbles trapped between capillary bridges (encircled in red) is less than (b) which results in a larger curvature for (b) than for (a). This observation is consistent with the Young-Laplace equation stating that the higher the curvature, the higher the pressure difference $\Delta P^*$.

When the first evaporation process is done, the granular system is condensed and evaporated again, the green and black curves of Figure (22), respectively. It is remarkable that the red and green curves do not match, even though they both represent the condensation phase. This is due to the fact that the initial distribution of capillary bridges
<figure><img src="image_21.png" /><figcaption>FIG. 20. Air bubbles trapped within the granular system are colored in blue with $S_r \approx 0.4$ during (a) condensation and (b) evaporation.</figcaption></figure>

is not the same, as shown in Figure (20). On the contrary, the evaporation curves (blue and black curves) are quite similar because they were in the same condition (capillary regime) when the evaporation process began—only one geometrical configuration for full saturation $S_r = 1$.

## VI. CONCLUSION

In this paper, we introduced a new interface-shifting method within a phase-field LBM framework to simulate condensation and evaporation processes in open systems with the purpose of investigating the hysteresis phenomenon within partially saturated granular assemblies in terms of capillary forces and matric suction. It is important to note that a limitation of our current condensation approach is its inability to capture the nucleation of capillary bridges. Instead, it only models the growth or shrinkage of pre-existing liquid bridges, without accounting for their initial formation during condensation. Capillary forces hysteresis is inspected on several configurations of small granular assemblies: triplets and quadruplets of particles for different wetting angles and different separation distances. Multiple crucial conclusions are

<figure><img src="image_22.png" /><figcaption>FIG. 21. Evolution of suction $\Delta P^*$ in terms of degrees of saturation $S_r$ for condensation and evaporation processes.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>FIG. 22. Evolution of suction $\Delta P^*$ in terms of degrees of saturation $S_r$ for condensation and evaporation processes for multiple cycles.</figcaption></figure>

drawn as follows:

1. We have recovered the compelling theoretical and experimental result of the sudden jump of 30% in capillary forces when isolated capillary bridges coalesce during condensation for equilateral triplets of particles for $\theta = 50^\circ$.

2. When the same granular assembly is subjected to evaporation, we observe that the steep drop in cap-
illary forces does not occur at the same volume for
which the coalescence took place.

3. We highlighted the effect of wetting angles and separation distances on the jumps or steep drops of capillary forces. According to our research, for equilateral triangle configuration, we have found that capillary force jumps and drops decrease when increasing the separation distances between particles. However, generally speaking, the capillary force jumps and drops increase with the wetting angle.

4. Regarding the regular quadruplets of particles, we have found that the steep jump and drop range between 35 and 55%. This range depends on the separation distances and whether the system is experiencing a condensation or evaporation process. The same trend as the triplet case has been reported meaning that the relative difference in capillary forces decreases when the separation distance increases.

5. The new approach has been extended to a large granular assembly composed of 1,000 fixed spherical particles to investigate the hysteresis phenomenon regarding matric suction. The proposed model demonstrates that hysteresis arises in the absence of solid particle movements, and even and under identical wetting angles, due to topological variations of trapped air bubbles.

The phase-field LBM model, along with the new simple
method of condensation/evaporation, has demonstrated
its capacity to simulate wet granular assemblies and re-
produce the *hysteresis* phenomenon.

Our model, although simplified, successfully captures es-
sential features of capillary bridge dynamics and hystere-
sis in wet granular media. By operating at the meso-
scopic scale, it directly resolves interface geometries and
does not rely on macroscopic closure relations. This en-
ables the simulation of path-dependent behavior that
would otherwise be lost in averaged continuum mod-
els. While more elaborate thermodynamic models exist,
our approach offers a tractable and physically grounded
framework for studying capillary-induced cohesion and
mechanical transitions. In this study, we focused on
the mesoscopic capillary geometry and the associated
interface transitions, leaving aside several deeper pro-
cesses, such as liquid-vapor equilibrium and diffusion-
driven mass transfer (leading to pressure equalization
over time and Ostwald ripening), nucleation phenomena,
and roughness-induced contact-line effects. We note that
the hysteresis observed here is expressed in the practical
reduced variables commonly used for wet granular mate-
rials under unsaturated conditions, whereas introducing
extended geometric and topological representations (be-
yond the volume of water) may describe the same inter-
face evolution without exhibiting closed loops, e.g., same

volume but different surface area. This may help clar-
ify whether the hysteresis observed in practical variables
are purely apparent due to representational limitations or
instead reveal fundamental properties of the underlying
interface dynamics.

This study will be extended to address more complex problems by coupling the LBM with the Discrete Element Method (DEM), a combination that would enable a comprehensive evaluation of capillary stress within a granular assembly—including capillary-induced grain rearrangement—and thus deepen our understanding of wet granular materials at the microscale.
# ACKNOWLEDGMENTS

The authors would like to express their sincere gratitude to the Nouvelle-Aquitaine Region for funding the StabDigue project, the French National Centre for Space Studies (CNES) for having supported this work, and the IRN Geomech network for its valuable discussions.

## Appendix A: Condensation and evaporation processes

### 1. Sensitivity analysis of $\Delta\xi$ on the interface-shifting method of $\phi^{eq}$

The choice of $\Delta\xi$ affects the distribution of whether $\phi^c(\xi)$ or $\phi^e(\xi)$. It is, therefore, reasonable to carry out a sensitivity analysis of the rate of condensation-evaporation $\Delta\xi$ on the profiles of condensed and evaporated profiles. Figure (23) describes the evolution of condensed as well evaporated liquid-gas capillary interfaces when $\Delta\phi$ is added or subtracted. According to Figure (23), for $\Delta\xi \in \{0.5; 1\}$ lu (dashed and dotted curves), the condensed and the evaporated capillary interfaces comply with the fact that the minimum, as well as the maximum values of the evaporated and the condensed curves, are higher than 0 and lower than 1, respectively. When $\Delta\xi$ is increased to $\Delta\xi \in \{1.5; 2\}$ lu (pointed-dashed and short-dotted curves) the minimum and the maximum values of $\phi$ after evaporation or condensation are lower than 0 and higher than 1, respectively. Consequently, the choice of $\Delta\xi = 1$ lu can be adopted. It roughly moves the interface over a distance corresponding to the lattice discretization.

<figure><img src="image_24.png" /><figcaption>FIG. 23. Profiles of condensed and evaporated liquid-gas capillary interfaces for multiple values of $\Delta\xi$.</figcaption></figure>
[1] D. G. Fredlund and H. Rahardjo, *Soil mechanics for unsaturated soils* (John Wiley & Sons, 1993).

[2] F. Xiao, J. Jing, S. Kuang, L. Yang, and A. Yu, Capillary forces on wet particles with a liquid bridge transition from convex to concave, Powder Technology **363**, 59 (2020).

[3] Y. Tang, S. Cheng, et al., Capillary forces on a small particle at a liquid-vapor interface: Theory and simulation, Physical Review E **98**, 032802 (2018).

[4] M. Dörmann and H.-J. Schmid, Simulation of capillary bridges between particles, Procedia Engineering **102**, 14 (2015).

[5] B. Saint-Cyr, F. Radjai, J.-Y. Delenne, and P. Sornay, Cohesive granular materials composed of nonconvex particles, Physical Review E **87**, 052207 (2013).

[6] C. Semprebon, M. Scheel, S. Herminghaus, R. Seemann, and M. Brinkmann, Liquid morphologies and capillary forces between three spherical beads, Physical Review E **94**, 012907 (2016).

[7] M. Pakpour, M. Habibi, P. Møller, and D. Bonn, How to construct the perfect sandcastle, Scientific reports **2**, 1 (2012).

[8] A. Gans, O. Pouliquen, and M. Nicolas, Cohesion-controlled granular material, Physical Review E **101**, 032904 (2020).

[9] V. Richefeu, M. S. El Youssoufi, and F. Radjai, Shear strength properties of wet granular materials, Physical Review E **73**, 051304 (2006).

[10] J. Duriez, M. Eghbalian, R. Wan, and F. Darve, The micromechanical nature of stresses in triphasic granular media with interfaces, Journal of the Mechanics and Physics of Solids **99**, 495 (2017).

[11] W. J. Likos and N. Lu, Hysteresis of capillary stress in unsaturated granular soil, Journal of Engineering mechanics **130**, 646 (2004).

[12] R. Hilfer, Macroscopic capillarity and hysteresis for flow in porous media, Physical Review E **73**, 016307 (2006).

[13] H. Liu, Y. Ju, N. Wang, G. Xi, and Y. Zhang, Lattice boltzmann modeling of contact angle and its hysteresis in two-phase flow with large viscosity difference, Physical Review E **92**, 033306 (2015).

[14] A. Poulovassilis, The effect of the entrapped air on the hysteresis curves of a porous body and on its hydraulic conductivity, Soil Science **109**, 154 (1970).

[15] H. Pham, D. Fredlund, and S. L. Barbour, A practical hysteresis model for the soil-water characteristic curve for soils with negligible volume change, Geotechnique **53**, 293 (2003).

[16] D. Gallipoli, A hysteretic soil-water retention model accounting for cyclic variations of suction and void ratio, Geotechnique **62**, 605 (2012).

[17] H. Q. Pham, D. G. Fredlund, and S. L. Barbour, A study of hysteresis models for soil-water characteristic curves, Canadian Geotechnical Journal **42**, 1548 (2005).

[18] B. Mielniczuk, T. Hueckel, and M. S. El Youssoufi, Laplace pressure evolution and four instabilities in evaporating two-grain liquid bridges, Powder Technology **283**, 137 (2015).

[19] O. Pitois, P. Moucheront, and X. Chateau, Liquid bridge between two moving spheres: An experimental study of viscosity effects, Journal of Colloid and Interface Science **231**, 26 (2000).

[20] B. Mielniczuk, O. Millet, G. Gagneux, and M. S. El Youssoufi, Characterisation of pendular capillary bridges derived from experimental data using inverse problem method, Granular Matter **20**, 1 (2018).

[21] H. N. G. Nguyen, O. Millet, and G. Gagneux, On the capillary bridge between spherical particles of unequal size: analytical and experimental approaches, Continuum Mechanics and Thermodynamics **31**, 225 (2019).

[22] H. N. G. Nguyen, O. Millet, and G. Gagneux, Liquid bridges between a sphere and a plane-classification of meniscus profiles for unknown capillary pressure, Mathematics and Mechanics of Solids **24**, 3042 (2019).

[23] H. N. G. Nguyen, O. Millet, C.-F. Zhao, and G. Gagneux, Theoretical and experimental study of capillary bridges between two parallel planes, European Journal of Environmental and Civil Engineering, 1 (2020).

[24] J.-P. Wang, E. Gallo, B. François, F. Gabrieli, and P. Lambert, Capillary force and rupture of funicular liquid bridges between three spherical bodies, Powder Technology **305**, 89 (2017).

[25] A. Di Renzo, G. Picarelli, and F. P. Di Maio, Numerical investigation of funicular liquid bridge interactions between spherical particles, Chemical Engineering & Technology **43**, 830 (2020).

[26] M. Miot, G. Veylon, A. Wautier, P. Philippe, F. Nicot, and F. Jamin, Numerical analysis of capillary bridges and coalescence in a triplet of spheres, Granular Matter **23**, 1 (2021).

[27] Z. Benseghier, O. Millet, P. Philippe, A. Wautier, N. Younes, and E. Liberge, Relevance of capillary interfaces simulation with the shan–chen multiphase lb model, Granular Matter **24**, 1 (2022).

[28] H. Liang, J. Xu, J. Chen, H. Wang, Z. Chai, and B. Shi, Phase-field-based lattice boltzmann modeling of large-density-ratio two-phase flows, Phys. Rev. E **97**, 033309 (2018).

[29] H. Liang, H. Liu, Z. Chai, and B. Shi, Lattice boltzmann method for contact-line motion of binary fluids with high density ratio, Phys. Rev. E **99**, 063306 (2019).

[30] A. Fakhari and M. H. Rahimian, Phase-field modeling by the method of lattice boltzmann equations, Physical Review E **81**, 036707 (2010).

[31] A. Fakhari and T. Lee, Multiple-relaxation-time lattice boltzmann method for immiscible fluids at high reynolds numbers, Physical Review E **87**, 023304 (2013).

[32] A. Fakhari, T. Mitchell, C. Leonard, and D. Bolster, Improved locality of the phase-field lattice-boltzmann model for immiscible fluids at high density ratios, Physical Review E **96**, 053301 (2017).

[33] N. Younes, Z. Benseghier, O. Millet, A. Wautier, F. Nicot, and R. Wan, Phase-field lattice boltzmann model for liquid bridges and coalescence in wet granular media, Powder Technology, 117942 (2022).

[34] J.-P. Gras, *Approche micromécanique de la capillarité dans les milieux granulaires: rétention d'eau et comportement mécanique*, Ph.D. thesis, Université Montpellier II-Sciences et Techniques du Languedoc (2011).

[35] G. Gagneux and O. Millet, An analytical framework for evaluating the cohesion effects of coalescence between capillary bridges, Granular Matter **18**, 1 (2016).
[36] S. Wang, F. Liu, J. Cui, M. Miao, and C. Pu, Experimental study on the rupture behavior of the liquid bridge between three rigid spheres, Langmuir **38**, 13857 (2022).

[37] H. Safari, M. H. Rahimian, and M. Krafczyk, Extended lattice boltzmann method for numerical simulation of thermal phase change in two-phase fluid flow, Physical Review E **88**, 013304 (2013).

[38] H. Safari, M. H. Rahimian, and M. Krafczyk, Consistent simulation of droplet evaporation based on the phase-field multiphase lattice boltzmann method, Physical Review E **90**, 033305 (2014).

[39] R. Ledesma-Aguilar, D. Vella, and J. M. Yeomans, Lattice-boltzmann simulations of droplet evaporation, Soft Matter **10**, 8267 (2014).

[40] Q. Li, P. Zhou, and H. Yan, Improved thermal lattice boltzmann model for simulation of liquid-vapor phase change, Physical Review E **96**, 063303 (2017).

[41] M. Sugimoto, Y. Sawada, M. Kaneda, and K. Suga, Consistent evaporation formulation for the phase-field lattice boltzmann method, Physical Review E **103**, 053307 (2021).

[42] S. Schlüter, S. Berg, M. Rücker, R. Armstrong, H.-J. Vogel, R. Hilfer, and D. Wildenschild, Pore-scale displacement mechanisms as a source of hysteresis for two-phase flow in porous media, Water Resources Research **52**, 2194 (2016).

[43] J. E. McClure, R. T. Armstrong, M. A. Berrill, S. Schlüter, S. Berg, W. G. Gray, and C. T. Miller, Geometric state function for two-fluid flow in porous media, Physical Review Fluids **3**, 084306 (2018).

[44] A. M. Boelens and H. A. Tchelepi, The effect of topology on phase behavior under confinement, Processes **9**, 1220 (2021).

[45] N. Younes, A. Wautier, R. Wan, O. Millet, F. Nicot, and R. Bouchard, Dem-lbm coupling for partially saturated granular assemblies, Computers and Geotechnics **162**, 105677 (2023).

[46] T. Hueckel, B. Mielniczuk, and M. S. El Youssoufi, Adhesion-force micro-scale study of desiccating granular material, Géotechnique **70**, 1133 (2020).