# Particulate transport in porous media at pore-scale.
## Part 2: CFD-DEM and colloidal forces

Laurez Maya Fogouang, Laurent André, Philippe Leroy, Cyprien Soulaine

► To cite this version:

Laurez Maya Fogouang, Laurent André, Philippe Leroy, Cyprien Soulaine. Particulate transport in porous media at pore-scale. Part 2: CFD-DEM and colloidal forces. Journal of Computational Physics, 2024, 519, pp.113439. 10.1016/j.jcp.2024.113439. insu-04729488

HAL Id: insu-04729488

https://insu.hal.science/insu-04729488v1

Submitted on 10 Oct 2024

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
# Particulate transport in porous media at pore-scale.
## Part 2: CFD-DEM and colloidal forces

Laurez Maya Fogouang<sup>a,b,*</sup>, Laurent André<sup>a,b</sup>, Philippe Leroy<sup>b</sup>, Cyprien Soulaine<sup>a,**</sup>

ARTICLE INFO

Keywords:

DLVO theory

JKR theory

Colloid deposition

Pore-scale modeling

CFD-DEM model

Pore-clogging

A B S T R A C T

Pore-clogging by aggregation of fine particles is one of the key mechanisms in particulate transport in porous media. In this work, the unresolved-resolved four-way coupling CFD-DEM (Computational Fluid Dynamics - Discrete Element Method) proposed in Part 1 is coupled with colloidal forces (long-range interactions) to model the transport of charged particles and retention by aggregation at the pore-scale. The model includes hydro-mechanical interactions (e.g. collision, drag, buoyancy, gravity) and electrochemical interactions (e.g. Van der Waals attraction, electrostatic double layer repulsion) between the particles, the fluid, and the porous formation. An adhesive contact force based on the Johnson-Kendall-Roberts theory allows for realistic particle adhesion on the walls. The model robustness is verified using reference semi-analytical solutions of the particle dynamics including long-range interactions. Finally, our CFD-DEM for particulate transport including DLVO and JKR adhesive contact forces is used to investigate the effect of fluid salinity on pore-clogging and permeability reduction. Importantly and unlike other approaches, our CFD-DEM model is not constrained by the size of the particle relative to the cell size. Our pore-scale model offers new possibilities to explore the impact of various parameters including particle size distribution, particle concentration, flow rates, and pore geometry structure on the particulate transport and retention in porous media.

1. Introduction

The transport of fine particles and colloids in porous media eventually reduces permeability due to particle deposition and pore-clogging. In most industrial processes involving particulate flow in porous media, clogging is synonymous with a drop in process efficiency and operating times [1–5]. Previous studies have identified three main mechanisms involved in particle clogging: (i) sieving, (ii) bridging, and (iii) aggregation of particles [6,7]. Sieving refers to size exclusion, and bridging refers to arch formation (bridge-like structures of particles). Both mechanisms depend mainly on hydromechanical interactions and their modeling has been discussed in Part 1 of this two-part paper [8]. Aggregation of particles, on the other hand, relies on electrochemical long-range interactions (repulsive or attractive) between the porous matrix surfaces, the particles, and the carrier fluids. Aggregation of particles refers to the successive deposition of particles which reduces the effective pore diameters until clogging [9–11]. In the two last decades,

<footnote><sup>a</sup> Institut des Sciences de la Terre d'Orléans (ISTO), UMR7327, Université d'Orléans, CNRS, BRGM, F-45071, Orléans, France</footnote>

<footnote><sup>b</sup> BRGM, F-45060, Orléans, France</footnote>

<footnote>* Corresponding author at: Institut des Sciences de la Terre d'Orléans (ISTO), UMR7327, Université d'Orléans, CNRS, BRGM, F-45071, Orléans, France.</footnote>

<footnote>** Principal corresponding author.</footnote>

<footnote>E-mail addresses: laurez.maya@cnrs-orleans.fr (L. Maya Fogouang), cyprien.soulaine@cnrs-orleans.fr (C. Soulaine).</footnote>
<table><caption>**Table 1**
Features of the resolved, semi-resolved [54], unresolved and unresolved-resolved CFD-DEM. DNS = Direct Numerical Simulation, FVM = Finite Volume Method.</caption><thead><tr><th></th><th>Resolved<br>CFD-DEM</th><th>Semi-resolved<br>CFD-DEM</th><th>Unresolved<br>CFD-DEM</th><th>Unresolved-resolved<br>CFD-DEM</th></tr></thead><tbody><tr><td>Particle surface</td><td>Resolved</td><td>Not resolved</td><td>Not resolved</td><td>Resolved for Δx/D<sub>p</sub> ≤ 1 and<br>unresolved for Δx/D<sub>p</sub> > 1</td></tr><tr><td>Particle-fluid force</td><td>Particle-resolved DNS</td><td>Drag force model</td><td>Drag force model</td><td>Combined particle-resolved DNS<br>and drag force model</td></tr><tr><td>Background velocity</td><td>Fluid velocities resolved<br>in FVM cells</td><td>Fluid velocities resolved<br>in neighboring FVM cells</td><td>Fluid velocities in the<br>local FVM cell only</td><td>Fluid velocities calculated in FVM<br>cells around resolved particles and<br>in local FVM cells for unresolved<br>particles</td></tr><tr><td>Cell to particle ratio Δx/D<sub>p</sub></td><td>&lt; 1/10</td><td>≈ 1</td><td>&gt; 3</td><td>No restriction</td></tr></tbody></table>

laboratory experiments (e.g. Bigna et al. [12], Oyeneyin et al. [13], Agbangla et al. [14], Ahfir et al. [15], Fetzer et al. [16], Jung et al. [17], Delouche et al. [18]) and numerical simulations (e.g. Ryan and Elimelech [19], Sefrioui et al. [20], Samari-Kermani et al. [21], Ramezanpour et al. [22], Li et al. [23]) have been made at the core-scale and pore-scale to investigate the migration of fine particles and colloids, and look at their impact on clogging in porous media. However, an efficient pore-scale model that captures sieving, bridging, and aggregation for various flow and water chemistry conditions is still missing.

Pore-scale studies are key to investigate clogging phenomena in which fluid-solid interface plays a major role. Indeed, at this scale, the porous media geometry is known, which enables an accurate description of fluid-solid interactions. Several methods are used to simulate particulate flow at the pore-scale using Navier-Stokes-based approaches. For example, the Euler-Euler approach models the cloud of particles as a continuum phase behaving as a non-miscible fluid with the carrier fluid [24,25]. In this approach, however, fluid-particle, particle-wall, and particle-particle interactions are not resolved but modeled using constitutive laws. Computational Fluid Dynamics-Discrete Element Method (CFD-DEM) [26,27] belongs to the Eulerian-Lagrangian family in which particles are seen as a Lagrangian discrete phase transported by fluids described by Navier-Stokes equations solved on an Eulerian grid [28,29]. A review of the different Euler-Lagrange methods including Lattice-Boltzmann-Method-DEM and Smooth Particle Hydrodynamics is found in Part 1.

Four-way CFD-DEM coupling accounts for the feedback of particles on the flow as well as particle-particle and particle-wall interactions. This enables the simulation of dense suspensions and overcomes a standard limitation due to the assumption of diluted suspension. If CFD-DEM simulations can capture both sieving and bridging mechanisms [30–33], standard resolved (i.e. particles larger than the size of the cells) and unresolved (i.e. particles smaller than the size of the cells) coupling approaches are limited for aggregation due to particle size relative to computational cells size [34,35]. Indeed, on the one hand, unresolved approaches miscapture the flow dynamics around the particles that guide their aggregation or remobilization. On the other hand, resolved approaches are computationally intensive with the increase of particle number. The unresolved-resolved CFD-DEM approach introduced in Part 1 [8] combined with an efficient cell search strategy overcomes these limitations and allows for simulations with a larger number of particles whose size is independent of the CFD grid refinement (see Table 1). Modeling colloidal transport leading to clogging by aggregation in porous media, however, requires additional forces and appropriate treatment of the colloid adhesion [36] that are not included in the unresolved-resolved approach.

Colloidal forces arise if particles are microns and submicrometer in size and electrochemical interactions become dominant [37]. These forces are due to long-range interactions whose influence applies from 0.1 µm down to nearly surface-to-surface contact. They are classically described using the Derjaguin-Landau-Verwey-Overbeek (DLVO) theory [38,39] that combines London-Van der Waals attraction and electrical double-layer (EDL) repulsion. DLVO forces are very sensitive to salinity and pH due to electrostatic forces depending on the surface electrical potential commonly assumed to be similar to the zeta potential [40]. Crystal structure and related electrochemical reactions at the surface of the particle in contact with water are responsible for the dependency of surface electrical potential on water chemistry [41].

DLVO forces are implemented in CFD-DEM by modifying the Lagrangian force balance [42,21]. Because the DLVO theory is based on the direct collision hypothesis, the force and torque caused by collisions are missing. Therefore when the separation distance between the particles approaches zero, there is a numerical singularity problem due to the infinitely attractive DLVO potential. This is a major concern for modeling pore-clogging by aggregation of particles because adhesive contact is not handled properly which leads to unphysical particle-particle and particle-wall interpenetrations. To circumvent this issue, Abbasfard et al. [43] introduces a cut-off distance between the surfaces below which DLVO forces vanish. The selection of an adequate cut-off length is challenging because it is system-dependent, and the particle displacement between time steps should be less than the cut-off distance to avoid numerical errors [44]. Another approach uses the Born repulsion – a highly repulsive short-range force – introduced by Ruckenstein and Prieve [45] to avoid the interpenetration of particle-particle and particle-wall [46–48]. Parametrization of the Born repulsion force requires the precise knowledge of the system’s atomic collision diameter (the finite intermolecular separation distance at which the Lennard-Jones potential associated with Van der Waals interactions are equal to zero). Moreover, the Born repulsion force becomes dominant if the separation distance is less than 1 nm [49]. With such an approach, however, small enough time steps are required to capture particle displacement under 1 nm of separation distance to avoid highly non-physical repulsion of the particles in contact. Alternatively, the Johnson-Kendall-Roberts (JKR) theory is a physically-rooted adhesion contact model known to accurately predict
the contact area [50]. JKR theory considers the interaction between surface energy on particles and the material surface. It does not rely on fitting parameters. It has been used already in some Euler-Lagrange approaches combined with DLVO theory for simulating fine migration [51–53], but not in the context of resolved and unresolved coupling approaches, which is the main purpose of this paper.

In this work, we implement an unresolved-resolved four-way coupling CFD-DEM combined with DLVO and the adhesive JKR contact forces to achieve pore-scale simulations of particulate flow through porous media including particle migration, deposition, and retention for various salinity conditions. The paper is organized as follows. In section 2, we introduce the governing equations of the four-way coupling CFD-DEM including a description of the DLVO forces and adhesive contact force and torque applied on particles. In section 3, we show the consistency and accuracy of the numerical implementation of the DLVO forces on colloid dynamics. In section 4, we apply our new modeling technique to the pore-scale simulation of particle-clogging by aggregation in a porous medium. Finally (section 5), we end with a summary and concluding remarks.

## 2. Computational model

In this section, we present the implementation of DLVO interactions within the unresolved-resolved four-way coupling CFD-DEM proposed in Maya et al. [8]. First, we briefly describe the CFD approach to model fluid flow (Section 2.1). Then, we describe our DEM approach including long-range interaction forces to simulate the displacement of colloidal particles (Section 2.2).

### 2.1. Fluid motion in the CFD Eulerian grid

Fluid flow in a fixed grid is modeled by solving incompressible Volume-Averaged Navier-Stokes (VANS) equations [55,56]. In this approach, the particle presence on the computational grid is described by a porosity field, $\epsilon$, defined as,

$$ \epsilon = \begin{cases} 1, & \text{if the cell is occupied by fluids only,} \\ 0, & \text{if the cell contains a fluid-solid aggregate or a fluid-solid interface,} \\ 0, & \text{if the cell is occupied by solids only.} \end{cases} \tag{1} $$

The mass balance equation for the fluid phase reads,

$$ \frac{\partial(\epsilon\rho_f)}{\partial t} + \nabla \cdot (\epsilon\rho_f \mathbf{v}^f) = 0, \tag{2} $$

where $\rho_f$ is the fluid density, and $\mathbf{v}^f$ is the cell-averaged fluid velocity.

The fluid momentum balance equation is:

$$ \frac{\partial(\epsilon\rho_f \mathbf{v}^f)}{\partial t} + \nabla \cdot (\epsilon\rho_f \mathbf{v}^f \mathbf{v}^f) = -\epsilon\nabla p + \epsilon\rho_f \mathbf{g} + \epsilon\nabla \cdot (\mu_f (\nabla \mathbf{v}^f + (\nabla \mathbf{v}^f)^T) - \epsilon^2 \mu_f (\mathbf{v}^f - \bar{\mathbf{v}}^f)/K, \tag{3} $$

where $p$ is the fluid pressure, $\mu_f$ is the fluid dynamic viscosity, $\mathbf{g}$ is the gravitational acceleration, and $\bar{\mathbf{v}}^f$ is the averaged particle velocity on the Eulerian grid, and $K$ is the cell-permeability. The last term of the right-hand side of Eq. (3) corresponds to the flow resistance due to the presence of particles within a computational cell. The description of this term whether the particle is larger (resolved) or smaller (unresolved) than the cell size is found in Part 1 [8].

### 2.2. Particle motion including colloidal interactions

The motion of a cloud of particles is solved in a Lagrangian frame using a four-way coupling Discrete Element Method (DEM) that accounts for fluid-particle, particle-particle, and particle-wall interactions. The total velocity, $\mathbf{v}_i^p$, of particle i reads,

$$ \mathbf{v}_i^p = \mathbf{U}_i^p + \boldsymbol{\omega}_i \times \mathbf{r}_i, \tag{4} $$

where $\mathbf{U}_i^p$, $\boldsymbol{\omega}_i$, and $\mathbf{r}_i$ are the translational, the angular velocity, and the position vector of particle i, respectively.

A spherical particle i of mass $m_i$, and moment of inertia $I_i$ is in contact with $n_i^c$ other particles and walls. Because electrochemical forces are long-range interaction forces, they apply to $n_i^{nc}$ objects (particles and walls) that are not necessarily in contact with particle i but in its close neighborhood. The force balance equation reads,

$$ m_i \frac{d\mathbf{U}_i^p}{dt} = \sum_j^{n_i^c} \mathbf{F}_{ij}^c + \sum_j^{n_i^{nc}} \mathbf{F}_{ij}^{nc} + \mathbf{F}_i^f + \mathbf{F}_i^g, \tag{5} $$

and,

$$ I_i \frac{d\boldsymbol{\omega}_i}{dt} = \sum_j^{n_i^c} \mathbf{M}_{ij}^c + \sum_j^{n_i^{nc}} \mathbf{M}_{ij}^{adh} + \mathbf{M}_i^{hyd}, \tag{6} $$
<figure><img src="image_2.png" /><figcaption>Fig. 1. The range of interaction of a particle (in green) is determined by a kernel (in purple) centered on the particle centroid and whose radius is equal to the particle diameter. The particles (in red) overlapping cells within this orbit are eligible for collision and long-range interactions. (For interpretation of the colors in the figure(s), the reader is referred to the web version of this article.)</figcaption></figure>

where $M_{ij}^c$, $M_i^{hyd}$, $F_i^f$, $F_i^g$ and $F_{ij}^c$, are the contact torque, the hydrodynamic torque, the particle-fluid interaction force, the gravitational force, and the contact force acting on particle $i$ at time $t$, respectively. These forces and torques are thoroughly described in Part 1 [8]. Long-range interaction forces are considered through $F_{ij}^{fc}$, the non-contact forces, and $M_{ij}^{adh}$, the adhesive torque acting on particle $i$ in interaction with $j = 1, 2, ..., n_i^{nc}$ objects (particles and walls) at time $t$, respectively. To reduce the computational costs, we consider that only the objects $j$ (particles and walls) within a kernel centered on the centroid of particle $i$ and with a radius equal to the particle diameter are candidates for particle-particle and particle-wall interactions. If an object $j$ covers a cell that is overlapping with that radius (see Fig. 1), the object $j$ is added to the interaction list of particle $i$, and long-range interaction forces are applied until they collide with each other or move away.

For particulate flow in a charged solution, adsorption of high concentrations of ions and water molecules at the particle surfaces creates a viscous thin layer known as the Stern layer (see Fig. 2a), where the ion concentration gradient is constant [41]. Therefore, the long-range interactions behave differently whether they occur inside or outside the Stern layer. The non-contact forces are, therefore, defined as,

$$F_{ij}^{nc} = \begin{cases} F_{ij}^{DLVO}, & \text{if the interactions take place outside the Stern layer,} \\ F_{ij}^{JKR}, & \text{if the interactions take place within the Stern layer,} \end{cases} \qquad (7)$$

where $F_{ij}^{DLVO}$ and $F_{ij}^{JKR}$ are the DLVO forces and the adhesive JKR force, respectively. These forces and torque are described in the following sections.

### 2.2.1. Long-range interactions

The long-range interactions outside the Stern layer are modeled using the well-established DLVO theory that consists of a combination of an attractive and a repulsive force (see Fig. 2b). The first, known as the London-Van der Waals attraction force, is caused by correlations in the fluctuating polarizations of the electron clouds around nearby atoms and molecules. This force is usually described as a combination of the London dispersion force between instantaneously induced dipoles, the Debye force between permanent dipoles and induced dipoles, and the Keesom force between permanent molecular dipoles whose rotational orientations are dynamically averaged over time [57]. The second, referred to as the electrostatic repulsion force, is related to the interactions between charged surfaces with the generation of a repulsive EDL (when the surface charges of the two interacting surfaces have the same sign) composed typically of a diffuse layer and a Stern layer (see Fig. 2a). The diffusive layer – also called the Gouy-Chapman layer – is the region where the ions are distributed under the action of electrical forces and thermal motions [58]. The electrostatic repulsion force occurs between charged objects across liquids and its strength increases with the magnitude of the electrical surface potential commonly considered to be equal to the zeta potential (the electrical potential located at the shear or slipping plane [59].).

DLVO forces applied on particle $i$ derive from the Van der Waals and electrostatic potentials,

$$F_{ij}^{DLVO} = -\frac{d}{dh_{ij}} (V_{ij}^{VDW} + V_{ij}^{EDL}) \mathbf{n}_{ij}, \qquad (8)$$

where $h_{ij}$ is the smallest surface-to-surface separation distance, $\mathbf{n}_{ij}$ is the unit vector normal to object $j$ surface and pointing to particle $i$ centroid. Particle-particle and particle-wall potentials have different formulations denoted with the superscript “pp” and “pw” in the following.
<figure><img src="image_3.png" /><figcaption>Fig. 2. (a) Illustration of the electrical double layer (EDL) around a negatively charged particle immersed in a 1:1 electrolyte like NaCl containing Na⁺ and Cl⁻ ions. (b) Typical profile of DLVO potential.</figcaption></figure>

The Van der Waals attraction potential between two spherical particles with smooth surfaces was approximated in 1937 by Hamaker [60], using London equation for the dispersion interaction energy between atoms/molecules as the starting point [61]. It reads,

$$V_{ij}^{\text{EDL}} = -\frac{A_{iLj}}{6} \left( \frac{2R_i R_j}{d_{ij}^2 - (R_i + R_j)^2} + \frac{2R_i R_j}{d_{ij}^2 - (R_i - R_j)^2} + \ln \left( \frac{d_{ij}^2 - (R_i + R_j)^2}{d_{ij}^2 - (R_i - R_j)^2} \right) \right), \quad (9)$$

where $A_{iLj}$ is the Hamaker constant of the particle i, the particle j, and the liquid medium L, expressed as $A_{iLj} = A_{ij} + A_{LL} - A_{iL} - A_{jL}$. $d_{ij} = h_{ij} + R_i + R_j$ is the center-to-center separation distance between particles, $R_i$ and $R_j$ are the radius of particle i and particle j, respectively. If the spheres are very close to each other, $h_{ij} \ll R_i, R_j$, and using the equivalent radius, $r_{ij} = R_i R_j / (R_i + R_j)$, the attraction potential becomes,

$$V_{ij}^{\text{EDL}} = -\frac{A_{iLj} r_{ij}}{6h_{ij}}. \quad (10)$$

The electrostatic repulsion potential between two spherical charged particles dispersed in a solvent is calculated by [39],

$$V_{ij}^{\text{EDL}} = 64\pi\epsilon_0\epsilon_r r_{ij} \left(\frac{k_B T}{Ze}\right)^2 \tanh\left(\frac{Ze\psi_i}{4k_B T}\right) \tanh\left(\frac{Ze\psi_j}{4k_B T}\right) e^{-\kappa h_{ij}}, \quad (11)$$

where $\epsilon_0\epsilon_r$ is the dielectric constant of the solvent, $\psi_i$ is the surface potential of the particle i, $\kappa^{-1} = \sqrt{\frac{\epsilon_0\epsilon_r k_B T}{2 \times 10^3 N_A e^2 I_S}}$ is the Debye screening length corresponding to the EDL thickness for a given electrolyte, $k_B$ is the Boltzmann's constant, $T$ is the absolute temperature, $Z$ is the valence of the electrolyte, $e$ is the charge of the electron, $N_A$ is the Avogadro's number, and $I_S$ is the electrolyte ionic strength. Eq. (11) is based on the constant surface potential approximation (no ion adsorption or condensation on the surface), with the surface potential taken as the potential at the slipping plane. In this paper, we consider particles of identical material, which leads to $\psi_i = \psi_j = \psi$. Therefore, the electrostatic potential becomes,

$$V_{ij}^{\text{EDL}} = 64\pi\epsilon_0\epsilon_r r_{ij} \left(\frac{k_B T}{Ze}\right)^2 \tanh\left(\frac{Ze\psi}{4k_B T}\right) e^{-\kappa h_{ij}}. \quad (12)$$

Finally, the particle-particle DLVO force reads,

$$F_{ij}^{\text{DLVO}} = \left( -\frac{A_{iLj} r_{ij}}{6h_{ij}^2} + 64\pi\epsilon_0\epsilon_r r_{ij}\kappa \left(\frac{k_B T}{Ze}\right)^2 \tanh\left(\frac{Ze\psi}{4k_B T}\right)^2 e^{-\kappa h_{ij}} \right) \mathbf{n}_{ij}. \quad (13)$$

To model the long-range interactions between a sphere ($R_i, \psi_i$) and a wall ($R_j \to \infty, \psi_j$), we use the retarded London-Van der Waals attraction potential [62] and the electrostatic repulsion potential proposed by Hogg et al. [63]. The attraction potential reads,
$$V_{ij}^{DW_{pw}} = - \frac{A_{iLj} R_i}{6h_{ij} \left(1 + 14 \frac{h_{ij}}{\lambda}\right)}, \qquad (14)$$

where $\lambda$ is the characteristic wavelength of the interaction (retardation length assumed to be around 100 nm). The retarded London-Van der Waals potential is an adequate approximation for separation distance, $h_{ij}$, up to $R_i/5$ [64].

The particle-wall electrostatic repulsion potential based on constant surface potential approximation is,

$$V_{ij}^{EDL_{pw}} = \pi \epsilon_0 \epsilon_r R_i (\psi_i^2 + \psi_j^2) \left( \frac{2\psi_i \psi_j}{\psi_i^2 + \psi_j^2} \ln \left( \frac{1+e^{-\kappa h_{ij}}}{1-e^{-\kappa h_{ij}}} \right) + \ln(1-e^{-2\kappa h_{ij}}) \right). \qquad (15)$$

Finally, the DLVO force for particle-wall interactions is,

$$\mathbf{F}_{ij}^{DLVO_{pw}} = \left( - \frac{A_{iLj} R_i \left( 1 + 28 \frac{h_{ij}}{\lambda} \right)}{6h_{ij}^2 \left( 1 + 14 \frac{h_{ij}}{\lambda} \right)^2} + 2\pi \epsilon_0 \epsilon_r R_i K \left( \frac{2\psi_i \psi_j e^{-\kappa h_{ij}} - \left( \psi_i^2 + \psi_j^2 \right) e^{-2\kappa h_{ij}}}{1 - e^{-2\kappa h_{ij}}} \right) \right) \mathbf{n}_{ij}. \qquad (16)$$

### 2.2.2. Adhesive force

DLVO potential has an infinite primary well due to the Van der Waals contribution if the distance between two particle surfaces approaches zero (see Fig. 2b). In CFD-DEM, this results in a large inter-penetration of the particles that is unphysical and unrealistic. To avoid the latter and describe accurately adhesive particle contacts, we modified the Hertzian contact [65,26] by adding an elastic solid-body behavior based on the Johnson-Kendall-Roberts (JKR) theory [50]. This was achieved by modifying the normal overlapping distance $|\delta_n|$ of the spring-sliding-dashpot model described in Part 1 [8].

In the JKR theory, a sphere in contact with an object (particle or wall) under the action of an external force deforms irreversibly (soft sphere model). It means that a finite contact area remains even if the external force vanishes. The JKR method consists of applying a constant attractive force (see Fig. 3b) based on the Derjaguin approximation [38] when the surface separation distance is smaller than the Stern layer thickness, $\sigma_{St}$ (a few nanometers maximum). Within the Stern layer ($h \le \sigma_{St}$), the JKR adhesive force of a particle $i$ in interaction with an object $j$ reads [66],

$$\mathbf{F}_{ij}^{JKR} = -\frac{3}{2}\pi r_{ij} W_{iLj} \mathbf{n}_{ij}, \qquad (17)$$

where $W_{iLj} = \gamma_{ij}^L - \gamma_{iL} - \gamma_{jL}$ is the surface energy or adhesion energy, and $\gamma$ is the interfacial energy.

The JKR adhesive potential is obtained by a linear interpolation between the potential at the wall ($h_{ij} = 0$) and the potential at the slipping plane ($h_{ij} = \sigma_{St}$),

$$V_{ij}^{JKR} = V_{St}^{DLVO} + |\mathbf{F}_{ij}^{JKR}| (\sigma_{St} - h_{ij}), \qquad (18)$$

where $V_{St}^{DLVO}$ is the DLVO potential at $h_{ij} = \sigma_{St}$.

The adhesive JKR force acting on particle $i$ is an attractive force that produces a slight overlap, $\zeta$, with object $j$, and, therefore, a finite contact radius $a$ (see Fig. 3a). The adhesive overlapping distance and the equilibrium contact radius are given by [50]

$$\zeta = \frac{a^2}{r_{ij}}, \quad \text{and,} \quad a = \sqrt[3]{\frac{6\pi r_{ij}^2 W_{iLj}}{k_{ij}}}, \qquad (19)$$

where $k_{ij} = \frac{4}{3\pi E_{ij}}$ and $E_{ij} = \left( \frac{1-\nu_i^2}{E_i} + \frac{1-\nu_j^2}{E_j} \right)^{-1}$ are the equivalent elastic constant and the equivalent Young modulus respectively, and $\nu$ refers to Poisson ratio. Note that, if object $j$ is a wall ($R_j \to \infty$), then the equivalent radius is $r_{ij} = R_i$.

When the contact appears (i.e. the contact overlap, $|\delta_n| > 0$), two configurations exist in the DEM calculation: (i) if $|\delta_n| \le \zeta$, the Hertzian contact force is not considered which leads to an elastic attachment of the particle and the object through the adhesive contact. (ii) If $|\delta_n| > \zeta$, then the rigid Hertzian contact force is applied with a new overlap distance as $|\delta_n^*| = |\delta_n| - \zeta$.

### 2.2.3. Adhesive torque

Adhesion of a particle to a wall or another particle occurs if the distance is within the range of the primary or secondary minimum for the interaction potential. The adhesive physicochemical forces, $\mathbf{F}_{ij}^{JKR}$ and $\mathbf{F}_{ij}^{DLVO}$, involved in the attachment process generate a torque that resists the particle detachment when hydrodynamic forces are applied. The adhesive – or resisting – torque is [67],

$$\mathbf{M}_{ij}^{adh} = \begin{cases} a\mathbf{t}_{ij} \times \mathbf{F}_{ij}^{JKR}, & \text{for adhesion at the primary minimum,} \\ l_x \mathbf{t}_{ij} \times \mathbf{F}_{ij}^{DLVO}, & \text{for adhesion at the secondary minimum,} \end{cases} \qquad (20)$$

where $\mathbf{t}_{ij} = \mathbf{v}_{ij}^p / |\mathbf{v}_{ij}^p|$ is the tangential unit vector, $\mathbf{v}_{ij}^p$ is the tangential slip velocity (see Maya et al. [8]). At the secondary minimum, the DLVO force characterizes the force acting on a lever arm $l_x$ that must be overcome to detach the particle. Since there is no overlapping area at the secondary minimum, $l_x$ is expressed as, [67]
<figure><img src="image_4.png" /><figcaption>Fig. 3. (a) Illustration of the JKR adhesion seen as an overlapping of two elastic spheres. (b) Profile of the potential function and forces including JKR adhesion for distance $0 \le h_{ij} < \sigma_{st}$.</figcaption></figure>

$$l_x = \sqrt[3]{\frac{r_{ij} |\mathbf{F}_{ij}^{\mathrm{DLVO}}|}{4k_{ij}}}. \qquad (21)$$

## 2.3. Numerical implementation

The electrochemical forces described in this section are implemented in the OpenFOAM-based unresolved-resolved four-way coupling CFD-DEM. The algorithmic procedure is illustrated in Fig. 4. It relies on OpenFOAM version 9 (https://www.openfoam.org). In terms of the coupling procedure, we use distinct time-steps, $\Delta t_{DEM}$ and $\Delta t_{CFD}$, for the particle dynamics and the fluid flow, respectively. As the constraints on the particle dynamics are stronger than on the fluid flow, we use sub-cycling of $\Delta t_{DEM}$ within a CFD time-step to capture particle collision events accurately. Additional details about the unresolved-resolved CFD-DEM implementation are found in Maya et al. [8].

## 3. Model validations

In this section, we verify the robustness of our CFD-DEM including DLVO forces using test cases for which we derived semi-analytical solutions. First, we verify the implementation of DLVO forces for particle-particle (section 3.1) and particle-wall (section 3.2) interactions separately. Then, we compare simulations of the formation of particle aggregates with reference experimental data to verify the coupling between the unresolved-resolved CFD-DEM model and the DLVO forces (section 3.3).

### 3.1. Particle-particle DLVO interactions

We consider two identical particles of density $\rho_p = 1050 \text{ kg/m}^3$ (lower than the average density of minerals around 2500 – 3000 kg/m}^3$), radius $R_1 = R_2 = R = 1 \text{ µm}$, and mass $m = 4.4 \times 10^{-15} \text{ kg}$ initially separated from each other by $h_0$ as presented in Fig. 5. Both particles are initially immobile and immersed in a brine (water + NaCl) solution. Table 2 summarizes particle and fluid properties. There is no flow, and the drag and sedimentation forces are neglected so that only the DLVO forces are involved in the particle displacement. The evolution of the surface-to-surface separation distance, $h(t)$, is computed using our CFD-DEM approach in three cases: pure attraction, pure repulsion, and mixed attraction-repulsion. Results are compared with semi-analytical solutions that are derived in the following.

The CFD-DEM computational domain is a $50 \text{ µm} \times 50 \text{ µm}$ square with a regular grid cell size of $1 \text{ µm} \times 1 \text{ µm}$. We use $\Delta t = 2.5 \times 10^{-6} \text{ s}$ and $h_0 = 2 \text{ µm}$ in the attraction and mixed attraction-repulsion cases, and $\Delta t = 10^{-8} \text{ s}$ and $h_0 = 0.04 \text{ µm}$ for the repulsive case.

Alternately, the evolution of the separation distance can also be solved using Ordinary Differential Equations (ODE). Indeed, from Eq. (5) and Eq. (13), the particle motion is governed by
<figure><img src="image_5.png" /><figcaption>Fig. 4. Flowchart for the unresolved-resolved four-way coupling CFD-DEM numerical procedure including colloidal forces.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. Representation of the two particles distant from h.</figcaption></figure>

<caption>Table 2</caption>

Fluid, wall, and particle properties used in the model verification simulations based on Yu et al. [68].

<table><thead><tr><th>Parameters</th><th>Particle-Particle</th><th>Particle-Wall</th></tr></thead><tbody><tr><td>Hamaker constant A</td><td>6.3 × 10<sup>-20</sup> J</td><td>0.1 × 10<sup>-20</sup> J</td></tr><tr><td>Surface electrical potential ψ</td><td>-58 mV</td><td>-52.3 mV</td></tr><tr><td>Inverse Debye length κ</td><td>1.05 × 10<sup>8</sup> m<sup>-1</sup></td><td></td></tr><tr><td>Relative permittivity ε<sub>r</sub></td><td>78.2</td><td></td></tr><tr><td>Salinity [NaCl]</td><td>1 mM</td><td></td></tr><tr><td>Temperature T</td><td>298 K</td><td></td></tr></tbody></table>

$$ \left\{ \begin{aligned} m \frac{d^2 x_1}{dt^2} \vec{x} &= \frac{AR}{12h^2} \vec{x} - \frac{R}{2} \kappa B e^{-\kappa h} \vec{x}, \\ m \frac{d^2 x_2}{dt^2} \vec{x} &= -\frac{AR}{12h^2} \vec{x} + \frac{R}{2} \kappa B e^{-\kappa h} \vec{x}, \end{aligned} \right. \qquad (22) $$
<figure><img src="image_7.png" /><figcaption>Fig. 6. Particle-particle separation distance $h(t)$ and the relative velocity of particles $h'(t)$: (a) for the attraction case with $h_0 = 2 \ \mu m$, (b) for the repulsion case with $h_0 = 0.04 \ \mu m$, and (c)-(d) for mixed attraction-repulsion with $h_0 = 2 \ \mu m$.</figcaption></figure>

where $x_1(t)$ and $x_2(t)$ denote the position of the two particles, and $B = 64\pi\epsilon_0\epsilon_r(\frac{k_B T}{Ze})^2 \tanh(\frac{Ze\psi}{4k_B T})$. Therefore, the evolution of the surface-to-surface separation distance, $h(t) = x_2(t) - x_1(t) - 2R$, reads

$$h^2 e^{kh} h'' = \frac{R}{6m} (6\kappa B h^2 - Ae^{kh}) \quad (23)$$

where $h''(t) = \frac{d^2h(t)}{dt^2}$. This ODE can be simplified for pure attraction and pure repulsion. On the one hand, for pure attraction, $B = 0$ and Eq. (23) becomes,

$$h^2 h'' = -\frac{AR}{6m}, \quad \text{for pure attraction.} \qquad (24)$$

On the other hand, for pure repulsion, $A = 0$ and we have,

$$e^{kh} h'' = \frac{\kappa B R}{m}, \quad \text{for pure repulsion.} \qquad (25)$$

These three ODEs are solved numerically using MATLAB high-order method for solving non-stiff differential equations – ode89 –, with absolute and relative tolerance of $10^{-10}$ and $10^{-8}$ respectively. The ODE results are used as reference solutions to compare with CFD-DEM predictions.

We see in Fig. 6 that the CFD-DEM including DLVO forces is in very good agreement with the reference ODE solutions. In the case of pure attraction, the surface-to-surface separation distance decreases from its initial value $h_0$ to a near-zero value, while the relative
<figure><img src="image_8.png" /><figcaption>Fig. 7. Representation of the particle and the wall distanced from h.</figcaption></figure>

velocity, h', increases exponentially until the contact. For the repulsion case, the surface-to-surface separation distance increases exponentially from its initial value and the particle relative velocity increases up to a threshold value as the repulsive force decreases with the increase of separation distance. If both repulsion and attraction are considered, the separation distance describes a periodic oscillatory motion ranging from its initial value $h_0$ down to a minimal value $h_{min} = 0.0856 \ \mu m$. The latter corresponds to the secondary minimum at which the repulsive force takes over attraction pushing back the particles to their initial positions. The three test cases confirm the robustness of the DLVO implementation in our CFD-DEM package.

## 3.2. Particle-wall DLVO interactions

In this test case, we verify the implementation of the particle-wall DLVO interactions. We consider a particle of density, $\rho = 1050 \ kg/m^3$, radius, $R = 1 \ \mu m$, and mass, $m = 4.4 \times 10^{-15} \ kg$, located at a distance, $h(t)$, from an immobile plane wall (see Fig. 7). Table 2 summarizes the particle, wall, and fluid properties. The particle is suspended in brine. The only force that applies are the attractive and repulsive DLVO interactions: there is no flow and the sedimentation and drag forces are neglected. Along the same line as the particle-particle verification cases, we consider three sets of simulations: pure attraction, pure repulsion, and mixed attraction-repulsion. The CFD-DEM predictions of the surface-to-wall separation distance $h(t)$ are compared with semi-analytical solutions.

The computational domain is a $50 \ \mu m \times 50 \ \mu m$ square with a regular grid cell size of $1 \ \mu m \times 1 \ \mu m$. We use $\Delta t = 10^{-5} \ s$ in the attraction case, $\Delta t = 10^{-8} \ s$ for the repulsion case and $\Delta t = 4 \times 10^{-5} \ s$ in the mixed attraction-repulsion case. $h_0 = 2 \ \mu m$ is used for the attraction and the mixed attraction-repulsion cases and $h_0 = 0.04 \ \mu m$ for the repulsion case.

Reference solutions for this configuration are obtained through Ordinary Differential Equations (ODE) that govern the evolution of the surface-to-wall distance, $h(t) = x_p(t) - x_w - R$, where $x_p(t)$ and $x_w$ are the particle and wall positions, respectively. The particle motion is obtained by combining Eq. (5) and Eq. (16),

$$m \frac{d^2 x_p}{dt^2} \bar{x} = - \frac{AR \left(1 + 28 \frac{h}{\lambda}\right)}{6h^2 \left(1 + 14 \frac{h}{\lambda}\right)^2} \bar{x} + B \kappa R \left(\psi_p^2 + \psi_w^2\right) \left[ \frac{2\psi_p\psi_w}{\psi_p^2 + \psi_w^2} - e^{-\kappa h} \right] \bar{x}. \quad (26)$$

where $B = 2\pi\epsilon_0\epsilon_r$. This equation is recast into an ODE solved numerically for mixed attraction-repulsion:

$$\frac{6h^2 \left(1 + h \frac{14}{\lambda}\right)^2 (e^{\kappa h} - e^{-\kappa h}) h''}{6\kappa B R h^2 \left(1 + h \frac{14}{\lambda}\right)^2 \left(2\psi_p \psi_w - \left(\psi_p^2 + \psi_w^2\right) e^{-\kappa h}\right) - A \left(1 + 28 \frac{h}{\lambda}\right) (e^{\kappa h} - e^{-\kappa h})} = \frac{R}{m}. \quad (27)$$

An ODE for pure attraction is obtained if $B = 0$,

$$6h^2 \left(1 + 14 \frac{h}{\lambda}\right)^2 h'' = - \frac{AR \left(1 + 28 \frac{h}{\lambda}\right)}{m}, \quad \text{for pure attraction.} \quad (28)$$

If $A = 0$, we obtain an ODE for pure repulsive interaction,

$$\frac{(e^{\kappa h} - e^{-\kappa h}) h''}{2\psi_p \psi_w - (\psi_p^2 + \psi_w^2) e^{-\kappa h}} = \frac{\kappa B R}{m}, \quad \text{for pure repulsion.} \quad (29)$$

Simulation results for pure attraction, pure repulsion, and mixed attraction-repulsion are shown in Fig. 8. The CFD-DEM results for the three cases are in very good agreement with the ODE reference solutions. We recover similar behaviors to the particle-particle
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Particle-wall separation distance $h(t)$ and the velocity of the particle $h'(t)$, (a) for the attraction case with $h_0 = 2~\mu\text{m}$, (b) for the repulsion case with $h_0 = 0.04~\mu\text{m}$, and (c)-(d) for mixed attraction-repulsion with $h_0 = 2~\mu\text{m}$.</figcaption></figure>

interactions including the periodic oscillatory motion already observed for mixed attractive-repulsive interactions with a different secondary minimum $h_{min} = 0.12~\mu\text{m}$. These three test cases ensure that the DLVO forces are properly calculated when a particle arrives at the vicinity of a solid wall described by a boundary condition in CFD-DEM.

### 3.3. Colloidal aggregation in a suspension

To verify the numerical implementation of the coupling between the unresolved-resolved CFD-DEM method and the DLVO + JKR theories, we simulate the aggregation kinetics of particles initially dispersed in a fluid at rest on a 2D plan. We compare the evolution of the total number of separated objects, $N(t)$ (i.e. monomers: unattached particles and clusters: attached particles), with the experimental data obtained by Earnshaw et al. [69]. Their two-dimensional colloidal systems provide an experimentally convenient approach to studying the impact of DLVO forces on aggregate dynamics.

The simulation domain consists in a two-dimensional $88.6\mu\text{m} \times 88.6\mu\text{m}$ square box discretized with a 205 × 205 regular grid. The box is closed and the lateral boundaries are impermeable walls. It is filled with a brine aqueous solution ($\text{CaCl}_2 + \text{H}_2\text{O}$) containing 0.73 M of salt concentration. Initially, $N_0 = 1000$ separated polystyrene latex particles ($\rho = 1050~\text{Kg/m}^3$, $E = 3~\text{GPa}$, $v = 0.34$) of 1 µm diameter are randomly placed in the box (see Fig. 9a). They occupy about 10% of the box surface area. The particle surface charge density is $0.4~\mu\text{C/cm}^2$ which corresponds to a surface potential of $\psi_p = 12~\text{mV}$, the particle-liquid-particle Hamaker constant
<figure><img src="image_10.png" /><figcaption>Fig. 9. (a) Image sequence of the aggregation of the colloidal layer over time obtained by simulation. (b) Comparison of the evolution of the number of separated objects (monomers and clusters) relative to the initial value in a 0.73 M CaCl₂ brine aqueous solution, obtained with simulation and experimentally (obtained with $N_0 \approx 6000$ objects). Time is non-dimensionalized by the duration of the experiment.</figcaption></figure>

is $3.4 \times 10^{-21}$ J. The particle motions are only due to long-range electrostatic double layer (repulsive) and Van der Waals (attractive) interactions. Gravity is neglected, and particle-wall interactions are not considered.

Fig. 9a shows the typical sequence of the formation of aggregates. We see that the number of objects decreases with time as the number of aggregates grows (see Fig. 9b). The growth kinetics demonstrates a ramp-up from slow aggregation early to fast aggregation until it reaches a steady state. We see that the prediction of the evolution of the number of separated objects, $N(t)$, is in good agreement with the experimental measurements. We observe, however, a slight deviation in the fast aggregation phase that can be explained by the different composition of the initial objects in the 2 cases. Indeed, in the simulation, initial objects are just composed of monomers, while in the experiment, initial objects are composed of both monomers and clusters which enhances the aggregation mechanism. This validation case gives strong confidence in our numerical model to capture the clusterization mechanisms of colloids under DLVO interactions.

## 4. Application: pore-scale modeling of pore-clogging in colloidal flows

In this section, we use the unresolved-resolved four-way coupling CFD-DEM including DLVO and JKR interactions to investigate the retention capacity and permeability reduction in a network of pores for various salinity conditions.

We use our CFD-DEM package to simulate pore-clogging by aggregation of particles in a porous medium which is a representation of a PDMS (Polydimethylsiloxane) microfluidic device [70]. The pore geometry and the meshing procedure are found in Soulaine [71]. We obtain a two-dimensional $1050 \ \mu m \times 310 \ \mu m$ pore-scale domain. The typical pore size is $d_{50} = 26 \ \mu m$, porosity is $\epsilon_i = 0.62$, and permeability is $K_i = 1.1 \times 10^{-11} \ m^2$. The domain is discretized into an unstructured mesh with 49330 cells using `snappyHexMesh`, the OpenFOAM automatic gridder. The left and right sides are inlet and outlet boundaries, respectively. The top, bottom, and grain surfaces are walls described with no-slip conditions. We consider a cloud of monodispersed polystyrene particles whose diameter is $D_p = 4 \ \mu m$. The mean mesh resolution is $D_p/\Delta x \approx 5$, which allows for an accurate calculation of the drag forces [8]. Fluid, particles, PDMS properties, and simulation parameters are listed in Table 3.

We consider three different brines with salinities [NaCl] = 1, 10, 100 mM. Both polystyrene and PDMS immersed in the fluid have a negative surface charge at $T = 293 \ K$ [72,73]. An increase in salinity increases the surface potential as well, which favors deposition [74]. In each case, the interactions between particles and the porous medium are different because the DLVO potential has different values. Fig. 10 shows that the repulsive barrier is more significant for 1 mM of salinity as the surface potential is the lowest. Overall, particle-particle and particle-wall interactions are weakly repulsive for [NaCl] = 1, 10, 100 mM with an attractive secondary minimum that will favor aggregation.

Simulations are run in two consecutive steps. First, a constant pressure difference, $\Delta p = 0.2 \ mbar$, is applied between the inlet and the outlet. The fluid flows through the porous formation until the flow field is stable. At that point, pore volume (PV) is set to $PV = 0$ (PV is a convenient measure of time in porous media defined as the ratio of the injected volume of water to the pore-space volume). Then, a concentration of particles, $C_0 = 0.3\%$, is continuously injected randomly from the inlet until the porous medium is completely
<table><caption>**Table 3**
Parameters of the CFD-DEM simulations used to investigate pore-clogging by aggregation of particles.</caption><thead><tr><th>CFD and fluid parameters</th><th></th><th>DEM and particle parameters</th><th></th></tr></thead><tbody><tr><td>CFD time-step Δt<sub>CFD</sub></td><td>10<sup>-5</sup> s</td><td>DEM time-step Δt<sub>DEM</sub></td><td>2 × 10<sup>-9</sup> s</td></tr><tr><td>Fluid density ρ<sub>f</sub></td><td>10<sup>3</sup> kg/m<sup>3</sup></td><td>Particle density ρ<sub>p</sub></td><td>1050 kg/m<sup>3</sup></td></tr><tr><td>Fluid viscosity μ<sub>f</sub></td><td>10<sup>-3</sup> Pa.s</td><td>Poisson ratio v<sub>p</sub> / v<sub>w</sub></td><td>0.34 / 0.5</td></tr><tr><td>Temperature T</td><td>293 K</td><td>Young modulus E<sub>p</sub> / E<sub>w</sub></td><td>3 GPa / 2 M Pa</td></tr><tr><td>Salinity [NaCl]</td><td>1 ; 10 ; 100 mM</td><td>Hamaker constant A<sub>pLp</sub> / A<sub>pLw</sub></td><td>3.4 × 10<sup>-21</sup> J / 10<sup>-21</sup> J</td></tr><tr><td>Ionic valence Z</td><td>1</td><td>Surface potential ψ<sub>p</sub></td><td>-58 ; -26 ; -20 mV</td></tr><tr><td>Model depth b</td><td>20 μm</td><td>Surface potential ψ<sub>w</sub></td><td>-52.3 ; -45 ; -26 mV</td></tr><tr><td></td><td></td><td>Surface energy W<sub>pLp</sub> / W<sub>pLw</sub></td><td>41 / 20 mJ/m<sup>2</sup></td></tr></tbody></table>

<figure><img src="image_11.png" /><figcaption>**Fig. 10.** DLVO potential (normalized by $k_B T$) for (a) particle-particle, and (b) particle-wall interactions for various NaCl salt concentrations.</figcaption></figure>

clogged (all entry pores of the porous medium are blocked). Because the simulations are 2D depth-integrated, a Hele-Shaw correction term, $12\mu_f c v^f / b^2$, is added in Eq. (3) to account for the hydrodynamic effects in the thickness, $b$, of the microfluidic device [75].

Snapshots of the particle migration and retention are shown in Fig. 11. For the three simulations, we observe that some particles deposit rapidly at the entrance of the porous medium. It means that long-range interactions are dominant over the hydrodynamic forces. Layers of deposited particles grow until they clog pores. Thus, The clogs change the system porosity and reroute the local flow lines, affecting the system permeability. When all entrance pores of the porous system are clogged, the mean flow velocity drops to near zero value and no more particles percolate through the domain.

Large-scale properties such as permeability and porosity are strongly impacted by pore-clogging. The impact of salinity on the evolution of porosity damage severity and the permeability damage severity is shown in Fig. 12a. The first is defined as $\epsilon_{rem} = \langle\epsilon\rangle / \epsilon_i$ where $\langle\epsilon\rangle$ is the domain porosity at time $t$. The second is calculated using $K_r = K_f / K_i$, where $K_f$ is the system permeability at time $t$ obtained using Darcy's law, $K_f = \frac{\langle v_x^f \rangle \mu_f L}{\Delta\rho}$, where $\langle v_x^f \rangle$ is the domain-averaged fluid velocity, and $L$ is the length of the domain. Both $K_r$ and $\epsilon_{rem}$ decrease at the beginning of the injection when $PV \in [0, 0.5]$. The same evolution is observed for the three salinities. This reduction corresponds to the flow resistance related to the suspension itself. For $PV > 0.5$, pores at the medium's entrance are clogged, resulting in a sudden drop in permeability damage. We observe fluctuations in the permeability damage evolution. The remobilization of clusters of deposited particles explains them. These remobilizations happen mostly for 1 mM and 10 mM salt concentrations for which the particle-wall potential barrier is more important. For identical flow conditions, the permeability impairment by aggregation of particles increases with salinity, as observed experimentally [19,14]. It can be seen that the $K_r$ plateau values – that indicate the clogging of the porous medium – are reached earlier for higher salinity.

Pore-scale simulations with CFD-DEM directly map the velocity and particle concentration. This constitutes a major advantage compared with other approaches classically used to investigate flow and transport in porous media. Better than a global insight into permeability impairment with average properties (e.g. $K_r$ and $\epsilon_{rem}$), the Probability Density Function (PDF) gives information about the redistribution of local flow rates during pore-clogging. From Fig. 12b, we see that for the various salinities used, the particulate flow with particle deposition inside the porous medium causes an increase in the densities of low fluid velocities and a reduction of the highest value of the fluid velocities compared to the state with no particles. If Fig. 12a characterizes the mobility reduction due to pore-clogging, Fig. 12b characterizes the flow velocities redistribution inside the porous medium. The results presented in Fig. 12b
<figure><img src="image_12.png" /><figcaption>Fig. 11. Snapshots of the particle migration and deposition for various salinity conditions: a) 1 mM, b) 10 mM, c) 100 mM NaCl. The background colors correspond to the fluid velocity field magnitude.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 12. (a) Evolution of permeability damage $K_r$ and porosity damage $\epsilon_{rem}$ over the injected PV for various salt concentrations. (b) Probability density function (PDF) of the longitudinal flow velocity $v_x^f$ with no particles, and with particles for various salt concentrations at PV = 2.2.</figcaption></figure>

support the conclusion that the porous medium permeability damage increases with increasing salinity (for the same concentration of injected particles).

This illustration highlights the ability of the CFD-DEM model for colloidal transport to simulate pore-clogging due to the aggregation of charged particles under various flow and chemistry conditions.
## 5. Conclusions

Understanding migration, deposition, and remobilization mechanisms of fine particles in porous media is key to deciphering pore-clogging formation and the subsequent permeability impairment. We implemented long-range interactions including London-Van der Waals attraction and electrical double layer repulsion forces in CFD-DEM to simulate particle clogging by aggregation at the pore-scale. Our generic model can simulate a wide range of particle size distribution, flow rates, particle concentrations, and pH/salinity conditions. High-resolution CFD-DEM simulations directly map the flow lines and clogs within the porous structure. They are, therefore, a real asset to probe local permeability reduction due to the injection and retention of charged particles.

We have implemented DLVO forces combined with the JKR theory in our unresolved-resolved four-way coupling CFD-DEM to model colloidal flows at the pore-scale in the approximation of constant surface potential. We also consider that an adhesive torque that resists the hydrodynamic torque applies to particle dynamics. The numerical implementation has been verified using benchmark cases solved with ODE solvers. These test cases demonstrate that the model accurately describes the particle-particle and particle-wall attraction and repulsion.

Combined with the unresolved-resolved feature described in Part 1 which overcomes the constraint of the particle size relative to the cell size, and maintains a decent computational cost, our CFD-DEM approach for colloidal transport in porous media can capture the three main clogging mechanisms: sieving (clogging by size exclusion), bridging (hydrodynamic formation of an arch of particles), and pore-clogging by aggregation (formation of multiple layers of deposited particles). The simulations presented in this paper illustrate the ability of CFD-DEM combined with DLVO forces to investigate the complex interplay between particle migration, deposition, and remobilization in porous media. Furthermore, our model offers new possibilities to explore deposition regimes around cylindrical collectors through hydrodynamic and electrochemical interaction effects and therefore, gives better insights into colloidal retention in porous media.

## Nomenclature

### List of superscripts and subscripts

adh Adhesive

c Contact

f Fluid

g Gravitational

hyd Hydrodynamic

l Liquid

nc Non-contact

p Particle

w Wall

### List of symbols

$\delta_n$ Normal overlap distance

$\omega$ Angular velocity vector

$\Delta p$ Pressure difference

$\Delta t$ Time-step

$\Delta x$ Cell's characteristic size

$\epsilon$ Porosity

$\epsilon_0 \epsilon_r$ Dielectric constant of the solvent

$\gamma$ Interfacial energy

$\kappa$ Inverse of the Debye screening length

$\lambda$ Characteristic wavelength

$\langle \mathbf{v}_x^f \rangle$ the Domain-averaged velocity

$\mu$ Dynamic viscosity

$\nu$ Poisson ratio

$\bar{\mathbf{v}}$ Averaged velocity vector

$\psi$ Surface potential

$\rho$ Density

$\sigma_{LJ}$ Lennard-Jones atomic collision diameter

$\sigma_{st}$ Stern layer thickness

$\mathbf{F}$ Force vector

$\mathbf{M}$ Torque vector

$\mathbf{n}$ Unit normal vector

$\mathbf{r}_i$ Position vector

$\mathbf{t}$ Tangential unit vector

$\mathbf{U}$ Translational velocity vector

$\mathbf{v}$ Velocity vector

$\mathbf{v}_t$ Particle tangential slip velocity vector

$\zeta$ Adhesive overlapping distance

$\mathbf{A}$ Hamaker constant

$\mathbf{a}$ Contact radius

$\mathbf{b}$ Geometry thickness

$\mathbf{C}$ Particle concentration

$\mathbf{D}$ Diameter

$\mathbf{d}$ Particle centers separation distance

$\mathbf{E}$ Young modulus

$\mathbf{e}$ Charge of the electron

$\mathbf{h}$ Smallest surface separation distance

$\mathbf{I}$ Moment of inertia

$\mathbf{I}_s$ Ionic strength

$\mathbf{K}$ Cell-permeability

$\mathbf{k}$ Equivalent elastic constant

$\mathbf{k}_B$ Boltzmann's constant

$\mathbf{L}$ Geometry length

$\mathbf{l}_x$ Lever arm

$\mathbf{m}$ Mass

$\mathbf{N}(t)$ Current number of objects (single particles and clusters)

$\mathbf{N}_0$ Initial number of objects

$\mathbf{N}_A$ Avogadro's number

$\mathbf{p}$ Fluid pressure

$\mathbf{R}$ Particle radius

$\mathbf{r}$ Equivalent radius

$\mathbf{T}$ Temperature

$\mathbf{t}$ Time

$\mathbf{V}$ Potential energy

$\mathbf{W}$ Adhesion energy

$\mathbf{Z}$ Valence of the solvent
CRediT authorship contribution statement

**Laurez Maya Fogouang:** Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Software, Validation, Visualization, Writing – original draft, Writing – review & editing. **Laurent André:** Conceptualization, Funding acquisition, Investigation, Methodology, Project administration, Supervision, Writing – review & editing. **Philippe Leroy:** Conceptualization, Funding acquisition, Investigation, Methodology, Project administration, Supervision, Writing – review & editing. **Cyprien Soulaine:** Conceptualization, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Software, Supervision, Validation, Writing – original draft, Writing – review & editing.

Declaration of competing interest

The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: Cyprien Soulaine reports financial support was provided by European Research Council. Laurez Maya Fogouang reports financial support was provided by Bureau for Geological and Mining Research. If there are other authors, they declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability

Data will be made available on request.

Acknowledgements

This research is funded by the European Union (ERC, COCONUT, 101043288). LM received a PhD fellowship from BRGM and Région Centre-Val de Loire. The authors benefitted from the Centre de Calcul Scientifique en région Centre-Val de Loire's (CaSciMoDOT) high-performance computational resources. The authors would like to thank Laurent Oxarango, Olivier Liot, and Noura Eddaoui for their interesting discussions and suggestions.

References

[1] S. Vigneswaran, R.B. Suazo, A detailed investigation of physical and biological clogging during artificial recharge, Water Air Soil Pollut. 35 (1987) 119-140, https://doi.org/10.1007/bf00183848.

[2] S. Rinck-Pfeiffer, Interrelationships between biological, chemical, and physical processes as an analog to clogging in aquifer storage and recovery (ASR) wells, Water Res. 34 (2000) 2110-2118, https://doi.org/10.1016/s0043-1354(99)00356-5.

[3] G. Kampel, G.H. Goldszein, J.C. Santamarina, Plugging of porous media and filters: maximum clogged porosity, Appl. Phys. Lett. 92 (2008), https://doi.org/10.1063/1.2883947.

[4] W. Song, X. Liu, T. Zheng, J. Yang, A review of recharge and clogging in sandstone aquifer, Geothermics 87 (2020) 101857, https://doi.org/10.1016/j.geothermics.2020.101857.

[5] R. Ines, K. Eric, F. Marc, P. Nicolas, L. Béatrice, Formation damage induced by clay colloids deposition in Triassic clastic geothermal fields: insights from an experimental approach, Appl. Clay Sci. 234 (2023) 106868, https://doi.org/10.1016/j.clay.2023.106868.

[6] L.M. McDowell-Boyer, J.R. Hunt, N. Sitar, Particle transport through porous media, Water Resour. Res. 22 (1986) 1901-1921, https://doi.org/10.1029/wr022i013p01901.

[7] E. Dressaire, A. Sauret, Clogging of microfluidic systems, Soft Matter 13 (2017) 37-48, https://doi.org/10.1039/c6sm01879c.

[8] L.F. Maya, L. André, C. Soulaine, Particulate transport in porous media at pore-scale. Part 1: unresolved-resolved four-way coupling CFD-DEM, 2024.

[9] G. Hua, W. Zhu, L. Zhao, J. Huang, Clogging pattern in vertical-flow constructed wetlands: insight from a laboratory study, J. Hazard. Mater. 180 (2010) 668-674, https://doi.org/10.1016/j.jhazmat.2010.04.088.

[10] B. Dersoir, M.R. de Saint Vincent, M. Abkarian, H. Tabuteau, Clogging of a single pore by colloidal particles, Microfluid. Nanofluid. 19 (2015) 953-961, https://doi.org/10.1007/s10404-015-1624-y.

[11] N. Delouche, A.B. Schofield, H. Tabuteau, Dynamics of progressive pore clogging by colloidal aggregates, Soft Matter 16 (2020) 9899-9907, https://doi.org/10.1039/d0sm01403f.

[12] Y. Bigna, M.S. Oyeneyin, J.M. Peden, Investigation of pore-blocking mechanism in gravel packs in the management and control of fines migration, https://doi.org/10.2118/27342-ms, 1994.

[13] M.B. Oyeneyin, J.M. Peden, A. Hosseini, G. Ren, Factors to consider in the effective management and control of fines migration in high permeability sands, https://doi.org/10.2118/30112-ms, 1995.

[14] G.C. Agbangla, É. Climent, P. Bacchin, Experimental investigation of pore clogging by microparticles: evidence for a critical flux density of particle yielding arches and deposits, Sep. Purif. Technol. 101 (2012) 42-48, https://doi.org/10.1016/j.seppur.2012.09.011.

[15] N.-D. Ahfir, A. Hammadi, A. Alem, H. Wang, G.L. Bras, T. Ouahbi, Porous media grain size distribution and hydrodynamic forces effects on transport and deposition of suspended particles, J. Environ. Sci. 53 (2017) 161-172, https://doi.org/10.1016/j.jes.2016.01.032.

[16] J. Fetzer, M. Holzner, M. Plötze, G. Furrer, Clogging of an Alpine streambed by silt-sized particles – insights from laboratory and field experiments, Water Res. 126 (2017) 60-69, https://doi.org/10.1016/j.watres.2017.09.015.

[17] J. Jung, H. Kang, S.C. Cao, R.I. Al-Raoush, K. Alshibli, J.Y. Lee, Effects of fine-grained particles' migration and clogging in porous media on gas production from hydrate-bearing sediments, Geofluids 2019 (2019) 1-11, https://doi.org/10.1155/2019/5061216.

[18] N. Delouche, J. van Doorn, T. Kodger, A. Schofield, J. Sprakel, H. Tabuteau, The contribution of colloidal aggregates to the clogging dynamics at the pore scale, J. Membr. Sci. 635 (2021) 119509, https://doi.org/10.1016/j.memsci.2021.119509.

[19] J.N. Ryan, M. Elimelech, Colloid mobilization and transport in groundwater, Colloids Surf. A, Physicochem. Eng. Asp. 107 (1996) 1-56, https://doi.org/10.1016/0927-7757(95)03384-x.

[20] N. Sefrioui, A. Ahmadi, A. Omari, H. Bertin, Numerical simulation of retention and release of colloids in porous media at the pore scale, Colloids Surf. A, Physicochem. Eng. Asp. 427 (2013) 33-40, https://doi.org/10.1016/j.colsurfa.2013.03.005.
[21] M. Samari-Kermani, S. Jafari, M. Rahnama, A. Raooft, Ionic strength and zeta potential effects on colloid transport and retention processes, Colloid Interface Sci. Commun. 42 (2021) 100389, https://doi.org/10.1016/j.colcom.2021.100389.

[22] M. Ramezanpour, M. Siavashi, A.Q. Raeini, M.J. Blunt, Pore-scale simulation of nanoparticle transport and deposition in a microchannel using a Lagrangian approach, J. Mol. Liq. 355 (2022) 118948, https://doi.org/10.1016/j.molliq.2022.118948.

[23] H. Li, S. Wang, X. Chen, L. Xie, B. Shao, Y. Ma, Cfd-dem simulation of aggregation and growth behaviors of fluid-flow-driven migrating particle in porous media, Geoenergy Sci. Eng. 231 (2023) 212343, https://doi.org/10.1016/j.geoen.2023.212343.

[24] P.-Y. Lagrée, L. Staron, S. Popinet, The granular column collapse as a continuum: validity of a two-dimensional Navier-Stokes model with a $\mu(i)$-rheology, J. Fluid Mech. 686 (2011) 378-408, https://doi.org/10.1017/jfm.2011.335.

[25] S. Cloete, S.T. Johansen, S. Amini, Grid independence behavior of fluidized bed reactor simulations using the two-fluid model: effect of particle size, Powder Technol. 269 (2015) 153-165, https://doi.org/10.1016/j.powtec.2014.08.055.

[26] Y. Tsuji, T. Kawaguchi, T. Tanaka, Discrete particle simulation of two-dimensional fluidized bed, Powder Technol. 77 (1993) 79-87, https://doi.org/10.1016/0032-5910(93)85010-7.

[27] C. O'Sullivan, Particulate Discrete Element Modelling: A Geomechanics Perspective, Applied Geotechnics, vol. 4, Spon Press, 2011.

[28] A.J.C. Ladd, R. Verberg, Lattice-Boltzmann simulations of particle-fluid suspensions, J. Stat. Phys. 104 (2001) 1191-1251, https://doi.org/10.1023/a:1010414013942.

[29] A. Parvan, S. Jafari, M. Rahnama, S.N. apourvari, A. Raooft, Insight into particle retention and clogging in porous media; a pore scale study using lattice Boltzmann method, Adv. Water Resour. 138 (2021) 103530, https://doi.org/10.1016/j.adwatres.2020.103530.

[30] G.C. Agbangla, É. Climent, P. Bacchin, Numerical investigation of channel blockage by flowing microparticles, Comput. Fluids 94 (2014) 69-83, https://doi.org/10.1016/j.compfluid.2014.01.018.

[31] J. Su, G. Chai, L. Wang, W. Cao, Z. Gu, C. Chen, X.Y. Xu, Pore-scale direct numerical simulation of particle transport in porous media, Chem. Eng. Sci. 199 (2019) 613-627, https://doi.org/10.1016/j.ces.2019.01.033.

[32] H. Sun, S. Xu, X. Pan, L. Shi, X. Geng, Y. Cai, Investigating the jamming of particles in a three-dimensional fluid-driven flow via coupled cfd-dem simulations, Int. J. Multiph. Flow 114 (2019) 140-153, https://doi.org/10.1016/j.ijmultiphaseflow.2019.01.017.

[33] A. Elrahmani, R.I. Al-Raoush, H. Abugazia, T. Seers, Pore-scale simulation of fine particles migration in porous media using coupled cfd-dem, Powder Technol. 398 (2022) 117130, https://doi.org/10.1016/j.powtec.2022.117130.

[34] S. Pirker, D. Kahrimanovic, C. Goniva, Improving the applicability of discrete phase simulations by smoothening their exchange fields, Appl. Math. Model. 35 (2011) 2479-2488, https://doi.org/10.1016/j.apm.2010.11.066.

[35] Z. Peng, E. Doroodchi, C. Luo, B. Moghtaderi, Influence of void fraction calculation on fidelity of CFD-DEM simulation of gas-solid bubbling fluidized beds, AIChE J. 60 (2014) 2000-2018, https://doi.org/10.1002/aic.14421.

[36] W. Poon, M. Haw, Mesoscopic structure formation in colloidal aggregation and gelation, Adv. Colloid Interface Sci. 73 (1997) 71-126, https://doi.org/10.1016/s0001-8686(97)90003-8.

[37] Y. Liang, N. Hilal, P. Langston, V. Starov, Interaction forces between colloidal particles in liquid: theory and experiment, Adv. Colloid Interface Sci. 134-135 (2007) 151-166, https://doi.org/10.1016/j.cis.2007.04.003.

[38] B. Derjaguin, L. Landau, Theory of the stability of strongly charged lyophobic sols and of the adhesion of strongly charged particles in solutions of electrolytes, Acta Physicochim. URSS 14 (1941) 633-662.

[39] E.J.W. Verwey, J.T.G. Overbeek, Theory of the Stability of Lyophobic Colloids; the Interaction of Sol Particles Having an Electric Double Layer, Elsevier Pub. Co., New York, 1948.

[40] I. Sameut Bouhaik, P. Leroy, P. Ollivier, M. Azaroual, L. Mercury, Influence of surface conductivity on the apparent zeta potential of tio2 nanoparticles: application to the modeling of their aggregation kinetics, J. Colloid Interface Sci. 406 (2013) 75-85, https://doi.org/10.1016/j.cjis.2013.05.034.

[41] P. Leroy, A. Maineult, S. Li, J. Vinogradov, The zeta potential of quartz. Surface complexation modelling to elucidate high salinity measurements, Colloids Surf. A, Physicochem. Eng. Asp. 650 (2022) 129507, https://doi.org/10.1016/j.colsurfa.2022.129507.

[42] G.C. Agbangla, P. Bacchin, E. Climent, Collective dynamics of flowing colloids during pore clogging, Soft Matter 10 (2014) 6303-6315, https://doi.org/10.1039/c4sm00869c.

[43] H. Abbasford, G. Evans, R. Moreno-Atanasio, Effect of van der Waals force cut-off distance on adhesive collision parameters in DEM simulation, Powder Technol. 299 (2016) 9-18, https://doi.org/10.1016/j.powtec.2016.05.020.

[44] M. Mihajlovic, I. Roghair, M.V.S. Annaland, On the numerical implementation of the van der Waals force in soft-sphere discrete element models for gas-solid fluidization, Chem. Eng. Sci. 226 (2020) 115794, https://doi.org/10.1016/j.ces.2020.115794.

[45] E. Ruckenstein, D.C. Prieve, Adsorption and desorption of particles and their chromatographic separation, AIChE J. 22 (1976) 276-283, https://doi.org/10.1002/aic.690220209.

[46] G.A. Schumacher, T.G.M. van de Ven, Brownian motion of charged colloidal particles surrounded by electric double layers, Faraday Discuss. Chem. Soc. 83 (1987) 75, https://doi.org/10.1039/dc9878300075.

[47] R. Muneer, M.R. Hashmet, P. Pourafshary, Fine migration control in sandstones: surface force analysis and application of dlv0 theory, ACS Omega 5 (2020) 31624-31639, https://doi.org/10.1021/acsomega.0c03943.

[48] S. Liu, I. Shikhov, Y. Cui, C.H. Arns, Pore-scale morphology effects on colloid deposition by trajectory tracking simulations, Geoenergy Sci. Eng. 227 (2023) 211772, https://doi.org/10.1016/j.geoen.2023.211772.

[49] T. Mahmood, A. Amirtharajah, T.W. Sturm, K.E. Dennett, A micromechanics approach for attachment and detachment of asymmetric colloidal particles, Colloids Surf. A, Physicochem. Eng. Asp. 177 (2001) 99-110, https://apisemanticscholar.org/CorpusID:95152653.

[50] K. Johnson, K. Kendall, A. Roberts, Surface energy and the contact of elastic solids, Proc. R. Soc. Lond. Ser. A, Math. Phys. Sci. 324 (1971) 301-313, https://doi.org/10.1098/rspa.1971.0141.

[51] M. Trofa, G. D'Avino, L. Sicignano, G. Tomaiuolo, F. Greco, P.L. Maffettone, S. Guido, CFD-DEM simulations of particulate fouling in microchannels, Chem. Eng. J. 358 (2019) 91-100, https://doi.org/10.1016/j.cej.2018.09.207.

[52] Y. Zhou, L. Chen, Y. Gong, S. Wang, Pore-scale simulations of particles migration and deposition in porous media using LBM-DEM coupling method, Processes 9 (2021) 465, https://doi.org/10.3390/pr9030465.

[53] K. Cheng, J. Zhu, F. Qian, B. Cao, J. Lu, Y. Han, Cfd-dem simulation of particle deposition characteristics of pleated air filter media based on porous media model, Particuology 72 (2023) 37-48, https://doi.org/10.1016/j.partic.2022.02.003.

[54] Z. Wang, Y. Teng, M. Liu, A semi-resolved CFD-DEM approach for particulate flows with kernel based approximation and Hilbert curve based searching strategy, J. Comput. Phys. 384 (2019) 151-169, https://doi.org/10.1016/j.jcp.2019.01.017.

[55] Z.Y. Zhou, S.B. Kuang, K.W. Chu, A.B. Yu, Discrete particle simulation of particle-fluid flow: model formulations and their applicability, J. Fluid Mech. 661 (2010) 482-510, https://doi.org/10.1017/s002211201000306x.

[56] S. Whitaker, The Forchheimer equation: a theoretical development, Transp. Porous Media 25 (1996) 27-61, https://doi.org/10.1007/bf00141261.

[57] C.V. Oss, R. Good, M. Chaudhury, The role of van der Waals forces and hydrogen bonds in "hydrophobic interactions" between biopolymers and low energy surfaces, J. Colloid Interface Sci. 111 (1986) 378-390, https://doi.org/10.1016/0021-9797(86)90041-x.

[58] P. Leroy, A. Maineult, Exploring the electrical potential inside cylinders beyond the Debye-Hückel approximation: a computer code to solve the Poisson-Boltzmann equation for multivalent electrolytes, Geophys. J. Int. (2018), https://doi.org/10.1093/gji/ggy124.
[59] P. Leroy, N. Devau, A. Revil, M. Bizi, Influence of surface conductivity on the apparent zeta potential of amorphous silica nanoparticles, J. Colloid Interface Sci. 410 (2013) 81–93, https://doi.org/10.1016/j.jcis.2013.08.012.

[60] H. Hamaker, The London—van der Waals attraction between spherical particles, Physica 4 (1937) 1058–1072, https://doi.org/10.1016/s0031-8914(37)80203-7.

[61] F. London, The general theory of molecular forces, Trans. Faraday Soc. 33 (1937) 8b, https://doi.org/10.1039/tf937330008b.

[62] J.H. Schenkel, J.A. Kitchener, A test of the Derjaguin-Verwey-Overbeek theory with a colloidal suspension, Trans. Faraday Soc. 56 (1960) 161–173, https://doi.org/10.1039/TF9605600161.

[63] R. Hogg, T.W. Healy, D.W. Fuerstenau, Mutual coagulation of colloidal dispersions, Trans. Faraday Soc. 62 (1966) 1638, https://doi.org/10.1039/tf9666201638.

[64] J. Gregory, Approximate expressions for retarded Van der Waals interaction, 1981.

[65] H. Hertz, Über die beruhrung fester elastischer korper (on the contact of elastic solids), J. Reine Angew. Math. (1882).

[66] C.-W. Hong, From long-range interaction to solid-body contact between colloidal surfaces during forming, J. Eur. Ceram. Soc. (1998).

[67] S. Torkzaban, S.A. Bradford, S.L. Walker, Resolving the coupled effects of hydrodynamics and DLVO forces on colloid attachment in porous media, Langmuir 23 (2007) 9652–9660, https://doi.org/10.1021/la700995e.

[68] Y.-S. Yu, M.-C. Wang, X. Huang, Evaporative deposition of polystyrene microparticles on PDMS surface, Sci. Rep. 7 (2017), https://doi.org/10.1038/s41598-017-14593-5.

[69] J.C. Earnshaw, M.B.J. Harrison, D.J. Robinson, Local order in two-dimensional colloidal aggregation, Phys. Rev. E 53 (1996) 6155–6163, https://doi.org/10.1103/physreve.53.6155.

[70] S. Roman, C. Soulaine, M.A. AlSaud, A. Kovscek, H. Tchelepi, Particle velocimetry analysis of immiscible two-phase flow in micromodels, Adv. Water Resour. 95 (2016) 199–211, https://doi.org/10.1016/j.advwatres.2015.08.015.

[71] C. Soulaine, Introduction to open-source computational fluid dynamics using openfoam® technology - part I, https://doi.org/10.5281/ZENODO.10647983, 2024.

[72] B.J. Kirby, E.F. Hasselbrink, Zeta potential of microfluidic substrates: 2. Data for polymers, Electrophoresis 25 (2004) 203–213, https://doi.org/10.1002/elps.200305755.

[73] Z. Liu, Y. Song, D. Li, Detecting zeta potential of polydimethylsiloxane (pdms) in electrolyte solutions with atomic force microscope, J. Colloid Interface Sci. 578 (2020) 116–123, https://doi.org/10.1016/j.jcis.2020.05.061.

[74] M.J. García-Salinas, M.S. Romero-Cano, F.J. Nieves, Zeta potential study of a polystyrene latex with variable surface charge: influence on the electroviscous coefficient, 2000, pp. 112–116, https://doi.org/10.1007/3-540-46545-6_23.

[75] C. Soulaine, J. Maes, S. Roman, Computational microfluidics for geosciences, Front. Water 3 (2021), https://doi.org/10.3389/frwa.2021.643714.