# Permeability of foam-filled granular packing: Numerical modeling

Vincent Langlois, Ali Salame, Olivier Pitois, Adrien Petit, Bartholomé Soltner

► To cite this version:

Vincent Langlois, Ali Salame, Olivier Pitois, Adrien Petit, Bartholomé Soltner. Permeability of foam-filled granular packing: Numerical modeling. Physical Review Fluids, 2025, 10 (5), pp.053604. 10.1103/PhysRevFluids.10.053604. hal-05243816

HAL Id: hal-05243816

https://enpc.hal.science/hal-05243816v1

Submitted on 7 Sep 2025

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 Attribution 4.0 International License
# Permeability of foam-filled granular packing: numerical modeling.

Vincent Langlois,* Ali Salame, and Olivier Pitois†

Université Gustave Eiffel, Ecole des Ponts et Chaussées,
CNRS, Laboratoire Navier, F-77454 Marne-la-Vallée, France

Adrien Petit and Bartholomé Soltner

Ecole des Ponts et Chaussées, 6 & 8 avenue Blaise Pascal, 77455 Marne-la-Vallée, France
Abstract

In recent experiments [Phys. Rev. Fluids 8, 024302 (2023) and Phys. Rev. Fluids 9, 074003 (2024)], it was shown that the permeability of foam-filled granular packings is strongly influenced by the ratio of bubble size to grain size. This effect appears to be primarily driven by geometric factors and is independent of interfacial mobility. It was suggested that the liquid from the foam introduced into a granular medium redistributes towards the grain surfaces and the contact points between grains, contributing to the drying of the foam and significantly reducing its hydraulic permeability. To investigate this phenomenon in detail, we developed a numerical model of foam within a granular medium, which describes the liquid distribution within the pores and allows to calculate the liquid permeability. This model offers improved predictions of liquid redistribution compared to existing analytical models and closely matches permeability measurements. Moreover, our results quantify the contributions to permeability of all liquid components formed within the pores, particularly bulk and wall foam liquid channels.

# I. INTRODUCTION

Liquid foam has a wide range of applications in porous media, including petroleum and natural gas production [1, 2], soil remediation [3, 4], and tunneling [5, 6]. The physical behavior of foam in such media has been shown to depend largely on the ratio of bubble size to pore size [7]. When this ratio is low, bubbles primarily interact with each other rather than with the pore walls, exhibiting behavior similar to bulk foam. Conversely, at high size ratios (such as those encountered in oil recovery, where bubbles are as large as or larger than the rock pores), the foam behaves differently and is described in terms of lamellae or films moving through pore channels [8]. In soil remediation, foam can help homogenize fluid transport [9], particularly in heterogeneous soils with significant permeability contrasts [10]. By forming a fluid network within the pores, foam facilitates the delivery of reagents to contaminated soils [11]. In tunnel drilling through granular soils, foams are used in Earth Pressure Balance (EPB) shields, where they are mixed with the soil at the cutting edge to stabilize the soil wall by balancing earth pressure and providing capillary cohesion to the

<footnote>* vincent.langlois@univ-eiffel.fr</footnote>

<footnote>† olivier.pitois@univ-eiffel.fr</footnote>
soil [5]. Whether for reagent delivery in soil remediation or tunnel drilling, the permeability of foam-filled soils is a critical factor. For decades, following the pioneering work of Bernard [12], the permeability of foam-filled sands was considered analogous to that of unsaturated soils [7, 13]. However, recent experimental studies reveal that surfactant properties and the ratio of bubble size to grain size, $r = D_b/D_g$, play key roles [14, 15]. These studies indicate that the permeability $k$ of foam-filled granular packings, which relates the apparent liquid's flow velocity to the pressure gradient and the liquid's viscosity, can be described as:

$$k = \phi_p f(r) k_{\text{bulk}}, \qquad (1)$$

where $\phi_p$ is the volume fraction of the pores and $k_{\text{bulk}}$ the bulk foam permeability, i.e. that of the foam when it is not confined in the pore space. $k_{\text{bulk}}$ can be measured by dedicated methods [16]. The effect of the surfactant and that of the liquid fraction are incorporated into $k_{\text{bulk}}$, and $\phi_p f(r)$ accounts for the geometric effects induced by foam confinement. Based on certain assumptions regarding the microstructure of the confined foam (assimilated to a set of parallel straight liquid channels with no grain-induced tortuosity and no liquid channel at the grain surface), Wang et al. [17] developed an analytical model for the permeability of foam-filled soils, expressed in Eq. (1) (see Ref. [15]). However, their expression for the function $f$ does not align with the experimental measurements conducted by Ceccaldi et al. [15] on model systems (monodisperse spherical grain packings filled with monodisperse foam). This discrepancy may stem from the assumptions Wang et al. made concerning liquid distribution within the pores. Following Pitois et al. [18], the liquid filling the pores can be divided into three main components [Fig. 1(a)]: core foam, where bubbles interact with each other [Fig. 1(b)], wall Plateau borders (wall PBs) where bubbles interact with the grains, and the liquid bridges at grain contacts [Fig. 1(c)]. The structure of the core foam is similar to that of bulk foam [Fig. 1(b)]: three films meet to form bulk Plateau borders (bulk PBs), while four bulk PBs converge at the vertices. By making a few assumptions about the curvature of gas-liquid interfaces and the volume of liquid bridges, Pitois et al. [18] developed an analytical model describing how fluid is distributed among the three microstructural elements, and how the curvature of these interfaces relates to the total liquid volume fraction within the pores, $\phi_L$. Finally, in terms of permeability, it is essential to consider how all PBs are interconnected and the specific properties of the wall PBs, as their flow conductance is known to be higher than that of bulk PBs [19]. Note that Eq. 1 does not explicitly account
<figure><img src="image_2.png" /><figcaption>FIG. 1. Liquid distribution for a foam confined in a grain packing: (a) Scheme describing the liquid distribution between the so-called *core foam*, which is similar to bulk foam, and specific components resulting from the interaction between the foam and the grains: wall Plateau Borders and liquid bridges, (b) Detail on the bulk foam microstructure, (c) Magnified view of a large bubble confined in a pore formed between three glass beads and the wall (through which we observe the pore): observe the liquid bridges at contact between grains, connected to both wall Plateau borders and thin liquid films.</figcaption></figure>

for the contribution of wall PBs to permeability. Furthermore, since the pioneering work of Leonard and Lemlich [20–22], it has been well established that flow properties through PBs depend on surfactants, as they can impose varying interfacial mobilities. Some surfactants, such as bovine serum albumin (BSA), saponin, or tetradecyltrimethylammonium bromide (TTAB) with dodecanol, result in non-moving (‘rigid’) interfaces, whereas others, such as sodium dodecyl sulfate (SDS) or TTAB, lead to mobile interfaces [14, 22, 23]. Although Eq. (1) accounts for the effects of gas-liquid interfaces through $k_{\text{bulk}}$, this simple contribution remains questionable, as the effects of interface mobility are not exactly the same for wall and volume PBs [21].

In this work, our goal is to develop a numerical model of foam-filled granular media that can describe both (I) the liquid distribution within the pores, without relying on certain assumptions made in [18]), and (II) the liquid flow through the PB network. The first section gives an overview of the methodology used to achieve this objective, and provides detailed explanations for each step of the numerical modeling process. Specific results for each step are provided in the appendices. The second section presents the predictions of our model, first applied to bulk foam, and then to foam confined within a granular medium. We
compare our model’s predictions with previous experimental measurements, as well as with
new experiments conducted as part of this study.

## II. METHODS

### A. Numerical model

*Modeling guidelines*

In this work, we aim to describe by numerical calculations the effect of foam confinement within pores on liquid flow. To do this, we adopt a macroscopic description of the permeability with a homogenization method based on periodicity [24]. To limit the computing cost, we use a packing of only four grains arranged with a face-centered cubic (FCC) structure within a cube $L_t \times L_t \times L_t$ [Fig. 2(a)]. Such a grain packing contains 12 pores which we consider sufficient to obtain a representative elementary volume for studying the effect of foam confinement on permeability. Moreover, since we want to model a confined foam with a low ratio $R_b/R_g < 0.05$, the number of bubbles that fill the pore space could be of the order of $10^4$ leading to a few $10^5$ of PBs and around $10^9$ to $10^{10}$ surface elements to model all the liquid/gas interfaces. Therefore, as the number of elements required to mesh liquid-gas interfaces is expected to become huge and to limit the computing possibilities, we choose to generate "dry foam" microstructures, where PBs are simple edges without volume. In our numerical model, this dry microstructure is considered as the skeleton of the wet foam confined within the pores. Liquid is taken into account by assigning volume to PBs depending on their lengths.

The generation of the numerical "dry foam" microstructure, detailed in Appendix A 1, consists of three main steps. First, by merging monodisperse polyhedra of size $D_c$ from a Voronoi tessellation generated using the Neper software [25], we obtain an initial structure with rough grains [Fig. 2(a)]. This structure is then refined [Fig. 2(b)] using Surface Evolver [26] to determine the final network of "dry" Plateau borders [Fig. 2(c)]. This PB network consists of wall PBs covering the grain surface and volume PBs within the pore volume. To achieve our main objective, we need to assign volume and flow properties to each PB. In this work, we assume that the foam is at equilibrium, which corresponds to a homogeneous capillary pressure, i.e., a homogeneous mean radius of curvature $\delta$ of the
gas-liquid interfaces. Consequently, the volume and flow conductance of a PB depend on $\delta$, its length $L$, and its type (wall or volume PB). From geometries computed using Surface Evolver, the PB volumes, $V_v$ for bulk PBs [Fig. 3(a)] and $V_w$ for wall PBs [Fig. 3(b)], were directly determined, while their flow conductances were calculated using Comsol software. Results concerning the volume of PBs and their conductances are detailed in Appendices A 2 and B 2 respectively. Finally, the liquid permeability of the foam-filled granular packing is determined by modeling the network of PBs as a network of conductances and calculating the equivalent conductance of the entire system. In doing so, we implicitly assume that the foam microstructure is not modified by either steady flow or drainage. Note that, at all stages of our modeling, the structures are 3D-periodic.

Before going any further, we need to take a closer look at how to deal with liquid bridges, as our “dry foam” numerical microstructure does not include them. To assess their contribution to the liquid distribution within the pores, we have calculated their volumes $V_m$ for various values of $\delta$, $R_g$, and the number of bubbles surrounding the grain-grain contact, $n_{mb}$, using Surface Evolver simulations [Fig. 3(c)]. Then, to account for their contribution to liquid flow, we associate each of them with the core PB located in the area where the bridge forms and that connects the two grains involved. The conductances of these selected PBs were modified by various factors, ranging from 1 to 100. As the value of the factor did not significantly alter the permeability of the liquid, we concluded that bridges have no significant effect on permeability.

Finally, we should mention that some simplifications have been made regarding the microstructural elements shown in Figs. 3(b)-(c). In Fig. 3(b), since the PB coming from the core foam, hereafter referred as transitional core PB, is connected to three wall PBs, it is not perfectly identical to a real core PB [Fig. 3(b)]. Similarly, the transitional wall PBs connecting to the liquid bridge do not have the same geometry as the wall PBs. However, for simplicity, we have considered transitional core PBs as core PBs and transitional wall PBs as wall PBs. Since the volume of transitional PBs is expected to be larger than that of the original PBs, this assumption leads to an underestimation of the overall liquid volume (e.g. for the distinction between transitional core PBs and core PBs, see Fig. 13 in Appendix A 2).
<figure><img src="image_3.png" /><figcaption>FIG. 2. Main steps to generate "dry foam" microstructure confined in a porous media. By merging polyhedra contained in a Voronoi tessellation generated by Neper software, we generate an initial structure with rough grains (a), which is then enhanced by using Surface Evolver (b) to find the final network of "dry" Plateau borders (c). The PBs network is used to calculate the liquid permeability. The microstructure shown in this figure corresponds to $r = 0.13$.</figcaption></figure>

*Foam liquid repartition within the pore space*

By considering the mean radius of curvature $\delta$ fixed and using the lengths of the PBs provided by our numerical "dry foam" structure, we can calculate the volume of fluid located inside the PBs. Moreover, by adding the contribution of bridges, we can obtain the relationship between $\delta$ and the total liquid fraction filling the pores, $\phi_L$, as it was done in Ref. [18]:

$$ \phi_L = \frac{1}{(1-c)L_t^3} \left( \sum_{i, \text{core PB}} V_{v,i} + \alpha_w \sum_{i, \text{wall PB}} V_{w,i} + \alpha_w N_m V_m \right), \quad (2) $$

where $N_m$ is the number of bridges included in the entire medium (= 24 in this work), and $c$ is the volume fraction of the grains (= 0.63 in this work, see Appendix A 1). The coefficient $\alpha_w$ should be equal to 1 if the model succeeds in reproducing the liquid distribution observed experimentally. But, as mentioned earlier, by ignoring the transitional microstructural
<figure><img src="image_4.png" /><figcaption>FIG. 3. Local description of the gas-liquid and solid-liquid interfaces: (a) four half bulk Plateau borders joining to a bulk vertex ($\delta/L = 0.2$), (b) three wall Plateau borders and a transitional bulk Plateau border joining to a wall vertex ($\delta/L = 0.19$, $R_g/L = 3$), (c) a liquid bridge at the contact of two grains with $n_{mb}$ membranes distributed around the bridge ($\delta/R_g = 0.1$, $n_{mb} = 10$). Parts highlighted in transparent red are duplicated to obtain the full shapes. In (a), for low $\delta/L$ ratios, the cross-section in the middle of the PB can be approximated by three circles of radius $\sim \delta$. In (b), the PB coming from core foam is not entirely similar to a core PB. The angle between core PBs is equal to 109.5°, while the transitional core PB is orthogonal to wall PBs. Similarly, in (c), the PBs joining to the bridge are not identical to wall PBs. Transitional core or wall PBs are in-between microstructural elements.</figcaption></figure>

elements, our model underestimates the liquid volume. The coefficient $\alpha_w$ should ultimately be considered as an adjustment parameter to account for additional volumes of liquid in the vicinity of the grains, which would otherwise have been ignored in the model. Unless otherwise stated, we will use $\alpha_w = 1$. For wet foam, we define the mean bubble size $D_b$ from $D_c$, as $D_b = (1 - \phi_L)^{1/3} D_c$.

As our numerical model of dry foam does not include the bridges, we must specify how much membranes distribute around the bridge, $n_{mb}$, as a function of the ratio $r$ and liquid fraction $\phi_L$. To do this, we used an empirical relationship based on the observation of bridges formed at grain contacts with the wall of a column filled with glass beads and foam (see section II B):

$$n_{mb}^S \approx (1.46\phi_L + 0.89) (r^{-0.83} + 4). \quad (3)$$

In our numerical model, we considered this relationship to be valid for bridges between two grains ($n_{mb} = n_{mb}^S$).
As shown by Cox *et al.* [27] when studying flow in the junction of four PBs, it is possible to associate a flow conductance *G* with each PB. Flow conductance is defined as the ratio between the flow rate passing through the PB to the pressure difference applied at its ends [Fig. 4(a)]. Note that the PBs are taken from one vertex center to another, i.e. a fraction of the vertices is included within the PB. In this way, we assume that the conductance is a property of the entire PB, and that fluid flow respects the geometric symmetry of the PBs. In this work, we perform a systematic determination of the conductance of PBs. Calculations were carried out by solving the Stokes equations with the finite element software Comsol (Creeping flow, "P3+P2" discretization and a "Normal" meshing). Boundary conditions applied in our simulations are described in Figure 4. Concerning the gas-liquid interface Γ<sub>lg</sub>, two conditions were assumed to model either mobile interface or rigid interface (no slip condition). Since the work of Leonard and Lemlich [20], mobile interfaces have often been modeled using viscous interfacial behavior, leading to the introduction of an interfacial shear viscosity μ<sub>s</sub>. In this model, the continuity of stresses at the interface imposes that bulk viscous stresses and interfacial shear forces are equal. By invoking the rheological behavior of the fluid, this continuity can be expressed as a condition on the fluid velocity field **u**: μ(∂**u**/∂**n**) = μ<sub>s</sub>Δ<sub>s</sub>**u** where μ is the fluid shear viscosity, **n** is the unity vector normal to Γ<sub>lg</sub> and directed outwards from the liquid, and Δ<sub>s</sub> is the surface Laplacian [28]. This model leads to introduce a dimensionless number *M* = μδ/μ<sub>s</sub> describing the interfacial mobility. By comparing experimental results to models of permeability based on this interfacial rheology, Lorenceau *et al.* [29] showed that the interfacial shear viscosity is not well appropriate to describe their measurements. A better description of the flow of liquid through foam would therefore require including other physicochemical phenomena at various scales and their interactions [22]. In this work, we use a simple classical model to consider the slip condition at the gas-liquid interface Γ<sub>lg</sub> [30, 31], known as Navier slip: ∂**u**/∂**n** = −**u**/L<sub>s</sub> where L<sub>s</sub> is called the slip length. We show in Appendix B1 that, under certain conditions, the two models predict similar conductance values for straight channels. Results concerning the conductances, *G*<sub>v</sub> for bulk PBs and *G*<sub>w</sub> for wall PBs, are detailed in Appendix B2. It is shown that conductances depend on the mean radius of curvature δ and the PB lengths. Therefore, at fixed δ, we can build a network of conductances associated with each simulated
<figure><img src="image_5.png" /><figcaption>FIG. 4. Boundary conditions applied for calculations of conductance: (a) core PB ($\delta/L = 0.2$), (b) wall PB ($\delta/L = 0.19$, $R_g/L = 2.5$). Only the red-highlighted parts of the whole PBs shown in Fig. 3 are used. A pressure difference is imposed between each end of the PB. The velocity field is imposed to be normal to the end faces. Different conditions were considered for the gas-liquid interface $\Gamma_{lg}$: either a no slip condition (rigid interface) or Navier slip (mobile interface). In (b), a no slip condition is imposed on the solid-liquid interface $\Gamma_{ls}$.</figcaption></figure>

microstructure of confined dry foam.

The determination of permeability is based on a “pore network” approach [32]. A macroscopic pressure difference $\Delta P_t$, imposed at the ends of the whole medium, leads to a fluid flow through the PB network [Fig. 2(c)]. After solving the pore-network problem, the permeability $k$ can be calculated from the global flow rate $Q_f$ passing through a cross-section of the sample: $k = \mu Q_f / (L_t^2 \lambda)$ where $\lambda = \Delta P_t / L_t$ is the macroscopic pressure gradient. Alternatively, the permeability can be determined from the pressure drops $\Delta p_i$ through each PB:

$$k = \frac{\mu}{(\Delta P_t)^2 L_t} \sum_{i, \text{PB}} G_i (\Delta p_i)^2 \quad (4)$$

where $\sum_{i, \text{PB}}$ is a sum over all PBs and $G_i$ is the conductance chosen according to the nature of the PB (wall or bulk). With this expression, the contributions of each type of PB can be determined by dividing the sum over PBs into two parts, one for each type of PB. Moreover, for each numerical microstructure, we calculated the permeability for each direction $(x, y$ and $z)$ and found very similar results (which is expected since no direction is favored in our model).

In this work, we calculate the permeability for two main hypotheses concerning the mean curvature radius $\delta$: (i) constant liquid volume fraction $\phi_L$ where $\delta$ is determined by inverting
<figure><img src="image_6.png" /><figcaption>FIG. 5. (a): Setup used to generate foam with controlled bubble sizes and liquid fractions, and to study liquid distribution in foam-filled bead packings. (b) Images of the liquid at the wall column captured through an optical prism for two size ratios $r = R_b/R_g$ (rows) and two liquid fractions $\phi_L$ (columns). We highlight the wall liquid bridges in green circles and the apparent core foam with transparent blue area. In the area of the core foam, we measured the surface liquid fraction $\phi_{L,core}^S$.</figcaption></figure>

Eq. (2), (ii) constant $\delta$ (equal to the value obtained for the bulk foam).

## B. Experiments

We used the same setup as in [14] to produce monodisperse foam of diameter $D_b \approx 425\mu m \pm 10\mu m$ [Fig. 5(a)]: a co-injection of the foaming solution at flow rate $Q_l = 10$ mL/min, with nitrogen gas at $Q_g = 10$ mL/min, into 2 inlets of a T-junction. The foaming solution was prepared with the nonionic surfactant Glucopon ® 225DK (BASF) at 20 g/L in water. The produced bubbles are accumulated in a vertical column filled with the same foaming solution. A syringe pump is connected to the column to withdraw the excess of the foaming solution at $Q_l$, so that the height of the foam inside the column is constant during the experiments [Fig. 5(a)]. Consequently, the foam is flowing out of the column at a
constant flow rate of $Q_g = 10$ mL/min and with a controlled liquid fraction $\phi_L = 8\%$. Then, foam can be directly transported to the measurement cell, or can pass by a mixing device where the foaming solution is injected to increase the liquid fraction to a targeted liquid fraction $\phi_L$ [Fig. 5(a)]. The measurement cell is a cuboid column of square cross section (27 x 27 x 50 mm), filled with monodisperse plastic beads of size $D_g$ (from 2 mm to 14 mm). The porosity of the bead pack is ranging between 0.37 and 0.47 (the largest value was measured for the largest bead size). To counteract the effect of gravity, the measurement cell is rotated horizontally during the filling step at a constant speed of 120°/s. Foam is injected into the cell until the outlet foam is similar to the inlet foam in terms of bubble size and liquid fraction (a sample of this foam is taken to ensure that the target liquid fraction is reached). After foam injection, pressure was released from both ends of the cell to limit the effect of gas compressibility. The rotation is then stopped and a picture is taken by a camera positioned on the upper side of the column through a right-angled prism (30 x 30 x 21.2 mm) fixed on the column in direct contact with foam and grains. The prism allows us to obtain images with sufficient contrast between the liquid phase and the gas phase [33] as shown in Fig. 5(b). From such an image, we determine the mean number of bubbles $n_{mb}$ located around a bridge, and the liquid fraction of the core foam $\phi_{L,core}$ by measuring the surface liquid fraction $\phi_{L,core}^S$, i.e. the proportion of liquid phase in the part of the image corresponding to the core foam. Indeed, as shown by Höhler et al. [34], the liquid fraction of the bulk foam is linked to the surface liquid fraction: $\phi_L^S = 1 - \tilde{\Pi}/(2 + \tilde{\Pi}/(1 - \phi_L))$. The ratio of the osmotic pressure to the capillary pressure of the bulk foam, $\tilde{\Pi}(\phi_L)$, is given by: $\tilde{\Pi}(\phi_L) = k(\phi_L^* - \phi_L)^2/\phi_L^{1/2}$ with $k \approx 7.3$ and $\phi_L^* \approx 0.26$ [35].

Note that, without any calculations, Fig. 5(b) exemplifies the effect of the confinement on the liquid fraction of the core foam. For the same total liquid fraction, foams confined between small grains appear drier than those confined between larger grains (compare the PB thickness).

### III. RESULTS AND DISCUSSION

In this section, we present the results of our numerical simulations supplemented by the analytical model of a single bubble per pore (Appendix C). We first apply the method to bulk foams, then detail the results obtained for confined foams.
<figure><img src="image_7.png" /><figcaption>FIG. 6. Bulk foam: (a) $\delta/D_b$ and (b) $k_{\text{bulk}}/D_b^2$ as a function of $\phi_L$ (with $L_s = 0.022D_b$ for mobile interfaces). The symbols (circles or triangles) represent the results of our numerical model. In (a), the dashed line corresponds to the relationship obtained for Kelvin foam [29]. In (b), the dashed lines '--' are plotted by using the approximate expressions provided by [23]. Stars are experimental data from [14]. The dashdotted '·--' and dotted '··' lines correspond to the results of our numerical model using the analogy between Navier slip and interfacial viscosity, i.e. $L_s = M\delta/9$, with $M = 1$ for the dashdotted line ('M cst'), or $M = 5\delta/D_b$ for the dotted line ('$\mu_s$ cst', $\mu_s = \mu D_b/5$).</figcaption></figure>

A. Bulk foam

The methods described in the section II can be applied in a simpler way to bulk foam. In this case, foam microstructures, containing only bulk PBs, are simply generated by surface minimization of Voronoi tessellations. Figure 6(a) shows that the relationship predicted by our numerical model using Eq. (2) [with $N_m = 0 = c$] between the ratio $\delta/D_b$ and the liquid fraction $\phi_L$ is in good agreement with the results obtained from Surface Evolver simulations on Kelvin foam (see [29]).

Regarding the permeability, Fig. 6(b) shows that our model also predicts quite well the values of literature as a function of the liquid fraction in the case of rigid interfaces. Indeed, our results are in good agreement with the semi-empirical relationship established by Rouyer *et al.* [23] and the experimental results obtained by Ceccaldi *et al.* [14]. It is worth recalling that with rigid interfaces, our model has no fitting parameter. The slight overestimation of the permeability for $\phi_L > 0.1$ is expected to be an effect of the limit of our "dry foam" model. For the case of mobile interfaces, the value of the slip length $L_s$ must be fixed before calculating the permeability. Using an approach similar to that used in
Ref. [29], we consider the slip $L_s$ as a macroscopic fitting parameter and fix it as a constant fraction of the bubble size. Using $L_s = 0.022D_b$, we obtain a pretty good estimate of the permeability compared to the experimental values [14, 23]. By applying the results obtained from straight channels (see Appendix B 1) to the bulk PBs, we were able to test the two hypotheses commonly found in the literature to describe the foam permeability [29, 36]: constant interfacial mobility $M$ or constant surface viscosity $\mu_s$. None of them achieves as satisfactory an agreement as that found with the length $L_s \propto D_b$. For this reason, we will proceed with this model to study the effect of confinement.

## B. Confined foam

We detail now the main results of this paper concerning foam confined in a granular medium. Figures 7 and 8 show the effects of foam confinement on mean radius of curvature $\delta$ and liquid distribution at fixed total liquid volume fraction $\phi_L$. The main effects, previously described in Ref. [18], are observed:

1. the mean curvature radius $\delta$ decreases as the size ratio $r$ increases, whatever the liquid fraction [Fig. 7];

2. foam confinement redistributes the liquid from the core foam to the bridges and the wall PBs, and for high confinement, from wall PBs to bridges [Fig. 8];

3. the proportions of liquid between core foam, wall PBs and bridges are mainly dependent on $r$ and only slightly on $\phi_L$ [Fig. 8];

Compared to the analytical model by Pitois et al. [18], our key improvements stem from the explicit consideration of the mean radius of curvature $\delta$ and the effect of membranes $n_{mb}$ on the volume of liquid bridges (see in Appendix A 2, how the bridge volume increases with the number of membranes). Thanks to the better description of the microstructure, the numerical model enables a more accurate prediction of the measured liquid fraction of the core foam [Fig. 8(a)]. The numerical model predicts that liquid bridges have a larger volume of liquid compared to the analytical model. However, our model does not describe the liquid fraction of the core foam as accurately as expected, based on our measurements. To better match the experimental data, it is necessary to increase the value of $\alpha_w$ in equation
<figure><img src="image_8.png" /></figure>

FIG. 7. Confined foam: mean curvature radius $\delta$ normalized by $D_b$ as a function of the size ratio $r$. Triangles represent the results obtained under the "constant liquid fraction" assumption, with $\phi_L = 0.05$, 0.1 or 0.15. Squares correspond to the results for foams with $\phi_L = 0.1$ and an additional volume of liquid located in the wall PB and the bridges [$\alpha_w = 1.6$ in Eq. (2)]. All data for $r \approx 0.55$ correspond to the results of the analytical model (Appendix C). The dashed lines correspond to the predictions of the analytical model of Pitois et al. [18].

2, i.e. to increase the volumes of wall PBs and bridges [$\alpha_w = 1.6$ in Fig. 8(a)]. It should be noted that such an increase in bridge volume corresponds to an increase in bridge size of around 12%. Nevertheless, the most notable effect of confinement is the drop in the liquid fraction of the core foam, associated with that of the ratio $\delta/D_b$. Therefore, as the flow conductances evolve roughly as $\delta^4/(\mu L)$ [see Appendix B 2], the liquid redistribution should have a significant impact on permeability.

Figure 9 shows that, at $\phi_L$ fixed, the permeability $k$ of the confined foam predicted by our numerical model decreases as the size ratio $r$ increases. Moreover, after normalized by the permeability of the bulk foam $k_{\text{bulk}}$, we observe that the permeability ratio $k/k_{\text{bulk}}$ is only slightly dependent on the liquid fraction $\phi_L$ and the interfacial mobility. These general observations are in perfect agreement with the experimental results of Ceccaldi et al. [14]. It should be emphasized that the decrease of the permeability with size ratio is a direct effect of the liquid redistribution at fixed $\phi_L$. Note that if, by assumption, we choose to keep $\delta$ constant, i.e. the liquid redistribution is now neglected, we arrive this time at an increase in the dimensionless permeability, which is completely contrary to the experimental results, as shown in Fig. 9(a). Finally, note that the analytical model of Wang et al. [17] does not reflect the trend observed with the experiments and our numerical simulations.
<figure><img src="image_9.png" /><figcaption>FIG. 8. Confined foam: liquid distribution between core foam (a), wall PBs (b) and bridges (c), as a function of the size ratio $r$. Triangles represent the results obtained under the "constant liquid fraction" assumption, with $\phi_L = 0.05$, 0.1 or 0.15. Squares correspond to the results for foams with $\phi_L = 0.1$ and an additional volume of liquid located in the wall PB and the bridges [$\alpha_w = 1.6$ in Eq. (2)]. All data for $r \approx 0.55$ correspond to the results of the analytical model (Appendix C). The dashed lines correspond to the predictions of the analytical model of Pitois et al. [18]. In (a), the red stars correspond to our experimental measurements. For each value of the size ratio $r$, we calculate the average value obtained with experiments having different $\phi_L = 8\%$, 13\%, 22\% and 28\%. Error bars correspond to minimal and maximal values. No systematic order was observed for the different $\phi_L$ values.</figcaption></figure>

The overall behavior described above shows that the ratio $k/(\phi_p k_{\text{bulk}})$ can be approximately expressed in terms of the parameter $r$ only, as proposed in Refs. [14, 15]: $k/(\phi_p k_{\text{bulk}}) = f(r)$. Although useful to describe the experimental data, this function $f$ conceals the contribution of the wall PBs. Numerical simulations show that the liquid flows at the surface of the grain are the main contributors to the whole permeability for $r > 0.06$ [Fig. 10]. As the ratio between the number of wall PBs to the number of core PBs decreases as the size ratio decreases, the wall PB contribution should become negligible for low size ratios. But, as the wall PBs are wider than bulk PBs at same $\delta$ [Fig. 15], this effect is only effective for very low size ratios. Moreover, by comparing the results for rigid interfaces [Fig. 10(a)] with those for mobile interfaces [Fig. 10(b)], it appears that the contribution of the core foam is only weakly dependent on the interfacial mobility, whereas the contribution of wall PBs is greater for non-mobile interfaces than mobile interfaces. Therefore, since the core foam and the bulk foam are made of similar PBs, it seems that Eq. 1 [or $f(r)$] can be firmly justified only in the case
<figure><img src="image_10.png" /><figcaption>FIG. 9. Confined foam: permeability as a function of the size ratio $r$ for (a) rigid interfaces and (b) mobile interfaces ($L_s = 0.022D_b$). Symbols '+' represent the experimental data from Ceccaldi et al. [14]. The shaded areas correspond to the error bars associated to these experimental results [including the systematic error due to $\phi_L$ effects obtained using Eq. (1) to analyze their experiments performed in the range 4% to 20% for $\phi_L$]. Triangles represent the results obtained under the "constant liquid fraction" assumption, with $\phi_L = 0.05$, 0.1 or 0.15. The circles correspond to the results obtained under the "constant $\delta$" hypothesis (where $\delta \approx 0.207D_b$, the value corresponding to a bulk foam with $\phi_L = 0.1$). Squares correspond to the results for foams with $\phi_L = 0.1$ and an additional volume of liquid located in the wall PB and the bridges [$\alpha_w = 1.6$ in Eq. (2)]. The black diamond at $r=0$ is the theoretical result expected for non conductive spheres embedded in a conductive medium [37]. The symbols for $r \approx 0.55$ correspond to the results of the analytical model (Appendix C). Dashed line corresponds to the prediction of the analytical model of Wang et al. [17].</figcaption></figure>

of low size ratios $r \le 0.02$. Additionally, it has been suggested that the function $f$ tends to one as the size ratio approaches zero. Our results clarify the expected behavior at very small values of $r$: the dimensionless permeability approaches a value of 0.74, which accounts for the tortuosity imposed by the presence of solid spheres. This limit can be understood by analogy with the effective conductivity of two-phase media. From [37], the effective conductivity $\sigma_e$ of a FCC packing of non conducting spheres embedded in a conducting matrix of conductivity $\sigma_m$ is given by: $\sigma_e/\sigma_m = (1+2(1-\phi_p)\beta - 2\phi_p\xi\beta^2)/(1-(1-\phi_p)\beta - 2\phi_p\xi\beta^2)$ with $\xi \approx 0.0824$ and $\beta = -0.5$. In our case, this expression leads to $\sigma_e/(\phi_p\sigma_m) \approx 0.74$. This value is more in agreement with the trend observed at low size ratios for both $k/(\phi_p k_{\text{bulk}})$ [Fig. 9] and the contribution of the core foam [Fig. 10].
<figure><img src="image_11.png" /><figcaption>FIG. 10. Details of the contributions of core PBs and wall PBs to total permeability under the "constant liquid fraction" assumption with $\phi_L = 0.1$ and (a) rigid interfaces or (b) mobile interfaces ($L_s = 0.022D_b$).</figcaption></figure>

Finally, although the general description is consistent with experimental observations, our model overestimates the permeability measured experimentally. Similarly to what we observed for liquid distribution, our predictions can be improved by considering larger wall PBs and liquid bridges [with $\alpha_w = 1.6$, see Fig. 9].

IV. CONCLUSION

We have developed a numerical model of foam-filled granular media aimed at describing both liquid distribution in the pores and liquid permeability. This model is constructed through several steps. Starting with a Voronoi tessellation generated by Neper software, we create an initial foam structure, which is subsequently modified using Surface Evolver software through a series of ad-hoc operations. The final structure comprises Plateau borders (PBs) located in the core of the pores and at the surface of the grains (wall PBs). To determine the fluid distribution within the pores, we enhance this "dry foam" model by assigning a volume and incorporating the contributions of liquid bridges at the grain contacts. For liquid flow determination, we assign a flow conductance to each Plateau border based on its length, the mean curvature radius, and grain size for wall PBs. Local properties (volume and conductance) are systematically determined through numerical simulations using Surface Evolver to establish the shape of the gas-liquid interface for each microstructural element (bulk PB, wall PB, and bridge), and Comsol calculations to
evaluate PB conductance. A Navier slip condition is used for PB conductances to account for interfacial mobility. The permeability is ultimately calculated by determining the global conductance of the PB network. This method has been applied to study both bulk foam and foam confined within a packing of spheres.

For bulk foam, the numerical model accurately predicts the evolution of the microstructure and permeability $k_{\text{bulk}}$ as the liquid volume fraction $\phi_L$ increases. We demonstrate that the inclusion of Navier slip provides a robust description of the relationship between $k_{\text{bulk}}$ and $\phi_L$ for foam with mobile interfaces. In the case of foam confined within granular media, we enhance the description of liquid distribution provided by Pitois et al. [18], and achieve a strong agreement between experimental permeability measurements and our calculations. We find that the observed decrease in permeability as a function of the size ratio at fixed $\phi_L$, observed experimentally by Ceccaldi et al. [14], results from the redistribution of liquid from core foam to bridges and wall PBs. This redistribution leads to a reduction in the mean curvature radius and the conductances. Furthermore, the model allows us to assess the respective contributions of bulk and wall PBs to liquid flow. We observe that wall PBs significantly contribute even at low size ratios $R_b/R_g \approx 0.04$. However, our model tends to overestimate both the liquid volume within the core foam and the permeability. The reasons for this discrepancy are not fully understood. But, we show that a corrective coefficient applied to increase the volumes of wall PBs and liquid bridges allows us to improve the predictions of core foam liquid fraction and permeability. Advanced imaging techniques such as X-ray microtomography could be useful to assess this assumption.

## ACKNOWLEDGMENTS

The authors gratefully acknowledge funding from the French National Research Agency (ANR 23-CE51-0030 BONDINGFOAM).

## Appendix A: Details on numerical microstructures

In this appendix, we detail the methods used to generate "dry foam" microstructures and the geometry of the "wet" microstructural elements (volume and wall PBs, liquid bridge), and we give some results concerning the volume of the microstructural elements.
# 1. Generation of "dry foam" structures

To produce numerical samples of confined "dry foam", we start by generating Voronoi tessellations with Neper software [25] and merge polyhedron located inside the expected grains. The original structure of foam is generated by Neper with small standard deviations for cell sizes and sphericities ($R_c/\bar{R}_c = 1 \pm 0.085$ for cell sizes, and $0.9 \pm 0.015$ for sphericities). Cell size is defined as the radius of the sphere having the same cell volume, $R_c$. Note that, to avoid complications due to pinched bubbles near sphere contacts, the final radius of the spherical grains, $R_g$, is chosen to prevent neighboring grains from touching in the final state. The final size of the spheres leads to a solid fraction $c$ (sphere volume fraction) equal to 0.63 similar to the experimental solid fraction of monodisperse grains. Structures based on Neper software are then modified by a series of ad-hoc operations in Surface Evolver software:

A. Set the surface tension of the surfaces covering the upcoming grains to 1, then apply surface minimization steps (gradient descent method and conjugate gradient method).

B. Calculate grain volume to determine grain expansion (or contraction) and contraction of pore-filling bubbles (or expansion).

C. In 20 steps, modify bubble and grain volumes up to their desired volumes by modifying target volumes and performing minimization steps.

D. In 10 steps, increase grain surface tension up to 10 with minimization steps in order to achieve quasi-spherical grains.

E. Gradually move the vertices of the grains to the desired positions with minimization steps;

F. Set the surface tension of grain to 0.5, and impose that grain vertices are located on spheres, then proceed with minimization steps.

Figure 2b shows an example of a final structure generated by the whole process, and the corresponding network of PBs is shown in Fig. 2c. Ten microstructures were generated for various ratios $r_c = D_c/D_g$, where $D_c$ is the mean diameter of the "dry" bubbles and $D_g$ is the diameter of the grains. To assess the suitability of the method, we calculated various geometrical parameters such as (i) the ratio between the number of "dry" bubbles $N_b$ and the
<figure><img src="image_12.png" /><figcaption>FIG. 12. Volume of PB, $V/(\delta^2L)$ as a function of $\delta/L$: (a) bulk PB, (b) wall PB. The symbols correspond to the Surface Evolver results. The diamonds in (b) highlight the values of $\delta$ and $V_w$ at which the films close (c). The dashed lines are plotted by using Eqs. A1 & A2.</figcaption></figure>

wall PBs for bridges (see Supplemental Material [38]). Results of this parametric study are detailed hereafter.

The gas-liquid interfaces of the various microstructural elements [Fig. 3] are determined by Surface Evolver. Starting from the initial raw geometries, the mesh is progressively refined, taking care to maintain a homogeneous mesh (equiangulation, removal of small triangles and edges). In the final mesh, the maximum edge length is around $\delta/9$. The calculation is stopped when a sequence of commands does not modify the surface of the liquid-gas interfaces by more than 0.001%. The sequence of commands involves the repetition 25 times of: a few iterations of gradient descent followed by the deletion of small edges and the refinement of longer ones to respect the maximum edge length.

Figure 12a shows the volume of a bulk PB as a function of $\delta/L$. Examples of PB microstructures obtained by Surface Evolver simulations are shown in Fig. 13. For low $\delta/L$ ratios, the volume could be approximate by a relationship of the type $A\delta^2L + B\delta^3$, where the term $\delta^2L$ would represent the contribution of the almost straight part of the Plateau border, and the term $\delta^3$ the one of the vertex. When $\delta/L$ increases, the PB widens and the area of the films forming the PBs decreases. For $\delta/L = \delta_w^*/L = 6^{0.5}/4$, the films close and the junction of the four PBs is formed by the packing of four spherical bubbles. For wall PBs, the volume depends on both $\delta/L$ and $R_g/L$ [Fig. 12(b)]. As for bulk PBs, the films are closing as $\delta/L$ increases. The PB volume and the value of $\delta/L$ for which the film close, $\delta_w^*/L$, depends on $R_g/L$ [Fig. 12(c)].
<figure><img src="image_13.png" /><figcaption>FIG. 13. Example of PB microstructures for various ratios $\delta/L$: (a) Bulk PBs, (b) Wall PBs ($R_g/L = 3$).</figcaption></figure>

Figure 14 shows the volume of a bridge as a function of $\delta/R_g$ for different number of surrounding membranes $n_{mb}$ [Fig. 3(c)]. Our numerical results show that, thanks to the surface tension in the membranes, the volume of bridge surrounded by membranes is higher than that of a capillary bridge without membrane. At a fixed ratio $\delta/R_g$, the more membranes there are, the greater the volume of the bridge. To this regard, note that the approximate relationship used by [18] to model the liquid distribution of a foam confined in a granular packing underestimates the volume of bridges.

Note that the comparison of the microstructure shown in Figs. 13(a)-13(b) for $\delta/L = 0.8$ clearly shows that the "core-like" PB connected to the three wall PBs is larger than a bulk PB. A more complete model should require to simulate transitional PBs linking the core foam to wall PBs, but also, linking the wall PB to the bridges. To avoid adding further complexity to our model, we have chosen not to take these transient microstructure elements into account. In our model, the transitional core PBs are considered as core PBs, and transitional wall PBs as wall PBs.

Below, we gather different approximate formula chosen to fit the numerical results. For convenience, we use the following notation: $x = \delta/L$, $x^* = \delta^*/L$, $\Delta x = x^* - x$ and $y = L/R_g$.

• for bulk PBs:

$$ \left\{ \begin{array}{ll} V_v = (3^{0.5} - \frac{\pi}{2}) (1 + 1.511x + 1.067x^3) x^2 L^3, & \text{if } x < 0.35, \\ V_v = \frac{V_v^*}{226.9\Delta x^3 - 39.41\Delta x^2 + 14.34\Delta x + 1}, & \text{if } 0.35 \le x \le x^*, \end{array} \right. \quad (A1) $$
<figure><img src="image_14.png" /><figcaption>FIG. 14. Volume of bridge $V_m/R_g^3$ as a function of $\delta/R_g$ for various number of membranes $n_{mb}$ surrounding the grain-grain contact. The symbols correspond to the Surface Evolver results. The continuous line is plotted by using a relationship established by Lian and Seville [39] for a capillary bridge without membrane. The dashdotted line '---' is determined from the approximate expression used in [18] assuming a cylindrical shape for the bridge. The dashed lines '---' are plotted by using Eq. (A3).</figcaption></figure>

with $V_v^* = (\frac{3}{2})^{0.5} (\frac{1}{2^{0.5}} + \pi - 3 \arccos(\frac{1}{3})) L^3$ and $x^* = \frac{\delta_v^*}{L} = \frac{6^{0.5}}{4}$.

• for wall PBs:

$$ \left\{ \begin{array}{ll} V_w = \frac{(2-\frac{\pi}{2})x^2}{1+(2.614y-1.616)x+(-1.364y^2-1.947y+2.276)x^3} L^3, & \text{if } x < 0.6x^*, \\ V_w = \frac{V_v^*}{1+(0.5603y^2-1.873y+5.447)(\frac{\Delta x}{x^*})+(21.32y^2-38.49y+44.39)(\frac{\Delta x}{x^*})^4} & \text{if } 0.6x^* \le x \le x^*, \end{array} \right. \quad (A2) $$

with $\frac{V_w^*}{L^3} = 0.01536y^3+0.01386y^2+0.02293y+0.1334$ and $x^* = \frac{\delta_w^*}{L} = 0.0678y^3+0.0439y^2+0.2205y+0.3677$.

• for bridges:

$$ V_m = \left[ \left( \frac{\delta}{R_g} \right)^{0.72} (0.0674n_{mb}^2 + 0.2436n_{mb}) + 1 \right] V_{m0}, \quad (A3) $$

where $V_{m0}$ is the volume of a capillary bridge between two spheres. From [39], we extract this approximate formula: $V_{m0} = \exp(-5.94 \cdot 10^{-3}z^4 - 0.1399z^3 - 1.037z^2 + 1.513z - 1.496)R_g^3$, with $z = \log(\frac{\delta}{R_g})$ and log the logarithm (base 10).
In our numerical model of "dry foam", some PBs should have small lengths such that $\delta_{v,w}^*(L) > \delta$. In a wet foam, these particular PBs should be integrated to neighboring vertices. Therefore, for these PBs and to be able to calculate the foam properties in our numerical model, we use $V_{v,w} = V_{v,w}^*$ when $\delta_{v,w}^*(L) > \delta$.

## Appendix B: Fluid properties at the scale of Plateau borders

In this appendix, we detail some results concerning fluid flow through straight Plateau borders and through the volume and wall Plateau borders studied in the present work (see Fig. 3 and Appendix A 2).

### 1. Straight Plateau borders

In this part, we consider the fluid flow through a straight channel of constant cross-section [Fig. 15]. Considering a slow (low Reynolds number) and incompressible flow, the steady-state is described by the Stokes equation. The fluid flow results from a driving force $\Lambda = -\lambda + \rho g \cos \theta$ where $\lambda$ is the pressure gradient and $\theta$ is the angle between the axis of the channel and the vertical axis. As detailed in [21] after [20], the equations followed by the fluid velocity field $u$ passing through straight channels are:

$$ \tilde{\nabla}^2 \tilde{u} = -1 \text{ over } \Omega_l, \qquad (B1) $$

$$ \tilde{u} = 0 \text{ over } \Gamma_{ls}, \qquad (B2) $$

$$ \tilde{u}_s(\tilde{s}) = 0 \text{ for } \tilde{s} = 0 \text{ and } \tilde{s} = \tilde{s}_{\max}, \qquad (B3) $$

$$ \frac{\partial \tilde{u}}{\partial \tilde{n}} = \frac{1}{M} \frac{\partial^2 \tilde{u}_s}{\partial \tilde{s}^2} \text{ over } \Gamma_{lg}, \qquad (B4) $$

where $\tilde{u} = \mu u/(\delta^2 \Lambda)$ and $M = \mu \delta/\mu_s$ with $\mu_s$ surface shear viscosity. Note that the lengths are normalized by $\delta$ ($\tilde{l} = l/\delta$), and $\tilde{s}_{\max} = \pi/3$ for bulk channels and $\pi/2$ for wall channels.

Equation B1 corresponds to the Stokes equation and Eq. (B2) to the no slip condition at the solid-liquid interface [involved in wall Plateau border, Fig. 15(b)]. Equation B3 was introduced by Leonard and Lemlich [20] to take into account the effects of surfactant counter-flow. Finally, Eq. (B4) derives from equilibrium of stresses at the liquid-gas interface. The term on the left-hand side of the equation shows the contribution of shear stresses from the bulk fluid, and the right-hand term corresponds to the contribution of forces per unit
<figure><img src="image_15.png" /><figcaption>FIG. 15. The cross-section of (a) a bulk channel and (b) a wall channel. The center lines show the symmetry. Purple lines correspond to liquid-gas interface $\Gamma_{lg}$, and dashed green line to liquid-solid interface $\Gamma_{ls}$.</figcaption></figure>

length acting on the interface. Leonard and Lemlich [20] considered an interfacial rheology governed by a surface shear viscosity $\mu_s$. In this work, we choose to use a Navier slip to model the interface mobility. With this simple model, Eq. (B4) must be replaced by:

$$ \frac{\partial \tilde{u}}{\partial \tilde{n}} = - \frac{\delta}{L_s} \tilde{u}_s \quad \text{over } \Gamma_{lg}, \qquad (B5) $$

where $L_s$ is the slip length.

Both problems can be solved numerically by various methods [19, 40, 41]. As Drenckhan et al. [41], we use the Comsol software to solve these equations (Coefficient Form PDE plus a weak contribution to take into account the effects of the interfacial forces on $\Gamma_{lg}$). By integrating the solution $\tilde{u}$ over $\Omega_l$, we calculate the flow conductance. Fig. 16 shows the evolution of flow rates as a function of $M$. For the Leonard and Lemlich model, we find the same evolution as previously published [19, 40]. Considering $L_s/\delta = M/9$ for bulk channels and $M/4$ for wall channels (justification follows), we show that the two models are almost equivalent over a wide range of interface mobility degrees, which gives us the opportunity to study the effects of different assumptions on $M$ through equivalent values for the slip length. For wall PBs, both solutions share the same asymptotic solutions at low and high interfacial mobilities, while a slight deviation is observed only in the transition regime.

As shown previously, both problems differ only on the way of modeling the slip condition. And, interestingly, the solution of each problem can be rewritten as a linear combination of
<figure><img src="image_16.png" /><figcaption>FIG. 16. Normalized flow rates through a straight channel considering either a Navier slip (circles) or the Leonard and Lemlich model (squares). The lines correspond to the solutions obtained by using the matrix approach [Eq. (B12) with $n_{max} = 200$]. Note in the case of the Navier slip, we consider $M = M_N = 9L_s/\delta$ for bulk channels and $4L_s/\delta$ for wall channels.</figcaption></figure>

functions $u_j^*$, solutions of the following problem:

$$ \tilde{\nabla}^2 u_j^* = -g(j) \quad \text{over } \Omega_l, \qquad (B6) $$

$$ u_j^* = 0 \quad \text{over } \Gamma_{ls}, \qquad (B7) $$

$$ u_j^*(\tilde{s}) = 0 \quad \text{for } \tilde{s} = 0 \text{ and } \tilde{s}_{\text{max}}, \qquad (B8) $$

$$ u_j^*(\tilde{s}) = \sin\left(j\pi\frac{\tilde{s}}{\tilde{s}_{\text{max}}}\right) \quad \text{over } \Gamma_{lg}, \qquad (B9) $$

where $g(0) = 1$ and $g(j) = 0$ for $j > 0$.

Few comments can be made: (i) The solution $u_0^*$ corresponds to the rigid interface solution obtained either when $M \to 0$ with the Leonard and Lemlich model or when $L_s \to 0$ with the Navier slip. (ii) $u_0^*$ and $u_1^*$ are positive everywhere over $\Gamma_{lg}$ on contrary to $u_j^*$ for $j > 1$ due to the boundary condition. Therefore, $u_0^*$ and $u_1^*$ are important contributions to the whole flow rate. (iii) The combination of $u_j^*$ must coincide with the Fourier series of $\tilde{u}(s) = \sum_{j=1}^{\infty} A_j \sin(j\pi\frac{\tilde{s}}{\tilde{s}_{\text{max}}})$ where $A_j$ are the Fourier coefficients [and $A_0 = 1$ to respect Eq. (B1)]. (iv) For bulk channels [Fig. 15(a)], the symmetry of the problem requires that the coefficients $A_j$ are zero for even $j > 0$.

By assuming $\tilde{u} = \sum A_j u_j^*$ and using Eq. (B4), we can show that the coefficients $A_j$ can be deduced by solving a linear problem:

$$ \frac{\partial u_{0,s}^*}{\partial \tilde{n}} = - \sum_{j>0}^{\infty} \left[ \frac{1}{M} \left( \frac{j\pi}{\tilde{s}_{\text{max}}} \right)^2 \sin \left( \frac{j\pi \tilde{s}}{\tilde{s}_{\text{max}}} \right) + \frac{\partial u_{j,s}^*}{\partial \tilde{n}} \right] A_j. \qquad (B10) $$
With the Navier slip, the same approach leads to:

$$ \frac{\partial u_{0,s}^*}{\partial \tilde{n}} = - \sum_{j>0}^{\infty} \left[ \frac{\delta}{L_s} \sin \left( \frac{j\pi \tilde{s}}{\tilde{s}_{\max}} \right) + \frac{\partial u_{j,s}^*}{\partial \tilde{n}} \right] A_j. \quad (B11) $$

Note that for low interface mobility (small $M$ or $L_s$), the two previous equations show that the coefficients $A_j$ should be proportional to $M$ or $L_s/\delta$. Consequently, the deviation to the rigid interface for the conductance is linearly dependent on $M$ or $L_s/\delta$ in the low interface mobility regime.

Now, if we limit the infinite sum to the terms ($j \le n_{\max}$) and discretize the boundary $\Gamma_{lg}$ with a number of vertices $n_v \ge n_{\max}$ of abscissa $s_i \in [4s_{\max}/(2n_{\max} - 1), s_{\max}/2]$, the linear problem can take a matrix form:

$$ \mathbf{K} \mathbf{A} = -\mathbf{S}, \quad (B12) $$

where $\mathbf{A}$ is the vector containing the coefficients $A_j$ (least-squares solution to the system of equations), $\mathbf{S}$ is the vector containing the coefficient $\frac{\partial u_{0,s}^*}{\partial \tilde{n}}(s_i)$ and $\mathbf{K}$ is the matrix defined as $[K_{LL}]_{i,j} = \frac{1}{M} \left( \frac{j\pi}{\tilde{s}_{\max}} \right)^2 \sin \left( \frac{j\pi s_i}{s_{\max}} \right) + \frac{\partial u_{j,s}^*}{\partial \tilde{n}}(s_i)$ for the Leonard and Lemlich model and $[K_N]_{i,j} = \frac{\delta}{L_s} \sin \left( \frac{j\pi s_i}{s_{\max}} \right) + \frac{\partial u_{j,s}^*}{\partial \tilde{n}}(s_i)$ for the Navier slip model.

Consider first the case of a bulk channel for which $A_{2k} = 0$. By assuming $M = M_N = \left(\frac{\pi}{\tilde{s}_{\max}}\right)^2 \frac{L_s}{\delta}$, we can write $\mathbf{K}_{LL} = \mathbf{K}_N + \mathbf{K}_c$ where the matrix $\mathbf{K}_c$ defined as $[K_c]_{i,k} = [4k(k+1)\sin\left(\frac{(2k+1)\pi s_i}{s_{\max}}\right)\frac{\delta}{L_s}]_{i,k}$ appears as a corrective matrix. For the case of wall channels, all coefficients $A_j$ are non zero, and the condition $M = M_N$ ensures that the two first terms ($j = 0, 1$) in the expansion of $u$ have the same form for both models. Consequently, the condition $M = M_N$ (leading to $L_s/\delta = M/9$ for bulk channels and $M/4$ for wall channels) allows both interfacial mobility models to predict similar conductance values.

Furthermore, Fig. 16 shows that this approximate approach gives a very good estimation of the conductance in the case of the Leonard and Lemlich model. For the Navier slip, the estimation is pretty good up to $M \approx 100$. Beyond $M > 200$, the approximate approach with $n_{\max} = 200$ fails to determine the flow conductance. But, increasing $n_{\max}$ improves the accuracy of the calculation for high interfacial mobilities (not shown here).

## 2. Curved Plateau borders

In this part, we detail the results concerning the fluid flow conductances of wall and volume calculated by Surface Evolver (see Appendix A 2). The results obtained for 3D
problems appeared to be highly sensitive to the choice of the discretization (either the default element order "P2+P1" or the higher-order elements "P3+P2"), particularly for high interfacial mobilities. We checked that the solutions obtained with the "P3+P2" discretization are only weakly sensitive to refinement (from a normal to a fine mesh). Moreover, we checked that the results of a 3D simulation of straight channels (with "P3+P2" discretization and the "Normal" meshing) are similar to those obtained with our accurate 2D simulations (see Appendix B 1). The estimated maximal error ≈ 5% is obtained for high interfacial mobilities.

Figure 17a shows the evolution of the conductance as a function of the ratio $\delta/L$ for bulk PBs. For low ratio $\delta/L$, the conductance tends to the conductance of the straight channel. As the ratio $\delta/L$ increases, the conductance deviates significantly from that of a straight channel. Figure 17b shows the effect of the interfacial mobility $M_{N,v}$ on the flow conductance. At low mobility $M_{N,v} < 0.1$, the flow conductance is up to 1.5 times greater than that of rigid interfaces, and the deviation from the rigid case depends linearly on $M_{N,v} = 9L_s/\delta$ (as mentioned in Appendix B 1). For $M_{N,v} > 1$, the flow conductance is significantly higher than that of rigid interfaces. We observe that this effect is more pronounced for small ratio $\delta/L$ (i.e. corresponding to fairly straight channels). Therefore, the curvature of the interfaces tends to limit the effect of the mobility.

For wall PBs, the conductance depends on both $\delta/L$ and $R_g/L$ [Figs. 18(a) & (b)]. For low ratios $\delta/L$, we find the conductance of straight channels. For higher ratios $\delta/L$, the conductance deviates from the case of straight channels, especially for large ratio $R_g/L$. However, similarly to straight channels, the conductance $G_w$ increases from that of a rigid interface, $G_{w,0}$, to that of a perfectly mobile interface, $G_{w,\infty}$ [Fig. 18(c)]. Numerical results show that the evolution can be described roughly by: $(G_w - G_{w,0})/(G_{w,\infty} - G_{w,0}) \approx (0.41 M_{N,w})^{0.92}/[1 + (0.41 M_{N,w})^{0.92}]$ with $M_{N,w} = 4L_s/\delta$.

As previously done for PB volume, we provide various approximate formulas chosen to fit the numerical results:

* for bulk PBs:

$$ \frac{\mu G_{v,0}}{\delta^3} = \frac{5.226 \cdot 10^{-4} \frac{\delta}{L}}{1 - 1.815 \frac{\delta}{L} - 0.6398 \left(\frac{\delta}{L}\right)^2 + 1.754 \left(\frac{\delta}{L}\right)^3}, \quad (B13) $$

$$ G_v = \left[ 1 + \frac{C_n M_{N,v}}{1 + (C_d M_{N,v})^{p_v}} \right] G_{v,0}, \quad (B14) $$
<figure><img src="image_17.png" /><figcaption>FIG. 17. Bulk Plateau borders: (a) Fluid flow conductances $G_{v,0}$ with rigid interface ($M_{N,v} = 0$), (b) deviation of the conductance $G_v$ to the rigid case $G_{v,0}$ as a function of the interfacial mobility ($M_{N,v} = 9L_s/\delta$). The crosses correspond to the results of our numerical simulations. The dashed lines are plotted by using Eqs. (B13) & (B14). In (b), the full line corresponds to the results for straight channels (see Appendix B1).</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>FIG. 18. Wall Plateau borders: (a) Fluid flow conductances $G_{w,0}$ with rigid interface ($M_{N,w} = 0$), (b) fluid flow conductances $G_{w,\infty}$ with perfectly mobile interface ($M_{N,w} \to \infty$), (c) evolution of the conductance $G_v$ as a function of the interfacial mobility ($M_{N,v} = 4L_s/\delta$). The crosses correspond to the results of our numerical simulations. In (a) and (b), the dashed lines are plotted by using Eqs. B15 & B16. In (c), the full line roughly describes all the results.</figcaption></figure>

with

$$C_n = -1.436 \left(\frac{\delta}{L}\right)^2 - 1.742 \frac{\delta}{L} + 3.253,$$

$$C_d = 0.2176 \left(\frac{\delta}{L}\right)^2 + 0.0725 \frac{\delta}{L} + 0.0134,$$

$$p_v = -1.397 \left(\frac{\delta}{L}\right)^2 + 1.667 \frac{\delta}{L} + 0.455.$$
• for wall PBs:

$$ \frac{\mu G_{w,0}}{\delta^3} = \frac{3.654 \cdot 10^{-3} \left(\frac{\delta}{L}\right)}{1 + D_1 \frac{\delta}{L} + D_2 \left(\frac{\delta}{L}\right)^2}, \qquad (B15) $$

$$ \frac{\mu G_{w,\infty}}{\delta^3} = \frac{2.731 \cdot 10^{-2} \left(\frac{\delta}{L}\right)}{1 + D_3 \frac{\delta}{L} + D_4 \left(\frac{\delta}{L}\right)^2}, \qquad (B16) $$

$$ G_w = \frac{(G_{w,\infty} - G_{w,0}) (C_w M_{N,w})^{p_w}}{1 + (C_w M_{N,w})^{p_w}} + G_{w,0}, \qquad (B17) $$

with

$$ D_1 = 1.293 \left(\frac{L}{R_g}\right)^2 + 5.987 \frac{L}{R_g} - 2.832, $$

$$ D_2 = 5.785 \left(\frac{L}{R_g}\right)^2 - 11.94 \frac{L}{R_g} + 1.101, $$

$$ D_3 = 6.488 \frac{L}{R_g} - 1.609, $$

$$ D_4 = 5.014 \left(\frac{L}{R_g}\right)^2 - 8.191 \frac{L}{R_g} - 1.811, $$

$$ C_w = 0.316 + 0.368 \frac{\delta}{L} $$

$$ p_w = 0.104 \frac{\delta}{L} - 0.068 \frac{\delta L}{R_g} + 0.90. $$

As for the calculation of the volume, we adapt the formulas to PBs having small lengths such that $\delta_{v,w}^*(L) > \delta$. For such PBs, we calculate the flow conductances by assuming $\delta/L = \delta^*/L$ to calculate the terms located on the right-hand side of Eqs B15-B17. Therefore, for such PBs, the conductance are proportional to $\delta^3/\mu$ as expected for a vertex contribution to flow conductance [36].

## Appendix C: Model of a single bubble per pore

Building an analytical model for all ratios $r$ is a task that lies beyond our capabilities, but the case of one single bubble per pore can be addressed. We consider the same packing as that used in our numerical model [Fig. 19(a)], i.e. four spherical grain of radius $R_g$ packed in a cube of side length $L_t = 2(\frac{2\pi}{3c})^{1/3} R_g$. In this compact configuration, the pore space is made of 12 pores: 8 tetrahedral sites and 4 octahedral sites. We assume that one single bubble fill each pore, i.e. a film is located in each constriction separating tetrahedral and octahedral sites [Figs. 19(b) & 19(c)]. The average "dry" bubble size is equal to $D_{cm} = (\frac{1-c}{2\pi})^{1/3} L_t$ and the size ratio to $r_c = ((1-c)/(3c))^{1/3} \approx 0.58$ with $c = 0.63$ as in our numerical model. Note that, in this configuration, no core foam fills the pore space and the liquid is distributed over 96 wall PBs and 24 bridges. Therefore, the liquid volume
<figure><img src="image_19.png" /><figcaption>FIG. 19. (a) view of the packing of the four spheres and the elemental cube, (b) clipped view of (a) with membranes in blue, (c) view after cutting with the median plane, (d) detail on films and on the PBs connecting the vertices located in the median plane to vertices below.</figcaption></figure>

$V_L$ included in the pore volume $V_L = \phi_L (1-c) L_t^3$ should respect the following equation: $\phi_L (1-c) L_t^3 = 96V_w (\delta, R_g, L) + 24V_m (\delta, R_g)$, where the volume of the wall PBs, $V_w$, is given by Eq. (A2) and that of bridges, $V_m$, by Eq. (A3). As in the numerical model, this equation can be solved to determine a relationship between $\delta/D_{cm}$ and $\phi_L$.

To determine the permeability of this configuration, we have to calculate the flow rate, $Q_t$, passing through the median plane [Fig. 19(d)]. Sixteen wall PBs of length $L = \pi R_g/3$ connect vertices located in the median plane to vertices below. We consider that the fluid flow is generated by imposing a pressure drop, $\Delta P_t$, to the elemental cube. Due to the symmetry of the configuration, the pressure drop through each PB emerging from the median plane [Fig. 19(d)] is equal to $\Delta P_t/4$. The flow rate, $q_w$, passing through each PB is therefore $q_w = G_w \Delta P_t/4$, leading to $Q_t = 4G_w \Delta P_t$. Finally, after some calculations, we obtain an expression relating the permeability $k = \mu Q_t/(L_t \Delta P_t)$ to the conductance of the wall PBs, $G_w$:

$$ \frac{k}{D_{cm}^2} = c \left( \frac{\pi^4}{108(1-c)} \right)^{2/3} \frac{\mu G_w}{L^3}, \quad (C1) $$

where $G_w$, is given by Eq. (B17).
[1] W. R. Rossen, Foams in Enhanced Oil Recovery, in *Foams: Theory, Measurements and Applications*, edited by R. K. Prud'homme and S. A. Khan (Marcel Dekker, New York, 1996) Chap. 11.

[2] R. Farajzadeh, A. Andrianov, R. Krastev, G. J. Hirasaki, and W. R. Rossen, Foam-oil interaction in porous media: Implications for foam assisted enhanced oil recovery, *Advances in Colloid and Interface Science* **183-184**, 1 (2012).

[3] R. Szafranski, J. B. Lawson, G. J. Hirasaki, C. A. Miller, N. Akiya, S. King, R. E. Jackson, H. Meinardus, and J. Londergan, Surfactant/foam process for improved efficiency of aquifer remediation, *Progress in Colloid and Polymer Science* **111**, 162 (1998).

[4] L. L. Schramm, *Emulsions, Foam and Suspensions: Fundamentals and Applications* (Wiley-VCH, Weinheim, Germany, 2005).

[5] S. Quebaud, M. Sibaï, and J. P. Henry, Use of Chemical Foam for Improvements in Drilling by Earth-Pressure Balanced Shields in Granular Soils, *Tunnelling and Underground Space Technology* **13**, 173 (1998).

[6] C. Budach and M. Thewes, Application ranges of EPB shields in coarse ground based on laboratory research, *Tunnelling and Underground Space Technology* **50**, 296 (2015).

[7] W. R. Rossen, Foam in Porous Media, *Foams and Emulsions*, 335 (1999).

[8] G. J. Hirasaki and J. B. Lawson, Mechanisms of Foam Flow in Porous Media: Apparent Viscosity in Smooth Capillaries., *Society of Petroleum Engineers journal* **25**, 176 (1985).

[9] L. Zhong, J. Szecsody, M. Oostrom, M. Truex, X. Shen, and X. Li, Enhanced remedial amendment delivery to subsurface using shear thinning fluid and aqueous foam, *Journal of Hazardous Materials* **191**, 249 (2011).

[10] H. J. Bertin, O. G. Apaydin, L. M. Castanier, and A. R. Kovscek, Foam flow in heterogeneous porous media: Effect of cross flow, *SPE Journal* **4**, 75 (1999).

[11] I. Bouzid, J. Maire, S. I. Ahmed, and N. Fatin-Rouge, Enhanced remedial reagents delivery in unsaturated anisotropic soils using surfactant foam, *Chemosphere* **210**, 977 (2018).

[12] G. G. Bernard, Effect of Foam on Trapped Gas Saturation and on Permeability of Porous Mediato Water, *Society of Petroleum Engineers journal* **5**, 295 (1965).

[13] A. R. Kovscek and C. J. Radke, *Foams: Fundamentals and Applications in the Petroleum*
Industry, in *Foams: Fundamentals and Applications in the Petroleum Industry*, edited by L. L. Schramm (1994) Chap. 3, pp. 115–163.

[14] M. Ceccaldi, V. Langlois, M. Guéguen, D. Grande, S. Vincent-Bonnieu, and O. Pitois, Liquid relative permeability through foam-filled porous media: Experiments, Phys. Rev. Fluids 8, 024302 (2023).

[15] M. Ceccaldi, V. Langlois, O. Pitois, M. Guéguen, D. Grande, and S. Vincent-Bonnieu, Coarsening effects on the liquid permeability in foam-filled porous media, Physical Review Fluids 9, 074003 (2024).

[16] I. Cantat, S. Cohen-Addad, F. Elias, F. Graner, R. Höhler, O. Pitois, F. Rouyer, and A. Saint-Jalmes, *Foams: structure and dynamics* (Oxford University Press, 2013).

[17] S. Wang, S. Huang, T. Qiu, J. Yang, J. Zhong, and Q. Hu, Analytical Study of the Permeability of a Foam-Conditioned Soil, International Journal of Geomechanics 20, 06020019 (2020).

[18] O. Pitois, A. Salamé, Y. Khidas, M. Ceccaldi, V. Langlois, and S. Vincent-Bonnieu, Daisy-shaped liquid bridges in foam-filled granular packings, Journal of Colloid and Interface Science 638, 552 (2023).

[19] S. Koehler, S. Hilgenfeldt, and H. Stone, Foam drainage on the microscale: I. Modeling flow through single Plateau borders, Journal of Colloid and Interface Science 276, 420 (2004).

[20] R. A. Leonard and R. Lemlich, A study of interstitial liquid flow in foam. Part I. Theoretical model and application to foam fractionation, AIChE Journal 11, 18 (1965).

[21] S. A. Koehler, S. Hilgenfeldt, E. R. Weeks, and H. A. Stone, Foam drainage on the microscale II. Imaging flow through single Plateau borders, Journal of Colloid and Interface Science 276, 439 (2004).

[22] A. Anazadehsayed, N. Rezaee, J. Naser, and A. V. Nguyen, A review of aqueous foam in microscale, Advances in Colloid and Interface Science 256, 203 (2018).

[23] F. Rouyer, O. Pitois, E. Lorenceau, and N. Louvet, Permeability of a bubble assembly: From the very dry to the wet limit, Physics of Fluids 22, 1 (2010).

[24] A. Bensoussan, J. Lions, and G. Papanicolaou, *Asymptotic Analysis for Periodic Structures*, north-holland ed. (Amsterdam, 1978).

[25] R. Quey and L. Renversade, Optimal polyhedral description of 3D polycrystals: Method and application to statistical and synchrotron X-ray diffraction data, Computer Methods in Applied Mechanics and Engineering 330, 308 (2018).
[26] K. A. Brakke, The surface evolver, Experimental Mathematics **1**, 141 (1992).

[27] S. J. Cox, G. Bradley, S. Hutzler, and D. Weaire, Vertex corrections in the theory of foam drainage, *Journal of Physics Condensed Matter* **13**, 4863 (2001).

[28] F. Rouyer, B. Haffner, N. Louvet, Y. Khidas, and O. Pitois, Foam clogging, Soft Matter **10**, 6990 (2014).

[29] E. Lorenceau, N. Louvet, R. F., and O. Pitois, Permeability of aqueous foams, Eur. Phys. J. E **28**, 293 (2009).

[30] F. Brochard and P. G. de Gennes, Shear-Dependent Slippage at a Polymer/Solid Interface, Langmuir **8**, 3033 (1992).

[31] S. G. Hatzikiriakos, Wall slip of molten polymers, Progress in Polymer Science (Oxford) **37**, 624 (2012).

[32] V. Langlois, C. T. Nguyen, F. Detrez, J. Guilleminot, and C. Perrot, Permeability of polydisperse solid foams, Physical Review E **105**, 015101 (2022).

[33] E. Forel, E. Rio, M. Schneider, S. Beguin, D. Weaire, S. Hutzler, and W. Drenckhan, The surface tells it all: relationship between volume and surface fraction of liquid dispersions, Soft Matter **12**, 8025 (2016).

[34] R. Höhler, J. Seknagi, and A. Kraynik, Capillary pressure, osmotic pressure and bubble contact areas in foams, Soft Matter **17**, 6995 (2021).

[35] M. Pasquet, N. Galvani, A. Requier, S. Cohen-Addad, R. Höhler, O. Pitois, E. Rio, A. Salonen, and D. Langevin, Coarsening transitions of wet liquid foams under microgravity conditions, Soft Matter **19**, 6267 (2023).

[36] I. Cantat, S. Cohen-Addad, F. Elias, F. Graner, R. Höhler, O. Pitois, F. Rouyer, A. Saint-Jalmes, R. Flatman, and S. Cox, *Foams : Structure and Dynamics*. (Oxford University Press, Incorporated, 2013) p. 288.

[37] S. Torquato, Effective electrical conductivity of two-phase disordered composite media, Journal of Applied Physics **58**, 3790 (1985).

[38] See Supplemental Material at [].

[39] G. Lian and J. Seville, The capillary bridge between two spheres: New closed-form equations in a two century old problem, Advances in Colloid and Interface Science **227**, 53 (2016).

[40] A. V. Nguyen, Liquid drainage in single Plateau borders of foam, Journal of Colloid and Interface Science **249**, 194 (2002).
[41] W. Drenckhan, H. Ritacco, A. Saint-Jalmes, A. Saugey, P. McGuinness, A. van der Net, D. Langevin, and D. Weaire, Fluid dynamics of rivulet flow between plates, Physics of Fluids 19, 102101 (2007).