The impact of diffusiophoresis on hydrodynamic dispersion and filtration in porous media

Abstract

It is known that the dispersion of colloidal particles in porous media is determined by medium structure, pore-scale flow variability and diffusion. However, much less is known about how diffusiophoresis, that is, the motion of colloidal particles along salt gradients, impacts large-scale particle dispersion in porous media. To shed light on this question, we perform detailed pore-scale simulations of fluid flow, solute transport and diffusiophoretic particle transport in a two-dimensional hyper-uniform porous medium. Particles and solute are initially uniformly distributed throughout the medium. The medium is flushed at constant flow rate, and particle breakthrough curves are recorded at the outlet to assess the macroscopic effects of diffusiophoresis. Particle breakthrough curves show non-Fickian behaviour manifested by strong tailing that is controlled by the diffusiophoretic mobility. Although diffusiophoresis is a short-time, microscopic phenomenon owing to the fast attenuation of salt gradients, it governs macroscopic colloid dispersion through the partitioning of particles into transmitting and dead-end pores. We quantify these behaviours by an upscaled analytical model that describes both the retention and release of colloids in dead-end pores and the observed long-time tailings. Our results suggest that diffusiophoresis is an efficient tool to control particle dispersion and filtration through porous media.

1. Introduction

Natural and engineered porous media are typically characterised by complex porous structures with stagnant flow regions that are not accessible by advection (Bear 1988). Understanding how nano-sized and micro-sized colloidal particles, for example fine powders, contaminants or biological materials, are transported along with nutrients or dissolved solutes through porous media is not only of fundamental interest but also important for technological applications such as groundwater remediation (Kahler & Kabala 2019), the design of water filtration systems (Shin et al. 2017; Miele, de Anna & Dentz 2019) and microfluidics for biomedical applications (Rasmussen, Pedersen & Marie 2020).

Diffusiophoresis (DP), that is, the motion of microscopic colloidal particles driven by local gradients of salt concentration (Derjaguin et al. 1947), has been demonstrated both theoretically (Prieve et al. 1984; Anderson 1989; Velegol et al. 2016; Ault, Shin & Stone 2018) and experimentally (Abécassis et al. 2008; Palacci et al. 2010; Kar et al. 2015; Shin et al. 2016; Battat et al. 2019; Singh et al. 2020) as a particle manipulation tool in simple microfluidic set-ups. The physical mechanisms that drive this physicochemical phenomenon can be broken down to two components: chemiphoresis that occurs due the osmotic pressure gradient along the surface of a charged particle (at the scale of the particle) and electrophoresis arising due to the difference in the diffusivities between the cation and the anion in the electrolyte solution (Prieve et al. 1984; Anderson 1989; Velegol et al. 2016). Note that DP for uncharged colloidal particles occurs via chemiphoresis alone (Derjaguin et al. 1947). Diffusiophoresis has been used to induce particle focusing (Shi et al. 2016), separation (Shin et al. 2017; Shin, Warren & Stone 2018; Rasmussen et al. 2020; Jotkar & Cueto-Felgueroso 2021), banding (Staffeld & Quinn 1989) and trapping (Singh et al. 2020).

Recent works have studied the impact of DP on colloid transport in microfluidic channels connected to cavities filled with hydrogel (Doan et al. 2021; Sambamoorthy & Chu 2023) or biofilm (Somasundar et al. 2023), or in a channel connected to a nano-porous medium (Lee et al. 2020). In this paper, we focus on more complex geometries in order to assess the impact of DP on hydrodynamic dispersion and particle filtration in disordered porous media. We consider an intricate porous structure characterised by dead-end pores (DEPs) that are connected via a network of percolating channels or transmitting pores (TPs). Such DEP–TP structures lead to anomalous transport of passive tracers (Bordoloi et al. 2022). To date, a fundamental understanding of the impact of microscopic interactions of DP and the complexity and disorder of the porous medium on particle transport is missing. We address this question using detailed numerical pore-scale simulations and analytical modelling to elucidate how DP couples with the medium structure to alter macroscopic particle transport. We provide a direct link between the colloidal diffusiophoretic mobility and salt Péclet number to predict, without any fitting parameter, the macroscopic colloidal transport at all times.

2. Numerical simulations

We study a fluid-saturated porous system where a particle suspension gets displaced by a continuously injected salt solution. This particular system mimics a cleanup scenario of an initially contaminated geological or biological porous medium. Nevertheless, the observed transport mechanisms, and the results on particle dispersion, can be transferred to other scenarios. The velocity experienced by each transported particle results from advection in the flow field $\boldsymbol u$, which is controlled by the medium structure and imposed flow rate (Dentz et al. 2011; de Anna et al. 2017; Dentz, Icardi & Hidalgo 2018), and from DP.

2.1. Governing equations

In the thin Debye layer limit, for dilute solutions with valence symmetric mono-valent solutes, that is, for $Z:Z$ electrolytes with $Z = 1$, like LiCl and NaCl, the diffusiophoretic drift is proportional to the gradient of the logarithm of the solute concentration s

(2.1)\begin{equation} \boldsymbol{u}_{dp} = \varGamma_p \boldsymbol{\nabla} \ln s \end{equation}

and the diffusiophoretic mobility $\varGamma _p$ is approximately a constant. For completeness, we mention that in the finite Debye layer limit $\varGamma _p$ is proportional to the logarithm of the solute concentration (Kirby & Hasselbrink 2004). Here, we consider dilute solutions in the thin Debye layer limit, for which $\varGamma _p$ is given by (Prieve et al. 1984; Anderson 1989)

(2.2)\begin{equation} \varGamma_p = \frac{\varepsilon (k_B T)^2}{ \mu (Z e)^2 } \left\{ \beta \frac{Z e \zeta}{k_B T} + 4 \ln \left[ \cosh \left(\frac{Z e \zeta}{4k_B T}\right) \right]\right\}, \end{equation}

where $\varepsilon$ is the dielectric permittivity of the medium, $\mu$ is the dynamic viscosity of the solution, $k_B$ is the Boltzmann constant, $T$ is the absolute temperature, $Z$ is the valence of the constituent ions of the solute, $\zeta$ is the zeta potential of the particle, $e$ is the proton charge and $\beta =(D_+ - D_-)/(D_+ + D_-)$ measures the difference in diffusivity $D_+$ of the cation and $D_-$ of the anion. The prefactor in (2.2) at room temperature is equal to $\varepsilon (k_B T)^2/\mu (Z e)^2 = 470$ $\mathrm {\mu }$m$^2$ s$^{-1}$. Expression (2.2) is valid for valence symmetric $Z:Z$ electrolytes. Note that $\varGamma _p$ depends implicitly on the particle size because $\zeta$ is affected by the absolute surface charge and size of the particle. The logarithmic dependence of the diffusiophoretic drift $\boldsymbol {u}_{dp}$ on the solute concentration gradients $\boldsymbol {\nabla } s$ allows for rapid and efficient particle motion, even in low concentration areas. Positive values of $\varGamma _p$ refer to upward particle migration along the direction of the concentration gradient, whereas negative values indicate downward migration along the concentration gradient.

In the low Reynolds number limit ($Re\ll 1$), the fluid–solute–particle dynamics is governed by the Stokes equation and the equation of continuity for fluid flow assuming an incompressible fluid

(2.3a)\begin{gather} - \boldsymbol{\nabla} p + \mu{\nabla}^2 \boldsymbol{u} = 0, \end{gather}

(2.3b)\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} =0, \end{gather}

where $\mu$ is the dynamic viscosity of the fluid. The advection-diffusion equations for solute and particle transport are (Ault et al. 2018)

(2.3c)\begin{gather} \frac{\partial s}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} s =D_s{\nabla}^2 s , \end{gather}

(2.3d)\begin{gather}\frac{\partial c}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} (\boldsymbol{u} +\boldsymbol{u}_{dp})c = D_p{\nabla}^2 c , \end{gather}

where $\boldsymbol {u}$ is the two-dimensional velocity field, $p$ is the pressure, $s$ is the solute concentration and $c$ is the particle concentration, $D_s$ and $D_p$ are the diffusion coefficients of solute and particles. Typically, solutes diffuse much faster than particles such that $D_s \gg D_p$. The volumetric flux through the porous medium is denoted by $U$ such that the average pore velocity is $\bar u = U/\phi$. The advection time over the mean pore length $\lambda$ is defined by $\tau _v = \lambda / \bar u$. We define the dimensionless position vector, pressure, time, velocity, volumetric flux and diffusiophoretic mobility as

(2.4a–f)\begin{equation} \boldsymbol x = \boldsymbol x' \lambda,\quad p = p' \frac{\mu \bar u}{\lambda} ,\quad t = t' \tau_v, \quad \boldsymbol u = \boldsymbol u' \bar u, \quad U = U' \bar u, \quad \varGamma_{p} = \varGamma_p^\ast D_s. \end{equation}

With these definitions, (2.3a)–(2.3d) can be written as

(2.5a)\begin{gather} - \nabla' p' + {\nabla}^2 \boldsymbol{u}' = 0, \end{gather}

(2.5b)\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}' =0, \end{gather}

(2.5c)\begin{gather}\frac{\partial s}{\partial t'} + \nabla' \boldsymbol{\cdot} \boldsymbol{u}' s = \frac{1}{Pe_s}{\nabla'}^2 s , \end{gather}

(2.5d)\begin{gather}\frac{\partial c}{\partial t'} + \nabla' \boldsymbol{\cdot} \left(\boldsymbol{u}' + \frac{\varGamma_p^\ast}{Pe_s} \nabla' s\right) c = \frac{1}{Pe_p}{\nabla'}^2 c , \end{gather}

where we used (2.1). The Péclet numbers, $Pe_s = \bar u \lambda / D_s$ and $Pe_p = \bar u \lambda / D_p$ compare the characteristic diffusion time scales $\tau _{D_s} = \lambda ^2 / D_s$ and $\tau _{D_p} = \lambda ^2 /D_p$ with the advection time $\tau _v$. Our study focuses on the interaction of pore-scale heterogeneity and DP while we ignore the interactions between the particles and the solid grains.

2.2. Porous medium

We consider a two-dimensional porous medium of length $L = 167 \lambda$ and width $W = 100 \lambda$ shown in figure 1. The medium characterised by DEPs of different size connected to a network of TPs, was generated using a solid-state de-wetting process (Salvalaglio et al. 2020; Bordoloi et al. 2022). It exhibits a complex pore structure interspersed among disordered solid grains. The medium is statistically homogeneous with a porosity of $\phi = 0.39$, and a narrow pore-size distribution of mean $\lambda =30\ {\rm \mu}$m. Solutes typically diffuse and dissipate gradients within pores over time scales shorter than the time $\tau _L = L/\bar u$ needed to elute a pore volume (also referred to as the mean breakthrough time). The junctures of TPs and DEPs serve as excellent candidates to retain gradients of solute concentration, which in turn trigger DP.

Figure 1. Hyper-uniform porous structure characterised by DEPs and TPs. Computational domain with white spaces indicating solid grains (bottom).

2.3. Numerical set-up

The domain is initially saturated with solute at uniform concentration $s_i=0.1$ mM and particles at uniform concentration $c_i=0.1$ mM. While $c_i$ may be unrealistically large, the specific value is not important here because particle concentration in the following is normalised by $c_i$. Thus, the results can be scaled to any desired (more realistic) initial value. Interactions of the particles with the solid matrix are not considered. At time $t > 0$, a sharp front of solute concentration $s_H=10$ mM is injected such that the ratio $\chi = s_i/s_H=10^{-2}$. Fluid flow is solved for no-slip boundary conditions at the fluid–solid boundary on the surface of the solid grains (see figure 1), by imposing the constant volumetric flux $U' = \phi$ at the left vertical boundary and constant pressure $p' = 0$ at the right boundary. For the solution of the solute and particle transport equations (2.3c) and (2.3d), constant flux is imposed at the inlet boundary, and zero gradient at the outlet. We use a finite element scheme to solve for the system of (2.3). For all the simulations shown here, the unstructured mesh is built using 197 723 triangular elements, the maximum size of which is $1\times 10^{-4}$ and minimum is $9\times 10^{-6}$. This corresponds to maximum and minimum resolutions of $3.33\lambda$ and $0.3\lambda$, respectively. An implicit second-order backward differentiation formulation with adaptive time stepping is used. We have checked that these results are numerically converged by performing simulations with up to 2 695 109 triangular elements with maximum and minimum element sizes of $1\times 10^{-4}$ and $6\times 10^{-6}$, respectively.

The dimensionless simulation parameters are listed in table 1. For a realistic solute diffusion coefficient of $D_s = 10^3$ $\mathrm {\mu }$m$^2$ s$^{-1}$, this simulation scenario corresponds to $D_p = 1$ $\mathrm {\mu }$m$^2$ s$^{-1}$, $U = 25$ $\mathrm {\mu }$m s$^{-1}$, $\varGamma _p \in [-1.2,0.6] \times 10^3$ $\mathrm {\mu }$m$^2$ s$^{-1}$ and valence of the $Z:Z$ electrolyte of $Z = 1$ with a prefactor in (2.2) of $\varepsilon (k_B T)^2/\mu (Z e)^2 = 470$ $\mathrm {\mu }$m$^2$ s$^{-1}$. Note that salt diffusion coefficients are typically of the order of $10^{3}$ $\mathrm {\mu }$m$^2$ s$^{-1}$, for instance, for LiCl and NaCl, $D_s = 1600$ $\mathrm {\mu }$m$^2$ s$^{-1}$ and $D_s = 1400$ $\mathrm {\mu }$m$^2$ s$^{-1}$, respectively. The value $D_p=1$ $\mathrm {\mu }$m$^2$ s$^{-1}$ corresponds to spherical colloidal particles with diameters of approximately $500$ nm at standard room temperature ($T = 300$ K), according to the Stokes–Einstein relation. To give a few examples, in some of the recent experiments conducted in microfluidic channels (Battat et al. 2019; Singh et al. 2020), polystyrene colloidal particles were driven by DP arising due to gradients in monovalent solutes like NaCl and LiCl. Typical values for the diffusiophoretic mobility are $|\varGamma _p| \lessapprox 1000$ $\mathrm {\mu }$m$^2$ s$^{-1}$ (Williams et al. 2024).

Table 1. Domain dimensions and parameters used in the numerical simulations.

3. Results and discussion

3.1. Diffusiophoretic particle trapping/extraction

Figure 2 shows the temporal evolution of the particle concentration field without DP ($\varGamma _p^\ast =0$), and for diffusiophoretic trapping ($\varGamma _p^\ast <0$), and extraction ($\varGamma _p^\ast >0$). In the absence of DP, the majority of the particles get dispersed through the TPs, leaving behind a small fraction of particles that accumulate within the DEPs, where the flow is stagnant, organised in convection rolls and described by closed streamlines (Bordoloi et al. 2022). The only mechanism through which these localised particles can escape into the TPs is diffusion, the time scale of which is typically orders of magnitude larger than $\tau _L$.

Figure 2. Temporal evolution of the dimensionless particle distributions $\log (c/c_i)$ for (a) trapping, (b) no DP and (c) extraction cases. White spaces indicate solid grains. Flow is from left to right. Note that $\tau _L/\tau _v \approx 167$ such that these plots correspond to times $0.5 \tau _L$, $1\tau _L$, $3 \tau _L$, and $10\tau _L$, approximately.

The injection of a sharp front of solute at a higher concentration results in local gradients of solute concentration that drive DP within the DEPs. For $\varGamma _p^\ast < 0$, that is, when particles migrate from high to low solute concentrations, DP leads to the trapping of particles inside the DEPs, as shown in figure 2(a).

When $\varGamma _p^\ast >0$, DP leads to particle mobilisation out of the DEPs because the particles move towards the higher solute concentration zones in the TPs. This is seen in the particle distributions shown in figure 2(c), where the particles within the DEPs rapidly escape the DEPs and leave the domain from the right outlet. Figure 3 shows the time evolution of the solute concentration at a point within a DEP close to the bottom of the DEP. The increase or decrease of concentration due to DP occurs on the time scale $\tau _{D_s}$, which here is smaller than $\tau _{D_p}$ by a factor of $10^3$. Thus, the action of DP on mass transfer between TPs and DEPs is limited to relatively short initial times.

Figure 3. Temporal evolution of dimensionless particle concentration ($c/c_i$) at an arbitrary point (indicated by the arrow) within a DEP shown in the inset for $\varGamma _p^\ast <0$ (blue, trapping), $\varGamma _p^\ast =0$ (black, no DP) and $\varGamma _p^\ast >0$ (red, extraction). The particle concentrations shown in the inset are the same as in figure 2. Note that time here is non-dimensionalised by $\tau _v$, where $\tau _L/\tau _v\approx 167$.

3.2. Breakthrough curves

Figure 4(a) shows the distribution of arrival times of particles at the outlet, also known as breakthrough curves. Similar to Bordoloi et al. (2022), we observe two distinct transport regimes for all the cases. At times of the order of $\tau _L$, particles at the outlet are produced by advection and dispersion from the TPs. For times $t \gg \tau _L$, the arrival time distribution deviates from the exponential decay predicted under the classical dispersion framework (denoted by the dotted curve in figure 4a) and displays a power law tailing. This is attributed to the particles that are initially trapped within the DEPs which can only escape the closed streamlines via diffusion, the time scale of which is typically large of the order of $\tau _{D_p} \gg \tau _v$ (here, $Pe = 1923$). For $\varGamma _p^\ast < 0$, a larger number of particles are trapped initially in the DEP, which manifests in a stronger tailing than the case without DP (indicated by black). For $\varGamma _p^\ast > 0$, particles are extracted from the DEP at initial times, and thus the tailing is weaker than for $\varGamma _p^\ast \leq 0$. Note that, albeit unphysical, extreme values of $|\varGamma _p^\ast |$ are chosen for the breakthrough curves in figure 4(a) to demonstrate the exaggerated effects of DP. The range of physically relevant values was discussed earlier in § 2.

Figure 4. (a) Arrival time distribution $F(t)$ at the outlet for (blue circles) $\varGamma _p^*=\varGamma _p/D_s = -1.2$, (black diamonds) $0$ and (orange circles) $0.6$. The dotted grey line indicates classical prediction under the advection-dispersion framework in (3.2), where the hydrodynamic dispersion is fitted to be $D_h = 7.43$. Solid lines correspond to the travel-time model in (3.1). Note that time here is made dimensionless using $\tau _v$, where $\tau _L/\tau _v\approx 167$. (b) Effective trapped particle fraction $\alpha$ in the DEPs from the (symbols) numerical data, and (solid line) the analytical expression (3.7) for $\ell _0^\ast = 1.6/Pe_p$.

The arrival time distribution computed using numerical simulations is predicted analytically by modelling this distribution as the superposition of the residence time distributions in the TPs and DEPs (Bordoloi et al. 2022)

(3.1)\begin{equation} F(t') = (1-\alpha) F_0(t')+ \alpha \int_0^\infty \text{d}\tau \frac{g(t'/\tau)}{\tau} f_D(\tau), \end{equation}

where $\alpha$ is the effective fraction of particles in the DEPs after particle redistribution due to DP at the short initial time interval, as discussed in the previous section. The first term on the right side of (3.1) describes transport in the TP by the solution of a one-dimensional advection-dispersion equation characterised by the mean flow velocity $\bar u = U/\phi$ and a constant hydrodynamic dispersion coefficient $D'_h$. The breakthrough curve at $x' = L'$ for a uniform distribution of particles in $0< x' < \infty$ is then given by (Kreft & Zuber 1978)

(3.2)\begin{equation} F_0(t') = 1 - \text{erfc}\left(\frac{L' - U't'/\phi}{\sqrt{4 D'_h t'}} \right) - \frac{D'_h \phi}{U'} \frac{\exp\left[\dfrac{(L' -U' t')^2}{4 D'_h t'} \right]}{\sqrt{4 D'_h t'}}. \end{equation}

The dimensionless hydrodynamic dispersion coefficient is obtained by fitting (3.2) to the numerical data. We find $D'_h = 7.43$. The dimensional dispersion coefficient is given by $D_h = U \lambda D_h'/\phi$. The second term on the right side of (3.1) corresponds to the breakthrough curve of the effective fraction $\alpha$ of particles in DEPs. The function $g(t')$ is given by the Gamma distribution

(3.3)\begin{equation} g(t') = \frac{t'^{{-}2/3}\exp ({-}t')}{\varGamma (1/3)}, \end{equation}

which denotes the residence time distribution in a DEP whose characteristic diffusion time is $\tau _\varLambda = 1$. The function $f_D(\tau )$ denotes the distribution of (dimensionless) diffusion times $\tau _\varLambda = \varLambda ^2 Pe_p$ within DEPs of (dimensionless) depth $\varLambda$, which is obtained from the distribution $f_\varLambda$ of depths through the map $\varLambda \to \varLambda ^2 Pe_p$. Equation (3.1) is used to model the numerical breakthrough curves by adjusting the effective trapped particle fraction $\alpha$ for different values of $\varGamma _p^\ast$.

As we see in figure 4(a), expression (3.1) is able to capture the early and late behaviours of the numerical breakthrough curves, but underestimates the arrival time distribution at intermediate times. This is due to the fact that expression (3.2) accounts for heterogeneity in the TPs by a constant hydrodynamic dispersion coefficient. This assumes that the characteristic advection times along TPs are much smaller than the mean arrival time at the outlet. The intermediate tailing observed here indicates that the heterogeneity of the TPs gives rise to larger advection times. While these effects can be quantified systematically in the framework of continuous time random walks (Dentz et al. 2018), we do not account for them in the current model because we focus here on the impact of DP on the long-time tailing of colloids, which is caused by mass transfer between TPs and DEPs. Figure 4(b) shows $\alpha$ vs $\varGamma _p^\ast$. For $\varGamma _p^\ast < 0$ particles are trapped in the DEPs by DP as the medium is flushed. Thus, for $\varGamma _p^\ast < 0$, $\alpha$ is larger than the trapped particle fraction $\alpha _0$ for the no-DP case ($\varGamma _p^\ast =0$). The tailing of the breakthrough, which is caused by particle retention in DEPs, is stronger than for $\varGamma _p^\ast = 0$. For $\varGamma _p>0$, particles are extracted from the DEPs when the medium is flushed. Thus, the effective trapped particle fraction is $\alpha < \alpha _0$. As a consequence, the tailing is weaker than for the no-DP scenario. The physical role of DP is to alter the effective trapped particle fraction $\alpha$ in the DEPs, and thus the tails of the particle breakthrough curves. This behaviour is correctly captured by (3.1).

3.3. Analytical model to understand how DP controls particle transport

To understand how DP controls the macroscopic fate of the suspended particles, we now estimate analytically the effective trapped particle fraction $\alpha$. To this end, we consider a single DEP of length $\ell _p$ connected to a TP with a width $\lambda$. We provide here an outline of the derivation of the analytical model. Further details can be found in the Appendix A. Note that, for clarity, the derivation is in dimensional terms, while the final result for the effective trapped fraction $\alpha$ is again in dimensionless terms.

We assume that particle motion in the DEPs is dominated by the diffusiophoretic drift when the salt gradients are large and represent particle transport by an advection equation. This allows us, by integration, to determine the total mass $m_{dp}$ of particles in the DEPs as a function of the particle concentration $c_0(t)$ and diffusiophoretic drift $u_0(t)$ at the interface between TPs and DEPs

(3.4)\begin{equation} m_{dp} = m_i + w \int_0^\infty {\rm d}t \,u_{0}(t) c_0(t), \end{equation}

where $m_i$ is the initial particle mass and $w$ the width of the interface between TP and DEP. Note that $\alpha = \alpha _0 m_{dp}/m_i$, where $\alpha _0$ is the effective trapped particle fraction for the no-DP case. In order to determine the diffusiophoretic drift at the interface, we estimate the salt gradient at the interface. We find an analytical solution by using the fact that salt transport in the DEPs is dominated by diffusion. From that and the definition (2.1) of the diffusiophoretic drift, we obtain the expression

(3.5)\begin{equation} u_0(t) ={-} \frac{\varGamma_p (1 - \chi) }{\sqrt{4 {\rm \pi}D_s t}} H(\tau_{D_s} {\rm \pi}- t), \end{equation}

where $H(t)$ denotes the Heaviside step function and $\tau _{D_s}$ is the characteristic salt diffusion time across the DEP. In order to estimate the particle concentration $c_0(t)$ at the interface, we distinguish between particle extraction ($\varGamma _p > 0$) and particle trapping ($\varGamma _p < 0$). In the former case, we set $c_0(t) = c_i$ equal to the resident particle concentration in the TP because the concentration of particles that are advected past the DEPs is constant. In the case of trapping ($\varGamma _p<0$), this is different. In this case, the particle concentration $c_0(t)$ at the interface is determined by the balance of the diffusive flux of particles from the TP toward the interface, and the diffusiophoretic flux from the interface into the DEP

(3.6)\begin{equation} - D_p \frac{c_0(t) - c_i}{\ell_0} = u_{0}(t) c_0(t). \end{equation}

The characteristic gradient scale $\ell _0$ can be estimated by equating the advective flux past the DEP and the diffusive flux toward it, which gives $\ell _0 \sim D_p/U$. Combining these relations, we obtain the following expression for $\alpha$ as a function of $\varGamma _p^\ast = \varGamma _p / D_s$:

(3.7)\begin{align} \alpha &= \alpha _0 [1 - \varGamma_p^\ast (1-\chi)] \nonumber\\ & \quad +H(-\varGamma_p^\ast) \left\{ \frac{{2 \varGamma_p^\ast}^2 (1 - \chi)^2 Pe_{p} \ell_0^\ast}{{\rm \pi} Pe_s} \ln\left[\frac{2\varGamma_p^\ast (1 - \chi) Pe_{p} \ell_0^\ast}{2\varGamma_p^\ast (1 - \chi) Pe_{p} \ell_0^\ast - {\rm \pi}Pe_s } \right] \right\}. \end{align}

The dimensionless gradient scale $\ell _0^\ast = \ell _0/\lambda$ is inversely proportional to the particle Péclet number such that we set $\ell _0^\ast = \varSigma /Pe_p$ with $\varSigma$ a number of order one. The latter is the only fitting parameter of the derived model, which fits the data perfectly for $\varSigma = 1.6$, as shown in figure 4(b).

As mentioned above, for positive $\varGamma _p^\ast > 0$ particles are extracted from the DEPs. In fact, by setting $\alpha = 0$ in expression (3.7), we find that particles are depleted from the DEPs for diffusiophoretic mobilities that are larger than the critical value

(3.8)\begin{equation} \varGamma_c^\ast = 1/(1 - \chi), \end{equation}

as shown in figure 4(b), where $\varGamma _c^\ast \sim 1$. This value corresponds to the extreme case of extraction, where the particles in the DEP will be completely depleted, i.e. when $\alpha =0$. In the opposite case of negative $\varGamma _p^\ast < 0$, DP leads to the trapping of particles in the DEPs. As $\varGamma _p^\ast$ decreases toward more negative values, the fraction $\alpha$ converges toward the asymptotic value

(3.9)\begin{equation} \alpha_{\infty} = \alpha_0 \left(1 + \frac{{\rm \pi} Pe_s}{\varSigma}\right), \end{equation}

which is obtained by taking the limit of $\varGamma _p^\ast \to \infty$. This result shows that the fraction of trapped particles cannot grow without bound. This is due to the fact that, for strongly negative diffusiophoretic mobilities, the rate at which particles are transferred to the interface between TPs and DEPs becomes smaller than the rate at which particles can be transferred into the DEPs by DP. In other words, the supply of particles to the interface is finite and does not increase as $\varGamma _p^\ast$ decreases. The analytical model describes the full dependence of $\alpha$ on $\varGamma _p^\ast$, as shown in figure 4(b) for both extraction from and trapping in DEPs, and thus seems to correctly capture the controls of DP on the macro-scale dispersion of particles.

4. Conclusions

To summarise, we investigate the effects of diffusiophoresis on particle dispersion in complex porous media. Using pore-scale simulations and analytical modelling, we have quantified the microscopic interactions between DP and flow and transport through porous media to show how this interplay impacts and governs the macroscopic fate of the colloidal particles. More precisely, we have demonstrated how DP occurring locally and microscopically with typically small time scales $\tau _{D_s}$, impacts the macroscopic particle breakthrough curves. We show that the physical role of DP is to alter the effective trapped particle fraction within the DEPs, referred to as $\alpha$. This is quantified using a novel analytical model derived in (3.7) to give a relation between $\alpha$ and the diffusiophoretic mobility $\varGamma _p$. Despite being a localised and short-term phenomenon, DP affects the long-term distribution of the particle arrival times by reorganising the local partitioning of particles between the DEPs and TPs. Depending on the sign of $\varGamma _p$, DP can facilitate particle mobilisation out of the DEPs or it may lead to particle entrapment within. The main implication of our results is that it is possible to access regions in the complex pore space that were initially inaccessible by merely exploiting inherent salt concentration landscapes. This suggests that DP can be used as a tool for controlling particle dispersion and filtration in complex porous media.

The potential of DP to control particle migration in simple microfluidic channels and cavities has been demonstrated and well established in the literature (Shin et al. 2016; Lee et al. 2020; Singh et al. 2020; Somasundar et al. 2023). However, despite the ubiquity of engineered and natural porous media in environmental and industrial applications, only little has been known about the impact of DP on macro-scale particle dispersion in complex disordered media. Our work is a first step towards closing this gap by shedding light on the impact of the interaction of DP, medium structure and flow heterogeneity on large-scale particle migration. Therefore, this work not only advances our fundamental understanding but also opens avenues for developing solutions to various technological problems of socio-economic relevance including groundwater remediation and microfluidics for biomedical applications, where DEPs are quite prevalent.

Lastly, it is important to note that this study considers the thin Debye layer limit, and DP that arises from monovalent salts e.g. LiCl and NaCl. In practical scenarios it is likely that the fluid composition is more complex. More investigations are required to understand DP under the influence of multivalent (Wilson et al. 2020) multiple salts (Alessio et al. 2021) and for finite Debye layer (Kirby & Hasselbrink 2004; Shin et al. 2016). In such scenarios, the diffusiophoretic mobility takes a more complex form than the one studied here. Furthermore, the present study does not consider stirring effects arising from diffusioosmosis, that is, bulk flow adjacent to the solid walls of the porous domain, which may be expected if the solid walls are themselves charged. Also, it needs to be pointed out that the model medium under consideration is two-dimensional with explicit dead-end regions. While DEPs in three-dimensional granular media and sphere packs are not likely to occur, they are an important transport-relevant medium property in many three-dimensional porous rocks (Coats & Smith 1964; Fatt, Maleki & Upadhyay 1966) and aggregated media (Philip 1968), and in general in porous materials characterised by a bi-porous structure, such as woven fabrics (Shin et al. 2018). Thus, we expect the reported results to be relevant also for particle transport in three-dimensional porous media.