# Spontaneous imbibition dynamics in two-dimensional porous media: a generalized interacting multi-capillary model

Shabina Ashraf, Yves Meheust, Jyoti Phirani

► To cite this version:

Shabina Ashraf, Yves Meheust, Jyoti Phirani. Spontaneous imbibition dynamics in two-dimensional porous media: a generalized interacting multi-capillary model. Physics of Fluids, 2023, 35, pp.012005. 10.1063/5.0123229. insu-03895688

HAL Id: insu-03895688

https://insu.hal.science/insu-03895688

Submitted on 13 Dec 2022

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
# Spontaneous imbibition dynamics in two-dimensional porous media: a generalized interacting multi-capillary model

**Accepted Manuscript:** This article has been accepted for publication and undergone full peer review but has not been through the copyediting, typesetting, pagination, and proofreading process, which may lead to differences between this version and the Version of Record.

Cite as: Physics of Fluids (in press) (2022); https://doi.org/10.1063/5.0123229

Submitted: 29 August 2022 • Accepted: 09 December 2022 • Accepted Manuscript Online: 10 December 2022

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

Shabina Ashraf,

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

Yves Meheust and

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

Jyoti Phirani

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

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

View Online

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

Export Citation

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

CrossMark

## ARTICLES YOU MAY BE INTERESTED IN

**Numerical study on the pulsating energy evolution in the cavitating flow around a mini Cascade**

Physics of Fluids **34**, 123308 (2022); https://doi.org/10.1063/5.0122844

**The effect of ribs on the sound radiation directivity of rectangular plates**

Physics of Fluids **34**, 122109 (2022); https://doi.org/10.1063/5.0127918

**Multiphysics modeling of femtosecond laser-copper interaction: from electron dynamics to plasma eruption**

Physics of Fluids (2022); https://doi.org/10.1063/5.0131001

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

Physics of Fluids

Special Topic: Food Physics

Submit Today!

<figure><img src="image_10.png" /></figure>
# Spontaneous imbibition dynamics in two-dimensional porous media: a generalized interacting multi-capillary model

Shabina Ashraf,¹ Yves Méheust,¹ and Jyoti Phirani², a)

¹) Univ. Rennes, CNRS, GosciencesRennes – UMR6118, F – 35000Rennes, France

²) Department of Civil and Environmental Engineering, University of Strathclyde, Glasgow.

(Dated: 7 December 2022)

The capillary bundle model, wherein the flow dynamics of a porous medium is predicted from that of a bundle of independent cylindrical tubes/capillaries whose radii are distributed according to the medium's pore size distribution, has been used extensively. But, as it lacks interaction between the flow channels, this model fails at predicting complex flow configuration, including those involving two-phase flow. We propose here to predict spontaneous imbibition in quasi-two-dimensional (quasi-2D) porous media from a model based on a planar bundle of interacting capillaries. The imbibition flow dynamics, and in particular, the breakthrough time, the global wetting fluid saturation at breakthrough, and which capillary carries the leading meniscus, are governed by the distribution of the capillaries' radii and their spatial arrangement. For an interacting capillary system consisting of 20 capillaries, the breakthrough time can be 39% smaller than that predicted by the classic, non-interacting, capillary bundle model of identical capillary radii distribution, depending on the spatial arrangement of the capillaries. We propose a stochastic approach to use this model of interacting capillaries for quantitative predictions. Comparing bundles of interacting capillaries with the same capillary diameter distribution as that of the pore sizes in the target porous medium, and computing the average behavior of a randomly-chosen samples of such interacting capillary bundles with different spatial arrangements, we obtain predictions of the position in time of the bulk saturating front, and of that of the leading visible leading front, that agree well with measurements taken from the literature. This semi-analytical model is very quick to run and could be useful to provide fast predictions on one-dimensional spontaneous imbibition in porous media whose porosity structure can reasonably be considered two-dimensional, e.g., paper, thin porous media in general, or layered aquifers.

## I. INTRODUCTION

When a wetting fluid is placed in contact with a porous medium, the fluid spontaneously imbibes into the pore space due to capillary suction. Such spontaneous imbibition in the porous matrix is crucial for applications such as oil recovery from reservoirs¹⁻³, Paper Analytic Devices (μPADs)⁴,⁵, textiles⁶, inkjet printing⁷,⁸, microfluidics⁹⁻¹³, lab-on-chip devices¹⁴,¹⁵, point-of-care diagnostics¹⁶,¹⁷, Polymer Electrolyte Membrane Fuel Cell (PEMFC)¹⁸,¹⁹, micro heat pipes²⁰,²¹, in understanding the motion of blood cells²² and in the design of bio-inspired drainage and ventilation systems²³. Capillary driven imbibition in a homogeneous porous medium follows diffusive dynamics, where the imbibition length is proportional to the square root of time²⁴⁻²⁶. This kind of dynamics was first characterized by Lucas²⁷ and³⁰ Washburn²⁸ for a horizontal cylindrical capillary tube: during the spontaneous imbibition of a wetting fluid of viscosity μ in a tube of radius r, the imbibition length (which here is simply the longitudinal position of the meniscus along the tube) is given by

$$l = \sqrt{\frac{r\sigma \cos \theta_w}{2\mu}} t, \quad (1)$$

where σ is the surface tension coefficient and θw is the wetting angle of the invading fluid on the tube's wall. In Eq. (1), the prefactor of the $\sqrt{t}$ law is proportional to $\sqrt{r}$, which implies that at any given time the meniscus will have advanced

more along a capillary of larger radius than along one of smaller radius. Later, the phenomenon of imbibition in a single pore/tube was observed to be strongly dependent on the geometries of the capillaries²⁹⁻⁴⁰.

Due to the similarity in the macroscopic laws describing the time evolution of the imbibition length between imbibition in a capillary tube and imbibition in a homogeneous porous medium, the capillary bundle model, considering a bundle of non-interacting capillaries of different radii, is classically considered as a proxy for porous media, in particular, soils⁴¹⁻⁴⁴. However, in a naturally occurring porous medium, the pores are of various shapes and sizes, and are interconnected⁴⁵,⁴⁶. In a quasi-two-dimensional (2D) porous medium such as paper, Bico and Quéré⁴⁷ showed that there are two imbibing fronts, a leading front in the small pores and a bulk saturating front which lags behind, which is contradictory to the predictions of the classic bundle of (non-interacting) capillaries, where the pores with larger radii have the leading front during imbibition.

The model geometry consisting of interacting capillaries (i.e., a capillary bundle where an opening allowing fluid exchange exists between adjacent capillaries, see e.g. Ref. 55) accounts for the effect of the interaction between pores on the pore scale flow dynamics, which in turn affects the Darcy scale flows in porous media⁴⁸⁻⁵⁵. In a system of two interacting capillaries, the imbibition in the capillary of smaller radius is found to be faster than that in the one of larger radii, unlike the behavior suggested by Eq. (1). However, a majority of these models were limited to predicting the imbibition dynamics in an ordered arrangement of pores or in two and three interacting capillary systems. For a system consisting of three interacting non-cylindrical capillaries, Unsal

<footnote>a) Electronic mail: jyoti.phirani@strath.ac.uk</footnote>
et al.⁵⁶⁻⁵⁸ showed experimentally that the imbibition speed is⁵⁸ fastest in the capillary of least effective radius. On the contrary, Ashraf et al.⁵⁵, using a one-dimensional lubrication approximation model and considering a system of three interacting cylindrical capillaries, showed that imbibition is not always fastest in the capillary of smallest radius. Furthermore,⁴³ both these studies⁵⁵,⁵⁶ showed that, for three capillary systems, the random positioning of the capillaries strongly impacts the invasion behaviour. But how the interconnected capillaries,⁴⁶ between capillaries impacts the overall imbibition dynamics⁴⁷ is far from being fully understood in the general case of a⁴⁸ larger number of tubes. Consequently, interacting capillary⁴⁹ systems, despite having a complexity which is intermediate between that of the classical bundle of non-interacting capillaries,⁵¹ have so far not been used to predict the generalized imbibition phenomenon observed in porous media consisting of⁵³ several pores of irregular sizes and varying connectivity. To⁵⁴ this aim, more complex models have been introduced since⁵⁵ based on pore-network geometries inferred from a geometry⁵⁶ cal analysis of the porous medium in which imbibition is investigated⁵⁹⁻⁶¹. We will present here a model of interaction between complexity between those early interacting-capillary⁵⁹ models and pore network models. Note that in many practical cases, the detailed porous structure is not known, and only⁶¹ an estimate of the pore size distribution is available; in such cases a pore network model cannot be applied without making assumptions on the unknown structure, whereas the model⁶³ presented here can be applied directly.⁶⁵

We thus propose a generalized one-dimensional model to predict spontaneous imbibition in a capillary bundle consisting of any number of randomly arranged cylindrical tubes that interact with each other, with any arbitrary distribution of the capillaries' radii. The model generalizes the study by Ashraf et al.,⁵⁵ for systems of two and three interacting capillaries,¹⁶⁸ to an arbitrary number of interacting capillaries. It is meant to model spontaneous imbibition in quasi-2D porous media⁶⁹ for which the pore size distribution is known. The model is⁷⁰ inspired from a model developed to tackle spontaneous imbibition in stratified geological porous media⁶². The two models are formally very similar to each other, but, due to the difference in geometries (flat layers for the stratified geology,⁷⁴ ical formation, cylindrical tubes in the present model), the equations are not identical. More importantly, the two studies differ widely in that the relative positioning of the layers in a geological medium is given, whereas, for a quasi-2D porous medium whose pore size distribution is known, the relative positioning of connected capillaries of different diameters within the 2D bundle that can predict the medium's behavior is not known a priori. Here, we explain the underlying physical phenomena causing the meniscus to advance at different rates in the different capillaries, and demonstrate that both the spatial arrangement of the interacting capillaries, and for a given arrangement, the contrasts in the capillaries' radii (i.e., their ratios), are crucial in predicting the imbibition dynamics. In contrast to the standard (non-interacting) capillary,⁸⁸ bundle, this model provides predictions that are qualitatively consistent with the phenomenology of spontaneous imbibition in real (quasi-)two-dimensional (2D) porous media. In

particular, this model correctly predicts that the smaller pores carry the leading front, while the larger pores carry the lagging saturating front responsible for the mass uptake of fluid in the porous medium, as measured in a paper-based porous medium⁴⁷. Furthermore, we provide a successful quantitative comparison between the measurements of Bico & Quéré on the leading and lagging imbibition fronts to predictions of the model obtained using a stochastic approach: the predicted behavior is the average of those obtained for all possible spatial organizations of the capillaries' diameter distribution. Though less accurate than fully numerical (and much more complicated) pore network models, this semi-analytical model has the advantage of running within seconds on any computer.

The presentation is organized as follows. We first review the model by Ashraf et al.,⁵⁵ (section II A). We then proceed to extend it to a system consisting of 4 interacting capillaries (section II B), before presenting the generalized one-dimensional model predicting spontaneous imbibition in an interacting multi-capillary system (section II C). We then examine the imbibition dynamics in a system of four interacting capillaries (section III A) and in a similar system consisting of 20 capillaries (section III B). In the discussion, we first compare the predictions of our model to those of the classic, non-interacting, capillary bundle (section III C 1), and, finally, confront its predictions of the leading and lagging fronts in a real quasi-2D porous medium from the literature to the published experimental measurements (section III C 2). Section IV contains a summary of the work and conclusive remarks, and discusses prospects to this study.

## II. MODELS

### A. Capillary imbibition in interacting capillaries

Using the capillary system shown in Fig. 1, Ashraf et al.,⁵⁵ used volume of fluid⁶³ (VOF) two-phase flow simulations to study spontaneous imbibition in a bundle of two or three interacting capillaries. These CFD (computational fluid dynamics) calculations provided the entire pressure and velocity fields inside the connected capillaries. They showed that (1) the invading wetting fluid transfers between two adjacent capillaries from the capillary of larger radius to that of smaller radius, but this transfer occurs only in the immediate vicinity of the (less advanced) meniscus of the capillary of larger radius; (2) that everywhere else (that is, everywhere except in the vicinity of that meniscus), the flow in the capillaries is not perturbed by the transfer of fluid between the capillaries; and (3) that, consequently, the pressure can be considered uniform over all transverse sections of the capillary system where both capillaries are filled with the same fluid, since no flow occurs along the transverse direction (if one neglects the small regions in the vicinity of the less advanced meniscus). These findings (1-3) served as basic assumptions to develop a reduced order, Washburn-like one dimensional model for a bundle of two and three interacting capillaries that can interact hydrodynamically with the neighbouring capillaries along their touching sides. The model predicted that in a bundle of two interacting capillaries,⁴⁶
ing capillaries the meniscus in the capillary of smaller radius moves ahead of the other one during the spontaneous imbibition, in consistency with the results of the VOF simulations. In this study we shall generalize the reduced order model of Ashraf et al.,⁵⁵ to an arbitrary number of capillaries positioned in the same plane and interacting with their neighbours.

<figure><img src="image_11.png" /><figcaption>FIG. 1. Spontaneous imbibition in two interacting capillaries, (a) cross-sectional view, (b) lateral view showing the contact angle $\theta_w$.</figcaption></figure>

For a flat bundle of three interacting capillaries, the model of Ashraf et al.,⁵⁵ showed that the distribution of radii and the spatial arrangement of the capillaries impact the imbibition behavior in the capillary system significantly. The meniscus in the capillary of smallest radius does not always move ahead of the others.

In the following sections, we examine the dynamics of the menisci during spontaneous imbibition in a flat bundle containing an arbitrary number of interacting capillaries. This generalization of the interacting capillaries' model follows the model development formulations from the study of Ashraf et al.,⁶² for imbibition in stratified porous media. In a stratified porous medium, the contrasts in layer transmissivities and the relative positioning of the layers control the imbibition dynamics, whereas in the present interacting capillaries the model, the positioning of the capillaries also plays a crucial role, but the role played by the transmissivities in the stratified medium is played by the product of the capillaries' permeabilities by their cross-sectional areas, both of which are controlled by the contrasts in the capillaries' radii.

We first describe below the one-dimensional model formulation for a system of four interacting capillaries to understand the underlying equations, before generalizing the model to multiple-interacting capillary system.

## B. Model development for four interacting capillaries

To predict the dynamics of spontaneous imbibition in porous medium using a system of interacting capillaries, we need to take the arrangement of capillaries into account, unlike for the classic capillary bundle (sometimes called bundle-of-tubes) model. For a porous medium made of $n$ interacting capillaries, there are $n!/2$ different arrangements. Fig. 2 shows a bundle of four interacting capillaries that are ordered spatially according to their radii $r_\alpha > r_\beta > r_\gamma > r_\delta$; we call this arrangement $\alpha\beta\gamma\delta$. The capillary pressure in tube (i = $\alpha, \beta, \gamma, \delta$) is given by the Young-Laplace equation as

$$Pc_i = \frac{2\sigma \cos \theta_w}{r_i}, \qquad (2)$$

where $\sigma$ is the surface tension and $\theta_w$ the contact angle. Hence, $Pc_\alpha < Pc_\beta < Pc_\gamma < Pc_\delta$. The corresponding imbibition lengths in the tubes at any time $t$ are denoted respectively.

<figure><img src="image_12.png" /><figcaption>FIG. 2. Schematic showing the spontaneous imbibition in an ordered system of four interacting capillaries. The imbibition lengths in capillaries $\alpha, \beta, \gamma, \delta$ of radii $r_\alpha, r_\beta, r_\gamma, r_\delta$ are denoted by $z_\alpha, z_\beta, z_\gamma, z_\delta$, respectively. The cross section of the system of interacting capillaries is also shown.</figcaption></figure>

We consider the assumptions from Ashraf et al.⁵⁵, according to which (1) the pressure equilibrates over the sections of the capillary system that are entirely filled with the invading fluid, and (2) fluid transfers from a capillary having a larger radius to an adjacent capillary having a smaller radius just before the meniscus, which in the model we assume to occur at the position of the meniscus. We show this fluid transfer between adjacent capillaries in the vicinity of the meniscus by vertical arrows in Fig. 2. We consider the interaction between the capillaries to be sufficiently low for the Poiseuille flow in each of the capillaries to be maintained. At any given time $t$, the less advanced meniscus (i.e., that for which the imbibition length is the smallest) will be in the capillary for which the driving capillary pressure jump across the meniscus is the smallest, hence it will be the meniscus in the $\alpha$ capillary. For $z < z_\alpha(t)$, the pressure field must be identical in all capillaries. Similarly, the next-less-advanced meniscus is necessarily the $\beta$ capillary driven by the capillary pressure $Pc_\beta$, so at any time $t$ the pressure field is identical in capillaries $\beta, \gamma$ and $\delta$ for $z_\alpha(t) < z < z_\beta(t)$, and so forth: the pressure field is identical in the $\delta$ and $\gamma$ capillaries for $z_\beta(t) < z < z_\gamma(t)$. The imbibition length in capillary $\delta, z_\delta(t)$ is the largest at any time $t$.

We now consider one of the random arrangements as shown in the schematic of Fig. 3, where the order of arrangement of the capillaries is $\beta\gamma\alpha\delta$. It was explained by Ashraf et al.,⁵⁴ that, for a randomly-arranged interacting capillary system, the meniscus in the smallest radius capillary does not always lead. For this arrangement, depending upon the contrasts in the radii, three different positionings of the menisci are possible as shown in Fig. 3 (a), (b) and (c). At any given time $t$, for $0 < z_\alpha(t)$, the pressure field is identical in all capillaries, and the pressure drop from the inlet to $z_\alpha(t)$ is $Pc_\alpha$. For $z > z_\alpha(t)$, the imbibing fluid is continuous in the capillaries $\beta$ and $\gamma$, since they are connected. Therefore, the pressure field is the same in the capillaries $\beta$ and $\gamma$ for $z_\alpha(t) < z < z_\beta(t)$. As $r_\beta > r_\gamma$ (meaning that the capillary suction in $\beta$ is less than that in $\gamma$), during the spontaneous imbibition, $z_\beta(t) < z_\gamma(t)$, at all times. Although the capillary $\delta$ is filled with the imbibing phase, the non-wetting fluid in $\alpha$ disconnects it from capillaries $\beta, \gamma$ for $z > z_\alpha(t)$. Therefore, for $z > z_\alpha(t)$ the pressure field in $\delta$ can be different from that in $\beta, \gamma$. For the arrangement $\beta\gamma\alpha\delta$ shown in the schematic of Fig. 3, $z_\alpha < z_\beta < z_\gamma$ and $z_\alpha < z_\delta$ during the imbibition process and the position of
<figure><img src="image_13.png" /><figcaption>FIG. 3. Spontaneous imbibition in a system of four interacting capillaries with a spatial arrangement of βγαδ of the capillaries. The imbibition lengths in capillaries α, β, γ, δ of radii $r_α$, $r_β$, $r_γ$, $r_δ$ are $z_α(t)$, $z_β(t)$, $z_γ(t)$, $z_δ(t)$, respectively. The schematics of the imbibition phenomenon show the fluid transfer at menisci locations with arrows. For this spatial arrangement, depending upon the contrasts in the capillaries' radii, the possible orders in the invasion lengths can be (a) $z_α < z_β < z_γ < z_δ$, (b) $z_α < z_β < z_δ < z_γ$ and (c) $z_α < z_δ < z_β < z_γ$. The cross section of the system of interacting capillaries is also shown for (a).</figcaption></figure>

283 $z_δ(t)$ relative to $z_β(t)$ and $z_γ(t)$ depends on the contrasts in the capillaries' radii.

284

285 The detailed development of the generalized one-dimensional model for this system of four interacting capillaries with arrangement βγαδ is described in Appendix A.

287 The pressure drop across each of the sections is determined

288 individually, i.e., for sections (I) $0 < z < z_α$, (II) $z_α < z < z_β$,

289 (III) $z_β < z < z_γ$, and (IV) $z_α < z < z_δ$. As spontaneous imbibition is driven by capillary forces, the sum of the pressure

290 drops across all the sections of a capillary is equal to the capillary pressure of that capillary.

291

292

293

$$Pc_i = \left( \sum_j P_{i(j)} \right), \qquad (3)$$

294 where $P_{i(j)}$ is the pressure drop across the section of index

295 $j=(I), (II), (III), (IV)$ of the capillary of index $i=\alpha, \beta, \gamma, \delta$

296 By solving the system of equations expressing (i) Darcy's law

297 in each of the capillaries, and (ii) the relations between the

298 meniscii's advancement and the fluid velocities and fluid ex

299 change between the capillaries, we obtain the equations governing the flow in the interacting capillaries, which are,

300

$$Pc_{\alpha} = \frac{8\mu z_{\alpha}(t)}{r_{\alpha}^{4} + r_{\beta}^{4} + r_{\gamma}^{4} + r_{\delta}^{4}} \left( r_{\alpha}^{2} \frac{dz_{\alpha}}{dt} + r_{\beta}^{2} \frac{dz_{\beta}}{dt} + r_{\gamma}^{2} \frac{dz_{\gamma}}{dt} + r_{\delta}^{2} \frac{dz_{\delta}}{dt} \right) \quad (4)$$

$$Pc_{\delta} - Pc_{\alpha} = \frac{8\mu(z_{\delta}(t) - z_{\alpha}(t))}{r_{\delta}^{2}} \left( \frac{dz_{\delta}}{dt} \right), \qquad (5)$$

$$Pc_{\beta} - Pc_{\alpha} = \frac{8\mu(z_{\beta}(t) - z_{\alpha}(t))}{r_{\beta}^{4} + r_{\gamma}^{4}} \left( r_{\beta}^{2} \frac{dz_{\beta}}{dt} + r_{\gamma}^{2} \frac{dz_{\gamma}}{dt} \right), \qquad (6)$$

$$Pc_{\gamma} - Pc_{\beta} = \frac{8\mu(z_{\gamma}(t) - z_{\beta}(t))}{r_{\gamma}^{2}} \left( \frac{dz_{\gamma}}{dt} \right). \qquad (7)$$

Eqs. (4) to (7) are rendered non-dimensional by normalizing the positions by the total capillary system's length, L, and time by $[8\mu L^2/(Pc_\alpha r_\alpha^2)]$, thus defining the non-dimensional positions and times

$$Z_i = \frac{z_i}{L}, \quad i = \alpha, \beta, \gamma, \delta \quad \text{and} \quad T = \frac{Pc_\alpha r_\alpha^2}{8\mu L^2 t}. \qquad (8)$$

Introducing the contrasts in radii, $\lambda_i = r_i/r_\alpha$, and in capillary pressures, $\epsilon_i = Pc_i/Pc_\alpha$, for $i = \beta, \gamma, \delta$, we then obtain the non-dimensional form of Eqs. (4) to (7) as

$$1 = \frac{Z_{\alpha}}{1 + \lambda_{\beta}^{4} + \lambda_{\gamma}^{4} + \lambda_{\delta}^{4}} \left( \frac{dZ_{\alpha}}{dT} + \lambda_{\beta}^{2} \frac{dZ_{\beta}}{dT} + \lambda_{\gamma}^{2} \frac{dZ_{\gamma}}{dT} + \lambda_{\delta}^{2} \frac{dZ_{\delta}}{dT} \right), \qquad (9)$$

$$\epsilon_{\delta} - 1 = \frac{Z_{\delta} - Z_{\alpha}}{\lambda_{\delta}^{2}} \left( \frac{dZ_{\delta}}{dT} \right), \qquad (10)$$

$$\epsilon_{\beta} - 1 = \frac{Z_{\beta} - Z_{\alpha}}{\lambda_{\beta}^{4} + \lambda_{\gamma}^{4}} \left( \lambda_{\beta}^{2} \frac{dZ_{\beta}}{dT} + \lambda_{\gamma}^{2} \frac{dZ_{\gamma}}{dT} \right). \qquad (11)$$

$$\epsilon_{\gamma} - \epsilon_{\beta} = \frac{Z_{\gamma} - Z_{\beta}}{\lambda_{\gamma}^{2}} \left( \frac{dZ_{\gamma}}{dT} \right), \qquad (12)$$

Further assuming that the contact angle $\theta_w$ is the same in all capillaries, we have $\epsilon_i = 1/\lambda_i$, and upon rearranging the governing Eqs. (9) to (12) and adding them, we obtain,

$$2 \left( 1 + \sum_{i=\beta,\gamma,\delta} \epsilon_i \lambda_i^4 \right) T = Z_{\alpha}^2 + Z_{\beta}^2 \lambda_{\beta}^2 + Z_{\gamma}^2 \lambda_{\gamma}^2 + Z_{\delta}^2 \lambda_{\delta}^2. \qquad (13)$$

Eq. (13) expresses that, in a system of interacting capillaries, the sum of the squares of the product of the non-dimensional radius with the non-dimensional distance invaded in all the capillaries is proportional to the invasion time $T$. For different arrangements of a system of 4 interacting capillaries having the same contrasts in capillary radii, the total capillary suction of the system remains the same. Therefore, for any of the $4!/2 = 12$ possible arrangements, rearranging the equations governing the imbibition process, and adding them, leads to Eq. (13). However, the velocity at which the individual meniscii travels in each of the tubes depends on the particular arrangement of the capillaries.
### C. Generalizing the one-dimensional spontaneous imbibition model in the interacting capillary system

Equation (13) is readily generalized to a system of $n$ interacting capillaries, in the form

$$2 \left( \sum_{i=1}^{n} \varepsilon_i \lambda_i^4 \right) T = \sum_{i=1}^{n} \psi_i Z_i \quad (14)$$

where $\psi_i = \pi r_i^2 z_i / (\pi r_\alpha^2 L)$ ($j = 1, 2, ..., n$) is the non-dimensional volume imbibed in the capillary of index $i$.$^{370}$ Eq. (14) expresses that the sum over all capillaries of the non-dimensional volumes times the corresponding non-dimensional imbibition lengths, is proportional to time. This$^{373}$ can be compared to the dynamics in a bundle of non-interacting capillaries, for which we know that the dynamics are diffusive, i.e., for each of the capillaries, the imbibed length square is proportional to time.$^{377}$

We note from the derivation of Eq. (13) for the system consisting of four capillaries, that each arrangement of the capillaries will have a different set of governing equations for menisci positions with time. This is because the knowledge of the arrangement is required to determine the regions of the capillaries across which the pressure equilibrates and the locations of fluid transfers. Therefore, for a system of $n$ interacting capillaries, we now propose an algorithm which can determine the imbibition behaviour in the bundle of interacting capillaries and form the governing equations for a generalized model of such systems of $n$ interacting capillaries. A MATLAB program has been written to implement this algorithm and obtain the advancement of the menisci, $z_l(t)$, where $l = 1, 2, 3, ..., n$, as a function of time. The step-by-step procedure is described in detail in Appendix B, but its principles can be described in the following manner.$^{384}$

<figure><img src="image_14.png" /><figcaption>FIG. 4. Schematic of spontaneous imbibition in an $n$-capillary system where the capillaries are positioned randomly. The capillaries are arranged in the capillary radii are denoted as $r_a, r_b, r_c, ..., $ and the corresponding imbibition distances at time $t$ are denoted by $z_a(t), z_b(t), z_c(t), ...$.</figcaption></figure>

First, the algorithm searches for the capillary of largest radius in the arbitrary arrangement, whose meniscus position is

$z_a$ at a given time; it is denoted $C_i$ in Fig. 4, where the capillaries in the order of arrangement are denoted from $C_1$ to $C_n$. The pressure drop in the region $0 < z < z_a$ is determined for all the capillaries and the algorithm then considers two regions: the 'top region' consisting of the capillaries $C_1$ to $C_{i-1}$ and the 'bottom region' consisting of the capillaries $C_{i+1}$ to $C_n$ (see Fig. 4). The largest radius capillaries in each of these two regions are determined and the pressure drop in the respective regions are determined for sections $z_a < z < z_b$ and $z_a < z < z_c$. Now, each of these two regions is further divided into two subregions each, i.e., containing the capillaries $C_1$ to $C_{j-1}$ on the one hand and $C_{j+1}$ to $C_{i-1}$ on the other hand in the 'top region', and $C_{i+1}$ to $C_{k-1}$ on the one hand and $C_{k+1}$ to $C_n$ on the other hand in the 'bottom region'. The pressure drops are determined in each of the subregions. This procedure is then performed recursively until the algorithm has identified the pressure drop in each of the sections for every capillary. It can then formulate the governing equations, which are consequently solved to obtain the advancement of all menisci as a function of time.

## III. RESULTS AND DISCUSSIONS

We first explore the imbibition of a system of four interacting capillaries, followed by the imbibition in a system consisting of 20 capillaries.

### A. Interacting four-capillary system

In section II A, we have anticipated that, in an ordered arrangement, the meniscus in the capillary of smallest radius, $\delta$, will always lead, followed by the capillary of second smallest radius, $\gamma$, as shown in Fig. 2, while the meniscus in the capillary $\alpha$ always lags behind. Solving the governing equations for this arrangement, we always get the same trend, i.e., $z_\alpha(t) < z_\beta(t) < z_\gamma(t) < z_\delta(t)$ for the imbibed lengths in the capillaries at any given time during the imbibition process. However, $4! / 2 = 12$ arrangements are possible for an interacting four-capillary system, for any given 4 radii of the capillaries. In section II B we chose one arrangement $\beta\gamma\alpha\delta$ and anticipated 3 cases of different relative positioning of menisci. The possibility of occurrence of these 3 cases depends upon the radii contrast in the capillaries. A change in radii contrast changes the pressure fields in the capillaries, which governs the menisci positions. Each of the 3 cases shown in Fig. 3 are shown in Fig. 5 (a), (c), (e). Solving Eqs. (9) to (12) over non-dimensional times, we show in Fig. 5 (b), (d), (f), how the relative positions of the plots of $Z_\beta, Z_\gamma, Z_\delta$ as a function of time change when the contrast in the radii of capillaries are changed according to the three configurations addressed in Fig. 5 (a), (c), (e).

We now consider two other random arrangements $\gamma\delta\alpha\beta$ and $\gamma\alpha\beta\delta$, which are illustrated in Figs. 6 and 7, respectively. In these figures, we show the schematic of the menisci locations at a given time during imbibition in (a), (c), (e). The corresponding time evolution of the positions of menisci in
<figure><img src="image_15.png" /></figure>

FIG. 5. Spontaneous imbibition in a system of four interacting capillaries which are positioned with respect to each other according to the arrangement $\beta\gamma\alpha\delta$, for three different contrasts in capillary radii. (a), (c), (e) represent the schematics of possible imbibition behavior at a given time $t$. The distribution of radii predicting the imbibition phenomenon are indicated in the plots (b), (d) and (f). The non-dimensional times at which the leading meniscus reaches the outlet end of the interacting capillary system ($T_{bt}$) for the cases (b), (d) and (f) are 0.43, 0.40 and 0.39, respectively.

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

FIG. 6. Spontaneous imbibition in a system of four interacting capillaries, spatially arranged as $\gamma\delta\alpha\beta$. Depending upon the contrasts in capillary radii, at a given time, the relative positions of the menisci vary. (a), (c), (e) represent the schematics of possible imbibition behavior. The non-dimensional meniscus positions and the radii contrasts corresponding to the schematics of (a), (c), (e) are shown in (b), (d), (f) respectively, as a function of the non-dimensional time. The times at which the invading fluid reaches the outlet end ($T_{bt}$) for the cases (b), (d) and (f) are 0.38, 0.40 and 0.39, respectively.

the four capillaries are shown in (b), (d), (f). Each of these figures shows that the contrast in the capillary radii, for a given arrangement, impacts the relative positions of the menisci at any given time. Conversely, in Figs. 5(f), 6(d), and 7(b), the radii of the capillaries in the interacting capillary system are identical but the arrangements of the capillaries are different. For the arrangement $\beta\gamma\alpha\delta$ shown in Fig. 5(f), the menisci's positions are ordered according to $Z_{\gamma} > Z_{\beta} > Z_{\delta} > Z_{\alpha}$ while for the arrangement $\gamma\delta\alpha\beta$ shown in Fig. 6(d) the menisci positions are ordered according to $Z_{\delta} > Z_{\beta} > Z_{\gamma} > Z_{\alpha}$ and so on. The arrangement $\gamma\alpha\beta\delta$ shown in Fig. 7(b), the menisci positions are ordered according to $Z_{\delta} > Z_{\gamma} > Z_{\beta} > Z_{\alpha}$. Hence, for an interacting multi-capillary system, both the contrasts in capillary radii and their arrangement are crucial in determining the imbibition behavior. The non-dimensional time

which the imbibing fluid first breaks through or reaches the non-dimensional length 1 in one of the interacting capillaries, and the radius of the capillary through which the breakthrough occurs, are impacted accordingly, as reported in the captions of Fig. 5, 6 and 7. Note that in Figs. 5, 6, 7, the schematics presented in (a), (c) and (e) are not necessarily to scale, either for the capillaries' radii (indicated in the legends of (b), (d) and (f)) or for the imbibition lengths.

We further illustrate the imbibition phenomenon in a system of four interacting capillaries for three arrangements out of the 12 possible arrangements in Fig. 8. The radii of the capillaries are $r_{\alpha} = 80$ m, $r_{\beta} = 60$ m, $r_{\gamma} = 40$ m, and $r_{\delta} = 20$ m for all the arrangements. In Fig. 8(a), where the capillaries are in the ordered arrangement ($\alpha\beta\gamma\delta$), the leading meniscus is in the capillary with the smallest radius ($\delta$). For the
<figure><img src="image_17.png" /><figcaption>FIG. 7. Spontaneous imbibition in a system of four interacting capillaries, spatially arranged as $\gamma\alpha\beta\delta$. Depending upon the contrasts in capillary radii, at a given time, the relative positions of the menisci vary. (a), (c), (e) represent the schematics of possible imbibition behavior. The non-dimensional meniscus positions and the radii contrasts corresponding to the schematics of (a), (c), (e) are shown in (b), (d), (f) respectively, as a function of the non-dimensional time. The times ($T_{bt}$) at which the invading fluid first reaches the outlet in any of the capillaries are 0.38, 0.42 and 0.38 for the cases (b), (d) and (f), respectively.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>FIG. 8. Spontaneous imbibition in a system of four interacting capillaries of radii $r_\alpha = 80$ m, $r_\beta = 60$ m, $r_\gamma = 40$ m and $r_\delta = 20$ m. The non-dimensional positions of the four meniscii are shown as functions of the non-dimensional time for three of the 12 possible arrangements (a) $\alpha\beta\gamma\delta$, (b) $\gamma\beta\alpha\delta$, and (c) $\alpha\delta\beta\gamma$ are shown. The relative position of the menisci with time and the breakthrough time depend upon the arrangement of the capillaries, for a given contrast in the radii.</figcaption></figure>

same contrast in radii and the arrangement $\gamma\beta\alpha\delta$ (Fig. 8(b))$^{464}$ the leading meniscus is in capillary $\gamma$. For arrangement $\alpha\delta\beta\gamma\delta$ shown in Fig. 8(c), the menisci in capillaries $\gamma$ and $\delta$ travel at the same velocity at all times. It can also be observed from Fig. 8 that the breakthrough times change with the arrangement of the capillaries; while the breakthrough for the ordered arrangement (Fig. 8(a)) occurs at $T = 0.33$, for the other two other arrangements shown in Fig. 8(b) and (c), the breakthrough occurs at $T = 0.40$. Similar plots are shown for all 12 possible arrangements in Fig. C.1 of Appendix C; all the arrangements are found to have breakthrough times in the range $T = 0.33$ to $T = 0.40$. For a wetting fluid of viscosity $10^{-3}$ Pa·s and surface tension of $73 \times 10^{-3}$ N/m impregnating the empty capillary system of length 1 m and with a maximum capillary radius of 80 m, the non-dimensional time corresponding to $T = 0.01$ is 6.84 s, so the breakthrough for the arrangements shown in Fig. 8 occurs between 225.7 s and 273.6 s. Hence, for the four-capillary system, we can summarize that the arrangement of the capillaries and the contrasts in

capillary radii significantly affect the breakthrough time and the index of the capillary through which breakthrough occurs.

## B. System consisting of 20 interacting capillaries

From the above analysis, we see that for any interacting multi-capillary system, the capillary having the leading meniscus and the breakthrough time both depend on the contrast in the capillary radii and on the spatial arrangement of capillaries. We now use the generalized model to predict imbibition in a system consisting of $n = 20$ interacting capillaries, focusing on the impact of the arrangement. We assume no spatial correlations in the capillaries' radii. The number of different arrangements for $n = 20$ is $20! / 2 = 1.216 \times 10^{18}$. We run the generalized model on 1000 random arrangements for capillaries whose radius distribution is uniform between 10 m (minimum radius) and 200 m (maximum radius).

We show in Fig. 9(a), the imbibition length in the capillary-
<figure><img src="image_19.png" /><figcaption>FIG. 9. Spontaneous imbibition in 7 systems of twenty interacting capillaries with identical radii but different spatial arrangements: 6 random arrangements and one ordered arrangement. (a) radii vs. imbibition length at $T = 0.2$; (b) radii vs. imbibition length at breakthrough time, $T = T_{bt}$, (c) saturation vs. longitudinal position at $T = 0.2$ and $T = 0.3$, (d) saturation vs. longitudinal position at breakthrough time, $T = T_{bt}$.</figcaption></figure>

ies vs the radii of the capillaries at the non-dimensional time $T = 0.2$ for 6 random arrangements (denoted arr1, arr2, arr3, arr4, arr5 and arr6 in the figure), and the ordered arrangement $T = 0.1$ through time and the capillary radius through which the break-513 through occurs can be observed for the given radii contrast of 514 the capillaries. We see from Fig. 9(a) that, at $T = 0.2$, the 515 capillary having the leading meniscus is different for different 516 arrangements and the menisci positions in the capillaries are 517 also dependent on the arrangement. For instance, at $T = 0.2$ 518 the meniscus in the capillary of radius 10 m (smallest radius) 519 has traveled a non-dimensional length of 0.79 for the ordered arrangement, whereas for random arrangement number 521 1, the non-dimensional length invaded in the smallest capillary is 0.51. In Fig. 9(b), we illustrate the relationship between 522 the radii and the imbibition length in all capillaries at break-523 through time. The breakthrough time for different arrangements is given in the legend of the arrangement in Fig. 9(b). 524 Breakthrough in the systems of 20 interacting capillaries occurs through different capillaries and at different times for the 525 6 random arrangements and the ordered arrangement.

The saturation at a given imbibition length $Z$ can be defined as the ratio of the cross-sectional area occupied by the imbibing fluid at $Z$ to the total cross-sectional area of the capillary 531 system, i.e., $((\sum_{j=1}^{n_f} r_j^2)/(\sum_{i=1}^n r_i^2))$, where $n_f(Z)$ is the number of capillaries filled by the imbibing fluid at $Z$, and the indices $j=34$ refer to all such capillaries. The plot of saturation vs. longitudinal position is shown in Fig. 9(c) at $T = 0.2$ and $T = 0.3$, for all the 7 spatial arrangements. These saturation profiles of the interacting capillary system depend significantly on the arrangement of the capillaries. For example, at $T = 0.3$, the saturation at $Z = 0.7$ is 0.43 for the random arrangement number 3, and 0.35 for the ordered arrangement as indicated in Fig. 9(c).

In Fig. 9(d) we show how saturation varies with the longitudinal position at breakthrough time for the 7 arrangements. The amount of non-wetting fluid displaced at the time of breakthrough is different between the different arrangements. We also observe from Fig. 9(a) that the random arrangements where the leading meniscus is in a capillary of larger radius, will have a longer breakthrough time as shown in Fig. 9(b). This will also cause the saturation of the random arrangement to be larger at the breakthrough time, which can be observed in Fig. 9(d).

However, since the contrast in the radii of the capillaries is identical for all arrangements, the effective capillary suction causing the imbibition phenomenon is also identical in all cases. Therefore, at a given time $T$, the global wetting fluid saturation in the interacting capillary system will be the same for all arrangements, which is determined as $S = \sum_{i=1}^{n_f} r_i^2 Z_i / \sum_{i=1}^{n_f} r_i^2$. The fraction of the interacting capillary system occupied with the imbibing phase at $T = 0.2$ is 0.55 and at $T = 0.3$, $S$ is 0.67 for all the 7 arrangements. But this is only applicable until breakthrough occurs in one of the arrangements.
<figure><img src="image_20.png" /><figcaption>FIG. 10. Radii of the capillaries in which breakthrough occurs vs. breakthrough time in 1000 randomly-chosen arrangements of a system of 20 interacting capillaries with radii uniformly distributed between 10 and 200 m (the upper boundary of the vertical scale is thus chosen to 200 m). The shortest breakthrough time is observed in the ordered arrangement, at $T = 0.31$, and the maximum observed breakthrough time is $T = 0.42$. The largest radius of a capillary through which breakthrough occurs is 100 m while the smallest one is 10 m.</figcaption></figure>

## C. Discussion

We now compare the predictions of our analytical model of interacting capillaries to those of the standard capillary bundle model, and discuss how the predictions of our model compare to experimental measurements in quasi-2D porous media. We use our model within a stochastic approach, that is, for a given number $n$ of capillaries of known radii we consider the average behavior of all $m = n!/2$ different spatial arrangements of the capillaries. When $m$ is too large to be tractable even for our very fast semi-analytical model (for example for $n = 20$, $m > 1.2110^{18}$), we consider the average behavior of a sufficiently large subsample of $R < m$ randomly-chosen spatial arrangements.

### 1. Confronting predictions from the classic (non-interacting) capillary bundle to our model

We show the spatial saturation profile for the classic capillary bundle model with $n = 20$ capillaries at three different times ($T = 0.1$, $T = 0.3$ and $T = T_{bt} = 0.5$) in Fig 11(a), and the average spatial saturation profile for 1000 randomly-chosen different spatial arrangements, for a system of 20 interacting capillaries (as predicted by our model) at the same three times in Fig 11(b). Note that the number of spatial arrangements was chosen after a convergence study which we present in Appendix D (see in particular Fig. D.1).

The capillary radii are identical in the two cases. For non-interacting capillaries, by non dimensionalizing the Washburn's law, $z_i^2 = (Pc_i r_i^2 / 4\mu)t$, we obtain

$$Z_i^2 = 2\epsilon_i \lambda_i^2 T, \qquad (15)$$

where $Z_i = z_i/L$ is the non-dimensional length imbibed in the capillary of radius $r_i$ and $L$ is the total length of the capillary system. The time is non-dimensionalised as $T = t(Pc_\alpha r_\alpha^2) / (8\mu L^2)$. In Eq. (15), $\epsilon = Pc_i / Pc_\alpha$ and $\lambda_i = r_i / r_\alpha$, where $Pc_\alpha$ and $r_\alpha$ are respectively the capillary pressure and radius of the widest capillary (200 m). The maximum value of $\epsilon_i$ and $\lambda_i$ are 1, which occurs for the largest radius capillary. For all other capillaries $\epsilon_i$ and $\lambda_i$ are always smaller than 1.

As discussed previously, in the classic capillary bundle model, imbibition follows Washburn's diffusive dynamics and therefore the invaded length is the largest in the capillary of largest radius. As illustrated in Fig. 11(a), due to the large cross-section area of that widest capillary, it contributes to a large fraction of the cross-sectional saturation for the bundle-of tubes model. On the contrary, in our interacting-capillary system, the largest radius capillary always has the least advanced meniscus, at any time. Consequently, the breakthrough time for the capillary bundle model is 136.5 s (at $T = 0.5$), at which the fractional volume occupied by the invading fluid is 0.86. This is considerably larger than the breakthrough time for interacting capillary systems, which occurs between 84.63 s and 114.66 s (between $T = 0.31$ and $T = 0.42$), depending on the configuration, and the fractional volumes occupied by the imbibing fluid across the 1000 arrangements lie between 0.69 and 0.79. In Fig. 11(b), we show

In Fig. 10 we have plotted the radius of the capillary having the leading meniscus vs. the breakthrough time for the 1000 randomly chosen arrangements, assumed to be representative of the entire statistics. We see that when a wetting fluid of viscosity of $10^{-3}$ Pa·s and surface tension of $73 \times 10^{-3}$ N/m² imbibes a twenty-capillary system of length 1 m and maximum capillary radius of 200 m, the non-dimensional time of $T = 0.01$ corresponds to 2.73 s. If such a wetting fluid were considered to imbibe into this interacting capillary system, the breakthrough which occurs between $T = 0.31$ and $T = 0.42$ corresponds to the dimensional times of 84.63 s and 114.66 s. Therefore, for the same contrast in capillary radii, the maximum mum and minimum breakthrough time are approximately 300 s apart, indicating that the breakthrough time significantly depends on the arrangement of the capillaries. It can also be observed from Fig. 10 that breakthrough in an ordered multi-capillary system occurs through the capillary of smallest radius at $T = 0.31$, which is the smallest breakthrough time compared to other arrangements. Fig 10 also shows that the largest radius of a capillary through which breakthrough occurs is as large as 100 m, while the minimum radius of the capillary through which breakthrough occurs is 10 m. For the arrangement number 6 (+ symbols), the leading meniscus is in the 100 m radius capillary and breakthrough occurs between $T_{bt} = 0.42$ as shown in Fig. 9(b). From Fig. 10, we also see that, when breakthrough occurs through the smallest radius, the capillary, the breakthrough time may vary between $T = 0.31$ and $T = 0.41$, and the total volume fraction of the interacting capillary system occupied by the invading phase can lie between 0.69 and 0.79. In contrast, if breakthrough occurs through the capillary of radius 70 m, the breakthrough time lies between $T = 0.38$ and $T = 0.42$ and the total volume fraction imbibed by the wetting phase lies between 0.76 and 0.8.
<figure><img src="image_21.png" /><figcaption>FIG. 11. (a) Spatial saturation profile during spontaneous imbibition in a bundle-of-tubes consisting of twenty non-interacting capillaries at $T = 0.1$, $T = 0.3$, and $T = T_{bt} = 0.5$. (b) Average spatial saturation profile for 1000 different spatial arrangements of the system consisting of twenty interacting capillaries of identical radii as in (a), at $T = 0.1$, $T = 0.3$, $T = T_{bt}$.</figcaption></figure>

the averaged saturation values along the length of the capillary system for all the 1000 arrangements of the twenty interacting capillary system at non-dimensional times $T = 0.1, 0.2$ and at breakthrough, i.e., $T_{bt}$. We see from fig. 11(b) that the standard deviation across the arrangements is due to the difference in the relative positioning of the menisci resulting from the spatial arrangement of the capillaries.

For instance, the leading meniscus for an orderly arranged interacting capillary system is in the smallest radius capillary and we know that the fraction of saturation contributed by the smallest radius capillary is small. For the arrangement 2 shown in Fig. 9, the leading meniscus is in the capillary of radius 100 m. In the capillary bundle model, the cross-section area of the leading capillary (200m) is 13.93% of the total cross-section area, whereas for the ordered arrangement 2 and the arrangement number 2, the respective cross-section area of the leading meniscus capillaries are 0.03% and 3.43% of the classic capillary bundle model, while in the case of interacting capillaries a steep decrease is observed already at small longitudinal positions. Fig. 11(b) also shows that the standard deviation in saturation from the average across the 1000 arrangements at $T = 0.1$ and $T = 0.2$, which is as high as 0.2 at $Z = 0.59$ and 0.69, respectively; whereas for $T = T_{bt}$, it is 0.18 at $Z = 0.76$. In real two-dimensional porous media, where the spatial arrangement of pores may vary, the interacting capillaries model will be more helpful in predicting the curvature of capillaries model. The saturation of the porous medium with length and the breakthrough time significantly differ for the classic (non-interacting) capillary bundle and for the different arrangements of the interacting multi-capillary system, although the contrast in the radii of the capillaries is the same.

## 2. Confronting predictions from the model to experimental measurements from previous studies

The spatial profiles of saturation for the interacting multi-capillary system are consistent with observations of imbibition phenomena in quasi-2D porous media described by Dong et al., Ding et al., Debbabi et al., and Akbari et al.,$^{48,66-68}$. In real porous media, the imbibing fluid saturation decreases gradually with longitudinal position, similarly to the trend shown by the interacting multi-capillary system. It was also previously described that the lagging macroscopic front is mostly responsible for the saturation of a porous medium$^{47}$, which is in good agreement with the saturation profile anticipated by the interacting multi-capillary system, as shown in Fig. 11(b). The saturation profile for the (classic) non-interacting capillary bundle (Fig. 11a) predicts that the large pores are responsible for the leading macroscopic front and the saturation of the porous medium, which is contrary to the interacting capillaries model (shown in Fig. 11(a)) and the experimental observations in real porous media.$^{47,54,62,69,70}$

Furthermore, in the following we compare the predictions of our model to two data sets from the literature, both taken from Ref. $^{47}$.

### a. Two capillary system:
We first compare our model predictions to measurements performed on a system of two capillaries consisting of a thread positioned inside a cylindrical tube. The time evolution of the menisci position squared, as predicted by our model, compares well with the experimental observations for both capillaries (Fig. 12). The radius of the large capillary was $r_α = 300$ m, that of the thread $r_β = 170$ m. From the experimental data$^{47}$, the value of $(Pc_αr_α^2)/(8μL^2)$ is 0.0108 s$^{-1}$, which is used to non-dimensionalize time in Fig. 12. The predictions from the classic (non-interacting) capillary bundle model (Eq. 15) are also shown in the inset of Fig. 12 for comparison. The imbibition in the wider capillary is little impacted by the imbibition in the (much) narrower capillary, so that the prediction of the non-interacting capillary bundle for the wider capillary are similar
692 to the experimental data; however the non-interacting capillary bundle underestimates the advancement of the meniscus in the narrower capillary (the thread) by a factor 5. 700

693 in the narrower capillary (the thread) by a factor 5. 701

694 702 Zα exp; Ref. 47 Zβ exp; Ref. 47 Zα sim; n = 2 Zβ sim; n = 2 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718

<caption>FIG. 12. Imbibition in a system of two interacting capillaries having radii $r_α = 300$ m and $r_β = 170$ m. Predictions from our semi-analytical model (solid lines) compare well to the data (symbols) Bico and Quéré$^{47}$. The inset of the figure shows the same comparison for predictions of the classic (non-interacting) capillary bundle model (solid lines), obtained through Eq. 15, which underestimates the advancement $Z_β$ of the meniscus in the narrower capillary (red line) by a factor 5. 719 720 721 722 723 724 725</caption>

695 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754

<caption>FIG. 13. Dependence of the square of the non-dimensional imbibition length on non-dimensional time. The experimental findings of Bico and Quéré$^{47}$ are shown with red squares (for leading front) and black triangles (for lagging front). The predictions of our model for two different samplings (6 and 12 interacting capillaries) of the uniform pore size distribution are shown with lines, respectively orange and dashed (for $n = 6$) or purple and solid (for $n = 11$) for the leading front, and thick, gray and dashed (for $n = 6$) or black and dashed for the lagging front. The results from the classic, non-interacting capillary bundle are presented for comparison for the leading front (green long-dashed line) and lagging fronts (green dashed line). 745 746 747 748 749 750 751 752 753 754</caption>

697 b. Imbibition in a paper filter: Bico and Quéré$^{47}$ also performed experiments in which a silicone oil of viscosity 755

698 one can assume that the Darcy law holds at any time across

16 · 10⁻³ Pa·s and surface tension 20.6 · 10⁻³ N/m spontaneously imbibes into a Whatman grade 4 filter paper, which has pore diameters in the range 20 to 25 m. They observed that the microscopic front propagating in small pores travels ahead of the saturating macroscopic front in large pores, again in contradiction to the predictions of the classic non-interacting capillary bundle model. In Fig. 13, we show a comparison of the experimental observations from these authors$^{47}$ (shown as symbols in the figure) with predictions of our model (shown as lines in the figure). Two capillary systems were simulated with our model, corresponding to two ways of sampling the pore size PDF (probability density function) of the paper filter: having no information on the functional form of that PDF, we assumed that it was uniform and sampled it first with $n = 6$ interacting capillaries of radii 10, 10.5, 11, 11.5, 12, and 12.5 m; we then performed a second calculation with a sampling twice finer, i.e., with $n = 11$ interacting capillaries of radii 10, 10.25, 10.5, 10.75, 11, 11.25, 11.5, 11.75, 12, 12.25 and 12.5 m. For $n = 6$ the non-dimensional leading front position was defined as the average of the positions of the two more advanced menisci, whereas that of the lagging front was defined as the average of the two less advanced menisci. For $n = 11$, a similar method was used, but involving the average of the 3 more advanced menisci positions for the leading front and that of the 3 less advanced menisci positions for the lagging front. A statistics of $R = 360$ arrangements (i.e., all possible arrangements) was chosen for $n = 6$, whereas for $n = 11$ we used $R = 1000$ randomly-chosen arrangements within more than 19.9 millions of different possible arrangements. The confidence interval defined from the standard deviations over the statistics is also shown in Fig. 13 as thin orange lines for the leading front computed with $n = 6$; for the lagging front the standard deviations are so small that they would be hardly visible, so we did not plot the corresponding confidence interval.

The predictions of our model for $n = 6$ and $n = 11$ are very similar to each other, especially for the leading front, which is a good test of consistency for the method. Indeed, it means that changing the sampling resolution for a given pore size distribution does not impact the predictions. Furthermore, these predictions appear to be quite consistent with the experimental data, for both the leading and lagging front. In other words, they exhibit the same Washburn-like dynamics as both the experimental leading front (at all times) and lagging front (for $T \le 0.3$ at least), with the same proportionality factors between $Z^2$ and $T$ (i.e., the slope in the plots). On the contrary, the predictions of the classic (non-interacting) capillary bundle, also shown in Fig 13 (as green dashed lines) are shown to be much less efficient at predicting the proportionality factor, especially the leading front; in addition they predict a leading front occupying the largest capillaries and a lagging front occupying the smallest ones, in contradiction to the experimental observations and to the predictions from our model.

Note that to non-dimensionalize the time in Fig. 13 we have relied on the observation by Bico and Quéré that most of the wetting fluid is carried by the lagging front (which they term macroscopic front). Adopting a macroscopic point of view, one can assume that the Darcy law holds at any time across
the porous medium's length, with a pressure gradient that is 809 $Pc_{eff}/z$, $Pc_{eff}$ being a constant effective capillary pressure defined for the entire medium. Then the Darcy law reads 810

$$ \frac{dz}{dt} = \frac{K}{\mu} \frac{Pc_{eff}}{z}, \quad \text{leading to} \quad z^2 = \frac{2Pc_{eff}K}{\mu} t, \quad (16) \quad 811 $$

where $K$ is the medium's permeability and we have assumed that at time $t=0$ no wetting fluid has yet invaded the medium. 810 If we choose to non-dimensionalize time by the characteristic time $(\mu L^2)/(Pc_{eff}K)$, we obtain from Eq. (16) the non-dimensional equation $Z^2 = 2T$. Since, according to Bico and Quéré's observation mentioned above, it is the lagging (macroscopic) front that carries most of the interface between the two fluids, Eq. 16, and therefore its non-dimensional counterpart, can be assumed to describe the behavior of the lagging front. From the experimental data for the lagging front, $(Pc_{eff}K)/(\mu L^2)$ is measured to be $9.7 \cdot 10^{-5} \text{ s}^{-1}$, which we thus use to non-dimensionalize all plots in Fig. 13. The dependence of $Z^2$ on $T$ for the lagging (macroscopic) front then has a slope 2 (as shown by the dotted gray line in Fig. 13), while that for the leading (microscopic) front exhibits a larger imbibition rate, with a slope 2.67 (as shown by the orange dotted line in Fig. 13). 820 821 822 823 824 825 826 827 828 829 830 831

## IV. CONCLUSIONS

In conclusion, we investigated spontaneous imbibition of a 836 wetting fluid in a randomly arranged planar system of interacting capillaries. This generalized model can predict the imbibition behavior for all the $n!/2$ possible arrangements of an 837 interacting $n$-capillary system. It is inspired from a previous 840 work on stratified geological formations, with planar layers 841 instead of cylindrical capillaries. 842

Using an interacting capillary system containing 4 capillaries, we showed that the imbibition dynamics depends significantly on the arrangement of the capillaries within the capillary system, for a given distribution of the capillary radii. Similarly, the dynamics are affected by that distribution for a given 843 844 845 846 847 arrangement of the capillaries. Furthermore, we showed that 848 the arrangement and radii distribution of the capillaries jointly 849 control the relative menisci's locations, the breakthrough time, 850 851 and which capillary carries the leading meniscus. The cross-sectional saturation of the impregnating fluid along the length 852 of the capillary system also changes with a change in the arrangement of the capillaries. However, the total capillary pressure driving the flow is identical for all arrangements, there 853 854 855 856 fore, the overall volume fraction occupied by the invading fluid (i.e., the global saturation of the wetting fluid) at a given 857 time remains the same across all arrangements, until breakthrough occurs in one of the arrangements. 858 859

Similarly, considering 1000 randomly-chosen different arrangements of an interacting twenty-capillary system having 860 uniform distribution of radii between 10 m and 200 m, we observed that, depending on the arrangement of the capillaries, 861 862 863 the leading meniscus can be in any of the capillaries whose radii are between 10 m and 100 m, and the non-dimensional breakthrough time lies between $T_{bt} = 0.31$ and $T_{bt} = 0.42$. 864 865 866

The dynamics of spontaneous imbibition as predicted by this new model is significantly different from that predicted by the classic bundle of non-interactive capillaries (or tubes), for which the leading meniscus is always in the largest radius capillary. For the interacting multi-capillary system mentioned above, on the contrary, the leading meniscus can be in any of the capillaries having radii between 10 m and 100 m. We observed that the breakthrough occurs earlier than in the classic capillary bundle, where it occurs at non-dimensional time $T_{bt} = 0.5$ for the aforementioned 20-capillary-system, to be compared to the 0.31–0.42 range for the 20-capillary-system mentioned above. Furthermore, for this system the saturation at breakthrough time falls in the range 0.69–0.79, whereas for the classic capillary bundle it is equal to 0.86. The dependence of the saturation as a function of the longitudinal position are also shows a stark contrast between the predictions of the classic capillary bundle and the average behavior of the 1000 arrangements of interacting capillaries. Indeed, the interacting capillary system shows a steep decrease in the saturation with length as compared to the classic capillary bundle. Additionally, the interacting multi-capillary system shows that the spatial arrangement of the capillaries may cause significantly different saturation values at a given longitudinal position.

So, how is this model consisting of a planar bundle of interacting capillaries to be used to predict spontaneous imbibition in quasi-two-dimensional porous media whose pore size distribution is known? We propose to use a stochastic approach, i.e., to consider the average behavior between a large number of randomly-picked spatial arrangements of the capillary diameters, the distribution of these diameters being equal to the pore size distribution of the real porous medium. We tested that method against data from the literature. Firstly, qualitative observations relative to which ranges of pore sizes mainly contribute to the leading and lagging fronts of the imbibition interface, and to the longitudinal saturation profile, are consistent between experiments from the literature and the predictions of our model. Secondly, to validate the model's quantitative predictive capacity, we compared its predictions to imbibition measurements in filter paper, performed by Bico and Quéré$^{47}$. The model predicts that the visible leading front is carried by smaller pores and that the bulk saturating front responsible for most of the fluid mass invasion is the lagging front carried by larger pores, which agrees very well with the experimental findings. The quantitative predictions for the positions in time of these two fronts, obtained from averaging over the statistics of randomly-chosen arrangements, agree well with the measurements.

This generalized model for spontaneous imbibition in a planar bundle of interacting capillaries, which is semi-analytical and runs extremely quickly, could be useful for fast assessment of one-dimensional imbibition dynamics in design-based porous media such as loop heat pipes, diagnostic devices and microfluidic devices, or in real porous media whose porosity structure can reasonably be considered two-dimensional, e.g., paper, thin porous media in general, or layered aquifers.

Prospects to this work include extending this approach to three-dimensional models by considering parallel capillaries,
867 the positions of whose axes in a transverse plane would be the nodes of a triangular grid.

868 900

901

902 CONFLICTS OF INTEREST

903

904

905

906

907 Defining $\omega$ and $(1-\omega)$ as the fractions of $dq_\alpha$ transferred

908 respectively to $\beta$ and $\gamma$, we can write an equation similar to

909 Eq. (A4) for both $\beta$ and $\alpha$, where $\omega dq_\alpha$ and $(1-\omega) dq_\alpha$ appear respectively as a differential velocity term arising from

910 fluid transfer. Considering that the pressure field is the same

911 SA and YM acknowledge the financial support granted by the French Agence Nationale de la Recherche, under project CO2-3D with project number ANR-16-CE06-0001.

912

913 CO2-3D with project number ANR-16-CE06-0001.

914

915

## 875 Appendix A: Mathematical formulation for the system of four interacting capillaires

877 In capillary $\alpha$, for $0 < z < z_\alpha(t)$, the pressure drop is given by the Hagen-Poiseuille law as,

878 916

$$P(z_{\alpha}(t), t) - P_0 = -\frac{8\mu z_{\alpha}(t)}{r_{\alpha}^2} v_{\alpha}(t), \quad (A1)^{916}$$

879 where $\mu$ is the imbibing fluid's viscosity, $v_\alpha(t)$ is the instantaneous velocity of the wetting fluid in the capillary $\alpha$, $P_0$ is the inlet pressure and $P(z_\alpha(t), t)$ is the pressure in the imbibing fluid at $z_\alpha(t)$, as shown in Fig. 3. Since the pressure fields are identical in all capillaries for $z < z_\alpha(t)$, the pressure gradient is the same in all capillaries, which from Eq. (A1) implies

880 917

$$\frac{v_{\alpha}(t)}{r_{\alpha}^2} = \frac{v_{\beta}(t)}{r_{\beta}^2} = \frac{v_{\gamma}(t)}{r_{\gamma}^2} = \frac{v_{\delta}(t)}{r_{\delta}^2}, \quad (A2)^{918}$$

881 920

882 921

883 where the index $i$ ($i = \alpha, \beta, \gamma, \delta$) indicates quantities relative to the capillary of radius $r_i$ and $v_i(t)$ ($i = \alpha, \beta, \gamma, \delta$) is the velocity of the imbibing fluid for $z < z_\alpha(t)$.

884 922 The capillary pressure jump through the fluid-fluid interface is $Pc_\alpha$ at $z_\alpha(t)$, where some of the imbibing fluid transfers from the capillary $\alpha$ to other capillaries. The volumetric fluid transfer from the capillary $\alpha$ to the capillaries $\beta$ and $\gamma$ is $dq_\alpha$, whereas the fluid transfer from the capillary $\alpha$ to the capillary $\delta$ is $dq'_\alpha$. The velocity of the advancing meniscus in capillary $\alpha$, $dz_\alpha/dt$, is thus given by

885 923

$$\frac{dz_{\alpha}}{dt} = v_{\alpha}(t) - \frac{dq_{\alpha} + dq'_{\alpha}}{\pi r_{\alpha}^2}. \quad (A3)$$

886 924 For $z_\alpha(t) < z < z_\delta(t)$, the velocity of the fluid in capillary $\delta$ is similarly given by

887 925

$$\frac{dz_{\delta}}{dt} = v_{\delta}(t) + \frac{dq'_{\alpha}}{\pi r_{\delta}^2}, \quad (A4)$$

888 926 so the pressure drop in the capillary $\delta$ between $z = z_\alpha(t)$ and $z = z_\delta(t)$ is

889 927

$$P(z_{\delta}(t), t) - P(z_{\alpha}(t), t) = -\frac{8\mu(z_{\delta}(t) - z_{\alpha}(t))}{r_{\delta}^2} \left( v_{\delta}(t) + \frac{dq'_{\alpha}}{\pi r_{\delta}^2} \right). \quad (A5)$$

890 928

891 929

892 930

893 931

894 932

At $z = z_\delta(t)$, the pressure jump across the meniscus is $Pc_\delta$, since the pressure in the non-wetting fluid is the atmospheric pressure.

The capillaries $\beta$ and $\gamma$ are on the other side of the capillary $\alpha$ with respect to the capillary $\delta$. As the capillary pressure jump of the capillary $\beta$ is smaller than that in the capillary $\gamma$, the meniscus in $\beta$ lags behind that in $\gamma$. Hence, the imbibing fluid in these capillaries is continuous for $z_\alpha(t) < z < z_\beta(t)$. Defining $\omega$ and $(1-\omega)$ as the fractions of $dq_\alpha$ transferred respectively to $\beta$ and $\gamma$, we can write an equation similar to Eq. (A4) for both $\beta$ and $\alpha$, where $\omega dq_\alpha$ and $(1-\omega) dq_\alpha$ appear respectively as a differential velocity term arising from fluid transfer. Considering that the pressure field is the same in the capillaries $\beta$ and $\gamma$ for $z_\alpha(t) < z < z_\beta(t)$, we then obtain in that $z$ range:

$$\frac{v_{\beta}(t) + \frac{\omega dq_{\alpha}}{\pi r_{\beta}^2}}{\pi r_{\beta}^2} = \frac{v_{\gamma}(t) + \frac{(1-\omega) dq_{\alpha}}{\pi r_{\gamma}^2}}{\pi r_{\gamma}^2}. \quad (A6)$$

Combining Eq. (A2) and Eq. (A6), we then obtain the fraction $\omega$ from the capillaries' radii: $\omega = r_{\beta}^4/(r_{\beta}^4 + r_{\gamma}^4)$. Therefore, the pressure drop in capillaries $\beta$ and $\gamma$ for $z_\alpha(t) < z < z_\beta(t)$ is

$$\begin{aligned} P(z_{\beta}(t), t) - P(z_{\alpha}(t), t) \\ = -\frac{8\mu(z_{\beta}(t) - z_{\alpha}(t))}{r_{\beta}^2} \left( v_{\beta}(t) + \omega \frac{dq_{\alpha}}{A_{\beta}} \right). \end{aligned} \quad (A7)$$

At the meniscus in the capillary $\beta$, the capillary pressure jump is $Pc_\beta$ and some of the impregnating fluid transfers from $\beta$ to $\gamma$, which we assume to correspond to a differential flow rate $dq_\beta$. The velocity of the meniscus in the capillary $\beta$ for $z > z_\beta(t)$ is then

$$\frac{dz_{\beta}}{dt} = v_{\beta}(t) + \omega \frac{dq_{\alpha}}{\pi r_{\beta}^2} - \frac{dq_{\beta}}{\pi r_{\beta}^2}. \quad (A8)$$

Similarly, for $z > z_\beta(t)$, the meniscus in the capillary $\gamma$ travels with a velocity given by

$$\frac{dz_{\gamma}}{dt} = v_{\gamma}(t) + (1-\omega) \frac{dq_{\alpha}}{\pi r_{\gamma}^2} + \frac{dq_{\beta}}{\pi r_{\gamma}^2}. \quad (A9)$$

The pressure drop between $z = z_\beta(t)$ and $z = z_\gamma(t)$ in capillary $\gamma$ is then given by,

$$P(z_{\gamma}(t), t) - P(z_{\beta}(t), t) = -\frac{8\mu(z_{\gamma}(t) - z_{\beta}(t))}{r_{\gamma}^2} \left( v_{\gamma}(t) + (1-\omega) \frac{dq_{\alpha}}{\pi r_{\gamma}^2} + \frac{dq_{\beta}}{\pi r_{\gamma}^2} \right). \quad (A10)$$

The pressure jump across the meniscus in each of the capillaries is given by the Young-Laplace equation$^{64,65}$, i.e., Eq. (2), from which it follows that

$$P(z_i, t) - P_0 = -Pc_i = -\frac{2\sigma \cos \theta_w}{r_i}, \quad (A11)$$
for $i = \alpha, \beta, \gamma, \delta$. Note that the prefactor 2 is controlled by circular cross-section of the tube, another geometry (e.g., square) cross section) would yield a different prefactor. Eq. (A11) imposes the total pressure drop within the impregnating wetting fluid in each of the capillaries. Substituting Eqs. (A3), (A4), (A8), (A9) in Eqs. (A1), (A5), (A7), (A10) respectively, we obtain the equations governing the flow in the interacting capillary system:

$$Pc_{\alpha} = \frac{8\mu z_{\alpha}(t)}{r_{\alpha}^{4} + r_{\beta}^{4} + r_{\gamma}^{4} + r_{\delta}^{4}} \left( r_{\alpha}^{2} \frac{dz_{\alpha}}{dt} + r_{\beta}^{2} \frac{dz_{\beta}}{dt} + r_{\gamma}^{2} \frac{dz_{\gamma}}{dt} + r_{\delta}^{2} \frac{dz_{\delta}}{dt} \right), \quad (A12) \quad (A13)$$

$$Pc_{\delta} - Pc_{\alpha} = \frac{8\mu(z_{\delta}(t) - z_{\alpha}(t))}{r_{\delta}^{2}} \left( \frac{dz_{\delta}}{dt} \right), \quad (A13)$$

$$Pc_{\beta} - Pc_{\alpha} = \frac{8\mu(z_{\beta}(t) - z_{\alpha}(t))}{r_{\beta}^{4} + r_{\gamma}^{4}} \left( r_{\beta}^{2} \frac{dz_{\beta}}{dt} + r_{\gamma}^{2} \frac{dz_{\gamma}}{dt} \right), \quad (A14)$$

$$Pc_{\gamma} - Pc_{\beta} = \frac{8\mu(z_{\gamma}(t) - z_{\beta}(t))}{r_{\gamma}^{2}} \left( \frac{dz_{\gamma}}{dt} \right). \quad (A15)$$

## Appendix B: Generalization of the model for an arbitrary number of capillaries

The following step-by-step procedure must be followed:

1. We initiate the model formulation by finding the largest radius capillary, $C_i$. The pressure field is identical in all capillaries for $z < z_a(t)$, and the corresponding pressure gradient is related to the fluid velocity in each capillary by Hagen-Poiseuille's law. Some of the invading fluid from capillary $i$ transfers to other capillaries in the immediate vicinity of the meniscus position $z_a(t)$.

2. For $z > z_a(t)$, the imbibing fluid in the capillaries $C_1$ to $C_{i-1}$ is separated from the imbibing fluid in the capillaries $C_{i+1}$ to $C_n$. We thus classify the capillaries on either sides of the capillary $C_i$ in two regions, the capillaries $C_1$ to $C_{i-1}$ in the first one, the capillaries from $C_{i+1}$ to $C_n$ in another one. The fluid transfer from the capillary $C_i$ is divided among the other capillaries according to their radii. If the fluid transfer to the 'top' region' is $dq_t$, the fraction of $dq_t$ flowing from capillary $C_i$ to a capillary of radius $r_p$ would be $r_p^4 dq_t / \sum_{q=1}^{i-1} (r_q^4)$. Similarly, for the 'bottom region', if $dq_b$ is the fluid transfer from $C_i$, the fractional flow in a capillary of radius $r_r$ will be $r_r^4 dq_b / \sum_{s=i+1}^{n} (r_s^4)$. This fluid transfer causes the flow rates to increase in capillaries $C_1$ to $C_{i-1}$ and $C_{i+1}$ to $C_n$.

3. The widest capillary among the capillaries $C_1$ to $C_{i-1}$ is now identified. For $z_a(t) < z < z_b(t)$ the pressure field in the imbibing fluid is identical in capillaries $C_1$ to $C_{i-1}$, and is related to the fluid velocity in each capillary by Hagen-Poiseuille's law. In the vicinity of

$z = z_b(t)$, some of the invading fluid transfers from $C_j$ to the capillaries $C_1$ to $C_{j-1}$ and $C_{j+1}$ to $C_{i-1}$, which increases the flow rate in these capillaries.

4. Similarly, the widest capillary among capillary $C_{i+1}$ to $C_n$, which we denote $C_k$, is chosen. The pressure field is identical in the capillaries $C_{i+1}$ to $C_n$ for $z_a(t) < z < z_c(t)$, and the pressure gradient is related to the fluid velocity in each of these capillaries from the Hagen-Poiseuille law. A $z = z_c(t)$, some of the fluid invading $C_k$ transfers into the capillaries $C_{i+1}$ to $C_{k-1}$ and $C_{k+1}$ to $C_n$, which increases the flow rate in these capillaries.

5. The impregnating fluids in the regions encompassing capillaries $C_1$ to $C_{j-1}$ and $C_{j+1}$ to $C_{i-1}$ are separated by displaced fluid in capillary $C_j$ for $z > z_j$. Again, the capillary of largest radius among the capillaries $C_1$ to $C_{j-1}$ is identified, as well as the capillary of largest radius among the capillaries $C_{j+1}$ to $C_{i-1}$. The similar procedure previously explained for the pressure field and its relation to the fluid velocity is repeated for those two regions.

6. The same procedure as explained in step 5. is performed in the regions encompassing capillaries $C_{i+1}$ to $C_{k-1}$ and $C_{k+1}$ to $C_n$.

7. This is repeated in all the regions which have been defined in steps 1 to 5, and this in a recursive manner, until the entire bundle of interacting capillaries is divided into regions containing only one capillary each.

8. The pressure jump across the meniscus in each of the capillaries is the corresponding Young-Laplace capillary pressure of that capillary. The $n$ equations relating the pressure drops to the velocities of the fluid-fluid interfaces are then solved to obtain the lengths impregnated in each of the capillaries at the considered time $t$.

## Appendix C: Imbibition in all possible arrangements of a system of four interacting capillaries

A four capillary system has 12 possible arrangements. For a set of capillaries with radii $r_\alpha = 80 \text{ m}$, $r_\beta = 60 \text{ m}$, $r_\gamma = 40 \text{ m}$ and $r_\delta = 20 \text{ m}$, we present in Fig. C.1 the time evolution of the meniscus' positions in all four capillaries for all 12 arrangements.

We see from Fig. C.1 that the leading meniscus is in capillary $\delta$ for arrangements shown in Fig. C.1(a),(b),(f),(g),(i),(j),(k),(l). For the arrangements shown in Fig. C.1(c),(d), the leading meniscus is in $\gamma$. For arrangements shown in Fig. C.1(e),(h), the capillaries $\gamma$ and $\delta$ impregnate the same distance with time. But the breakthrough times are different for all the arrangements, varying from $T = 0.33$ to $T = 0.40$. The minimal breakthrough time is 0.33, observed in arrangements (a), (g), (k) and (l) of Fig. C.1.
<figure><img src="image_22.png" /><figcaption>FIG. C.1. Spontaneous imbibition in a system of four interacting capillaries of radii $r_\alpha = 80$ m, $r_\beta = 60$ m, $r_\gamma = 40$ m and $r_\delta = 20$ m. The non-dimensional positions of the four menisci are shown as a function of non-dimensional time for all the 12 possible arrangements in (a) to (l). The arrangement, the ordering of the menisci locations, and the breakthrough times for each of the cases (a) to (l) are provided as legends of the plots.</figcaption></figure>

1020

1021

The breakthrough for all the arrangements shown in Fig. C.1 occurs between 225.7 s and 273.6 s.

1032

1033

1034

1022

1035

## Appendix D: Convergence of the computations for a system of 20 interacting capillaries

1036

1023

1038

1039

1024

For the study of the bundle consisting of 20 interacting capillaries, the convergence of the results as a function of the number of randomly-chosen spatial arrangements was verified in the following manner.

1025

1040

1026

1041

1027

1042

1028

1043

Three sets of $R = 100, 1000$ and $2000$ randomly-chosen arrangements were simulated independently, and their results were compared with each other. Fig. D.1(a) shows the spatial profile of wetting phase saturation at three different times ($T = 0.2, T = 0.3, T = T_{bt}$), obtained as the average of the spatial profiles for all $R$ arrangements. Fig. D.1(b) shows the standard deviation over the statistics of the spatial wetting phase saturation profiles for the $R$ arrangement, also at times $T = 0.2, T = 0.3, T = T_{bt}$. Obviously the average behavior for 1000 arrangements (in contrast to the case $R = 100$) cannot be distinguished from that for 2000 arrangements, and even the spatial profiles of the standard deviation over the statistics are quite similar for the two cases. Therefore, we consider $R = 1000$ to be a sufficiently large number of randomly-chosen arrangements for the imbibition dynamics to be well predicted in a system of 20 interacting capillaries.

1029

1044

1030

<footnote n="1">¹B. Xiao, J. Fan, and F. Ding, “Prediction of relative permeability of unsaturated flow in a capillary network,” in *J. Petrol. Sci. Technol.* **7**(2) (June 1982).</footnote>
<figure><img src="image_23.png" /><figcaption>FIG. D.1. Convergence of the simulations for a system of 20 interacting capillaries, based on $R = 100, 1000$ and $2000$ arrangements at $T = 0.2$, $T = 0.3$, and $T = T_{bt}$ (breakthrough time): (a) Mean saturation as a function of the longitudinal coordinate. (b) Standard deviation (SD) of the statistics as a function of the longitudinal coordinate.</figcaption></figure>

urated porous media based on fractal theory and monte carlo simulation 1091
Energy & fuels 26, 6971-6978 (2012). 1092

2Y.-J. Lin, P. He, M. Tavakkoli, N. T. Mathew, Y. Y. Fatt, J. C. Chai, A. Gu 1093
harzadeh, F. M. Vargas, and S. L. Biswal, "Characterizing asphaltene dep 1094
sition in the presence of chemical dispersants in porous media micromod 1095
els," Energy & fuels 31, 11660-11668 (2017). 1096

3S. Saraji, L. Goual, and M. Piri, "Adsorption of asphaltenes in porou 1097
media under flow conditions," Energy & fuels 24, 6009-6017 (2010). 1098

4M. Taghizadeh-Behbahani, B. Hemmateenejad, M. Shamsipur, and 1099
A. Tavassoli, "A paper-based length of stain analytical device for naked cy 1100
(readout-free) detection of cystic fibrosis," Analytica Chimica Acta (2019) 1101
5Y. Soda, D. Citterio, and E. Bakker, "Equipment-free detection of k+ on mi 1102
crofluidic paper-based analytical devices based on exhaustive replacement 1103
with ionic dye in ion-selective capillary sensors," ACS sensors 4, 670-675 1104
(2019). 1105

6B. Dai, K. Li, L. Shi, X. Wan, X. Liu, F. Zhang, L. Jiang, and S. Wang 1106
"Bioinspired janus textile with conical micropores for human body moi 1107
ture and thermal management," Advanced Materials (2019). 1108

7M. Rosello, S. Sur, B. Barbet, and J. P. Rothstein, "Dripping-onto-substrat 1110
capillary breakup extensional rheometry of low-viscosity printing inks 1111
Journal of Non-Newtonian Fluid Mechanics 266, 160-170 (2019). 1111

8Y. Wang, R. Deng, L. Yang, and C. D. Bain, "Fabrication of monolayers on 1112
uniform polymeric particles by inkjet printing of monodisperse emulsions 1113
produced by microfluidics," Lab on a Chip 19, 3077-3085 (2019). 1114

9Y. Liu, J. Kaszuba, and J. Oakey, "Microfluidic investigations of crude oil 1116
brine interface elasticity modifications via brine chemistry to enhance 1116
recovery," Fuel 239, 338-346 (2019). 1117

10R. Gharibshahi, M. Omidkhah, A. Jafari, and Z. Fakhroueian, "Hybridiza 1118
tion of superparamagnetic fe3o4 nanoparticles with mwents and effect on 1119
surface modification on electromagnetic heating process efficiency: A mi 1120
crofluidics enhanced oil recovery study," Fuel 282, 118603 (2020). 1121

11C. Carrell, A. Kava, M. Nguyen, R. Menger, Z. Munshi, Z. Call, M. Nus 1122
baum, and C. Henry, "Beyond the lateral flow assay: A review of pap er 1123
based microfluidics," Microelectronic Engineering 206, 45-54 (2019). 1124

12F. Schaumburg and C. L. Berli, "Assessing the rapid flow in multilayer 1126
paper-based microfluidic devices," Microfluidics and Nanofluidics 23, 98-102 1127
(2019). 1127

13M. Rich, O. Mohd, F. S. Ligler, and G. M. Walker, "Characterization 1130
of glass frit capillary pumps for microfluidic devices," Microfluidics and 1131
Nanofluidics 23, 70 (2019). 1130

14J.-H. Lin, W.-H. Chen, Y.-J. Su, and T.-H. Ko, "Performance analysis of 1133
a proton-exchange membrane fuel cell (pemfc) with various hydrophobic 1133
agents in a gas diffusion layer," Energy & fuels 22, 1200-1203 (2008). 1133

15K. K. Lee, M.-O. Kim, and S. Choi, "A whole blood sample-to-answ er 1134
polymer lab-on-a-chip with superhydrophilic surface toward point-of-c 1135

technology," Journal of pharmaceutical and biomedical analysis 162, 28-33 162
(2019).

16C. Liang, Y. Liu, A. Niu, C. Liu, J. Li, and D. Ning, "Smartphone-app 169
based point-of-care testing for myocardial infarction biomarker ctni using 169
an autonomous capillary microfluidic chip with self-aligned on-chip focus- 169
ing (sof) lenses," Lab on a Chip 19, 1797-1807 (2019).

17H.-A. Joung, Z. S. Ballard, A. Ma, D. K. Tseng, H. Teshome, S. Bu- 171
rakowski, O. B. Garner, D. Di Carlo, and A. Ozcan, "based multiplexed 171
vertical flow assay for point-of-care testing," Lab on a Chip 19, 1027-1034 171
(2019).

18B. Xiao, W. Wang, X. Zhang, G. Long, H. Chen, H. Cai, and L. Deng, "A 184
novel fractal model for relative permeability of gas diffusion layer in proton 184
exchange membrane fuel cell with capillary pressure effect," Fractals 27, 1950012 184
(2019).

19P. Carrere and M. Prat, "Liquid water in cathode gas diffusion layers of pem 196
fuel cells: Identification of various pore filling regimes from pore network 196
simulations," International Journal of Heat and Mass Transfer 129, 1043- 196
1056 (2019).

20M. Singh, N. V. Datla, S. Kondaraju, and S. S. Bahga, "Enhanced ther- 201
mal performance of micro heat pipes through optimization of wettability 201
gradient," Applied Thermal Engineering 143, 350-357 (2018).

21M. Chernysheva and Y. Maydanik, "Simulation of heat and mass transfer in 204
a cylindrical evaporator of a loop heat pipe," International Journal of Heat 204
and Mass Transfer 131, 442-449 (2019).

22C. Pozrikidis, "Axisymmetric motion of a file of red blood cells through 208
capillaries," Physics of fluids 17, 031503 (2005).

23K. Singh, B. P. Muljadi, A. Q. Raeini, C. Jost, V. Vandeginste, M. J. Blunt, 205
G. Theraulaz, and P. Degond, "The architectural design of smart ventila- 205
tion and drainage systems in termite nests," Science advances 5, eaat8520 205
(2019).

24K. Li, D. Zhang, H. Bian, C. Meng, and Y. Yang, "Criteria for applying the 209
lucas-washburn law," Scientific reports 5, 14085 (2015).

25S. Gruener and P. Huber, "Capillarity-driven oil flow in nanopores: Darcy 212
scale analysis of lucas-washburn imbibition dynamics," Transport in Porous 212
Media 126, 599-614 (2019).

26J. Cai, Y. Chen, Y. Liu, S. Li, and C. Sun, "Capillary imbibition and flow 215
of wetting liquid in irregular capillaries: A 100-year review," Advances in 215
Colloid and Interface Science, 102654 (2022).

27R. Lucas, "Ueber das zeitgesetz des kapillaren aufstiegs von flüssigkeiten," 219
Colloid & Polymer Science 23, 15-22 (1918).

28E. Washburn, "The dynamics of capillary flow," Physical Review 17, 273 219
(1921).

29R. Lenormand, C. Zarcone, et al., "Role of roughness and edges during 222
imbibition in square capillaries," in *SPE annual technical conference and 222
exhibition* (Society of Petroleum Engineers, 1984).
1136 $^{30}$M. Dong and I. Chatzis, "The imbibition and flow of a wetting liquid along the corners of a square capillary tube," Journal of colloid and interface science 172, 278-288 (1995).

1137 31 M. Ramezanzadeh, S. Khasi, and M. H. Ghazanfari, "Simulating imbibition in capillary bundle model with corner flow: A study of capillary morphology," Journal of Petroleum Science and Engineering 176, 62-73 (2019).

1140 32 D. Zheng, W. Wang, and Z. Reza, "Integrated pore-scale characterization of mercury injection/imbibition and isothermal adsorption/desorption in porous media using dendroidal model for shales," Journal of Petroleum Science and Engineering 178, 751-765 (2019).

1141 33 M. Reyssat, L. Courbin, E. Reyssat, and H. A. Stone, "Imbibition in porous liquids with axial variations," Journal of Fluid Mechanics 615, 335-344 (2008).

1145 34 A. Budaraju, J. Phirani, S. Kondaraju, and S. S. Bahga, "Capillary displacement of viscous liquids in geometries with axial variations," Langmuir 34, 10513-10521 (2016).

1154 35 F. F. Ouali, G. McHale, H. Javed, C. Trabi, N. J. Shirtcliffe, and M. I. Newton, "Wetting considerations in capillary rise and imbibition in closed square tubes and open rectangular cross-section channels," Microfluidics and Nanofluidics 15, 309-326 (2013).

1156 36 U. Rosendahl, A. Grah, and M. E. Dreyer, "Convective dominated flows in open capillary channels," Physics of Fluids 22, 052102 (2010).

1159 37 M. M. Weislogel, "Capillary flow in interior corners: The infinite column," Physics of Fluids 13, 3101-3107 (2001).

1161 38 D. Dimitrov, L. Klushin, A. Milchev, and K. Binder, "Flow and transport in brush-coated capillaries: A molecular dynamics simulation," Physics of Fluids 20, 092102 (2008).

1164 39 J. Wang, A. Salama, and J. Kou, "Experimental and numerical analysis of capillary imbibition processes in a corrugated capillary tube," Capillarity 5, 83-90 (2022).

1167 40 A. Salama, "On the dynamics of a meniscus inside capillaries during imbibition and drainage processes: A generalized model, effect of inertia, and numerical algorithm," Physics of Fluids 33, 082104 (2021).

1170 41 H. K. Dahle, M. A. Celia, and S. M. Hassanizadeh, "Bundle-of-tubes model for calculating dynamic effects in the capillary-pressure-saturation relationship," Transport in Porous media 58, 5-22 (2005).

1173 42 R. Douglas and J. Bartley, "Capillary tube models with interaction between the tubes [a note on "immiscible displacement in the interacting capillary bundle model."]," Transport in Porous media 86, 479-482 (2011).

1178 43 J. Bartley and D. Ruth, "Relative permeability analysis of tube bundle models, including capillary pressure," Transport in porous media 45, 445-478 (2001).

1181 44 J. Bartley and D. Ruth, "Relative permeability analysis of tube bundle models," Transport in Porous Media 36, 161-188 (1999).

1183 45 Y. Shiri and S. M. J. Seyed Sabour, "Analytical, experimental, and numerical study of capillary rise dynamics from inertial to viscous flow," Physics of Fluids 34, 102105 (2022).

1186 46 J. Kim, M.-W. Moon, and H.-Y. Kim, "Capillary rise in superhydrophobic rough channels," Physics of Fluids 32, 032105 (2020).

1188 47 J. Bico and D. Quéré, "Precursors of impregnation," EPL (Europhysics Letters) 61, 348 (2003).

1190 48 M. Dong, J. Zhou, et al., "Characterization of waterflood saturation profiles in porous media," Transport in porous media 31, 213-237 (1998).

1193 49 M. Dong, F. A. Dullien, L. Dai, and D. Li, "Immiscible displacement in the interacting capillary bundle model part i. development of interacting capillary bundle model," Transport in Porous media 59, 1-18 (2005).

1196 50 M. Dong, F. A. Dullien, L. Dai, and D. Li, "Immiscible displacement in the interacting capillary bundle model part ii. applications of model and computation," Transport in Porous Media 63, 289-304 (2006).

51 S. Krishnamurthy and Y. Peles, "Gas-liquid two-phase flow across a bank of micropillars," Physics of fluids 19, 043302 (2007).

52 J. Wang, F. A. Dullien, and M. Dong, "Fluid transfer between tubes in interacting capillary bundle models," Transport in Porous Media 71, 115-131 (2008).

53 S. Li, M. Dong, and P. Luo, "A crossflow model for an interacting capillary bundle: Development and application for waterflooding in tight oil reservoirs," Chemical Engineering Science 164, 133-147 (2017).

54 S. Ashraf, G. Visavale, S. S. Bahga, and J. Phirani, "Spontaneous imbibition in parallel layers of packed beads," The European Physical Journal E 40, 39 (2017).

55 S. Ashraf, G. Visavale, and J. Phirani, "Spontaneous imbibition in randomly arranged interacting capillaries," Chemical Engineering Science 192, 218-234 (2018).

56 E. Unsal, G. Mason, N. Morrow, and D. Ruth, "Co-current and counter-current imbibition in independent tubes of non-axisymmetric geometry," Journal of Colloid and Interface Science 306, 105-117 (2007).

57 E. Unsal, G. Mason, D. Ruth, and N. Morrow, "Co-and counter-current spontaneous imbibition into groups of capillary tubes with lateral connections permitting cross-flow," Journal of Colloid and Interface Science 315, 200-209 (2007).

58 E. Unsal, G. Mason, N. R. Morrow, and D. W. Ruth, "Bubble snap-off and capillary-back pressure during counter-current spontaneous imbibition into model pores," Langmuir 25, 3387-3395 (2009).

59 T. Bultreys, K. Singh, A. Q. Raeini, L. C. Ruspini, P.-E. Øren, S. Berg, M. Rücker, B. Bijeljic, and M. J. Blunt, "Verifying pore network models of imbibition in rocks using time-resolved synchrotron imaging," Water Resources Research 56, e2019WR026587 (2020).

60 T. Bultreys, K. Singh, A. Q. Raeini, P.-E. Øren, S. Berg, B. Bijeljic, and M. J. Blunt, "Improving the description of two-phase flow in rocks by integrating pore scale models and experiments," in *InterPore 11th Annual Meeting and Jubilee* (2019) p. 87.

61 S. Foroughi, B. Bijeljic, and M. J. Blunt, "Pore-by-pore modelling, validation and prediction of waterflooding in oil-wet rocks using dynamic synchrotron data," Transport in Porous Media 138, 285-308 (2021).

62 S. Ashraf and J. Phirani, "A generalized model for spontaneous imbibition in a horizontal, multi-layered porous medium," Chemical Engineering Science 209, 115175 (2019).

63 C. W. Hirt and B. D. Nichols, "Volume of fluid (vof) method for the dynamics of free boundaries," Journal of computational physics 39, 201-225 (1981).

64 T. Young, "III. an essay on the cohesion of fluids," Philosophical Transactions of the Royal Society of London 95, 65-87 (1805).

65 P. S. de Laplace, *Supplément au dixième livre du Traité de mécanique céleste: sur l'action capillaire* (1806).

66 L. Ding, Q. Wu, L. Zhang, and D. Guérillot, "Application of fractional flow theory for analytical modeling of surfactant flooding, polymer flooding, and surfactant/polymer flooding for chemical enhanced oil recovery," Water 12, 2195 (2020).

67 Y. Debbabi, M. D. Jackson, G. J. Hampson, P. J. Fitch, and P. Salinas, "Viscous crossflow in layered porous media," Transport in Porous Media 117, 281-309 (2017).

68 S. Akbari, S. M. Mahmood, H. Ghaedi, and S. Al-Hajri, "A new empirical model for viscosity of sulfonated polyacrylamide polymers," Polymers 11, 1046 (2019).

69 S. Ashraf and J. Phirani, "Capillary displacement of viscous liquids in a multi-layered porous medium," Soft matter 15, 2057-2070 (2019).

70 S. Ashraf and J. Phirani, "Capillary impregnation of viscous fluids in a multi-layered porous medium," in *Fluids Engineering Division Summer Meeting*, Vol. 59087 (American Society of Mechanical Engineers, 2019) p. V005T05A057.
<figure><img src="image_24.png" /></figure>

(b)

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

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

$\theta_w$
<figure><img src="image_27.png" /></figure>
$$z_{\delta} > z_{\gamma} > z_{\beta} > z_{\alpha}$$

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

$$z_{\gamma} > z_{\delta} > z_{\beta} > z_{\alpha}$$

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

$$z_{\gamma} > z_{\beta} > z_{\delta} > z_{\alpha}$$

<figure><img src="image_30.png" /></figure>
<figure><img src="image_31.png" /></figure>
$$z_{\delta} > z_{\gamma} > z_{\beta} > z_{\alpha}$$

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

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

$$z_{\gamma} > z_{\delta} > z_{\beta} > z_{\alpha}$$

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

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

$$z_{\gamma} > z_{\beta} > z_{\delta} > z_{\alpha}$$

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

<figure><img src="image_37.png" /></figure>
$$z_{\delta} > z_{\gamma} > z_{\beta} > z_{\alpha}$$

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

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

$$z_{\delta} > z_{\beta} > z_{\gamma} > z_{\alpha}$$

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

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

$$z_{\beta} > z_{\delta} > z_{\gamma} > z_{\alpha}$$

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

<figure><img src="image_43.png" /></figure>
$z_{\delta} > z_{\gamma} > z_{\beta} > z_{\alpha}$

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

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

$z_{\delta} > z_{\beta} > z_{\gamma} > z_{\alpha}$

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

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

$z_{\gamma} > z_{\delta} > z_{\beta} > z_{\alpha}$

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

<figure><img src="image_49.png" /></figure>
<figure><img src="image_50.png" /></figure>

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

<figure><img src="image_52.png" /></figure>
<figure><img src="image_53.png" /></figure>

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

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

<figure><img src="image_56.png" /></figure>
<figure><img src="image_57.png" /></figure>
<figure><img src="image_58.png" /></figure>

<figure><img src="image_59.png" /></figure>
<figure><img src="image_60.png" /></figure>
<figure><img src="image_61.png" /></figure>

slope 2.67

slope 2

lead exp; Ref. 47

lead sim; n = 6

lead sim; n = 11

lag exp; Ref. 47

lag sim; n = 6

lag sim; n = 11

lead Eq. 15

lag Eq. 15
<figure><img src="image_62.png" /></figure>

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

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

<figure><img src="image_65.png" /></figure>
<figure><img src="image_66.png" /></figure>

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