# Experimental and simulation analysis of capillary impregnation and wettability in cotton fiber porous media

Shaolin Liu, Amine Ben-Abdelwahed

► To cite this version:

Shaolin Liu, Amine Ben-Abdelwahed. Experimental and simulation analysis of capillary impregnation and wettability in cotton fiber porous media. Colloids and Surfaces A: Physicochemical and Engineering Aspects, 2025, 720, pp.136991. 10.1016/j.colSURFA.2025.136991. hal-05133400

HAL Id: hal-05133400

https://artsetmetiers.hal.science/hal-05133400v1

Submitted on 27 Jun 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Archive Ouverte - Open Repository

Science Arts & Métiers (SAM)

is an open access repository that collects the work of Arts et Métiers Institute of
Technology researchers and makes it freely available over the web where possible.

This is an author-deposited version published in: https://sam.ensam.eu
Handle ID: http://hdl.handle.net/10985/26369

To cite this version :

Shaolin LIU, Amine BEN ABDELWAHED - Experimental and simulation analysis of capillary impregnation and wettability in cotton fiber porous media - Colloids and Surfaces A - Vol. 720, p.136991 - 2025
<figure><img src="image_3.png" /></figure>

# Experimental and simulation analysis of capillary impregnation and wettability in cotton fiber porous media

Shaolin Liu<sup>*</sup>, Amine Ben-Abdelwahed<sup>*</sup>

Univ. Bordeaux, CNRS, Bordeaux INP, I2M, UMR 5295, F-33400, Talence, France

Arts et Metiers Institute of Technology, CNRS, Bordeaux INP, I2M, UMR 5295, F-33400 Talence, France

## GRAPHICAL ABSTRACT

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

tensometer

sample holder dimensions

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

n-heptane saturation

Comparison of simulation and experimental results

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

## ARTICLE INFO

### Keywords:

Wettability

Numerical study

Fiber swelling

Permeability

Liquid absorption

## ABSTRACT

Capillary-driven flow is a critical phenomenon in fibrous porous media, influencing a wide range of industrial processes such as filtration, drying, and biomass conversion. Understanding the interplay between material properties and fluid behavior is essential for optimizing these processes. This study aims to analyze capillary flow behavior in cotton fibers by combining numerical simulations with experimental data. An inverse analysis approach was employed to determine several unknown parameters associated with three different liquids: n-heptane, water, and glycerol. Based on simulation data, the wicking velocity of each liquid in the porous medium was calculated, followed by the determination of the corresponding capillary numbers (Ca). The results indicate that the capillary numbers for n-heptane, water, and glycerol are 4.91 × 10⁻⁵, 1.61 × 10⁻⁵, and 1.42 × 10⁻⁴, respectively, highlighting the differences in infiltration behavior among these fluids. Subsequent simulations systematically examined the effects of porosity, surface tension, and dynamic contact angle on liquid transport. The capillary number was used as a quantitative measure to assess the influence of these parameters on infiltration behavior. Lower Ca values were generally associated with faster liquid absorption, particularly in cases with higher surface tension and smaller contact angles, where capillary forces dominated.

<footnote>* Correspondence to: Arts et Métiers Institute of Technology, 33400, Talence, France.</footnote>

<footnote>E-mail address: Amine.BEN-ABDELWAHED@ensam.eu (A. Ben-Abdelwahed).</footnote>
This study provides a modeling framework to evaluate parameter effects on capillary flow in fibrous porous media, offering guidance for optimizing materials in filtration and drying applications. It also contributes to a deeper understanding of interfacial transport phenomena in porous structures.

## Nomenclature

### Latin symbols

<table><tr><td>a</td><td>Coefficient related to fiber swelling, s<sup>-1</sup></td></tr><tr><td>b</td><td>Absorption coefficient, s<sup>-1</sup></td></tr><tr><td>c</td><td>A constant value accounting for the tortuous path of liquid</td></tr><tr><td>H</td><td>sample height, m</td></tr><tr><td>h</td><td>Distance traveled by the liquid front, m</td></tr><tr><td>K</td><td>Permeability, m<sup>2</sup></td></tr><tr><td>m</td><td>Mass, kg</td></tr><tr><td>P</td><td>Pressure, Pa</td></tr><tr><td>Q</td><td>Mass source, kg m<sup>-3</sup> s<sup>-1</sup></td></tr><tr><td>r</td><td>Radius of the capillary tube, m</td></tr><tr><td>r̄</td><td>Mean capillary radius of a porous medium, m</td></tr><tr><td>R</td><td>Inner radius of the sample holder, m</td></tr><tr><td>S</td><td>Sink term due to liquid absorption by the solid matrix, s<sup>-1</sup></td></tr><tr><td>s<sub>i</sub></td><td>Saturation or volume fraction of each phase i</td></tr><tr><td>t</td><td>Time of flow, s</td></tr><tr><td>u</td><td>Darcy velocity vector, m s<sup>-1</sup></td></tr><tr><td>u</td><td>Fluid velocity in porous medium, m s<sup>-1</sup></td></tr><tr><td>V<sub>f</sub></td><td>The original, untreated fiber volume, m<sup>3</sup></td></tr><tr><td>V<sub>final</sub></td><td>The final volume of the sample after compression, m<sup>3</sup></td></tr><tr><td>x</td><td>X-coordinate</td></tr></table>

### Greek symbols

<table><tr><td>γ<sub>L</sub></td><td>Liquid surface tension, mN m<sup>-1</sup></td></tr><tr><td>λ<sub>p</sub></td><td>Pore size distribution index</td></tr><tr><td>μ</td><td>Dynamic viscosity of the liquid, Pa s</td></tr><tr><td>∇</td><td>Gradient operator</td></tr><tr><td>ρ</td><td>Density of the liquid, kg m<sup>-3</sup></td></tr><tr><td>θ</td><td>Contact angle, rad</td></tr><tr><td>ε</td><td>Mean porosity</td></tr></table>

### Subscripts

<table><tr><td>atm</td><td>Atmospheric</td></tr><tr><td>a</td><td>Advancing</td></tr><tr><td>cap</td><td>Capillary</td></tr><tr><td>e</td><td>Experiment</td></tr><tr><td>l</td><td>Liquid</td></tr><tr><td>n</td><td>Numerical</td></tr><tr><td>s</td><td>Static</td></tr></table>

# 1. Introduction

Natural fibers, such as cotton, hemp, jute, and wood, have gained widespread attention as sustainable materials due to their biodegradability, cost-effectiveness, and potential for biogenic carbon storage [1, 2]. Their applications range from composite reinforcement to bioenergy and filtration systems [3]. However, their inherent hydrophilicity makes them highly sensitive to moisture, which affects their dispersion and interactions with liquid media [4]. The effectiveness of

fiber impregnation directly influences the efficiency and quality of industrial processes such as composite manufacturing and chemical treatments [5]. An inadequate understanding of liquid transport in fibrous media can result in uneven wetting, void formation [6,7], and compromised material performance [8].

Capillary-driven flow governs liquid infiltration in fibrous porous media, and its behavior is influenced by various factors, including fiber geometry, liquid properties, surface wetting characteristics, and the medium's overall permeability [9]. The arrangement and packing density of fibers determine the effective pore size distribution and affect the rate at which liquid spreads through the structure [10,11]. Additionally, the surface chemistry of the fibers influences wettability, while liquid viscosity and surface tension define the capillary forces driving infiltration [12].

Theoretical models, such as the Washburn equation [13,14], are widely used to describe capillary rise and liquid absorption in porous structures. These models provide a simplified mathematical framework by balancing viscous dissipation and capillary forces to estimate infiltration rates [15]. However, they often rely on idealized assumptions, such as uniform fluid geometry, negligible inertia, and a fixed contact angle. In practice, during the initial phase of infiltration, liquid acceleration cannot be ignored, and inertial effects play a role [16]. Dynamic wetting, including variations in the advancing contact angle, further complicates the process, particularly in fibrous porous media where surface roughness and chemical heterogeneity influence liquid spreading [6,17,18]. In later stages, gravity becomes more relevant, affecting the final equilibrium position of the liquid front. Fibrous porous media also exhibit a high degree of heterogeneity due to variations in fiber orientation, diameter, and packing density [19]. These structural irregularities lead to deviations from classical capillary models, resulting in variations in local infiltration dynamics. Unlike rigid porous materials with well-defined pore networks, fibrous systems can undergo structural changes upon wetting, which further influences permeability and capillary transport. As a result, liquid infiltration in these materials does not always follow a straightforward or predictable pattern. Experimental studies have provided insight into these infiltration behaviors by measuring parameters such as capillary pressure [9], dynamic contact angles [17,20], and liquid front propagation [21,22]. However, directly observing the infiltration process in dense fiber networks remains a challenge due to optical limitations and the complexity of tracking liquid movement in three-dimensional structures. Traditional measurement techniques, such as weight-based analysis and tensiometry [20,23], offer macroscopic data but lack spatial resolution, making it difficult to capture local variations in liquid distribution. These limitations highlight the need for complementary approaches, such as numerical simulations, to develop a more detailed understanding of capillary-driven flow in fibrous porous media.

Numerical simulations based on the Navier-Stokes equations are widely employed to model fluid flow, often relying on simplifying assumptions such as inviscid or incompressible flow [15,24]. To address these limitations, researchers have developed various computational methods. One such method is the Lattice Boltzmann Method (LBM), which models macroscopic fluid dynamics by simulating particle motion on discrete lattices. Raiskinmaki et al. [25] applied the Shan-Chen model to study capillary rise, demonstrating its applicability in larger capillary tubes and analyzing the relationship between dynamic contact angle and capillary number. Wolf et al. [26] further investigated the characteristics of capillary rise, focusing on the evolution of the dynamic contact angle. Additionally, Lu et al. [27] examined the initial stages of capillary rise and highlighted the limitations of the Lucas-Washburn (L-W) equation in describing this phase accurately. Another
well-established tool is COMSOL Multiphysics, which employs the finite element method (FEM) to solve continuum mechanics problems. It can also be coupled with wetting dynamics models — such as the Cahn–Hilliard equation [28] or Gibbs free energy models [29] — to track interface motion and contact angle variations. Patari et al. [21] used a numerical method to study liquid wicking in a porous medium (paper strips), integrating the L–W equation and evaporation effects to explain the underlying physics of capillary rise. Similarly, Masoodi et al. [30] applied Darcy’s law to simulate liquid infiltration into a swelling porous medium composed of cellulose and superabsorbent fibers. Their model coupled the flow equation with the mass conservation equation, incorporating additional source and sink terms to account for matrix swelling and liquid absorption. Furthermore, OpenFOAM, an open-source finite volume method (FVM)-based CFD solver, simulates free-surface flows and large-scale fluid absorption processes. Horgue et al. [31] developed a numerical model and integrated it into OpenFOAM to extend its capabilities for modeling absorption phenomena. Liu et al. [32] further utilized this model to investigate fluid transport in absorbing materials, demonstrating that capillary forces govern internal fluid movement within the medium. These studies highlight the advantages of different numerical approaches in modeling capillary-driven flow and wetting dynamics. However, further research is needed to refine the understanding of flow mechanisms in multiphase systems, particularly in scenarios involving heterogeneous porous structures, dynamic wetting behavior, and the coupling between capillary forces and other transport phenomena. Such advancements would improve simulation accuracy and provide deeper insights into experimental observations.

In our previous study, we conducted experimental investigations on the spontaneous capillary impregnation of cotton fibers, utilizing tensiometric methods to quantify key parameters such as fiber geometric products, apparent advancing contact angles, and capillary pressure [20]. While the experimental results provided an understanding of the impregnation dynamics, they had limitations in resolving the internal fluid distribution within the porous medium. To overcome these constraints, the present study incorporates numerical simulations using COMSOL Multiphysics to complement the experimental findings. A 2D axisymmetric model was developed to simulate capillary infiltration in compressed cotton fibers, accounting for porosity, density, permeability, and dynamic wetting behavior. This approach enables a detailed analysis of the transient evolution of liquid penetration, helping to identify the limitations of theoretical models and improve predictive accuracy. The numerical results are systematically compared with experimental data to assess model reliability and further refine key parameters such as permeability.

The structure of this paper is as follows: Section 2 introduces the theoretical framework for describing capillary-driven flow in fibrous porous media. Section 3 provides a brief overview of the previous experimental setup and details the numerical modeling approach based on COMSOL. Section 4 presents a comparative analysis between the numerical simulations and experimental results, examining the influence of various parameters on wetting behavior during fluid transport through porous media. Finally, Section 5 summarizes the key findings and outlines potential directions for future research.

## 2. Theoretical framework for capillary-driven flow in porous media

This section presents the theoretical basis for capillary-driven flow in fibrous porous media, divided into two parts. The first part discusses wicking in rigid porous media, while the second focuses on wicking in nonrigid, swelling porous media. For each case, two fundamental models are examined: the Washburn equation (also known as the Lucas–Washburn equation) and Darcy’s law, which describe fluid transport under different conditions.

### 2.1. Capillary flow in rigid porous media

#### 2.1.1. Washburn equation

In rigid fibrous porous media, fiber dimensions, and pore structure remain unchanged during the wicking process. The Washburn equation describes capillary rise in a tube, which can be extended to porous media [33]. To provide context, examining the capillary rise of a liquid in a tube is useful. Numerous studies have addressed this topic. Considering the balance between viscous and capillary forces and assuming initial conditions of $t = 0$ and $dh/dt = 0$, the height of the rising liquid front and its corresponding mass can be described as follows: [9,30,34].

$$h(t) = \sqrt{\frac{r\gamma_L \cos\theta_a}{2\mu}} t \quad (1)$$

$$m(t) = \rho \cdot h(t) \cdot \pi r^2 = \rho \sqrt{\frac{r\gamma_L \cos\theta_a}{2\mu}} t \pi r^2 \quad (2)$$

where $h(t)$ and $m(t)$ represent the height and mass of the rising liquid front, respectively. The parameters include $r$, the capillary tube radius; $\theta_a$, the apparent advancing contact angle; $\mu$, the liquid viscosity; and $t$, the flow time. The term $\gamma_L$ denotes the liquid surface tension.

Taking into account only capillary and viscous forces, an ideal porous medium, packed into a column and represented as a collection of individual capillary tubes with an average radius $\bar{r}$, follows the Washburn law, which can be written as [9,23]:

$$h(t) = \sqrt{\frac{c\bar{r}\gamma_L \cos\theta_a}{2\mu}} t \quad (3)$$

$$m(t) = \rho \cdot \varepsilon \cdot h(t) \cdot \pi R^2 = \rho \varepsilon \sqrt{\frac{c\bar{r}\gamma_L \cos\theta_a}{2\mu}} t \pi R^2 \quad (4)$$

where $c$ is a geometric constant that accounts for the liquid’s tortuous path and is inversely related to tortuosity. $\bar{r}$ represents the average capillary radius. $R$ denotes the inner radius of the sample holder (see Fig. 1), while $\varepsilon$ refers to the mean porosity of the medium.

#### 2.1.2. Darcy law

The flow of a Newtonian liquid through an isotropic and rigid porous medium is governed by Darcy’s law and the continuity equation [30]. Originally formulated to predict the pressure drop of water through sand beds, Darcy’s law has been extended to relate the Darcy velocity vector $\mathbf{u}$ to the pressure gradient $\nabla P$, incorporating permeability tensor $\mathbf{K}$ as a parameter that defines the medium’s ability to transmit a liquid of viscosity $\mu$.

$$\mathbf{u} = -\frac{1}{\mu} \mathbf{K} \cdot \nabla P \quad (5)$$

$$\nabla \cdot \mathbf{u} = 0 \quad (6)$$

Pucci et al. [9] demonstrated that, based on the one-dimensional form of Eqs. (5) and (6), capillary wicking in a porous medium within a cylindrical holder leads to an equation for the liquid-front position:

$$\frac{dh}{dt} = \frac{u}{\varepsilon} = -\frac{K}{\varepsilon\mu} \frac{dP}{dh} = -\frac{K}{\varepsilon\mu} \frac{(P_{\text{atm}} - P_{\text{cap}}) - P_{\text{atm}}}{dh} = \frac{K}{\varepsilon\mu} \frac{P_{\text{cap}}}{dh} \quad (7)$$

where the capillary pressure ($P_{\text{cap}}$) characterizes the wicking behavior in an equivalent homogeneous porous medium. Integrating the equation over time yields expressions for the liquid front height $h(t)$ and mass $m(t)$ as functions of time $t$, with initial conditions of $t = 0$, $h(t) = 0$, and $m(t) = 0$,

$$h(t) = \sqrt{\frac{2K P_{\text{cap}}}{\varepsilon\mu}} t \quad (8)$$

$$m(t) = \rho \cdot \varepsilon h(t) \cdot \pi R^2 = \rho \varepsilon \sqrt{\frac{2K P_{\text{cap}}}{\varepsilon\mu}} t \cdot \pi R^2 \quad (9)$$
<figure><img src="image_7.png" /><figcaption>Fig. 1. Experimental images and schematic diagram of the sample holder.</figcaption></figure>

By comparing the expressions for the absorbed mass $m(t)$ in Eqs. (4) and (9), a new relationship for the capillary pressure $P_{cap}$ can be derived.

$$P_{cap} = (c\bar{r})\epsilon \frac{\gamma_L \cos\theta_a}{4K} \qquad (10)$$

Unlike the conventional expression obtained through the Laplace equation $P_{cap} = (2\gamma_L \cos\theta_a)/r$ [30], this new formulation incorporates additional parameters, including the tortuosity factor $c$, mean capillary radius $\bar{r}$, porosity $\epsilon$, and permeability $K$, as well as fluid-porous medium interaction parameters such as liquid surface tension $\gamma_L$ and apparent advancing contact angle $\theta_a$. This relationship has also been validated by Pucci et al. [9,23].

## 2.2. Capillary flow in swelling porous media

Unlike rigid porous media, swelling porous media exhibits time- and space-dependent changes in pore radius, porosity, and permeability during capillary wicking. As these parameters evolve dynamically with the wetting process, modifications to the previously discussed models are necessary. The following discussion is structured around the Washburn equation and Darcy's law.

### 2.2.1. Washburn equation

Schuchardt and Berg [35] introduced a modified Washburn equation incorporating swelling dynamics to predict capillary wicking in

natural fiber bundles:

$$h(t) = \sqrt{\frac{r_0 \gamma_L \cos\theta_a}{2\mu}} \sqrt{t - \frac{a}{r_0} t^2 + \frac{a^2}{3r_0^2} t^3} \qquad (11)$$

In this model, fiber swelling reduces the hydraulic radius behind the flow front, expressed as: $r = r_0 - at$ where $r$ is the hydraulic pore radius, $r_0$ is its initial value, and $a$ represents the swelling effect. When swelling is neglected ($a = 0$, Eq. (11) simplifies to Eq. (1)), recovering the original Washburn equation for rigid porous media.

### 2.2.2. Darcy law

For swelling porous media, Masoodi et al. [30] have reviewed the modified forms of Darcy's law; thus, only a brief discussion is provided. The modified continuity equation is given as:

$$\nabla \cdot \mathbf{u} = -S - \frac{\partial \epsilon}{\partial t} \qquad (12)$$

where $-S$ represents the sink term due to liquid absorption by the solid matrix, while $\partial\epsilon/\partial t$ accounts for liquid generation within the pores. As swelling progresses, porosity gradually decreases, contributing to this effect. To simplify the analysis, the wicking flow can be assumed to be one-dimensional (x-direction), allowing Darcy's law and the continuity equation (Eq. (12)) to be reduced to the following form.

$$\text{Darcy's Law: } u = -\frac{K}{\mu} \frac{dP}{dx} \qquad (13)$$

$$\text{Continuity Equation: } \frac{du}{dx} = -S - \frac{\partial \epsilon}{\partial t} \qquad (14)$$
<table><caption>**Table 1**
Parameters used in simulations for different liquids based on experimental methods.</caption><thead><tr><th>Name</th><th>Value</th><th>Description</th></tr></thead><tbody><tr><td>H (mm)</td><td>31</td><td>Sample height</td></tr><tr><td>R (mm)</td><td>6</td><td>Inner radius of the sample holder</td></tr><tr><td>c&#773; (µm)</td><td>6.97</td><td>Geometric products</td></tr><tr><td>K (m²)</td><td>10<sup>-10</sup>-10<sup>-13</sup></td><td>Permeability</td></tr><tr><td>ρ<sub>g</sub> (kg/m³)</td><td>1</td><td>Air density</td></tr><tr><td>μ<sub>g</sub> (Pa s)</td><td>1.76 × 10<sup>-5</sup></td><td>Air viscosity</td></tr><tr><td>P<sub>cap</sub> (Pa)</td><td>(c<sup>2</sup>F)ε<sub>g</sub> cosθ<sub>a</sub>/4K</td><td>Entry capillary pressure</td></tr><tr><td>ε</td><td>0.53</td><td>Porosity</td></tr><tr><td>γ<sub>t</sub> (mN/m)</td><td>n-heptane</td><td>Water</td></tr><tr><td>γ<sub>a</sub> (°)</td><td>20.10</td><td>72.75</td></tr><tr><td>ρ<sub>i</sub> (kg/m³)</td><td>0</td><td>72.8</td></tr><tr><td>μ<sub>i</sub> (Pa s)</td><td>684</td><td>1000</td></tr><tr><td>μ<sub>t</sub> (mN/m)</td><td>0.402 × 10<sup>-3</sup></td><td>1.002 × 10<sup>-3</sup></td></tr><tr><td>ε<sub>t</sub></td><td>0.53</td><td>Glycerol</td></tr><tr><td>ε<sub>a</sub></td><td>0.53</td><td>Liquid surface tension [36]</td></tr><tr><td>γ<sub>t</sub></td><td>0</td><td>63</td></tr><tr><td>γ<sub>a</sub></td><td>0</td><td>69.55</td></tr><tr><td>ρ<sub>i</sub></td><td>684</td><td>1261</td></tr><tr><td>μ<sub>i</sub></td><td>0.402 × 10<sup>-3</sup></td><td>1.490</td></tr><tr><td>μ<sub>t</sub></td><td>0.402 × 10<sup>-3</sup></td><td>Liquid viscosity</td></tr></tbody></table>

Finally, under a set of assumptions — such as the liquid absorption rate being linearly related to the change in solid volume fraction and porosity and permeability depending only on time — Masoodi et al. [30] derived the following expression for the liquid front height:

$$h(t) = \sqrt{\frac{2P_{cap}}{\epsilon_0 \mu}} e^{(b-1)\frac{\epsilon_t}{\epsilon}} \int_0^t e^{(1-b)\frac{\epsilon_t}{\epsilon}} K(t') dt' \quad (15)$$

where b represents the absorption coefficient: b = 0 indicates no liquid absorption by the solid matrix, Eq. (15) simplifies to the form of Eq. (8). While b = 1 corresponds to the rate of solid matrix growth. ε denotes the initial porosity, and ε_t is the surface porosity at the liquid front.

## 3. Experimental setup and numerical modeling

This section provides a brief overview of the previous experimental setup [20] and outlines the numerical modeling approach using COMSOL. The objective is to validate the experimental results through numerical simulations and further analyze key parameters required in the simulations and their specific effects.

### 3.1. Overview on experiment setup

In our previous work, the behavior of natural fibers (cotton fibers) during capillary impregnation was experimentally investigated using a dynamic approach [20]. Three different liquids — n-heptane, water, and glycerol — were used. The experimental materials and setup are illustrated in Fig. 1, where a series of photographs and diagrams provide a visual overview of the sample preparation and loading process. Fig. 1(a) shows the rolled cotton fiber porous medium. Fig. 1(b) presents photographs of the experimental components used in the capillary test. Fig. 1(c) shows the process of inserting the prepared sample into the Washburn cell. Fig. 1(d) displays the Washburn cell mounted on the test instrument. Fig. 1(e) illustrates the structure and dimensions of the sample holder. Finally, Fig. 1(f) compares the sample volume before and after compression.

The sample preparation procedure consists of the following steps: (1) The cotton fiber bundle is first tightly rolled, with care taken to align the fiber orientation along the vertical axis of the Washburn cell (Fig. 1(a)). (2) A filter paper is placed at the bottom of the Washburn cell to separate the sample from the base (Fig. 1(b)). (3) The rolled sample is carefully inserted into the cell to ensure a snug fit (Fig. 1(c)). (4) The maximum recommended weight is applied five times to pre-compress the sample and eliminate air gaps. (5) A small metal cylinder is placed on top of the sample to ensure uniform pressure distribution. (6) The spacer is screwed down to its minimum height to hold the sample firmly during testing. (7) Finally, the spacer is fully secured, completing the preparation. The resulting compressed sample and its dimensions are shown in Fig. 1(f).

The surface tension of each liquid was measured first. Fiber samples were then placed in the sample holder for wicking tests. Test curves were obtained using a tensiometric method, and the modified Washburn equation was applied to determine parameters such as the geometric products c&#773; and apparent advancing contact angles θ_a, and surface free energy. Experimental results included the geometric products of cotton fibers within the sample holder, 6.97 ± 0.86 µm, and the apparent advancing contact angles for water and glycerol, 74.93 ± 2.20° and 69.55 ± 1.83° under a porosity of 0.53. Detailed results are presented in Table 1, and a more in-depth discussion is provided in the simulation part.

The purpose of reviewing the experimental setup here is to provide specific structural parameters and experimentally obtained values, which will be used in the following simulations. Detailed information can be found in our previous work [20].

### 3.2. Numerical modeling approach using COMSOL

#### 3.2.1. Numerical modeling

This part describes the numerical methods used in this work. Specifically, COMSOL Multiphysics, based on the finite element method (FEM), was employed to simulate wicking flow in a porous medium (cotton fibers). The simulation involved multiphysics coupling using the porous media and subsurface flow module, combining the Darcy's Law interface with the phase transport in the porous media interface [21]. Darcy's Law was applied to model liquid movement within the porous medium, while the transport of diluted species interface governed the independent equations for the volume fraction s_i of wetting (liquid) or non-wetting fluid i (air), as shown below [21]:

$$\frac{\partial}{\partial t}(\epsilon \rho_{s_i} s_i) + \nabla \cdot \left( -\rho_{s_i} K \frac{K_{rs_i}}{\mu} (\nabla P_{s_i} - \rho_{s_i} g) \right) = Q_m = 0 \quad (16)$$

The parameters used in the equations are defined as follows: Q_m represents the mass source. ε denotes the porosity of the medium. K is the permeability of the medium, measured in m², K<sub>rs_i</sub> refers to the relative permeability, which depends on the saturation of the given fluid. μ is the dynamic viscosity of the fluid, expressed in Pa·s. P<sub>s_i</sub> indicates the pressure of phase i, measured in Pa. ρ<sub>s_i</sub> represents the density of the fluid in phase i, given in kg/m³. The volume fractions of each phase add up to 1, expressed as: s<sub>1</sub> + s<sub>2</sub> = 1. The equations for the wetting phase in the Brooks and Corey model [37] are expressed as follows:

$$P_c = P_{cap} \cdot s_1^{-1/\lambda_p} \quad (17)$$

$$K_{rs_1} = s_1^{(3+2/\lambda_p)} \quad (18)$$

$$K_{rs_2} = s_2^2 \cdot (1 - s_1^{1+2/\lambda_p}) \quad (19)$$

where P_c is the capillary pressure, and P<sub>cap</sub> represents the entry capillary pressure, as shown in Eq. (10) ((c&#773;)ε<sub>g</sub> cosθ<sub>a</sub>/4K). λ_p denotes the pore size distribution index.

#### 3.2.2. Parameter selection and implementation

The two-dimensional axisymmetric model, developed to correspond to the experimental structure, is shown in Fig. 2. The geometry represents the sample holder with a height of H = 31 mm and a radius of R = 6 mm. This model is designed to simulate the wicking process in a porous medium under similar conditions to the experimental setup.

The porous medium (cotton fibers) is initially filled with air, with the water saturation set to 0.01. A time-dependent simulation is carried out using the Brooks and Corey model, which describes capillary pressure and relative permeabilities. In the two-dimensional axisymmetric structure, the right-hand side boundary is defined as impermeable, meaning no fluid flow is allowed across it. The Darcy's Law boundary conditions are as follows: At the bottom boundary, the water phase
<figure><img src="image_8.png" /><figcaption>Fig. 2. Model geometry and boundary conditions.</figcaption></figure>

is set to atmospheric pressure. At the top boundary, the pressure is specified as the hydrostatic atmospheric pressure minus the capillary pressure. For the phase transport boundary conditions: At the bottom boundary, no flux is assumed for the air phase, representing the absence of air movement through the bottom of the porous medium. At the top boundary, a mass flux is applied for the air phase. This flux is driven by the pressure gradient determined by Darcy's Law, ensuring consistency between the phase transport and the flow model.

The initial conditions specify that the cotton fibers are saturated only with air, and the liquid phase starts infiltrating from the bottom boundary. The boundary and initial conditions are carefully set up to simulate the capillary-driven wicking process realistically. The porosity and permeability of the cotton fibers are assumed constant during the simulation, while the pressure and saturation fields evolve dynamically over time to capture the interaction between the air and water phases within the porous medium.

The parameters used in simulations for different liquids, based on experimental methods, are summarized in Table 1. These parameters include fixed geometric properties of the model, such as the sample height H and radius R, as well as constant physical properties of the fluids and air, such as viscosity µ and density ρ. The geometric products of the sample c̄r, porosity ε, liquid surface tension γL, and advancing contact angle θa were obtained from previous experiments [20]. Permeability values K refer to the flow resistance along the axis of cotton fibers, assuming a regular hexagonal packing arrangement. The estimation is based on the following theoretical relationship [10,38]:

$$K = \frac{8r_f^2}{53} \frac{\epsilon^3}{(1-\epsilon)^2} \qquad (20)$$

where r_f is the diameter of an individual fiber and ε is the material porosity. The fiber diameter r_f typically ranges from 5 to 24 µm, as reported in the literature [10]. A porosity value of 0.53 was used, as shown in Table 1. Based on this relation, the estimated permeability falls within the range of 10⁻¹² to 10⁻¹⁰ m². To account for possible variations in fiber deformation and non-ideal arrangement, a more conservative order of magnitude range of 10⁻¹³ to 10⁻¹⁰ m² was considered for further analysis.

## 4. Results and discussion

In this section, the simulated results are compared with experimental data, and the specific values of solid permeability for different

<figure><img src="image_9.png" /><figcaption>Fig. 3. n-heptane saturation in the cotton fibers after 1, 5, 10, 20, and 40 s, (K = 6.3 × 10⁻¹² m², c̄r = 6.97 µm, ε = 0.53, γL = 20.10 mN/m, θa = 0°, ρl = 684 kg/m³, μl = 0.402 × 10⁻³ Pas).</figcaption></figure>

liquids are determined using an inverse analysis approach. Additionally, the fluid flow behavior within the cotton fibers is analyzed. Subsequently, the influence of various parameters listed in Table 1 on the flow behavior is examined. Specifically, the effects of different parameters on the distribution of liquid saturation and the variation of absorbed liquid mass over time are analyzed. Each aspect is discussed in detail below.

### 4.1. Comparison of simulated and experimental results

In the experiments, three different fluids — n-heptane, water, and glycerol — were used. This subsection focuses on analyzing the cases corresponding to each fluid.

#### 4.1.1. n-heptane - cotton fibers

First, the case with n-heptane as the fluid is examined. Since the exact permeability value is unknown, an estimated value derived from theoretical equations, 6.3 × 10⁻¹² m², is used. Other parameters correspond to the n-heptane case listed in Table 1.

A transient simulation was conducted for 100 s. Actually, after 40 s, the cotton fibers are fully saturated with n-heptane. Fig. 3 presents the spatial distribution of n-heptane saturation (the volume fraction of n-heptane) in the cotton fibers at selected time steps: 1, 5, 10, 20, and 40 s. The saturation variable is dimensionless and ranges from 0 to 1, divided into ten levels in the colormap. A value of 0 indicates that the liquid has not yet penetrated the cotton fibers, while a value of 1 corresponds to complete saturation. The saturation field over time demonstrates the progression of n-heptane infiltration into the cotton fibers. At t = 1 s, the liquid saturation front is localized near the bottom boundary, indicating the initial stage of infiltration. By t = 5 s, the liquid has advanced further, and a distinct transition zone between the saturated and unsaturated regions is observed. At t = 10 s and beyond, the saturation becomes more uniform, with the fibers nearing full saturation by t = 40 s. This progression highlights the efficiency of n-heptane in quickly saturating cotton fibers due to its low viscosity and surface tension. Although the contour plots at t = 10 s, 20 s, and 40 s appear visually similar due to the system approaching full saturation, this trend is quantitatively confirmed in Fig. 5(a), where the squared mass of absorbed liquid remains almost constant after t = 10 s (to be discussed in detail later). However, Fig. 5(b) shows that when material parameters such as permeability are varied, the absorption process may continue beyond 10 s. Therefore, to maintain consistency with the subsequent parametric analysis and to support comparative discussions, the result at t = 40 s is retained in Fig. 3.

To further improve clarity in the presentation of spatial saturation fields, a supplementary figure has been included in Appendix A. This
<figure><img src="image_10.png" /><figcaption>Fig. 4. n-heptane pressure field in the cotton fibers after 5 and 40 s.</figcaption></figure>

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

figure provides the saturation distribution of n-heptane at a single time step ($t = 1$ s) with fully labeled X and Y axes, which correspond to the width and height of the computational domain, respectively. While composite figures such as Figs. 3, 7, and others omit the horizontal axis labels to allow for side-by-side comparison across time steps, the Appendix figure serves as a reference to help readers interpret the spatial dimensions involved. More detailed information regarding the axis definition and colorbar labeling is provided in Appendix.

Fig. 4 displays the pressure field (in Pa) of n-heptane within the cotton fibers at two representative time steps: $t = 5$ s and $t = 40$ s. Streamlines of the velocity field are superimposed to illustrate the direction and distribution of the overall flow. The pressure values are visualized using a colormap, where the maximum and minimum values are indicated in the colorbar. At $t = 5$ s, a large pressure drop (around $1.98 \times 10^4$ Pa) is evident near the advancing front, driving the liquid into the porous medium. By $t = 40$ s, the pressure drop diminishes (around 212 Pa), and the pressure field stabilizes, reflecting the equilibrium condition as the cotton fibers become fully saturated.

In our experiments, the exact permeability of the cotton fibers was not directly known. By comparing simulation results with experimental data, the permeability was determined using an inverse analysis approach. Additionally, in the experiments, filter paper was placed at the bottom of the sample holder to prevent leakage, which introduced an additional effect. To account for this, the initial time in the experimental data was redefined. Fig. 5 shows the comparison of simulation and experimental results across different permeability values.

In Fig. 5(a), the simulation results under theoretical permeability conditions ($K = 6.3 \times 10^{-12}$ m$^2$) are shown. The simulation shows noticeable differences from the experimental results, especially during the early stages, suggesting that the theoretical value does not accurately represent the actual permeability. Fig. 5(b) presents simulation results for four different permeability values ($K = 1.0 \times 10^{-11}$, $5.0 \times 10^{-11}$, $8.0 \times 10^{-11}$, $9.0 \times 10^{-11}$), chosen to define the upper and lower bounds of the permeability range. The variation in these curves shows that permeability affects the mass increase. To refine the permeability estimation, a Python script was developed to automatically adjust the permeability parameter and compute the relationship between the square of the mass and time. An error function was defined to determine the permeability value corresponding to the minimum error. The error function is expressed as:

$$ \text{Error} = \sqrt{\frac{1}{n} \sum_{i=1}^{n} \left( \frac{m_{n,i}^2 - m_{e,i}^2}{m_{e,i}^2} \right)^2} \quad (21) $$

where $m_{e,i}^2$ is the squared mass measured in the experiment at time step $i$, $m_{n,i}^2$ is the squared mass from the simulation at the same time step, and $n$ is the total number of time steps (data points every 0.2 s were used for comparison). Fig. 5(c) shows the variation of the defined error as a function of permeability, with the minimum error occurring at around $K = 7.0 \times 10^{-11}$, with an uncertainty of $\pm 0.5 \times 10^{-11}$. Finally,

Fig. 5(d) shows the simulation results corresponding to the minimum error. The simulated curve aligns closely with the experimental data, confirming that the estimated permeability accurately describes the capillary wicking behavior of the cotton fibers.

To nondimensionalize the results, the time evolution of capillary rise height and velocity was further analyzed. Fig. 6(a) shows the variation of capillary rise height $h$ over time, exhibiting a rapid initial rise followed by a gradual approach to equilibrium. The data in this figure is based on the simulation results presented in Fig. 5(d). Fig. 6(b) shows the corresponding rise velocity $dh/dt$, which decreases over time as the driving capillary forces balance with viscous resistance. A threshold velocity of 0.0001 m/s was defined, and the time at which this threshold was reached was determined to be 13.50 s, marking the practical equilibrium of the capillary infiltration process. Based on this, the average velocity over the time interval from 0 to 13.5 s was calculated as 0.002453 m/s.

Using this average velocity, the capillary number $Ca$ was computed, as shown in Eq. (22).

$$ \begin{aligned} Ca &= \frac{\text{viscous forces}}{\text{surface tension}} = \frac{\mu(dh/dt)}{\gamma_L} = \frac{0.402 \times 10^{-3} \times 0.002453}{0.02010} \\ &= 4.91 \times 10^{-5} \end{aligned} \quad (22) $$

The resulting $Ca = 4.91 \times 10^{-5}$ represents the balance between viscous forces and surface tension, indicating that capillary forces dominate infiltration while viscosity affects the rate. The low $Ca$ suggests minimal viscous resistance, enabling rapid liquid uptake. This value serves as a reference for comparing capillary-driven infiltration and optimizing fiber materials for liquid transport [39].

### 4.1.2. Water - cotton fibers

The analysis now transitions to the case where water is used as the infiltrating fluid in cotton fibers. Similar to the simulation conducted for n-heptane, a transient simulation was performed to study the infiltration dynamics of water. Water, as a polar liquid with higher viscosity (0.001 Pa s) and surface tension (72.75 mN/m) compared to n-heptane, is expected to exhibit different infiltration behavior. Besides, the apparent contact angle ($\theta_a$) of water, experimentally measured at around 72.8°, is higher than that of n-heptane (0°), indicating reduced wettability on the fiber surface, which may slow down the infiltration rate. Other relevant parameters are listed in Table 1.

Fig. 7 shows the water saturation (the volume fraction of water) in cotton fibers at various time steps during the simulation, providing a comparative perspective on how fluid properties influence the infiltration process. The saturation variable is dimensionless, ranging from 0 to 1, and the colormap is divided into ten discrete levels. A value of 0 indicates regions that have not yet been wetted by the liquid, while a value of 1 corresponds to complete saturation. At $t = 1$ s, the water saturation front is localized near the bottom boundary, indicating the initial stage of infiltration where water starts to penetrate
<figure><img src="image_12.png" /><figcaption>(a) simulation results under theoretical permeability conditions</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(c) variation of defined error with permeability</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(b) simulation results under different permeabilities</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(d) simulation results at the moment of minimum error</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(a) wicking height *h* vs. time *t*</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>(b) velocity $dh/dt$ vs. time *t*</figcaption></figure>

<caption>**Fig. 6.** Evolution of wicking height and velocity over time (n-heptane).</caption>

the cotton fibers. By $t = 5$ s, the water has advanced further, and a noticeable transition zone between saturated and unsaturated regions appears, highlighting the gradual progression of the liquid front. At $t = 10$ s and $t = 20$ s, the saturation continues to increase, with the transition zone moving upward and becoming less distinct as the fibers absorb more water. By $t = 40$ s, the cotton fibers are nearly fully saturated, showing a uniform saturation distribution across the domain. This progression demonstrates the influence of water's higher viscosity and surface tension than n-heptane, resulting in slower infiltration. These results highlight the importance of fluid properties in influencing infiltration rate and saturation distribution.

For the case of water as the infiltrating liquid, the analysis was carried out following the same methodology as for n-heptane. Transient simulations were performed to compare experimental and simulated results and determine the optimal permeability for water. Fig. 8(a) shows the simulation results under the obtained permeability (in the case of n-heptane) value of $K = 7 \times 10^{-11}$, which was based on the initial assumptions. However, similar to the n-heptane case, a clear mismatch is observed between the simulation and experimental

<figure><img src="image_18.png" /><figcaption>**Fig. 7.** Water saturation in the cotton fibers after 1, 5, 10, 20, and 40 s, ($K = 7.0 \times 10^{-11}$ m$^2$, $c\bar{r} = 6.97$ µm, $\epsilon = 0.53$, $\gamma_L = 72.75$ mN/m, $\theta_a = 72.8^\circ$, $\rho_L = 1000$ kg/m$^3$, $\mu_L = 1.002 \times 10^{-3}$ Pa s).</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>Fig. 8. Comparison of simulation and experimental results across different permeabilities (water).</figcaption></figure>

results, particularly during the initial infiltration phase. This suggests that the input permeability does not align well with the experimental conditions. To address this, Fig. 8(b) presents the simulation results for different permeability values, ranging from $K = 1.0 \times 10^{-11}$ to $K = 9.0 \times 10^{-11}$. This analysis helps establish the lower and upper bounds of the permeability that best represent the experimental data. Next, Fig. 8(c) shows the variation of the defined error (Eq. (21)) with permeability. The error function, calculated as the squared difference between experimental and simulated mass at different time steps, identifies the minimum error at a permeability value of around $K = 6.5 \times 10^{-11}$. This value represents the most accurate permeability for the water infiltration within the cotton fibers. Finally, Fig. 8(d) demonstrates the simulation results at the optimal permeability ($K = 6.5 \times 10^{-11}$).

Similarly, as a reference, the capillary number ($Ca$) was determined based on the simulation data from Fig. 8(d) using Eq. (23). The average velocity was calculated as 0.00117 m/s (as shown in Fig. 9), resulting in a capillary number of $1.61 \times 10^{-5}$.

$$Ca = \frac{\mu(dh/dt)}{\gamma_L} = \frac{1.002 \times 10^{-3} \times 0.00117}{0.07275} = 1.61 \times 10^{-5} \quad (23)$$

Although Fig. 8(d) shows the best agreement between the simulation and experimental results, some discrepancies remain. While the simulated time required for water to fully saturate the cotton fibers aligns well with the experimental findings, the final absorbed liquid mass in the simulation is lower than that observed in the experiments. These discrepancies may come from uncertainties in the estimation of porosity during the experiment. In particular, ensuring that the porous material was uniformly compressed to the intended position may not have been feasible, resulting in a true porosity different from the assumed value of 0.53. To address this, we adjusted the porosity input in the simulation and conducted further analysis.

Fig. 10 compares the simulation and experimental results for different porosity values. The results indicate that when the porosity is set to 0.63, the simulated liquid absorption closely matches the experimental data under the same permeability conditions. However,

while increasing the porosity improves the agreement between simulation and experiment, some deviations persist. This suggests that adjusting a single parameter, such as porosity, may not be sufficient to fully reconcile the simulation with the experimental results. A more comprehensive analysis incorporating additional factors is necessary for a more accurate representation of the physical system.

To better understand these deviations, factors such as structural deformation and transient flow effects must also be taken into account. The deviation observed in the water infiltration results compared to simulation can be attributed to several factors. First, cotton fibers may swell upon absorbing water, which reduces local porosity and effective permeability. Prior studies (e.g., Testoni et al. [40]) have shown that natural fibers such as flax can swell by over 30% in diameter, absorbing up to 39% of their own weight. While swelling was not explicitly included in the model, it may lead to changes in the capillary pathways, particularly in compressed fiber networks. Second, the assumption of uniform pore structure with fixed porosity may not fully capture the local heterogeneity in real cotton samples. Structural irregularities, compounded by swelling, could contribute to non-uniform flow. Finally, although the simulation is time-dependent, it does not include early-stage transient effects such as entrance resistance or inertial acceleration [21], which may influence the dynamics for low-viscosity liquids like water.

### 4.1.3. Glycerol - cotton fibers

Glycerol, compared to n-heptane and water, exhibits a really higher dynamic viscosity of 1.490 Pa s, around 1490 times that of water and 3457 times that of n-heptane. Such a high viscosity creates unique challenges in modeling its infiltration behavior. Specifically, during simulations, the extremely slow flow rate can lead to singularity issues in numerical solutions, which occur due to abrupt gradients or instability in the computed values. To address this, an initial glycerol saturation condition of 0.9 was applied to the cotton fibers. This adjustment avoids sharp discontinuities in the saturation field near the initial time step, ensuring numerical stability.
<figure><img src="image_20.png" /><figcaption>Fig. 9. Evolution of wicking height and velocity over time (water).</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Fig. 10. Comparison of simulation and experimental results across different porosities (water).</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Fig. 11. Glycerol saturation in the cotton fibers after 0, 400, 800, 1200, and 1500 s, ($K = 6.5 \times 10^{-11}$ m², $c\bar{r} = 6.97$ µm, $\epsilon = 0.53$, $\gamma_L = 63$ mN/m, $\theta_a = 69.55^{\circ}$, $\rho_l = 1.490$ Pa s).</figcaption></figure>

Fig. 11 shows the saturation distribution of glycerol in the cotton fibers at five time intervals: 0, 400, 800, 1200, and 1500 s. The saturation variable is dimensionless and ranges from 0.9 to 1, divided into ten discrete levels in the colormap. Although the minimum value is set to 0.9 for numerical reasons, it corresponds to an initially unsaturated state and can be interpreted as a dry region in the simulation context. A value of 1 represents full saturation. At the initial time ($t = 0$ s), the saturation value is set to 0.9 throughout the domain, which serves as a numerical approximation rather than an indication of actual saturation. This initial condition ensures numerical stability during the simulation and can be interpreted as the fibers being unsaturated. Over time, the glycerol front progresses slowly upwards due to its high viscosity, as shown by the gradual color transition from blue (unsaturated) to pink (saturated). By 1500 s, the fibers have achieved near-complete saturation, though the saturation front advancement is much slower compared to n-heptane or water.

To further analyze the case of glycerol, the initial permeability value was estimated based on simulation results from n-heptane and water, around $K = 6.5 \times 10^{-11}$. However, as shown in Fig. 12(a), the simulated curve differs in shape from the experimental results. To explore further, three different permeability values ($K = 1 \times 10^{-11}$, $K = 5 \times 10^{-11}$, $K = 6.5 \times 10^{-11}$) were tested. Among these, the permeability $K = 1 \times 10^{-11}$ showed a curve shape closest to the experimental data, suggesting that the permeability of glycerol may be lower than that of other liquids. Nevertheless, differences remain in the total liquid mass absorbed. Experimentally, the advancing contact angles for glycerol were measured as 68.41°, 71.38°, and 67.71°. The variation in these measurements indicates potential experimental errors, likely due to the sensitivity of contact angle measurements to experimental conditions, such as surface preparation or measurement techniques. For this numerical analysis, an average contact angle of 69.55° was used. Based on these observations, the advancing contact angle was hypothesized to also contribute to the observed differences. Keeping the permeability fixed at $K = 1 \times 10^{-11}$,

the contact angle was varied between 60° and 80°, as shown in Fig. 12(b). The results indicate that the contact angle could influence the infiltration rate and saturation distribution. A smaller contact angle ($\theta_a = 60^{\circ}$) results in faster infiltration, whereas a larger contact angle ($\theta_a = 80^{\circ}$) slows down the liquid absorption process.

As a reference, the data corresponding to $\theta_a = 69^{\circ}$ in Fig. 12(b) was analyzed to determine the average velocity. A threshold velocity of $1 \times 10^{-6}$ m/s was set, and the results indicate that within a 1800 s timeframe, the computed velocity values remained above this threshold. The average velocity was calculated as $6 \times 10^{-6}$ m/s, corresponding to a capillary number of $1.42 \times 10^{-4}$ (Eq. (24)), which is nearly an order of magnitude higher than the values observed for n-heptane ($4.91 \times 10^{-5}$) and water ($1.61 \times 10^{-5}$).

$$C_a = \frac{\mu(dh/dt)}{\gamma_L} = \frac{1.490 \times 6 \times 10^{-6}}{0.063} = 1.42 \times 10^{-4} \quad (24)$$

Unlike the previous cases, detailed error analysis to determine the optimal permeability and contact angle parameters was not performed here. However, this approach remains feasible for future studies. By defining error metrics based on the squared mass curves from simulation and experiments over time, a systematic method for optimizing permeability and advancing contact angles could be developed. The current analysis highlights the combined effects of permeability and contact angle on the infiltration process, emphasizing the complexity of modeling high-viscosity liquids like glycerol in porous media.

## 4.2. Effects of parameters on liquid saturation distribution

In this subsection, the impact of various parameters listed in Table 1 on the flow behavior is thoroughly investigated. Specifically, the effects
<figure><img src="image_23.png" /><figcaption>**Fig. 12.** Comparison of simulation and experimental results across different permeabilities and advancing contact angle (glycerol).</figcaption></figure>

of these parameters on the spatial distribution of liquid saturation within the porous medium and the temporal evolution of absorbed liquid mass are analyzed. The analysis aims to isolate and evaluate the role of each parameter in influencing the infiltration process. To ensure clarity and focus, the baseline physical properties of n-heptane are used as a reference, with only one variable altered at a time for analysis.

### 4.2.1. Liquid surface tension

The influence of liquid surface tension ($\gamma_L$) on the infiltration process was analyzed by considering four different surface tension values: 20.10, 40, 63, and 72.75 mN/m, while keeping all other parameters consistent with the properties of n-heptane. It is important to note that these surface tension values may not reflect actual physical properties but were chosen specifically to isolate and investigate the effect of surface tension on the infiltration behavior.

Fig. 13(a) to (d) use a consistent colormap to represent the liquid saturation (the volume fraction of liquid) in cotton fibers under different surface tension conditions. The saturation variable is dimensionless and ranges from 0 to 1, divided into ten discrete levels. A value of 0 indicates that the region is completely dry, while a value of 1 corresponds to full saturation. The progression of the saturation front across the four subfigures highlights the influence of surface tension on infiltration behavior. Fig. 13(e) presents the simulated evolution of the squared liquid mass absorbed over time. Fig. 13(f) further shows the variation of wicking velocity as a function of time, along with the precise stopping time.

As shown in Fig. 13, the infiltration rate is primarily influenced by capillary pressure rather than solely determined by the capillary number (Ca). Higher surface tension ($\gamma_L$) enhances capillary pressure, significantly accelerating the infiltration process. Although the capillary number (Ca) represents the balance between viscous and capillary forces, numerical data indicate that its value does not exhibit a strictly monotonic relationship with infiltration time. This is because the primary driving force for liquid penetration arises from capillary pressure, which is directly proportional to surface tension. As $\gamma_L$ increases, capillary pressure correspondingly strengthens, leading to faster liquid propagation even when the variation in Ca is minimal. Consequently, the infiltration behavior should be analyzed mainly in terms of capillary pressure rather than relying solely on the capillary number to predict infiltration rates.

### 4.2.2. Advancing contact angle

To analyze the effect of advancing contact angle ($\theta_a$) on the infiltration process in cotton fibers, simulations were conducted under four different contact angles: $\theta_a = 0^\circ, 30^\circ, 60^\circ$, and $80^\circ$, while other parameters were kept consistent with the properties of n-heptane.

The results, presented in Fig. 14, display both the spatial distribution of liquid saturation at various time points and the evolution of the absorbed liquid mass squared over time. In Fig. 14(a) to (d), the liquid

saturation profiles show a clear dependence on the advancing contact angle. The colormap represents the dimensionless liquid saturation, ranging from 0 (dry) to 1 (fully saturated). At $\theta_a = 0^\circ$ (Fig. 14(a)), the liquid infiltrates the cotton fibers most efficiently, achieving nearly complete saturation by $t=10$ s. The absence of a significant contact angle promotes strong capillary forces, driving rapid infiltration. As the contact angle increases to $\theta_a = 30^\circ$ (Fig. 14(b)), the infiltration rate decreases slightly, with the saturation front advancing more gradually over time. At $\theta_a = 60^\circ$ (Fig. 14(c)), the infiltration slows further, as the increased contact angle reduces the effective capillary force acting within the fibers. By $\theta_a = 80^\circ$ (Fig. 14(d)), the liquid's progression is markedly slower, with incomplete saturation observed even after $t=40$ s. The higher contact angle diminishes the wettability of the fiber surface, which restricts the infiltration process. The absorbed liquid mass squared over time, shown in Fig. 14(e), further supports these observations. The black curve ($\theta_a = 0^\circ$) exhibits the fastest rise, reaching saturation quickly, followed by the blue ($\theta_a = 30^\circ$), red ($\theta_a = 60^\circ$), and green ($\theta_a = 80^\circ$) curves, which demonstrate progressively slower infiltration rates. The differences in the curves highlight the sensitivity of the infiltration dynamics to variations in the advancing contact angle.

Fig. 14(f) further shows the evolution of wicking velocity under different advancing contact angles. For $\theta_a = 0^\circ$ (Ca = $5.93 \times 10^{-5}$), infiltration occurs rapidly, with velocity dropping below 0.0001 m/s within 10.8 s. As the contact angle increases to $\theta_a = 30^\circ$ (Ca = $4.56 \times 10^{-5}$), the infiltration rate decreases, extending the cutoff time to $14.0$ s. At $\theta_a = 60^\circ$ (Ca = $1.36 \times 10^{-5}$), the infiltration rate slows significantly, with the cutoff time reaching 39.8 s. However, at $\theta_a = 80^\circ$ (Ca = $1.40 \times 10^{-5}$), while the infiltration rate is the lowest, the cutoff time shortens to 13.5 s due to the definition of the threshold velocity (<0.0001 m/s).

These findings confirm that the advancing contact angle directly influences both the rate and extent of liquid infiltration. Smaller contact angles facilitate rapid and complete saturation, whereas larger contact angles weaken capillary action, slowing infiltration and leaving portions of the fiber structure unsaturated. To further explore the sensitivity of infiltration dynamics to contact angle, a supplementary simulation study was performed across a wider range of $\theta_a$ values. The results, including quantitative trends and the non-linear response of absorbed mass, are presented and discussed in Appendix B.

### 4.2.3. Porosity

To analyze the effect of porosity ($\epsilon$) on the infiltration process, simulations were conducted using four different porosity values: $\epsilon = 0.1, 0.3, 0.7, \text{and } 0.9$. All other parameters were kept consistent with the baseline properties of n-heptane. The results are presented in Fig. 15, which shows both the spatial distribution of liquid saturation at different times and the squared liquid mass absorbed over time.
<figure><img src="image_24.png" /></figure>

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

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

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

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

<figure><img src="image_29.png" /><figcaption>**Fig. 13.** Simulation results of liquid saturation in cotton fibers under four different surface tension conditions, ($K = 7 \times 10^{-11}$ m², $c\bar{r} = 6.97$ µm, $\epsilon = 0.53$, $\gamma_L = 20.10-72.75$ mN/m, $\theta_a = 0°$, $\rho_l = 684$ kg/m³, $\mu_l = 0.402 \times 10^{-3}$ Pa s).</figcaption></figure>

**Fig. 15(a)–(d)** demonstrate the liquid saturation profiles for the selected porosity values. The colormap illustrates the liquid saturation (volume fraction of liquid), which is a dimensionless variable ranging from 0 to 1. A value of 0 indicates completely dry regions, while a value of 1 corresponds to full saturation. At a lower porosity ($\epsilon = 0.1$, the liquid front advances slowly, and the saturation remains relatively low across the fiber domain, even at later times. This slow infiltration is due to the reduced void space available for the liquid to fill, which limits the overall fluid transport. As the porosity increases, the liquid infiltration becomes faster. For $\epsilon = 0.3$ and $\epsilon = 0.7$, the saturation profiles show a more uniform distribution at earlier times, with the liquid front advancing more rapidly compared to $\epsilon = 0.1$. At the highest porosity ($\epsilon = 0.9$), the liquid quickly permeates the fiber domain, achieving near-complete saturation within a much shorter time. This is attributed to the larger void spaces that facilitate higher liquid mobility. **Fig. 15(e)** plots the squared liquid mass absorbed over time for all cases. The curve corresponding to $\epsilon = 0.9$ reaches the highest mass within the shortest time, highlighting the direct relationship between porosity and

liquid absorption capacity. Conversely, the curve for $\epsilon = 0.1$ remains lower throughout the time range, reflecting the limited liquid uptake due to low porosity.

**Fig. 15(f)** further shows the evolution of wicking velocity under different porosity conditions. For $\epsilon = 0.1$ ($Ca = 3.2 \times 10^{-5}$), infiltration occurs at a slower rate, with a cutoff time of 7.0 s, indicating that the limited pore space restricts fluid transport. As porosity increases to $\epsilon = 0.3$ ($Ca = 2.3 \times 10^{-5}$), despite a lower capillary number, the infiltration process is prolonged, with the velocity dropping below the threshold only after 24.8 s, suggesting a more complex interaction between capillary flow and pore structure. When porosity is further increased to $\epsilon = 0.7$ ($Ca = 6.5 \times 10^{-5}$), the infiltration rate improves, reducing the cutoff time to 9.5 s. At the highest porosity, $\epsilon = 0.9$ ($Ca = 7.2 \times 10^{-5}$), liquid transport is the most efficient, with the cutoff time further decreasing to 8.5 s.

These findings confirm that porosity plays a crucial role in governing both the infiltration rate and the overall saturation distribution. Higher porosity enhances liquid transport within the porous medium
<figure><img src="image_30.png" /><figcaption>**Fig. 14.** Simulation results of liquid saturation in cotton fibers under four different advancing contact angles, ($K = 7 \times 10^{-11}$ m², $c\bar{r} = 6.97$ µm, $\epsilon = 0.53$, $\gamma_L = 20.10$ mN/m, $\theta_a = 0°-80°$, $\rho_l = 684$ kg/m³, $\mu_l = 0.402 \times 10^{-3}$ Pa s).</figcaption></figure>

by increasing the available void space, whereas lower porosity restricts infiltration and reduces absorption capacity.

In this subsection, while other parameters such as fluid density and viscosity, have not been explicitly analyzed, a systematic approach has been proposed to evaluate the impact of various factors on the infiltration process. This method enables both qualitative and quantitative assessments of the influence of individual parameter changes on the results. By isolating and examining each parameter individually, the approach provides a clear framework for understanding their effects. This framework can be extended to study additional physical and material properties in future analyses.

## 5. Conclusions

Using a numerical framework combined with experimental data, this research explores capillary-driven flow behavior in fibrous porous media. Through an inverse analysis approach, it determines the permeability values for various liquids and analyzes the effects of multiple

parameters on liquid saturation and mass absorption dynamics. The findings are summarized as follows:

1. The simulations were compared with experimental data, and permeability values for different liquids were determined using an inverse analysis approach. For example, the permeability values for n-heptane and water were estimated as $K = (7.0 \pm 0.5) \times 10^{-11}$ m² and $K = (6.5 \pm 0.5) \times 10^{-11}$ m², respectively. These differences reflect variations in liquid transport within the fibrous porous medium.

2. Based on simulation data, the wicking velocity of different liquids in the porous medium was calculated, followed by the determination of the corresponding capillary numbers. The results show that the capillary numbers for n-heptane, water, and glycerol are $4.91 \times 10^{-5}$, $1.61 \times 10^{-5}$, and $1.42 \times 10^{-4}$, respectively.

3. The effects of liquid surface tension, advancing contact angle, and porosity on infiltration behavior were individually examined through simulations, with the capillary number used for quantitative analysis. Lower $Ca$ values generally correspond to faster
<figure><img src="image_31.png" /><figcaption>(a) liquid saturation, ε=0.1</figcaption></figure>

<figure><img src="image_32.png" /><figcaption>(b) liquid saturation, ε=0.3</figcaption></figure>

<figure><img src="image_33.png" /><figcaption>(c) liquid saturation, ε=0.7</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>(d) liquid saturation, ε=0.9</figcaption></figure>

<figure><img src="image_35.png" /><figcaption>(e) squared liquid mass absorbed over time</figcaption></figure>

<figure><img src="image_36.png" /><figcaption>(f) velocity dh/dt over time</figcaption></figure>

**Fig. 15.** Simulation results of liquid saturation in cotton fibers under four different porosities, (K = 7 × 10⁻¹¹ m², c̃r = 6.97 µm, ε = 0.1–0.9, γL = 20.10 mN/m, θa = 0°, ρl = 684 kg/m³, μl = 0.402 × 10⁻³ Pa s).

liquid absorption, particularly under higher surface tension and smaller contact angles.

In conclusion, this study presents a comprehensive analysis methodology that integrates experimental and simulation approaches, shedding light on the multi-factor interactions driving capillary flow in fibrous porous media. By quantitatively linking interfacial properties (e.g., surface tension and contact angle) to saturation dynamics via the capillary number, this work offers a physicochemical perspective on capillary transport phenomena. These findings contribute to a broader understanding of liquid-solid interactions at the microscale and provide theoretical support for future studies in colloid and interface science. Beyond its practical relevance to material design and process optimization, the work also advances fundamental insights into wetting-driven flow behavior in porous structures.

CRediT authorship contribution statement

**Shaolin Liu:** Writing – review & editing, Writing – original draft, Software, Investigation, Data curation. **Amine Ben-Abdelwahed:** Writing – review & editing, Supervision, Methodology, Funding acquisition, Formal analysis.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A. Example of saturation distribution at a single time step (with axis labels)

To provide a clear representation of the spatial distribution of saturation, Fig. 16 shows the n-heptane saturation (volume fraction of
<figure><img src="image_37.png" /><figcaption>Fig. 16. Spatial distribution of n-heptane saturation (volume fraction) in cotton fibers at t = 1 s.</figcaption></figure>

liquid) in cotton fibers at t = 1 s. This corresponds to one of the selected time points shown in Fig. 3 and is presented here separately with full axis labels for clarity. The horizontal axis (X) and vertical axis (Y) represent the width and height of the computational domain, respectively, with units in meters (m), consistent with the geometry shown in Fig. 2. The colorbar indicates the volume fraction of liquid, which is a dimensionless variable ranging from 0 (dry) to 1 (fully saturated), divided into multiple color levels for improved visual resolution.

In this study, Figs. 3, 4, 7, 11, 13, 14, and 15 display side-by-side subplots at multiple time steps to highlight the temporal evolution of saturation. As the horizontal axis in these composite figures does not correspond to a physical dimension, X-axis labels are omitted. Fig. 16 provides a representative field plot at a single time point to assist readers in understanding the actual spatial dimensions of the saturation field.

## Appendix B. Influence of contact angle variation on capillary infiltration

To further assess the sensitivity of the model to variations in contact angle, a series of numerical simulations were conducted, with the advancing contact angle $\theta_a$ systematically varied from 0° to 80°, while all other parameters were kept constant. The entry capillary pressure boundary condition used in the simulations is defined as:

$$P_{\text{cap}} = (c\bar{r})\epsilon \frac{\gamma_L \cos\theta_a}{4K} \qquad (25)$$

This expression shows that capillary pressure decreases as the contact angle increases, due to the $\cos\theta_a$ term. As $\theta_a$ increases from 0° to 90°, $\cos\theta_a$ decreases monotonically from 1 to 0, which directly reduces the driving force for infiltration.

The simulation results are presented in Fig. 17, showing the evolution of the square of absorbed mass (Mass²) over time for each contact angle. Several key trends can be observed: (1) When the contact angle is below 30°, infiltration behavior remains relatively consistent, with similar absorption rates and final liquid uptake; (2) For angles greater than 40°, the infiltration slows significantly, and a saturation plateau appears earlier in the simulation; (3) At $\theta_a \ge 70°$, the liquid barely penetrates the fiber structure, and the absorbed mass curve rapidly flattens, indicating near non-wetting behavior. Additionally, an increase of just 10° in $\theta_a$ results in noticeably different curve trajectories, highlighting a magnified dynamic effect of the capillary pressure change.

These results demonstrate a nonlinear sensitivity of the infiltration process to contact angle. Although the general trend follows the $\cos\theta_a$ dependence, the effect accumulates over time, leading to substantially inhibited capillary uptake once a threshold angle (around 50°–60°) is exceeded. This finding also suggests that even small local variations in contact angle (e.g., 2°–5°), as may occur in real fibrous materials, can result in differences in local wetting behavior and absorption efficiency.

<figure><img src="image_38.png" /><figcaption>Fig. 17. Simulated evolution of square of absorbed liquid mass for different advancing contact angles $\theta_a$ ranging from 0° to 80°.</figcaption></figure>

## Data availability

Data will be made available on request.

## References

[1] X. Li, L.G. Tabil, S. Panigrahi, Chemical treatments of natural fiber for use in natural fiber-reinforced composites: a review, J. Polym. Environ. 15 (2007) 25–33, http://dx.doi.org/10.1007/s10924-006-0042-3.

[2] A. Keller, Compounding and mechanical properties of biodegradable hemp fibre composites, Compos. Sci. Technol. 63 (9) (2003) 1307–1316, http://dx.doi.org/10.1016/S0266-3538(03)00102-7.

[3] J. Wang, G. Geng, X. Liu, F. Han, J. Xu, Magnetically superhydrophobic kapok fiber for selective sorption and continuous separation of oil from water, Chem. Eng. Res. Des. 115 (2016) 122–130, http://dx.doi.org/10.1016/j.cherd.2016.09.032.

[4] M. Mohammed, R. Rahman, A.M. Mohammed, T. Adam, B.O. Betar, A.F. Osman, O.S. Dahham, Surface treatment to improve water repellence and compatibility of natural fiber with polymer matrix: Recent advancement, Polym. Test. 115 (2022) 107707, http://dx.doi.org/10.1016/j.polymtesting.2022.107707.

[5] M. po Ho, H. Wang, J.-H. Lee, C. kit Ho, K. tak Lau, J. Leng, D. Hui, Critical factors on manufacturing processes of natural fibre composites, Compos. Part B: Eng. 43 (8) (2012) 3549–3562, http://dx.doi.org/10.1016/j.compositesb.2011.10.001.

[6] Y. Wielhorski, M.A. Ben Abdelwahed, L. Bizet, J. Bréard, Wetting effect on bubble shapes formed in a cylindrical T-junction, Chem. Eng. Sci. 84 (2012) 100–106, http://dx.doi.org/10.1016/j.ces.2012.08.008.

[7] M.A.B. Abdelwahed, Y. Wielhorski, L. Bizet, J. Bréard, Bubble formation and transport in T-junction for application to liquid composite molding: Wetting effect, J. Compos. Mater. 48 (1) (2014) 37–48, http://dx.doi.org/10.1177/0021998312467553.

[8] T. Mesogitis, A. Skordos, A. Long, Uncertainty in the manufacturing of fibrous thermosetting composites: A review, Compos. Part A: Appl. Sci. Manuf. 57 (2014) 67–75, http://dx.doi.org/10.1016/j.compositesa.2013.11.004.

[9] M.F. Pucci, P.-J. Liotier, S. Drapier, Capillary wicking in a fibrous reinforcement – orthotropic issues to determine the capillary pressure components, Compos. Part A: Appl. Sci. Manuf. 77 (2015) 133–141, http://dx.doi.org/10.1016/j.compositesa.2015.05.031.

[10] V. Rougier, Modélisation Multi-échelle de l'Imprégnation d'un Milieu Fibreux : Morphologie, Mouillage et Perméabilité (Ph.D. thesis), Normandie Université, 2021, 2021NORMLH02.

[11] S. Liu, A. Ahmadi-Senichault, C. Levet, J. Lachaud, Experimental investigation on the validity of the local thermal equilibrium assumption in ablative-material response models, Aerosp. Sci. Technol. 141 (2023) 108516, http://dx.doi.org/10.1016/j.ast.2023.108516.

[12] F.A. Ogunmokun, R. Wallach, Effect of surfactant surface and interfacial tension reduction on infiltration into hydrophobic porous media, Geoderma 441 (2024) 116735, http://dx.doi.org/10.1016/j.geoderma.2023.116735.

[13] L. Galet, S. Patry, J. Dodds, Determination of the wettability of powders by the Washburn capillary rise method with bed preparation by a centrifugal packing technique, J. Colloid Interface Sci. 346 (2) (2010) 470–475, http://dx.doi.org/10.1016/j.jcis.2010.02.051.

[14] E.W. Washburn, Note on a method of determining the distribution of pore sizes in a porous material, Proc. Natl. Acad. Sci. 7 (4) (1921) 115–116, http://dx.doi.org/10.1073/pnas.7.4.115.
[15] Y. Shiri, S.M.J. Seyed Sabour, Analytical, experimental, and numerical study of capillary rise dynamics from inertial to viscous flow, Phys. Fluids 34 (10) (2022) 102105, http://dx.doi.org/10.1063/5.0111688.

[16] N. Fries, M. Dreyer, An analytic solution of capillary rise restrained by gravity, J. Colloid Interface Sci. 320 (1) (2008) 259–263, http://dx.doi.org/10.1016/j.jcis.2008.01.009.

[17] A. Mohammad Karim, W.J. Suszynski, Physics of dynamic contact line: Hydrodynamics theory versus molecular kinetic theory, Fluids 7 (10) (2022).

[18] M.A. Ben Abdelwahed, Y. Wielhorski, L. Bizet, J. Bréard, Characterisation of bubbles formed in a cylindrical T-shaped junction device, Chem. Eng. Sci. 76 (2012) 206–215, http://dx.doi.org/10.1016/j.ces.2012.04.025.

[19] B. Koohestani, A. Darban, P. Mokhtari, E. Yilmaz, E. Darezereshki, Comparison of different natural fiber treatments: a literature review, Int. J. Environ. Sci. Technol. 16 (2019) 629–642, http://dx.doi.org/10.1007/s13762-018-1890-9.

[20] S. Liu, A. Ben-Abdelwahed, Experimental study of capillary impregnation and wettability effects in porous cotton fiber structures, Colloids Surf. A: Physicochem. Eng. Asp. 708 (2025) 135995, http://dx.doi.org/10.1016/j.colsurfa.2024.135995.

[21] S. Patari, P.S. Mahapatra, Liquid wicking in a paper strip: An experimental and numerical study, ACS Omega 5 (36) (2020) 22931–22939, http://dx.doi.org/10.1021/acsomega.0c02407.

[22] Y. Wielhorski, A.B. Abdelwahed, J. Bréard, Theoretical approach of bubble entrapment through interconnected pores: supplying principle, Transp. Porous Media 96 (1) (2013) 105–116.

[23] M.F. Pucci, P.-J. Liotier, S. Drapier, Tensiometric method to reliably assess wetting properties of single fibers with resins: Validation on cellulosic reinforcements for composites, Colloids Surf. A: Physicochem. Eng. Asp. 512 (2017) 26–33, http://dx.doi.org/10.1016/j.colsurfa.2016.09.047.

[24] S. Liu, A. Ahmadi-Senichault, C. Levet, J. Lachaud, Development and validation of a local thermal non-equilibrium model for high-temperature thermal energy storage in packed beds, J. Energy Storage 78 (2024) 109957, http://dx.doi.org/10.1016/j.est.2023.109957.

[25] P. Raiskinmäki, A. Shakib-Manesh, A. Jäsberg, A. Koponen, J. Merikoski, J. Timonen, Lattice–Boltzmann simulation of capillary rise dynamics, J. Stat. Phys. 107 (2002) 143–158, http://dx.doi.org/10.1023/A:1014506503793.

[26] F.G. Wolf, L.O. dos Santos, P.C. Philippi, Capillary rise between parallel plates under dynamic conditions, J. Colloid Interface Sci. 344 (1) (2010) 171–179, http://dx.doi.org/10.1016/j.jcis.2009.12.023.

[27] G. Lu, X.-D. Wang, Y.-Y. Duan, Study on initial stage of capillary rise dynamics, Colloids Surf. A: Physicochem. Eng. Asp. 433 (2013) 95–103, http://dx.doi.org/10.1016/j.colsurfa.2013.05.004.

[28] P. Zhang, Periodic Phase Separation: A Numerical Study Via a Modified Cahn-Hilliard Equation, Canadian Theses on Microfiche, Simon Fraser University, 2006, URL https://books.google.fr/books?id=HjyC0AEACAAJ.

[29] R. Pan, Y. Yang, B.G. Lougou, L. Wu, W. Wang, Y. Guo, Y. Shuai, Thermal performance evaluation of a novel solar-driven pyrolysis reactor, Energy 313 (2024) 134051, http://dx.doi.org/10.1016/j.energy.2024.134051.

[30] R. Masoodi, K. Pillai, Darcy's law-based model for wicking in paper-like swelling porous media, AIChE J. 56 (2010) 2257–2267, http://dx.doi.org/10.1002/aic.12163.

[31] P. Horgue, C. Soulaine, J. Franc, R. Guibert, G. Debenest, An open-source toolbox for multiphase flow in porous media, Comput. Phys. Comm. 187 (2015) 217–226, http://dx.doi.org/10.1016/j.cpc.2014.10.005.

[32] D. Liu, Multiphase fluid handling in absorbing materials-simulating porous media using OpenFoam, 2018, doi:https://hdl.handle.net/20.500.12380/256606.

[33] E.W. Washburn, The dynamics of capillary flow, Phys. Rev. 17 (1921) 273–283, http://dx.doi.org/10.1103/PhysRev.17.273.

[34] M.F. Pucci, P.-J. Liotier, S. Drapier, Capillary effects on flax fibers – modification and characterization of the wetting dynamics, Compos. Part A: Appl. Sci. Manuf. 77 (2015) 257–265, http://dx.doi.org/10.1016/j.compositesa.2015.03.010.

[35] D.R. Schuchard, J.C. Berg, Liquid transport in composite cellulose—superabsorbent fiber networks, Wood Fiber Sci. (1991) 342–357.

[36] M.A. Ben Abdelwahed, Mécanismes d'Imprégnation en Milieux Fibreux : Modélisation et Application à la Mise en Oeuvre des Matériaux Composites à Fibres Longues (Ph.D. thesis), Université du Havre, 2011, URL https://theses.hal.science/tel-00715952.

[37] R.H. Brooks, A.T. Corey, Properties of porous media affecting fluid flow, J. Irrig. Drain. Div. 92 (2) (1966) 61–88.

[38] B. Gebart, Permeability of unidirectional reinforcements for RTM, J. Compos. Mater. 26 (8) (1992) 1100–1133, http://dx.doi.org/10.1177/002199839202600802.

[39] B.E. Rapp, Chapter 9 - fluids, in: B.E. Rapp (Ed.), Microfluidics: Modelling, Mechanics and Mathematics, in: Micro and Nano Technologies, Elsevier, Oxford, 2017, pp. 243–263, http://dx.doi.org/10.1016/B978-1-4557-3141-1.50009-5.

[40] G.A. Testoni, S. Kim, A. Pisupati, C.H. Park, Modeling of the capillary wicking of flax fibers by considering the effects of fiber swelling and liquid absorption, J. Colloid Interface Sci. 525 (2018) 166–176, http://dx.doi.org/10.1016/j.jcis.2018.04.064.