1. Introduction
Predicting and controlling colloid transport in porous media is essential across diverse fields, from targeted drug delivery (Manzari et al. Reference Manzari, Shamay, Kiguchi, Rosen, Scaltriti and Heller2021; Mitchell et al. Reference Mitchell, Billingsley, Haley, Wechsler, Peppas and Langer2021), to managing the spread of microplastics and contaminants in subsurface environments and coastal aquifers (Alimi et al. Reference Alimi, Budarz, Jeffrey, Laura and Tufenkji2018; Brewer, Dror & Berkowitz Reference Brewer, Dror and Berkowitz2021; Pahlavan Reference Pahlavan2024; Spielman-Sun et al. Reference Spielman-Sun, Boye, Dwivedi, Engel, Thompson, Kumar and Noël2024). Chemical gradients, prevalent in natural and engineered porous systems (Dentz et al. Reference Dentz, Borgne, Englert and Bijeljic2011, Reference Dentz, Hidalgo and Lester2023; Berkowitz et al. Reference Berkowitz, Dror, Hansen and Scher2016; Rolle & Le Borgne Reference Rolle and Le Borgne2019), can induce colloidal motion through diffusiophoresis (Anderson Reference Anderson1989; Velegol et al. Reference Velegol, Garg, Guha, Kar and Kumar2016; Marbach & Bocquet Reference Marbach and Bocquet2019; Shim Reference Shim2022; Ault & Shin 2024). For example, ionic gradients resulting from agricultural practices, such as pesticide application and irrigation, can mobilise microplastics or natural colloids and clay particles (Müller et al. Reference Müller, Magesan and Bolan2007). These mobilised colloids can facilitate the transport of contaminants, impacting groundwater quality, ecosystem health and subsurface biodiversity (McCarthy & Zachara Reference McCarthy and Zachara1989; de Jonge, Kjaergaard & Moldrup Reference de Jonge, Kjaergaard and Moldrup2004; Sen & Khilar Reference Sen and Khilar2006; Weber et al. Reference Weber, Voegelin, Kaegi and Kretzschmar2009; Sheng et al. Reference Sheng, Jing, He, Klein, Köhler and Wanger2024). A critical open question is whether solute gradients significantly influence the macroscopic transport of colloids in these complex environments.
Our understanding of the fundamental aspects of solute-mediated migration of colloids in idealised scenarios characterised by sharp solute gradients and no background flows, e.g. in dead-end pores, has significantly improved over the past decade (Shin et al. Reference Shin, Um, Sabass, Ault, Rahimi, Warren and Stone2016; Battat et al. Reference Battat, Ault, Shin, Khodaparast and Stone2019; Gupta, Shim & Stone Reference Gupta, Shim and Stone2020; Wilson et al. Reference Wilson, Shim, Yu, Gupta and Stone2020; Alessio et al. Reference Alessio, Shim, Mintah, Gupta and Stone2021, Reference Alessio, Shim, Gupta and Stone2022; Shim Reference Shim2022; Akdeniz, Wood & Lammertink Reference Akdeniz, Wood and Lammertink2023; Lee, Lee & Ault Reference Lee, Lee and Ault2023). However, in porous medium flows, geometric disorder weakens the solute gradients (Koch & Brady Reference Koch and Brady1985; Sahimi Reference Sahimi1993; Dentz et al. Reference Dentz, Icardi and Hidalgo2018, Reference Dentz, Hidalgo and Lester2023). It is therefore not a priori clear if diffusiophoresis could play an important role in these complex environments.
Recent experimental and theoretical works offer evidence for the role of diffusiophoresis in porous environments. In the absence of background flows, colloidal diffusiophoresis has been experimentally demonstrated and theoretically rationalised in fibrous porous media such as collagen gels and biofilms, where pore sizes are comparable to the colloid size (Doan et al. Reference Doan, Chun, Feng and Shin2021; Sambamoorthy & Chu Reference Sambamoorthy and Chu2023; Somasundar et al. Reference Somasundar, Qin, Shim, Bassler and Stone2023; Sambamoorthy & Chu Reference Sambamoorthy and Chu2025). In the presence of flows, microfluidic experiments and numerical simulations have illustrated how diffusiophoresis influences colloidal transport into and out of dead-end pores (Park et al. Reference Park, Lee, Yoon and Shin2021; Jotkar et al. Reference Jotkar, de Anna, Dentz and Cueto-Felgueroso2024a ). More recently, we have shown how the interplay between flow disorder and solute gradients governs the transport of colloids in porous media, highlighting the role of diffusiophoresis in both preferential flow pathways and stagnant fluid pockets (Alipour et al. Reference Alipour, Li, Liu and Pahlavan2026).
Despite these advances, it remains unclear whether diffusiophoretic removal of colloids from stagnant fluid pockets remains effective over larger scales, as solute fronts become increasingly diffuse, weakening spatial gradients and diffusiophoretic velocities. Here, we address this central question: Does dispersion of the solute front weaken the phoretic effects? To address this question, here we combine numerical simulations and microfluidic experiments, demonstrating that solute front dispersion in disordered porous media does not diminish colloid removal efficiency. To rationalise this observation, we analyse diffusiophoretic colloid migration from idealised one-dimensional dead-end pores exposed to diffuse solute fronts, deriving an analytical model to describe the spatio-temporal evolution of solute and colloid fields. Our model shows excellent agreement with numerical simulations for linear and error-function solute profiles, showing that broadening of the solute front leads to a persistent spatial gradient beyond the solute-diffusion time scale along the dead-end pore, leading to a more efficient and uniform extraction of colloids. Our work therefore points to the potential relevance of diffusiophoresis in realistic porous medium scenarios.
2. Influence of solute dispersion on phoretic colloid removal from porous media
2.1. Experimental set-up
We probe the transport of colloids in a quasi-2-D microfluidic geometry with height
$h \approx 50\ {\unicode{x03BC}} \text{m}$
, patterned with a disordered array of cylindrical obstacles with diameter
$\approx 175\ {\unicode{x03BC}} \text{m}$
. The mean pore size is
$\approx 60\ {\unicode{x03BC}} \text{m}$
and the medium porosity is
$\approx 0.36$
. The medium is initially filled with an aqueous colloidal solution with density
$n_0=1$
and solute concentration
$c_0=1\,\textrm {mM}$
. We then inject a second aqueous solution with solute concentration
$c_1=100\,\textrm {mM}$
at constant flow rate to displace the colloids and monitor the removal process in different fields of view along the medium (figure 1
a). In the experiments, we use LiCl as the solute with diffusivity of
$D_{{s}} = 1.37 \times 10^{-9}\,\mathrm{m}^2\,\mathrm{s}^{- 1}$
, and 1 micron size colloids with diffusivity
$D_{{p}} = 4.37 \times 10^{-13}\,\mathrm{m}^2\,\mathrm{s}^{- 1}$
, and diffusiophoretic mobility of
$\varGamma _{{p}} \approx 0.7 \times 10^{-9}\,\mathrm{m}^2\,\mathrm{s}^{- 1}$
(Shin et al. Reference Shin, Um, Sabass, Ault, Rahimi, Warren and Stone2016; Velegol et al. Reference Velegol, Garg, Guha, Kar and Kumar2016; Gupta et al. Reference Gupta, Shim and Stone2020; Alessio et al. Reference Alessio, Shim, Gupta and Stone2022; Alipour et al. Reference Alipour, Li, Liu and Pahlavan2026; Liu & Pahlavan Reference Liu and Pahlavan2025). We note that the colloidal suspension is very dilute with a volume fraction of
$\approx 0.1$
%, leading to a packing fraction of
$\approx 10^{-3}$
, which together with the
$1/r^3$
decay of hydrodynamic disturbance velocity for diffusiophoretic colloids, implies that we can safely ignore particle–particle or particle–surface interactions (Shin et al. Reference Shin, Um, Sabass, Ault, Rahimi, Warren and Stone2016).
Figure 1. The dispersion of a solute front does not weaken the diffusiophoretic removal of colloids. (a–c) Numerical simulations; (d) simulations and experiments. (a) We displace an aqueous solution of colloids with salt concentration
$c_0$
with a colloid-free solution with salt concentration
$c_1$
at constant flow rate. Medium disorder leads to the heterogeneous distribution of flow velocity magnitude, solute concentration, diffusiophoretic velocity magnitude and particle density. Three fields of view (FOVs) are highlighted by white boxes. (b) Cross-sectionally averaged profiles of the normalised solute concentration (red dashed), particle density (black solid) and diffusiophoretic velocity magnitude (purple dash-dotted), shown at the same time as panel (a). (c) Temporal evolution of the average solute concentration and diffusiophoretic velocity magnitude within the three FOVs. The solute profiles are fitted using the solution to the one-dimensional (1-D) advection–diffusion equation, (2.3). The inset presents the fitted characteristic solute-transition time
$\tau _{\textit{erf}}$
for each FOV. (d) Temporal evolution of the average particle density across the three FOVs, in the presence (single lines), and absence of solute gradients (double solid lines). Symbols represent the corresponding experimental observations.
2.2. Numerical simulations
We use OpenFOAM to numerically solve for the evolution of solute and colloid fields (Ault et al. Reference Ault, Shin and Stone2018, Reference Ault, Shin and Stone2019; Alipour et al. Reference Alipour, Li, Liu and Pahlavan2026; Liu & Pahlavan Reference Liu and Pahlavan2025). We calculate the steady-state flow field using the simpleFoam solver, applying a pressure drop across the medium, which results in a mean velocity magnitude of
$\langle U_{\textit{mag}} \rangle = 184.3\ {\unicode{x03BC}} \rm {m\,s^{- 1}}$
. The geometric disorder of the medium leads to a heterogeneous velocity field, where most of the transport occurs through preferential flow pathways surrounding low permeability/stagnant fluid pockets (figure 1
a). With the steady-state flow field established, we then solve for the transient transport of solute,
$c$
, and particles,
$n$
, governed by the advection–diffusion equations
\begin{align} \frac {\partial c}{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }(\boldsymbol{u} c) & = D_{{s}} {\nabla} ^2 c, \end{align}
\begin{align} \frac {\partial n}{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }\left ( (\boldsymbol{u} + \boldsymbol{u}_{\textit{DP}}) n \right ) & = D_{{p}} {\nabla} ^2 n, \end{align}
where
$\boldsymbol{u}$
is the background velocity, and
$\boldsymbol{u}_{\textit{DP}} = \varGamma _{{p}} \boldsymbol{\nabla }\ln c$
represents the diffusiophoretic velocity. The recent work of Jotkar et al. (Reference Jotkar, de Anna, Dentz and Cueto-Felgueroso2024a
) also studies the diffusiophoretic transport of colloids in porous media using these equations, reporting on the migration of colloids into/from dead-end pores due to solute gradients.
2.3. Flow disorder enhances solute dispersion
Flow disorder leads to a non-uniform solute front, enhancing the dispersion (Saffman Reference Saffman1959; Koch & Brady Reference Koch and Brady1985; Koch et al. Reference Koch, Cox, Brenner and Brady1989; Sahimi Reference Sahimi1993; Dentz, Icardi & Hidalgo Reference Dentz, Icardi and Hidalgo2018; Puyguiraud, Gouze & Dentz Reference Puyguiraud, Gouze and Dentz2021; Liu et al. Reference Liu, Xiao, Aquino, Dentz and Wang2024; Woods Reference Woods2025) and weakening the diffusiophoretic velocities (figure 1
a). The broadening of the solute profile leads to the persistence of solute gradients around the stagnant fluid pockets (de Anna et al. Reference de Anna, Pahlavan, Yawata, Stocker and Juanes2021) as evidenced in the snapshot of the diffusiophoretic velocity field. The cross-sectionally averaged fields show the extent of dispersion along the medium (figure 1
b). We note that the peak of the diffusiophoretic velocity is ahead of the solute front due to the
$\boldsymbol{\nabla }c / c$
functional form.
As the solute front travels along the medium, it becomes more dispersed, as demonstrated by the evolution of mean solute concentration in the three FOVs along the medium (figure 1 c). To characterise this dispersion, we fit the solute profiles using the analytical solution of the 1-D advection–diffusion equation
\begin{equation} c(t) = c_0 + \frac {c_1 - c_0}{2} \left ( 1 + \text{Erf} \left ( \frac {U_{{x}} t - x_{\textit{FOV}}}{\sqrt {4D_{\textit{eff}} t}} \right ) \right )\!, \end{equation}
where
$U_{{x}} = 141\ {\unicode{x03BC}} \rm{m\,s^{-1}}$
is the average velocity in the mean flow direction. Here, the effective distance from the inlet,
$x_{\textit{FOV}}$
, and the effective dispersion coefficient,
$D_{\textit{eff}}$
, are determined as fitting parameters. We can then define the characteristic solute-transition time for the error-function profile as
\begin{equation} \tau _{\textit{erf}} = \frac {2\sqrt {4D_{\textit{eff}} (x_{\textit{FOV}} / U_{{x}})}}{U_{{x}}}. \end{equation}
We observe that
$\tau _{\textit{erf}}$
increases by
$60\,\%$
, from
$41.7\ \text{s}$
in FOV1 to
$68.6\ \text{s}$
in FOV3, over a distance of
$10\ \text{mm}$
, which corresponds to approximately
$160$
pore sizes (inset of figure 1
c). These results highlight the significant dispersion of the solute front as it travels through the porous medium, further emphasising the influence of the disordered structure on solute transport. In figure 1(c), the evolution of the average magnitude of the diffusiophoretic velocity in each FOV is plotted with semi-transparent dashed lines. As expected, the increasingly dispersed solute front generates a weaker but more persistent diffusiophoretic velocity field in downstream FOVs compared with upstream regions.
2.4. Solute dispersion does not weaken the phoretic removal
To examine the influence of diffusiophoresis on particle transport, we monitor the temporal evolution of the average particle density,
$n(t)/n_{0}$
, for each FOV (figure 1
d). The results are represented by thick black, red and blue lines. We further conduct simulations and experiments for the control case, where no solute gradients are present (
$c_0=c_1$
). The corresponding average particle-density profile in this control case is plotted as double lines. In the control case, the average particle density decreases rapidly from its initial value of
$1.0$
to approximately
$n_{f}/n_0 \approx 0.015$
, where
$n_0$
and
$n_f$
represent the initial and final particle densities, respectively. This initial drop is primarily due to the displacement of particles in main flow pathways. Subsequently, the particle density decreases more gradually, as diffusion becomes their only escape route from dead-end pores.
In the presence of solute gradients, the colloid removal becomes much more efficient as diffusiophoresis removes the colloids from dead-end pores and lead to their cross-streamline migration in the main flow pathways (figure 1
d) (Alipour et al. Reference Alipour, Li, Liu and Pahlavan2026). We might expect the sharper solute front in the upstream region (FOV 1) to lead to a more efficient particle removal than the downstream region (FOV 3). Our observations, however, clearly show this not to be the case as the final particle fraction after the passage of solute front is the same in all three FOVs (
$n_{f}/n_0 \approx 0.001$
, as shown in figure 1
d). Our numerical simulations agree reasonably well with our microfluidic experiments shown by the symbols in figure 1(d). To gain insight into this surprising observation, we next focus on the idealised case of 1-D dead-end pores, probing how solute dispersion influences the diffusiophoretic removal of colloids from stagnant fluid pockets.
3. Diffusiophoretic removal of colloids from dead-end pores
3.1. Sharp vs diffuse solute fronts: experiments and 2-D simulations
3.1.1. Experiments
To probe the role of solute dispersion on the removal of colloids from dead-end pores, we conduct two sets of experiments with sharp and diffuse solute fronts. In both, we first fill the chip with an aqueous colloidal solution with solute concentration
$c_0=1\,\textrm {mM}$
. We then displace this solution with a second colloid-free aqueous solution with salt concentration
$c_1=100\,\textrm {mM}$
(figure 2
a,b). In each chip, six dead-end pores of length
$600\,{\unicode{x03BC}} \textrm {m}$
and width
$100\,{\unicode{x03BC}} \textrm {m}$
are connected to the main channel. To create a sharp step-like solute-concentration gradient at the inlet of dead-end pores, we use an air bubble to separate the two solutions (Shin et al. Reference Shin, Um, Sabass, Ault, Rahimi, Warren and Stone2016; Wilson et al. Reference Wilson, Shim, Yu, Gupta and Stone2020; Alessio et al. Reference Alessio, Shim, Gupta and Stone2022). In the absence of this bubble, the interface between the two solutions becomes diffuse, leading to a slowly varying solute-concentration profile. We added a trace amount of fluorescein dye to the colloid-free solution to probe the evolution of the solute concentration at the inlet of dead-end pores, demonstrating the distinction between the sharp and diffuse solute profiles (figure 2
c). Here, we use LiCl as the solute, matching the porous medium experiments above. We also repeated the experiments with NaCl, obtaining similar results (Appendix A.3).
Figure 2. The evolution of colloid density in dead-end pores exposed to (a) sharp and (b) diffuse solute fronts. (c) The dispersion of the solute front is represented by the evolution of the fluorescein light intensity at the inlet of the dead-end pores. (d) The evolution of colloid density exposed to sharp versus diffuse solute fronts. For both cases, snapshots at
$t/\tau _{{s}}=0,\,2,\,4$
. (e) The evolution of fraction of residual particles over time showing the cross-over at ‘late times’. Symbols: experiments; solid lines: 2-D simulations. The blue dash–dot line indicates the analytical prediction of the final residual (4.22), consistent with both. In the presence of sharp solute front, colloids migrate toward the inlet of the dead-end pore, leading to a non-uniform density profile (f). However, in the diffuse-front case, colloidal density in the pore remains nearly uniform (g). For the sharp front, data points correspond to
$t/\tau _{{s}}=0,\,0.1,\,0.2,\,0.5,\,1,\,2$
, and for the diffuse front, we have
$t/\tau _{{s}}=0,\,0.25,\,0.5,\,1,\,2,\,4$
. Solid lines represent smoothed density profiles obtained via a five-point adjacent-averaging method. (h) Two-dimensional simulation domain. Colour map: particle density after the flushing stage; streamlines: background flow showing the primary recirculation in the pore. (i, j) Simulated particle-density profiles for the sharp (i) and diffuse (j) fronts at the same sampling times as in (f, g).
Surprisingly, displacement using a diffuse solute front leads to a more effective removal of the colloids (figure 2
d). In the sharp-front case, colloids migrate toward the inlet of the dead-end pore, but not all can escape before the solute gradients vanish (
$\tau _s=L^2/D_s \approx 250\,\textrm {s}$
). In the diffuse-front case, however, the removal of particles is more uniform. This is evident in 1-D colloid density profiles constructed using particle detection and averaged over the six dead-end pores (figure 2
f, g). The lines show an interpolation using the five-point adjacent-averaging method.
The evolution of particle density in dead-end pores clearly shows that the sharp solute front is more effective in removing the colloids at ‘early times’. However, we observe a cross-over at ‘late times’ to the diffuse solute front becoming the more efficient way to remove particles from dead-end pores (figure 2
e). In the sharp solute front case, solute gradients disappear beyond
$t=\tau _s$
, and particle diffusivity is negligible over the time scales probed here since
$L^2/D_p\gg L^2/D_s$
. Colloids can migrate in dead-end pores as long as solute gradients persist. Therefore, we expect the persistence of solute gradients in the diffuse-front case to lead to the cross-over at late times.
3.1.2. Numerical simulations
We solve the coupled solute and colloid transport in a 2-D dead-end pore attached to a straight main channel (figure 2
h) using OpenFOAM. The pore length and width are
$L=600\,{\unicode{x03BC}} \mathrm{m}$
and
$W=100\,{\unicode{x03BC}} \mathrm{m}$
, matching the experiments; the main channel has the same in-plane dimensions. A symmetry condition at the top wall halves the computational channel height relative to the experiments. The background flow is left to right with mean speed
$U_{ m}=2\,\mathrm{mm\,s^{-1}}$
(measured experimentally). We initialise the particle number density uniformly at
$n_0$
, apply a short flushing stage to clear high-velocity regions (the main channel and the first primary eddy in the pore), and then impose a time-dependent inlet solute history
$c_{\textit{in}}(t)$
at the left boundary. For a sharp front we impose a step from
$c_0=1\,\mathrm{mM}$
to
$c_1=100\,\mathrm{mM}$
. For a diffuse front we use the fluorescein-based inlet profile measured in the experiments (figure 2
c), but rescale its time axis by the factor
$\sqrt {D_{{f}}/D_{{s}}}$
to account for the different molecular diffusivities of salt and fluorescein, where
$D_{{s}}=1.37\times 10^{-9}\,\mathrm{m^2\,s^{-1}}$
and
$D_{{f}}\approx 4.6\times 10^{-10}\,\mathrm{m^2\,s^{-1}}$
. This scaling follows from high-Péclet-number dispersion, for which an effective coefficient of the form
$D_{\textit{eff}}/D\sim {Pe}^2$
yields a front width (and thus a transition time)
$\Delta t \sim (D_{\textit{eff}})^{1/2}/U_{\textit{m}} \propto D^{-1/2}$
; hence, the fluorescein profile is compressed by
$\sqrt {D_{{f}}/D_{{s}}}$
to represent the salt front. Here, we take the diffusiophoretic mobility to be constant,
$\varGamma _{p}=0.7\times 10^{-9}\,\mathrm{m^2\,s^{-1}}$
. We also repeat the simulations with the constant–surface–charge (CSC) model (Lee et al. Reference Lee, Lee and Ault2023), which yields similar results (Appendix A.4). Other transport parameters are as specified above; in particular, the particle diffusivity is
$D_{p}=4.37\times 10^{-13}\,\mathrm{m^2\,s^{-1}}$
.
As shown by the solid curves in figure 2(e), the simulations reproduce the experimental ‘cross-over’ in the residual fraction
$N/N_0$
: sharp fronts remove particles more rapidly at early times, whereas diffuse fronts ultimately leave fewer particles. In both cases, the agreement with the measurements is close, and for the diffuse front the experiment and simulation converge to the same final residual once the solute gradient has vanished. The profiles in figure 2(i–j) further illustrate this behaviour: the sharp front generates a non-uniform
$n/n_0$
with a peak near the inlet, while the diffuse front yields an almost uniform profile with the same final residual; both profiles agree reasonably with the experiments in figure 2(f–g).
3.2. Sharp vs diffuse solute fronts: 1-D model with a linear solute profile
To elucidate why a diffuse solute front can enhance particle removal, we further simplify our model to a 1-D dead-end pore of length
$L$
, with the inlet located at
$x = 0$
(figure 3
a). The influence of flow within the dead-end pore is assumed negligible, which is a reasonable approximation for pores with high aspect ratio. The initial solute concentration and particle density are set to
$c_0 = 1$
and
$n_0 = 1$
, respectively, with boundary conditions at the inlet
$c_1 = c_{\textit{in}}(t)$
and
$n_1 = 0$
.
Figure 3. Transport of solute and particles in a 1-D dead-end pore. (a) Schematic of the set-up: the initial solute concentration and particle density inside the pore are
$c_0 = 1$
and
$n_0 = 1$
, respectively. The inlet boundary conditions are given by
$c_{\textit{in}}(t)$
and
$n_1 = 0$
. As shown on the right,
$c_{\textit{in}}(t)$
increases linearly from 1 to 100 over the transition time
$T$
and remains at 100 thereafter. (b) Maximum diffusiophoretic velocity over time for different normalised solute-transition times
$T/\tau _{{s}} \in \{0,\,1,\,10\}$
(red: 0, yellow: 1, blue: 10). (c, d) Evolution of the solute field (c) and diffusiophoretic velocity field (d) for
$T/\tau _{{s}} = 0,\,1,\,10$
. Sampling times are
$t/\tau _{{s}} = \{0.001,\,0.02,\,0.05,\,0.1,\,1\}$
for
$T/\tau _{{s}} = 0$
,
$\{0.01,\,0.1,\,0.5,\,1,\,2\}$
for
$T/\tau _{{s}} = 1$
and
$\{0.1,\,0.5,\,1,\,2,\,10\}$
for
$T/\tau _{{s}} = 10$
. (e) Particle trajectories
$x(t,\xi _0)$
as functions of both time
$t$
and initial position
$\xi _0$
, using the same sampling times. (f) Particle-density fields: non-diffusive particles from particle tracking using (3.4) (solid lines) and diffusive particles from the continuum model in (3.5) (symbols), evaluated at the same sampling times.
For simplicity, we first impose a linearly varying inlet solute-concentration profile
$c_{\textit{in}}(t)$
that increases from
$c_0$
to
$c_1$
over a duration
$T$
, as shown in figure 3(a). The evolution of the solute-concentration field is governed by the diffusion equation
\begin{equation} \frac {\partial c\left (x,t\right )}{\partial t} - D_{{s}}\frac {\partial ^2 c\left (x,t\right )}{\partial x^2} = 0, \end{equation}
subject to the initial condition
$c(x,0) = c_0$
and boundary conditions
\begin{equation} c\left (0,t\right ) = \begin{cases} c_0 + \dfrac {c_1 - c_0}{T} \boldsymbol{\cdot }t, & t \in [0, T), \\ c_1, & t \in [T, \infty ), \end{cases} \qquad \frac {\partial c}{\partial x}(L,t) = 0, \end{equation}
where
$D_{{s}}$
is the solute diffusivity, and
$L$
is the pore length. Defining the characteristic solute-diffusion time as
$\tau _{{s}} = L^2 / D_{{s}}$
, we non-dimensionalise the system using
$\bar x = x / L$
,
$\bar t = t / \tau _{{s}}$
,
$\bar T = T / \tau _{{s}}$
and
$\bar c = (c - c_0) / (c_1 - c_0)$
. The resulting analytical solution for the non-dimensional solute field is
\begin{equation} \bar c\left (\bar x, \bar t\right ) = \begin{cases} \sum \limits _{n=0}^\infty \dfrac {2}{\lambda _n^3 \bar T} e^{-\lambda _n^2 \bar t} \sin \left (\lambda _n \bar x\right ) + \dfrac {1}{2\bar T} \left (\bar x^2 - 2\bar x\right ) + \dfrac {\bar t}{\bar T}, & \bar t \in [0, \bar T], \\[12pt] -\sum \limits _{n=0}^\infty \dfrac {2}{\lambda _n^3 \bar T} \left (1 - e^{-\lambda _n^2 \bar T}\right ) e^{-\lambda _n^2 (\bar t - \bar T)} \sin \left (\lambda _n \bar x\right ) + 1, & \bar t \in (\bar T, \infty ), \end{cases} \end{equation}
where
$\lambda _n = ( {(2n+1)\pi })/{2}$
.
Given that the particle diffusivity is significantly smaller than the solute diffusivity, we first assume particles to be non-diffusive. The trajectory
$\bar x(\bar \xi _0, \bar t)$
of a particle starting at
$\bar \xi _0 = \xi _0 / L$
is governed by
\begin{equation} \frac {\partial \bar x\left (\bar \xi _0, \bar t\right )}{\partial \bar t} = \bar U_{\textit{DP}}\left (\bar x, \bar t\right ) = \bar \varGamma _{{p}} \frac {\partial \ln \left (\bar c\left (\bar x, \bar t\right ) + \dfrac {1}{\beta - 1}\right )}{\partial \bar x}, \end{equation}
where
$\bar \varGamma _{{p}} = \varGamma _{{p}} / D_{{s}} \approx 0.7$
, and
$\beta = c_1 / c_0 = 100$
.
The evolution of the solute and diffusiophoretic velocity fields for cases with
$\bar T = 0,\ 1,\ 10$
are shown in figure 3(c,d), obtained from (3.3) and (3.4). For the
$\bar T = 0$
case, L’Hôpital’s rule is applied to evaluate the
$\bar T \rightarrow 0$
limit. In this case, the solute diffuses into the pore over time scale
$\approx \tau _s$
leading to a rapid decay of the diffusiophoretic velocity. In contrast, for
$\bar T \gg 1$
, the normalised solute-concentration profile approximately follows the form
$({1}/{\bar T})(\bar x^2 - 2\bar x) + ( {\bar t}/{\bar T})$
as the exponential terms in (3.3) decay rapidly (figure 8
a). This expression indicates that the solute concentration in the pore has a time-independent quadratic profile that increases at a constant rate to match the inlet boundary condition. This results in a weaker, but more persistent diffusiophoretic velocity compared with the sharp-gradient case (figure 3
b,d).
Using (3.3) and (3.4), we numerically determine the trajectories of non-diffusive particles,
$\bar {x}(\bar {\xi }_0, \bar {t})$
, for the three considered scenarios (figure 3
e). Over time, the
$x - \xi _0$
profiles shift downward, reflecting the diffusiophoretic migration of particles toward the inlet. The point where these profiles intersect the horizontal axis
$x/L=0$
at any given moment defines the critical depth (
$\xi _{0,c}$
): particles originating above this depth (
$\xi _0\lt \xi _{0,c}$
) have exited the pore while those below (
$\xi _0\gt \xi _{0,c}$
) remain trapped. The red line in each panel denotes the fraction of residual particles when the solute gradient approaches zero. Increasing
$T$
reduces the residual particle fraction, demonstrating that more dispersed solute fronts enhance particle extraction. Further, for larger
$T$
, the
$x - \xi _0$
profiles become more linear, indicating a more uniform particle distribution.
To assess the effect of finite particle diffusivity, we solve the 1-D advection–diffusion equation for particle density
\begin{equation} \frac {\partial n}{\partial \bar t} + \frac {\partial }{\partial \bar x}\left (\bar U_{\textit{DP}} n\right ) = \bar D_{{p}} \frac {\partial ^2 n}{\partial \bar x^2}, \end{equation}
where
$\bar D_{{p}} = D_{{p}} / D_{{s}} \approx 3.2 \times 10^{-4}$
is the normalised particle diffusivity. The initial condition is
$n(\bar x, 0) = n_0$
, and the boundary conditions are
$n(0, \bar t) = 0$
and
$( {\partial n}/{\partial \bar x})(1, \bar t) = 0$
. Figure 3(f) compares the non-diffusive particle-density fields (solid black lines) obtained from the Lagrangian trajectories
$n / n_0 = (\partial \bar {x} / \partial \bar {\xi }_0)^{-1}$
with the density fields obtained from (3.5), represented by the square symbols. The good agreement confirms that particle diffusivity has negligible influence on the time scales probed here, except in the vicinity of the pore inlet where a thin diffusion layer induces minor deviations.
For a step-like solute boundary condition (
$\bar {T} = 0$
), the particle-density profile is non-monotonic at early times, before the solute diffuses along the pore: particles in the middle region are rapidly drawn toward the inlet, whereas those near the bottom remain largely undisturbed. This results in a density peak near the inlet and a valley in the middle. As time progresses and the solute front dissipates, particles at the bottom also begin to move toward the inlet, ultimately producing a monotonic profile with a peak density of approximately
$n / n_0 \approx 0.8$
at
$\bar {t} = 1$
.
In contrast, when the solute concentration changes gradually
$(\bar {T} = 10)$
, the particle-density profile evolves smoothly, decreasing almost uniformly along the pore. By the time the diffusiophoretic velocity approaches zero, the profile is nearly flat, with a notably low density of
$n / n_0 \approx 0.03$
. This is consistent with our experimental observations (figure 2). Because the particle diffusion time
$\tau _{{p}} = L^2 / D_{{p}}$
is much larger than the solute-diffusion time
$\tau _{{s}}$
, once the solute gradients (and thus the diffusiophoretic velocity) vanish, the particle distribution remains effectively unchanged. Hence, the grey shaded areas in figure 3(f), represent the unrecoverable particle fraction.
3.2.1. Evolution of residual particle fraction
To gain insight into the effect of solute-concentration profile at the inlet, we use (3.3) and (3.4) to numerically solve for the evolution of residual particle fraction in the pores (figure 4
a). The residual fraction decreases initially before reaching a plateau at longer times when the solute gradients vanish. For larger
$T/\tau _{{s}}$
values, the extraction process is slower, but more particles are ultimately removed. Consequently, at any given normalised time
$t/\tau _{{s}}$
, there is an optimal
$T/\tau _{{s}}$
that maximises particle extraction efficiency and minimises the residual fraction. This optimum can be found numerically, and is represented by the cyan envelope in figure 4(a).
Figure 4. Longer solute-transition times slow the extraction of particles from dead-end pores but ultimately enhance their overall removal. (a) Temporal evolution of the fraction of residual particles in dead-end pores for different normalised solute-transition times (
$T/\tau _{{s}} = 0,\,1,\,10,\,100$
), derived from the non-diffusive particle trajectories (3.4), is shown by black lines. The cyan solid line is the envelope of these curves, indicating the minimal achievable residual particle fraction at any given time. (b) Final residual particle fraction (evaluated once the solute gradient has vanished) as a function of the solute-transition time. Results are obtained by solving continuous equations (3.5) for diffusive particles (diamond symbols) and via particle tracking (3.4) for non-diffusive particles (black solid lines). The magenta dash-dot line indicates the residual particle fraction at the moment when the inlet solute concentration reaches its final value, assuming particles are non-diffusive. The yellow solid line represents the prediction from our analytical model (4.26) that accounts for both diffusiophoresis and particle diffusion.
The final residual particle fraction,
$\alpha$
, decreases as
$T$
increases, from
$\alpha _{{H}}\approx 0.17$
to a lower plateau of
$\alpha _{{L}}\approx 0.03$
as shown by the solid black line in figure 4(b). This indicates that in the absence of particle diffusivity, a small fraction of particles remain unrecoverable regardless of how diffuse the solute front is. However, particle diffusion allows these particles to slowly diffuse out of the pores as shown by the diamond symbols obtained from the numerical solution of (3.5). Note that the diamond symbols and the solid line match very well for times much smaller than particle diffusion time scale
$T \ll \tau _{{p}} = L^2/D_{{p}}$
. (grey vertical dash-dot line).
4. Analytical model for particle extraction with slowly varying solute profiles
Consider a general, slowly varying inlet solute profile
$c_{\textit{in}}(\bar {t}) = C_{\textit{in}}(\bar {T}(\bar {t}))$
, where the stretched time variable
$\bar T = \delta \,\bar t$
with
$\delta \ll 1$
and the derivatives of
$C_{\textit{in}}$
with respect to
$\bar T$
are taken to be
$O(c_0)$
. Thus, over the characteristic solute-diffusion time
$\tau _{{s}}$
(i.e.
$\Delta \bar t=O(1)$
), the relative change of the inlet concentration is only
$O(\delta )$
, which defines the ‘slowly varying’ limit. The solute field inside the pore evolves according to
\begin{equation} \frac {\partial c(\bar {x},\bar {t})}{\partial \bar {t}} - \frac {\partial ^2 c(\bar {x},\bar {t})}{\partial \bar {x}^2} = 0, \end{equation}
with initial and boundary conditions
\begin{equation} c(\bar {x},0) = c_{\textit{in}}(0) = c_0,\qquad c(0,\bar {t}) = c_{\textit{in}}(\bar {t}),\qquad \frac {\partial c}{\partial \bar {x}}(1,\bar {t}) = 0. \end{equation}
To homogenise the boundary conditions, we define
\begin{equation} \phi (\bar {x}, \bar {t}) = c(\bar {x}, \bar {t}) - c_{\textit{in}}(\bar {t}) - \frac {1}{2}(\bar {x}^2 - 2\bar {x})c_{\textit{in}}^{\prime}(\bar {t}), \end{equation}
where
$c_{\textit{in}}^{\prime} = \mathrm{d}c_{\textit{in}}/\mathrm{d}\bar {t}$
. Substituting into (4.1) yields
\begin{equation} \frac {\partial \phi (\bar {x},\bar {t})}{\partial \bar {t}} - \frac {\partial ^2 \phi (\bar {x},\bar {t})}{\partial \bar {x}^2} = -\frac {\delta ^2}{2}(\bar {x}^2 - 2\bar {x})\,C_{\textit{in}}^{\prime\prime}(\bar {T}(\bar {t})), \end{equation}
with the transformed initial and boundary conditions
\begin{equation} \phi (\bar {x},0) = - \frac {\delta }{2}(\bar {x}^2 - 2\bar {x})\,C_{\textit{in}}^{\prime}(0),\qquad \phi (0,\bar {t}) = 0,\qquad \frac {\partial \phi }{\partial \bar {x}}(1,\bar {t}) = 0, \end{equation}
where
$C_{\textit{in}}^{\prime}=\mathrm d C_{\textit{in}}/\mathrm d\bar T$
and
$C_{\textit{in}}^{\prime\prime}= d^2 C_{\textit{in}}/\mathrm d\bar T^2$
.
We expand
$\phi =\phi _0+\delta \,\phi _1+O(\delta ^2)$
. Since the source term on the right of (4.4) is
$O(\delta ^2)$
and the initial condition is
$O(\delta )$
, the leading problems are
\begin{equation} \phi _0\equiv 0,\qquad \frac {\partial \phi _1}{\partial \bar t}-\frac {\partial ^2 \phi _1}{\partial \bar x^2}=0,\quad \phi _1(\bar x,0)=-\frac 12\bigl(\bar x^2-2\bar x\bigr)\,C_{\textit{in}}^{\prime}(0), \end{equation}
with
$\phi _1(0,\bar t)=0$
and
$\partial _{\bar x}\phi _1(1,\bar t)=0$
. Using separation of variables with eigenfunctions
$\sin (\lambda _n\bar x)$
,
$\lambda _n=( {(2n+1)\pi })/{2}$
$(n=0,1,\ldots )$
, yields
\begin{equation} \phi _1(\bar x,\bar t)=C_{\textit{in}}^{\prime}(0)\sum _{n=0}^{\infty }\frac {2}{\lambda _n^{3}}\,e^{-\lambda _n^{2}\bar t}\,\sin (\lambda _n \bar x). \end{equation}
Therefore, the solute field valid up to
$O(\delta )$
, is
\begin{equation} c(\bar x,\bar t)\approx c_{\textit{in}}(\bar t)+\frac 12\bigl(\bar x^2-2\bar x\bigr)\,c_{\textit{in}}^{\prime}(\bar t) +c_{\textit{in}}^{\prime}(0)\sum _{n=0}^{\infty }\frac {2}{\lambda _n^{3}}\,e^{-\lambda _n^{2}\bar t}\,\sin (\lambda _n \bar x), \end{equation}
where
$\lambda _n=( {(2n+1)\pi })/{2}$
.
Since we focus on slowly varying solute profiles, the overall process duration is much longer than
$\tau _{{s}}$
, ensuring that
$\bar {t} \gt 1$
for most of the extraction. Under these conditions, the exponentially decaying terms in (4.8) become negligible, effectively eliminating the influence of the initial condition. As a result, the solute field simplifies to
\begin{equation} c(\bar {x}, \bar {t}) \approx c_{\textit{in}}(\bar {t}) + \frac {1}{2}\bigl(\bar {x}^2 - 2\bar {x}\bigr)\,c_{\textit{in}}^{\prime}(\bar {t}). \end{equation}
To retain generality, we allow the (non-dimensional) diffusiophoretic mobility to depend on concentration,
$\bar \varGamma _{ p}=\bar \varGamma _{p}(c)$
, which captures the dependence of the particle
$\zeta$
-potential and the Debye length on the local solute concentration (Prieve et al. Reference Prieve, Anderson, Ebel and Lowell1984; Gupta et al. Reference Gupta, Shim and Stone2020; Lee et al. Reference Lee, Lee and Ault2023). The diffusiophoretic velocity then reads
\begin{equation} \bar U_{\textit{DP}}(\bar x,\bar t) \;\simeq \; \bar \varGamma _{p}\!\left (c_{\textit{in}}(\bar t)\right )\, \frac {(\bar x-1)\,c_{\textit{in}}^{\prime}(\bar t)}{c_{\textit{in}}(\bar t)}\,, \end{equation}
where we have ignored the slowly quadratic correction in (4.9) compared with the leading term for slowly varying inlet profiles. Therefore, the particle transport equation can be written as
\begin{equation} \frac {\partial n}{\partial \bar {t}} + \frac {\partial }{\partial \bar {x}}\left [\bar {\varGamma }_{{p}}(c_{\textit{in}})\frac {(\bar {x}-1)c_{\textit{in}}^{\prime}}{c_{\textit{in}}}n\right ] = \epsilon \frac {\partial ^2 n}{\partial \bar {x}^2}, \end{equation}
where
$\epsilon = \bar {D}_{{p}} = D_{{p}}/D_{{s}} \ll 1$
, and subject to the initial and boundary conditions
\begin{equation} n(\bar {x},0) = n_0,\quad n(0,\bar {t}) = 0,\quad \frac {\partial n}{\partial \bar {x}}(1,\bar {t}) = 0. \end{equation}
We employ a singular-perturbation analysis to solve (4.11). Far from the inlet, where diffusion is negligible, we set
$\epsilon = 0$
and balance only the advective terms. By employing the method of characteristics, the outer solution,
$n_{\mathrm{o}}$
, is obtained
\begin{equation} \frac {n_{\mathrm{o}}(\bar {x},\bar {t})}{n_0}=\frac {n_{\mathrm{o}}(\bar {t})}{n_0} = \exp \!\left (-\int _{c_{\textit{in}}(0)}^{\,c_{\textit{in}}(\bar {t})} \bar {\varGamma }_{p}(c)\,\frac {\mathrm{d}c}{c}\right )\!. \end{equation}
When
$\bar {\varGamma }_{p}$
is constant, this reduces to the power law
\begin{equation} \frac {n_{\mathrm{o}}(\bar {t})}{n_0} = \left (\frac {c_{\textit{in}}(\bar {t})}{c_{\textit{in}}(0)}\right )^{-\bar {\varGamma }_{p}}. \end{equation}
The outer solution
$n_{\mathrm{o}}(\bar {t})$
is independent of
$\bar {x}$
and therefore cannot satisfy the inlet boundary condition
$n(0,\bar {t})=0$
.
To resolve this mismatch, we consider a boundary layer near the inlet, introducing a stretched coordinate
$\bar {y}=\bar {x}/\epsilon$
. In the inner region, (4.11) can be written as
\begin{equation} \bar {\varGamma }_{{p}}\frac {c_{\textit{in}}^{\prime}(\bar {t})}{c_{\textit{in}}(\bar {t})}\frac {\partial n_{{i}}}{\partial \bar {y}} + \frac {\partial ^2 n_{{i}}}{\partial \bar {y}^2} = \epsilon \left [\frac {\partial n_{{i}}}{\partial \bar {t}} + \bar {\varGamma }_{{p}}\frac {c_{\textit{in}}^{\prime}(\bar {t})}{c_{\textit{in}}(\bar {t})}\left (n_{{i}} + \bar {y}\frac {\partial n_{{i}}}{\partial \bar {y}}\right )\right ]\!. \end{equation}
Neglecting the
$\mathcal{O}(\epsilon )$
terms and applying the boundary conditions
$n_{{i}}(0,\bar {t})=0$
and
$\lim _{\bar {y}\to \infty }n_{{i}}(\bar {y},\bar {t})=\lim _{\bar {x}\to 0}n_{\text{o}}(\bar {x},\bar {t})$
, we find
\begin{equation} \frac {n_{{i}}(\bar {y},\bar {t})}{n_0} = \frac {n_{\mathrm{o}}(\bar {t})}{n_0} \left [1 - \exp \left (-\bar {\varGamma }_{{p}}\frac {c_{\textit{in}}^{\prime}(\bar {t})}{c_{\textit{in}}(\bar {t})}\bar {y}\right )\right ]\!. \end{equation}
Transforming back to the original coordinate, we have
\begin{equation} \frac {n_{{i}}(\bar {x},\bar {t})}{n_0} = \frac {n_{\mathrm{o}}(\bar {t})}{n_0} \left [1 - \exp \left (-\mathcal{K}(\bar t)\bar {x}\right )\right ]\!, \end{equation}
where for simplicity we define the diffusion layer thickness,
$\mathcal{K}^{-1}(\bar t) \equiv (({\bar {\varGamma }_{p}}/{\bar {D}_{p}})( { c_{\textit{in}}^{\prime}(\bar t)}/ { c_{\textit{in}}(\bar t)}) )^{-1}$
.
To construct the composite solution that satisfies both the inlet boundary condition and the behaviour away from the inlet, we combine the inner and outer solutions as
\begin{equation} \frac {n(\bar {x},\bar {t})}{n_0} = \underbrace {\frac {n_{\mathrm{o}}(\bar {t})}{n_0}}_{\frac {n_{\text{o}}}{n_0}} + \underbrace {\frac {n_{\mathrm{o}}(\bar {t})}{n_0}\left [1 - \exp \left (-\mathcal{K}(\bar t)\bar {x}\right )\right ]}_{\frac {n_{{i}}}{n_0}} - \underbrace {\frac {n_{\mathrm{o}}(\bar {t})}{n_0}}_{\lim _{\bar {x}\to 0}\frac {n_{\text{o}}}{n_0}}, \end{equation}
which simplifies to
\begin{equation} \frac {n(\bar {x},\bar {t})}{n_0} = \exp \!\left (-\int _{c_{\textit{in}}(0)}^{\,c_{\textit{in}}(\bar {t})} \bar {\varGamma }_{p}(c)\,\frac {\mathrm{d}c}{c}\right ) \left [1 - \exp \left (-\mathcal{K}(\bar t)\bar {x}\right )\right ]\!, \end{equation}
that is just the inner solution as the outer solution is constant.
To ensure the validity of this perturbation solution, we substitute it back into (4.11) and verify that the diffusion term remains negligible compared with the advective term in the outer region. This leads to the thin diffusion layer condition
\begin{equation} \lvert \mathcal{K}^{-1}(\bar t)\rvert = \bigg \lvert \left (\frac {\bar {\varGamma }_{{p}}}{\bar {D}_{{p}}}\frac { c_{\textit{in}}^{\prime}(\bar {t}) }{ c_{\textit{in}}(\bar {t}) }\right )^{-1}\bigg \rvert \ll 1. \end{equation}
Integrating (4.19) over
$\bar {x}$
from 0 to 1 yields the average particle density
\begin{equation} \frac {n_{\textit{ave}}(\bar {t})}{n_0} = \exp \!\left (-\int _{c_{\textit{in}}(0)}^{\,c_{\textit{in}}(\bar {t})} \bar {\varGamma }_{p}(c)\,\frac {\mathrm{d}c}{c}\right ) \left [1 - \frac {1}{\mathcal{K}(\bar t)} \left (1 - \exp \left (-\mathcal{K}(\bar t)\right )\right )\right ]\!. \end{equation}
In the non-diffusive limit,
$D_{{p}}/\varGamma _{{p}} \rightarrow 0$
, and with constant diffusiophoretic mobility, the final residual particle fraction simplifies to
\begin{equation} \frac {N_{{F}}}{N_0} = \lim _{\bar {D}_{{p}} \rightarrow 0}\frac {n_{\textit{ave}}(\infty )}{n_{0}} = \left (\frac {c_{\textit{in}}(\infty )}{c_{\textit{in}}(0)}\right )^{-\bar {\varGamma }_{{p}}} = \beta ^{-\varGamma _{{p}}/D_{{s}}}, \end{equation}
where
$\beta = c_{\textit{in}}(\infty )/c_{\textit{in}}(0)$
. Here,
$N_0$
and
$N_{{F}}$
denote the initial and final total number of particles in the dead-end pore, while
$n_0$
denotes the initial particle number density in the pore. Notably, this expression is independent of the solute-transition time and depends solely on the final solute contrast and the ratio
$\varGamma _{{p}}/D_{{s}}$
. This result is fully consistent with that of Migacz, Castleberry & Ault (Reference Migacz, Castleberry and Ault2024), which focused on non-diffusive particles and linear solute profiles from a Lagrangian point of view. In contrast, we adopt an Eulerian perspective by solving the continuous advection–diffusion equation for the particle-density field, thereby incorporating finite particle diffusivity into the model. This approach captures the diffusion-driven corrections that become important for large
$T/\tau _{{s}}$
. Moreover, our framework is applicable to a wide range of slowly varying inlet solute-concentration profiles, including those of the error-function type (Appendix A.2).
4.1. Linear solute profile under slowly varying conditions with constant
$\bar {\varGamma }_{{p}}$
To apply the above framework to the previously examined linear solute profile, we must first ensure that the inlet solute concentration is slowly varying and the thin diffusion layer condition (4.20) is satisfied, which implies
\begin{equation} \bar {T} \gg 1 \ \text{and}\ \lvert \mathcal{K}_{\textit{linear}}(\bar t) \rvert \gg 1,\ \bar t\in [0,\bar T], \end{equation}
where
$\beta = c_1/c_0$
and
$\mathcal{K}_{\textit{linear}}(\bar t) = ({\bar \varGamma _{{p}}}/{\bar D_{{p}}})( {1}/{\bar t + ({\bar T}/{\beta -1})})$
. Under these constraints, the particle-density field and the average particle density simplify to
\begin{equation} \frac {n(\bar x,\bar t)}{n_{0}} = \left (1+\frac {\beta -1}{\bar T}\bar t\right )^{-\bar \varGamma _{{p}}}\left [1-\exp \left (-\mathcal{K}_{\textit{linear}}(\bar t)\bar x\right )\right ]\!,\ \bar t\in [0,\bar T], \end{equation}
and
\begin{equation} \frac {n_{\textit{ave}}(\bar t)}{n_{0}} = \left (1+\frac {\beta -1}{\bar T}\bar t\right )^{-\bar \varGamma _{{p}}}\left [1-\frac {1}{\mathcal{K}_{\textit{linear}}(\bar t)}\left (1-\exp {\left (-\mathcal{K}_{\textit{linear}}(\bar t)\right )}\right )\right ]\!,\ \bar t\in [0,\bar T]. \end{equation}
The above expressions apply only for
$\bar {t} \lt \bar {T}$
. However, when
$\bar {T} \gg 1$
, diffusiophoretic effects after
$\bar {t} = \bar {T}$
are negligible. As shown in figure 4, for
$\bar {T} \gg 1$
, the magenta dash-dot line, which indicates the residual particle fraction at
$\bar {t} = \bar {T}$
, converges to the black solid line, representing the final residual particle fraction derived under the non-diffusive assumption. Therefore, for
$\bar {T} \gg 1$
, we can approximate the final residual particle fraction by
\begin{equation} \frac {N_{{F}}}{N_{0}} \approx \frac {n_{\textit{ave}}(\bar T)}{n_{0}} = \beta ^{-\bar \varGamma _{{p}}}\left [1-\frac {\bar D_{{p}}}{\bar \varGamma _{{p}}}\frac {\beta \bar T}{\beta -1}\left (1-\exp \left (-\frac {\bar \varGamma _{{p}}}{\bar D_{{p}}}\frac {\beta -1}{\beta \bar T}\right )\right )\right ]\!. \end{equation}
Figure 4(b) compares this prediction (yellow solid line) against the simulation results (diamonds). As expected, when
$\bar T$
falls within
$[10, ( {\varGamma _{{p}}}/{D_{{p}}})\boldsymbol{\cdot }( {\beta -1}/{\beta })]$
, where condition (4.23) remains valid for most of the extraction process, the theoretical prediction closely matches the numerical simulations. The validity of the solute-field expression (4.9) and the particle-field expression (4.19) is also examined (Appendix A.1).
Once the solute gradients vanish (
$\bar {t} \gt \bar {T}$
), particle transport is governed solely by diffusion. For a pure diffusion process starting from a uniform particle distribution
$n_0$
, the evolution of the average particle density is given by
\begin{equation} \frac {n_{\textit{ave}}^{\text{C}}(\bar {t})}{n_0} = \sum _{m=0}^{\infty } \frac {2}{\lambda _m^2}\exp \left (-\lambda _m^2 \bar {D}_{{p}}\bar {t}\right )\!, \end{equation}
where
$\lambda _m = ({(2m+1)\pi })/{2}$
. Under conditions of large solute-transition times, the particle distribution remains nearly uniform; therefore, its concentration at
$\bar {t} = \bar {T}$
can be used as the initial condition for the subsequent pure diffusion. By combining (4.27) with (4.25), we derive an expression valid over the entire process
\begin{equation} \frac {n_{\textit{ave}}(\bar {t})}{n_0} = \begin{cases} \displaystyle \left (1+\frac {\beta -1}{\bar {T}}\bar {t}\right )^{-\bar {\varGamma }_{{p}}}\left [1-\frac {1}{\mathcal{K}_{\textit{linear}}(\bar t)}\left (1-\exp \left (-\mathcal{K}_{\textit{linear}}(\bar t)\right )\right )\right ]\!, & \bar {t}\in [0,\bar {T}], \\[8pt] \displaystyle \frac {n_{\textit{ave}}(\bar {T})}{n_0}\sum _{m=0}^{\infty }\frac {2}{\lambda _m^2}\exp \left (-\lambda _m^2\bar {D}_{{p}}(\bar {t}-\bar {T})\right ), & \bar {t}\gt \bar {T}. \end{cases} \end{equation}
This composite solution characterises the entire particle extraction process, encompassing the initial diffusiophoretic-driven phase (
$\bar {t} \leqslant \bar {T}$
) and the subsequent diffusion-only phase (
$\bar {t} \gt \bar {T}$
). Using this equation, we can also derive the particle escape time probability density function, or equivalently, the breakthrough curve (BTC), defined as
\begin{equation} p(\bar t) = - \frac {1}{n_0}\frac {{\rm d} n_{\textit{ave}}(\bar t)}{{\rm d}\bar t}. \end{equation}
Figures 5(a) and 5(b) respectively show the residual particle fraction and BTCs for linear solute profiles at various
$\bar {T}$
values, along with a control case for comparison. Symbols represent numerical solutions of the advection–diffusion equations, and solid lines are our theoretical predictions. For the cases with
$\bar {T}=10$
and
$\bar {T}=100$
, the BTCs display a distinct transition at
$\bar {t}=\bar {T}$
, reflecting the shift in escape mechanisms from combined diffusiophoretic and diffusive transport to pure diffusion after the solute gradient dissipates.
Figure 5. Analytical and numerical results for particle transport in a dead-end pore with large normalised solute-transition times,
$T/\tau _{{s}}$
. (a, b) Residual particle fraction (a) and breakthrough curves (b) for
$T/\tau _{{s}} = 10^1, 10^2$
and
$10^3$
, compared with a control case without solute gradients. Symbols indicate numerical results, and solid lines show analytical predictions, demonstrating good agreement across a range of solute-transition times.
5. Solute-concentration-dependent diffusiophoretic mobility (
$\varGamma _p(c)$
)
So far, we have assumed the diffusiophoretic mobility to be constant. In practice, the local solute concentration modulates the particle zeta potential
$\zeta$
and the Debye length, and thus the mobility. To test the robustness of our conclusions, we retain the slowly varying, linear inlet profile but allow the non-dimensional mobility
$\bar {\varGamma }_{p}(c)\equiv \varGamma _{p}(c)/D_{\mathrm s}$
to depend on
$c$
. The diffusiophoretic mobility is expressed as (Prieve et al. Reference Prieve, Anderson, Ebel and Lowell1984; Lee et al. Reference Lee, Lee and Ault2023)
\begin{equation} \varGamma _{p} = \frac {\varepsilon }{\mu }\!\left (\frac {k_{B}T}{Z e}\right )^{2} \!\left \{ \frac {Z e\,\zeta \,b}{k_{B}T}\,g(\lambda ) + 4\,\ln \!\left [\cosh \!\left (\frac {Z e\,\zeta }{4 k_{B}T}\right )\right ]\,h(\lambda ) \right \}\!, \end{equation}
with
\begin{align} g(\lambda ) & = \left [ 1 - \frac {\lambda k_{B}T}{2 Z e\,\zeta } \!\left ( F_{1}(\bar {\zeta }) + \frac {\varepsilon }{2\mu D_{ s}}\!\left (\frac {k_{B}T}{Z e}\right )^{2} \big (b F_{4}(\bar {\zeta })+F_{5}(\bar {\zeta })\big ) \right )\! \right ]^{-1}\!, \end{align}
\begin{align} h(\lambda ) & = \left [ 1 - \frac {\lambda }{8\,\ln \!\left (\cosh \!\left (\frac {Z e\,\zeta }{4 k_{B}T}\right )\right )} \!\left ( F_{0}(\bar {\zeta }) + \frac {\varepsilon }{2\mu D_{ s}}\!\left (\frac {k_{B}T}{Z e}\right )^{2} \big (F_{2}(\bar {\zeta })+b F_{3}(\bar {\zeta })\big ) \right )\! \right ]^{-1}. \end{align}
Here,
$\varepsilon$
is the medium permittivity,
$\mu$
the dynamic viscosity,
$k_B$
Boltzmann’s constant,
$T$
temperature,
$Z$
the electrolyte valence,
$e$
the elementary charge and
$\zeta$
the zeta potential;
$\bar {\zeta }=\zeta /\! (k_B T/(Z e) )$
is its non-dimensional form. The coefficient
$b=(D_{+}-D_{-})/(D_{+}+D_{-})$
(for a
$Z{:}Z$
electrolyte) involves the ionic diffusivities
$D_{\pm }$
. The size-dependent functions
$g(\lambda )$
and
$h(\lambda )$
depend on
$\lambda =(\kappa a)^{-1}$
, the ratio of Debye length to particle radius
$a$
, with
$\kappa ^{-1}=\sqrt { {\varepsilon k_B N_A T}/{2 Z^{2} e^{2} c}}$
, where
$N_A$
is Avogadro’s number and
$c$
the molar ion concentration (Prieve et al. Reference Prieve, Anderson, Ebel and Lowell1984; Lee et al. Reference Lee, Lee and Ault2023). Following Lee et al. (Reference Lee, Lee and Ault2023), we evaluate the auxiliary functions
$F_n(\bar {\zeta })$
using fitted forms
$F_n(\bar {\zeta })\approx -A_n \exp (\bar {\zeta }/T_n)$
with
$A_n=\{4.53,\,5.25,\,4.30,\,3.35,\,6.77,\,9.26\}$
and
$T_n=\{1.37,\,1.41,\,1.11,\,1.08,\,1.17,\,1.23\}$
for
$n=0,1,\ldots ,5$
.
To close the model, we relate the
$\zeta$
potential to the local salt concentration. Under a CSC approximation (Lee et al. Reference Lee, Lee and Ault2023),
$\zeta$
varies approximately logarithmically with concentration. For a
$Z{:}Z$
electrolyte we adopt the empirical fit (Kirby & Hasselbrink Jr. Reference Kirby and Hasselbrink2004a
,
Reference Kirby and Hasselbrinkb
)
\begin{equation} \zeta (c) \;=\; a_{0} \;+\; a_{1}\,\log _{10}\!\Big (\frac {Z^{2}\,c}{c_{0}}\Big ), \end{equation}
where
$c_{0}=1\,\mathrm{M}$
is a reference concentration,
$a_{0}$
(mV) is an offset, and
$a_{1}$
(mV/decade) is the slope. For latex particles in NaCl (
$Z=1$
) solutions, representative values are
$a_{0}\approx 0\,\mathrm{mV}$
and
$a_{1}\approx -43.7\,\mathrm{mV/decade}$
(Staffeld & Quinn Reference Staffeld and Quinn1989; Lee et al. Reference Lee, Lee and Ault2023).
Substituting (5.3) together with
$\lambda (c)= [\kappa (c)a ]^{-1}$
(with
$\kappa$
the inverse Debye length) into (5.1)–(5.2) yields
\begin{equation} \bar {\varGamma }_{p}(c) \;=\; \frac {1}{D_{ s}}\, \varGamma _{p}\!\big (\zeta (c),\,\lambda (c)\big ), \end{equation}
which we use below to quantify how local concentration modulates the particle mobility.
For a linear inlet ramp from
$0.1\,\mathrm{}$
to
$10\,\mathrm{mM}$
, we numerically solve (3.4) and (3.5) with the concentration-dependent mobility
$\varGamma _{{p}}(c)$
, and obtain the final residual particle fraction
$N_{{F}}/N_0$
as a function of the non-dimensional solute-transition time
$T/\tau _{{s}}$
(figure 6
a). The trend mirrors the constant-
$\varGamma _{{p}}$
case: as
$T/\tau _{{s}}$
increases from
$0$
to
$\infty$
, the non-diffusive model predicts a decrease in the final residual from
$\alpha _{{H}}=0.33$
to
$\alpha _{{L}}=0.06$
. Including particle diffusion further reduces the residual at larger
$T$
, as indicated by the symbols in figure 6(a).
Figure 6. Our analytical and numerical predictions using a solute-concentration-dependent diffusiophoretic mobility
$\varGamma _{{p}}(c)$
mirrors those obtained earlier using constant diffusiophoretic mobility: longer solute-transition times enhance particle extraction. (a) Final residual fraction
$N_{{F}}/N_0$
versus
$T/\tau _{{s}}$
. Diamonds: diffusive model from the continuum solution of (3.5). Black solid line: non-diffusive particle tracking using (3.4). Yellow line: analytical prediction (5.7) that accounts for both diffusiophoresis and particle diffusion. Inset: predicted dependence of the normalised diffusiophoretic mobility
$\varGamma _{{p}}/D_{{s}}$
on bulk NaCl concentration
$c$
for latex particles by (5.4); dashed magenta lines delimit the concentration range used in our simulations. (b) Temporal evolution of the residual fraction for
$T/\tau _{{s}}=10^{1},10^{2},10^{3}$
. Symbols: numerical results; solid lines: analytical predictions from (5.6). In both panels the agreement is good across a wide range of solute-transition times.
The analytical framework for slowly varying inlet profiles also applies when
$\varGamma _{{p}}$
depends on
$c$
. Following the procedure in § 4.1, but now with
$\varGamma _{{p}}=\varGamma _{{p}}(c)$
, the particle-density field and its cross-sectional average are obtained by combining (5.4), (4.19) and (4.21)
\begin{align} \frac {n(\bar x,\bar t)}{n_{0}} & = \exp \!\left (-\int _{c_{\textit{in}}(0)}^{\,c_{\textit{in}}(\bar {t})} \bar {\varGamma }_{p}(c)\,\frac {\mathrm{d}c}{c}\right ) \left [1-\exp \!\left (-\mathcal{K}_{\textit{linear}}(\bar t)\,\bar x\right )\right ]\!, \quad \bar t\in \bigl[0,\bar T\bigr], \end{align}
\begin{align} \frac {n_{\textit{ave}}(\bar t)}{n_{0}} & = \exp \!\left (-\int _{c_{\textit{in}}(0)}^{\,c_{\textit{in}}(\bar {t})} \bar {\varGamma }_{p}(c)\,\frac {\mathrm{d}c}{c}\right ) \!\left [\!1-\frac {1}{\mathcal{K}_{\textit{linear}}(\bar t)} \left (1-\exp {\left [-\mathcal{K}_{\textit{linear}}(\bar t)\right ]}\right )\!\right ]\!, \!\!\quad \bar t\in \bigl[0,\bar T\bigr], \end{align}
where
$\mathcal{K}_{\textit{linear}}(\bar t)=( {\bar \varGamma _{p}(c_{\textit{in}}(\bar t))}/{\bar D_{p}} )(\bar t+( {\bar T}/{\beta -1} ))^{-1}$
with
$\beta \equiv c_1/c_0$
.
An approximation for the final residual follows by evaluating (5.6) at
$\bar t=\bar T$
\begin{align} \frac {N_{{F}}}{N_{0}} \approx \frac {n_{\textit{ave}}(\bar T)}{n_{0}} & = \exp \!\left (-\int _{c_0}^{\,c_1} \bar {\varGamma }_{p}(c)\,\frac {\mathrm{d}c}{c}\right ) \notag\\&\quad\times\left [1-\frac {\bar D_{{p}}}{\bar \varGamma _{{p}}(c_1)}\, \frac {\beta \bar T}{\beta -1}\left (1-\exp \!\left (-\frac {\bar \varGamma _{{p}}(c_1)}{\bar D_{{p}}}\, \frac {\beta -1}{\beta \bar T}\right )\right )\right ]\!. \end{align}
We compare these predictions with simulations in figure 6. In panel (a), the approximation (5.7) (yellow dashed) is plotted against the numerical results (diamonds). In panel (b), the curves from (5.6) (solid lines) for
$T/\tau _{{s}}\in \{10, 100, 1000\}$
agree closely with the corresponding simulations (symbols).
Within the solute-concentration range we probed, variation in diffusiophoretic mobility did not alter our conclusions about the relative effects of diffuse versus sharp solute fronts. Nevertheless, a systematic analysis of regimes where the diffusiophoretic mobility approaches the particle diffusivity is needed to test the generality of these conclusions over a broader concentration range (Keh & Wei Reference Keh and Wei2000; Prieve et al. Reference Prieve, Malone, Khair, Stout and Kanj2019; Gupta et al. Reference Gupta, Shim and Stone2020; Ault & Shin 2024).
6. Electrolytes versus non-electrolytes
While our model in § 4 applies for either sign of the diffusiophoretic mobility
$\varGamma _{p}$
, our measurements and simulations above focused on
$\beta \gt 1$
with
$\varGamma _{p}\gt 0$
, for which a slowly varying (‘diffuse’) solute concentration at the inlet (larger
$T$
) removes more particles and leaves a smaller residual particle fraction than an abrupt change in solute concentration. This raises a question on whether our observations are limited to positive mobilities and stem from the ‘log-sensing’ effect (Palacci et al. Reference Palacci, Cottin-Bizonne, Ybert and Bocquet2012; Banerjee et al. Reference Banerjee, Williams, Azevedo, Helgeson and Squires2016; Banerjee & Squires Reference Banerjee and Squires2019; Raj, Shields & Gupta Reference Raj, Shields and Gupta2023).
To address this question, here, we first examine the case
$\beta \lt 1$
with constant
$\varGamma _{p}\lt 0$
for non-Brownian particles (negligible
$D_{p}$
) driven by linearly varying solute concentration at the inlet. The two data sets, shown together in figure 7(a), reveal three trends:
-
(i) at fixed
$|\ln \beta |$
and
$|\varGamma _{p}|$
, extraction is stronger for
$\varGamma _{p}\lt 0$
than for
$\varGamma _{p}\gt 0$
; -
(ii) in the
$\varGamma _{p}\lt 0$
regime, increasing the ramp duration
$T$
(i.e. a more diffuse front) slightly weakens the extraction; and -
(iii) for both signs of
$\varGamma _{p}$
, when
$T/\tau _{s}\gg 1$
the residual particle fraction approaches the asymptotic prediction
$\beta ^{-\varGamma _{p}/D_{ s}}$
from (4.22) in § 4.
Figure 7. Comparison of electrolyte (
$\varGamma _{{p}}\gt 0$
and
$\varGamma _{{p}}\lt 0$
) and non-electrolyte (
$\varGamma _{\textit{NE}}$
) cases. (a) Final residual fraction
$N_{{F}}/N_0$
versus the normalised solute-transition time
$T/\tau _{{s}}$
, obtained from particle tracking with non-diffusive particles. Red: electrolyte with
$\varGamma _{{p}}\gt 0$
; blue: electrolyte with
$\varGamma _{{p}}\lt 0$
; black: non-electrolyte model
$U_{\textit{DP}}=\varGamma _{\textit{NE}}\boldsymbol{\nabla }(c/c_0)$
. (b–c) Temporal evolution of the residual fraction
$N/N_0$
for the
$\varGamma _{{p}}\lt 0$
electrolyte (b) and the non-electrolyte (c) cases at
$T/\tau _{{s}}\in \{0,1,10,100\}$
, computed from non-diffusive particle trajectories. Across cases, larger
$T/\tau _{{s}}$
corresponds to more diffuse inlet ramps.
These trends are consistent with the observations of Migacz et al. (Reference Migacz, Castleberry and Ault2024). The temporal evolution for the
$\varGamma _{p}\lt 0$
cases is shown in figure 7(b). Compared with the
$\varGamma _{p}\gt 0$
case (figure 4
a), the residual
$N/N_0$
initially decreases more slowly but then drops rapidly near the end of the ramp. This acceleration follows from the functional form of the diffusiophoretic velocity
$U_{DP}=\varGamma _{ p}\boldsymbol{\nabla }\ln c=\varGamma _{ p}\,\boldsymbol{\nabla }c / c$
: for
$\varGamma _{p}\lt 0$
we have
$\beta \lt 1$
and the ramp drives
$c$
downward, so the factor
$1/c$
grows in time, amplifying
$|\boldsymbol{\nabla }c|/c$
and hence the late-time drift. Finally, for non-electrolytes, where the diffusiophoretic velocity takes the form
$U_{\textit{DP}}=\varGamma _{\textit{NE}}\boldsymbol{\nabla }(c/c_0)$
, we find that a diffuse solute front consistently removes more particles than a sharp step.
6.1. Sign asymmetry from the solute equation written in
$\ln c$
To rationalise these trends, recall that the diffusiophoretic velocity
$U_{\textit{DP}}=\varGamma _{p}\,\boldsymbol{\nabla }\ln (c/c_0)$
is controlled by the gradient of the effective concentration
$e\equiv \ln (c/c_0)$
. Dividing the solute-transport (3.1) by
$c$
gives
\begin{equation} \partial _t e - D_{ s}\,\partial _{xx}e = D_{s}\bigl (\partial _x e\bigr )^2, \end{equation}
with
$e(x,0)=0$
, the inlet ramp
\begin{equation} e(0,t)= \begin{cases} \ln \bigl [1+(\beta -1)t/T\bigr ], & 0\le t\lt T,\\[3pt] \ln \beta , & t\ge T, \end{cases} \end{equation}
and
$\partial _x e(L,t)=0$
.
Because
$|U_{\textit{DP}}|\propto |\partial _x e|$
, the maximum solute-concentration contrast along the pore
$|\Delta e_{\textit{max}}|=|\ln \beta |$
serves as the fixed available ‘driving budget’ for extraction. However, (6.1) is non-conservative: the positive source term
$D_{ s}(\partial _x e)^2$
generates effective concentration in the interior for any
$e(x,t)$
field. Its influence, however, depends on the sign of the target level
$\ln \beta$
.
-
(i) For
$\varGamma _{p}\gt 0$
and
$\beta \gt 1$
, effective concentration
$e$
increases from
$0$
to
$+\ln \beta$
: the source assists equilibration toward the inlet value, shortening the lifetime of gradients. -
(ii) For
$\varGamma _{p}\lt 0$
and
$\beta \lt 1$
, effective concentration
$e$
decreases from
$0$
to
$-\ln \beta$
: the same positive source now opposes the downward relaxation, keeping gradients alive for longer.
In short, the positive source in (6.1) accelerates convergence when the target
$e$
is positive and impedes it when the target is negative. This sign-asymmetric interplay could explain the stronger overall extraction for
$\varGamma _{p}\lt 0$
at fixed
$|\ln \beta |$
.
6.2. Non-electrolyte case
$ (U_{\textit{DP}}=\varGamma _{\textit{NE}}\,\boldsymbol{\nabla }(c/c_0) )$
For non-electrolytes the diffusiophoretic velocity is
$U_{\textit{DP}}=\varGamma _{\textit{NE}}\,\boldsymbol{\nabla }(c/c_0)$
, with
$c_0$
a reference concentration (Anderson, Lowell & Prieve Reference Anderson, Lowell and Prieve1982; Anderson Reference Anderson1989; Velegol et al. Reference Velegol, Garg, Guha, Kar and Kumar2016; Williams et al. Reference Williams, Lee, Apriceno, Sear and Battaglia2020; Shim et al. Reference Shim, Gouveia, Ramm, Valdez, Petry and Stone2024). Unlike the electrolyte case written in terms of
$e=\ln (c/c_0)$
, the solute equation expressed in
$c$
has no quadratic source term analogous to
$D_{{s}}(\partial _x e)^2$
; therefore the sign-asymmetric argument does not apply directly.
For slowly varying inlet histories, an exact counterpart of the analysis in § 4 can be developed. Proceeding as there, the particle density obeys
\begin{equation} \frac {\partial n}{\partial \bar {t}} +\frac {\partial }{\partial \bar {x}} \!\left [\, \bar {\varGamma }_{\textit{NE}} \big (c_{\textit{in}}(\bar t)\big )\,\frac {(\bar x-1)\,c_{\textit{in}}^{\prime}(\bar t)}{c_0}\,n \right ] =\epsilon \,\frac {\partial ^2 n}{\partial \bar {x}^2}, \end{equation}
where
$\bar {\varGamma }_{\textit{NE}}\equiv \varGamma _{\textit{NE}}/D_{{s}}$
and
$\epsilon =\bar D_{{p}}$
. Using the same singular-perturbation (inner/outer) construction yields
\begin{equation} \frac {n(\bar x,\bar t)}{n_0} =\exp \!\left ( -\!\int _{c_{\textit{in}}(0)}^{\,c_{\textit{in}}(\bar t)} \bar {\varGamma }_{\textit{NE}}(c)\,\frac {\mathrm{d}c}{c_0} \right ) \left [1-\exp \!\big (-\mathcal{K}_{\textit{NE}}(\bar t)\,\bar x\big )\right ]\!, \end{equation}
with diffusion-layer parameter
$\mathcal{K}_{\textit{NE}}(\bar t)=( {\bar {\varGamma }_{\textit{NE}} (c_{\textit{in}}(\bar t) )}/{\bar D_{p}})\,( {c_{\textit{in}}^{\prime}(\bar t)}/{c_0})$
. Averaging over the pore gives
\begin{equation} \frac {n_{{\textit{ave};\textit{NE}}}(\bar t)}{n_0} =\exp \!\left ( -\!\int _{c_{\textit{in}}(0)}^{\,c_{\textit{in}}(\bar t)} \bar {\varGamma }_{\textit{NE}}(c)\,\frac {\mathrm{d}c}{c_0} \right ) \left [ 1-\frac {1}{\mathcal{K}_{\textit{NE}}(\bar t)} \Big (1-e^{-\mathcal{K}_{\textit{NE}}(\bar t)}\Big ) \right ]\!. \end{equation}
In the non-diffusive limit
$\bar D_{{p}}\to 0$
and for constant mobility
$\bar {\varGamma }_{\textit{NE}}$
, the final residual fraction is
\begin{equation} \frac {N_{{F;\textit{NE}}}}{N_0} =\lim _{\bar D_{{p}}\to 0}\frac {n_{{\textit{ave};\textit{NE}}}(\infty )}{n_0} =\exp \!\left [-\,\bar {\varGamma }_{\textit{NE}}\, \frac {c_{\textit{in}}(\infty )-c_{\textit{in}}(0)}{c_0}\right ]\!. \end{equation}
We also simulated the non-electrolyte case with constant
$\varGamma _{\textit{NE}}$
, using non-Brownian particles subjected to linear inlet ramps of varying duration
$T$
. For ease of comparison, we set
$c_0=c_{\textit{in}}(0)$
and choose
$\varGamma _{\textit{NE}}=\varGamma _{p}\,( {\ln \beta }/{\beta -1})$
so that, in the limit
$T/\tau _{s}\to \infty$
, the asymptotic residual matches the electrolyte result:
$N_{{F}}/N_0\to \beta ^{-\varGamma _{p}/D_{ s}} =\exp [-\bar \varGamma _{\textit{NE}}(\beta -1)]$
with
$\bar \varGamma _{\textit{NE}}=\varGamma _{\textit{NE}}/D_{ s}$
. The outcomes are shown in figure 7(a). As
$T/\tau _{ s}$
increases,
$N_{{F}}/N_0$
decreases and approaches the same asymptote, mirroring the trend for
$\varGamma _{p}\gt 0$
. The corresponding time evolutions likewise resemble the
$\varGamma _{p}\gt 0$
case (compare figure 7
c with figure 4
a).
7. Discussion
Our work probes a fundamental question on the role of diffusiophoresis in the transport of colloids in porous media: Does dispersion of the solute front weaken the phoretic effects? Combining microfluidic experiments and numerical simulations, we showed that, while the solute front becomes more dispersed as it travels through the medium, the residual particle fraction plateaus to the same value in different fields of view along the medium. This observation suggests that perhaps surprisingly the dispersion of solute front does not weaken the overall diffusiophoretic migration from low permeability pockets of the medium.
To gain further insight, we then probed the phoretic migration of colloids from 1-D dead-end pores, varying the dispersion of the solute front. Our microfluidic experiments showed that dispersed fronts perhaps paradoxically enhance the removal efficiency of the particles compared with sharp fronts. Using numerical simulations and analytical modelling, we then showed that this observation is due to the persistence of solute gradients in the dispersed front case, leading to a more uniform removal of particles. Our analytical model shows that the final residual particle fraction is mainly a function of solute evolution time
$\bar {T} = T/\tau _{{s}}$
.
Applying this insight to the porous medium studied in § 2, where the area-weighted average pore length is approximately
$L \approx 130\,{\unicode{x03BC}} \text{m}$
and
$\tau _{{s}} \approx 12\,\text{s}$
, we find an equivalent
$\bar {T}$
of approximately
$4$
to
$6$
. The dead-pore analysis predicts a final residual particle fraction of
$\approx 0.04$
for this
$\bar {T}$
. Using the measured fraction of the dead-end pore areas in the porous medium (
$\eta = 0.03$
) and the mean particle density after flushing the porous medium (
$\theta _{{F}}^{{A}} \approx 0.001$
in figure 1
d), we can calculate the average final residual particle fraction
$\alpha = \theta _{{F}}^{{A}}/\eta \approx 0.03$
, which agrees well with our estimated value of
$\alpha \approx 0.04$
for a single dead-end pore.
Our work therefore shows that, contrary to the intuitive expectation that sharper solute fronts should be more effective in driving the diffusiophoretic migration of colloids, smoothing the solute-concentration gradient improves the extraction of particles from dead-end pores. Our work generalises the analysis of Migacz et al. (Reference Migacz, Castleberry and Ault2024), who considered non-Brownian particles under linear inlet ramps and obtained closed-form expressions for trajectories and residual fractions. We extend that framework in two directions. First, we incorporate finite particle diffusivity, capturing diffusion-driven corrections that become relevant at large
$T/\tau _{{s}}$
. Second, we present a general formulation that accommodates a broad class of slowly varying inlet solute-concentration profiles, including error-function-type transitions (Appendix A.2). In the limiting case of negligible
$D_{p}$
, our predictions reduce to those of Migacz et al. (Reference Migacz, Castleberry and Ault2024), ensuring consistency between the two approaches.
In the current descriptions of colloid transport in porous media, typically hydrodynamic and non-hydrodynamic Derjaguin–Landau–Verwey–Overbeek (DLVO)-type interactions between colloids and surrounding matrix are considered (Ryan & Elimelech Reference Ryan and Elimelech1996; Kretzschmar et al. Reference Kretzschmar, Borkovec, Grolimund and Elimelech1999; Bradford et al. Reference Bradford, Yates, Bettahar and Simunek2002; Torkzaban, Bradford & Walker Reference Torkzaban, Bradford and Walker2007; Bradford & Torkzaban Reference Bradford and Torkzaban2008; Molnar et al. Reference Molnar, Johnson, Gerhard, Willson and O’Carroll2015; Morales et al. Reference Morales, Dentz, Willmann and Holzner2017; Miele et al. Reference Miele, de Anna and Dentz2019; Molnar et al. Reference Molnar, Pensini, Asad, Mitchell, Nitsche, Pyrak-Nolte, Miño and Krol2019; Bizmark et al. Reference Bizmark, Schneider, Priestley and Datta2020; Perez et al. Reference Perez, Patiño, Soos and Morales2020; Bordoloi et al. Reference Bordoloi, Scheidweiler, Dentz, Bouabdellaoui, Abbarchi and de Anna2022; Fan et al. Reference Fan, Chapman, Khan, Iacoviello, Mikutis, Pini and Striolo2022; Patino, Johnson & Morales Reference Patino, Johnson and Morales2023; Al-Zghoul et al. Reference Al-Zghoul, Volponi, Johnson and Bolster2024; Volponi et al. Reference Volponi, Al-Zghoul, Porta, Bolster and Johnson2025). While the solute concentration is known to influence the range of equilibrium DLVO-type interactions (Liu, Johnson & Elimelech Reference Liu, Johnson and Elimelech1995; Roy & Dzombak Reference Roy and Dzombak1997; Elimelech et al. Reference Elimelech, Nagai, Ko and Ryan2000), the phoretic effects driven by solute gradients have been mostly ignored so far. Our work demonstrates that, while these phoretic velocities can be much weaker than the typical background flow velocities, they remain persistent due to solute dispersion. Our work therefore builds on the growing evidence that solute gradients can have a strong effect on the diffusiophoretic transport of colloids and emulsions in porous media (Doan et al. Reference Doan, Chun, Feng and Shin2021; Sambamoorthy & Chu Reference Sambamoorthy and Chu2023; Jotkar et al. Reference Jotkar, de Anna, Dentz and Cueto-Felgueroso2024a , Reference Jotkar, Ben-Noah, Hidalgo and Dentzb ; Alipour et al. Reference Alipour, Li, Liu and Pahlavan2026; Sambamoorthy & Chu Reference Sambamoorthy and Chu2025). Our findings could therefore have implications for controlling particle transport in a broad range of applications, from microfluidics and drug delivery to contaminant remediation and oil recovery in subsurface flows, where tuning the solute dispersion could be used to optimise particle extraction efficiency.
Declaration of interests
The authors report no conflicts of interest.
Appendix A. Validation of analytical model
A.1. Field expression validation
We evaluated the accuracy of our approximations for both the solute and particle fields. It is found that the approximate solute field becomes valid with a relatively small
$T/\tau _{{s}}$
. Figure 8(a) displays the normalised solute field
$(c - c_{\textit{in}})/c_{\textit{in}}^{\prime}$
for
$T/\tau _{{s}} = 10$
at various times. As predicted by (4.8), the field converges to the analytical form
$( {1}/{2})(2\bar {x} - \bar {x}^2)$
(yellow dashed line) once
$t/\tau _{{s}} \gt 1$
.
Figure 8. Analytical and numerical results for the evolution of solute and particle fields in a dead-end pore with large normalised solute-transition times,
$T/\tau _{{s}}$
. (a) Evolution of the solute-concentration profile for
$T/\tau _{{s}} = 10$
. The normalised solute concentration
$(c - c_{\textit{in}}) / c_{\textit{in}}^{\prime}$
(symbols) converges to the analytical solution
$( {1}/{2})(2\bar x - \bar x^2)$
(yellow dashed line) for
$t/\tau _{{s}} \gt 1$
. (b) Evolution of the particle-density profile for
$T/\tau _{{s}} = 10^3$
. Symbols show numerical solutions of the advection–diffusion equation, while solid black lines denote the analytical predictions (4.24). The red dashed lines represent the non-diffusive particle model, which diverges from numerical results at later times due to diffusion effects.
In contrast, the theoretically predicted particle-density expression requires a much greater
$T/\tau _{{s}}$
to be accurate. Figure 8(b) shows the particle-density field at different times for
$T/\tau _{{s}} = 10^3$
. Symbols represent numerical solutions to the advection–diffusion equation, while the solid black lines correspond to the predictions of (4.24). The good agreement confirms the validity of our analytical model under conditions of large solute-transition times. During the particle extraction process, the particle-density field remains nearly uniform, a hallmark of scenarios with slowly varying inlet solute profiles, except near the inlet, where a pronounced drop is induced by particle diffusion. For comparison, the red dashed lines show the non-diffusive particle model predictions, which deviate from the numerical results at later times due to the omission of diffusion effects.
A.2. Error-function-type slowly varying solute profile with constant
$\bar {\varGamma }_{{p}}$
Having verified our framework for linear solute profiles, we now consider more realistic inlet conditions. In many porous media, dispersion often yields inlet solute-concentration profiles that closely resemble an error function. Such profiles arise as solutions of the 1-D advection–diffusion equation. Here, we examine particle extraction under an error-function-type inlet solute profile
\begin{equation} c_{\textit{in}}(\bar t) = c_{0} + \frac {1}{2}(c_{1} - c_{0})\left (1 + {erf}\left (\frac {2\bar t}{\bar \tau _{\textit{erf}}} - \delta \right )\right )\!, \end{equation}
where
$\bar \tau _{\textit{erf}} = \tau _{\textit{erf}}/\tau _{{s}}$
sets the time scale of solute-concentration change, and
$\delta =3$
is a shift chosen so that
$c_{\textit{in}}(0)$
remains close to
$c_{0}$
.
To ensure that the inlet solute concentration is slowly varying and the thin diffusion layer condition (4.20) is satisfied, enabling the application of our analytical model, we require
$\bar \tau _{\textit{erf}} \gg 1$
and
\begin{equation} \frac {\bar \varGamma _{{p}}}{\bar D_{{p}}} \gg \frac {c_{\textit{in}}(\bar {t})}{c_{\textit{in}}^{\prime}(\bar {t})} = \frac {\bar \tau _{\textit{erf}}}{\beta - 1}\exp \left [\left (\frac {2\bar t}{\bar \tau _{\textit{erf}}} - \delta \right )^2\right ]\left (1 + \frac {1}{2}(\beta - 1)\left [1 + {erf}\left (\frac {2\bar t}{\bar \tau _{\textit{erf}}} - \delta \right )\right ]\right )\!, \end{equation}
where
$\beta = c_1/c_0$
. Although appropriate choices of
$\bar \tau _{\textit{erf}}$
and
$\beta$
can satisfy the second inequality, the first inequality may be violated at early and late times, when the inlet solute concentration is nearly constant. In these regimes, we approximate transport as pure diffusion (see (4.27)).
As a result, the extraction process is divided into three stages: pure diffusion (4.27) at early times, combined diffusiophoresis and diffusion (4.21) at intermediate times and pure diffusion (4.27) again once the solute gradient has effectively vanished. We determine the two transition points,
$\bar t_{{c1}}$
and
$\bar t_{{c2}}$
, between these stages using the criterion
\begin{equation} \frac {c_{\textit{in}}(\bar t_{\text{c}i})}{c_{\textit{in}}^{\prime}(\bar t_{\text{c}i})} = \frac {r_i}{\bar \tau _{\textit{erf}}} \boldsymbol{\cdot }\frac {\bar \varGamma _{{p}}}{\bar D_{{p}}}, \end{equation}
where
$r_i \ll 1$
(
$i=1,2$
) are fitting parameters. Here, we choose
$r_1 = 6\times 10^{-4}$
and
$r_2 = 8\times 10^{-4}$
.
Under these conditions, the average particle density evolves as
\begin{equation} \frac {n_{\textit{ave}}(\bar t)}{n_0} = \begin{cases} \displaystyle \sum _{m=0}^{\infty }\frac {2}{\lambda _m^2} \exp (-\lambda _m^2 \bar D_{{p}} \bar t), & \bar t \in [0,\bar t_{{c1}}], \\[10pt] \displaystyle \frac {n_{\textit{ave}}(\bar t_{{c1}})}{n_0}\left (\frac {c_{\textit{in}}(\bar t)}{c_{\textit{in}}(\bar t_{{c1}})}\right )^{-\bar \varGamma _{{p}}} \frac {\left [1 - \left (1 - 1/K(\bar t)\exp \bigl (-K(\bar t)\bigr )\right )\right ]}{\left [1 - \left (1 - 1/K(\bar t_{{c1}})\exp \bigl (-K(\bar t_{{c1}})\bigr )\right )\right ]}, & \bar t \in (\bar t_{{c1}}, \bar t_{{c2}}], \\[10pt] \displaystyle \frac {n_{\textit{ave}}(\bar t_{{c2}})}{n_0}\sum _{m=0}^{\infty }\frac {2}{\lambda _m^2}\exp [-\lambda _m^2 \bar D_{{p}} (\bar t - \bar t_{{c2}})], & \bar t \in (\bar t_{{c2}}, \infty ), \end{cases} \end{equation}
where
$\lambda _m = ({(2m+1)\pi })/{2}$
and
$K(\bar t) = ({\bar \varGamma _{{p}}}/{\bar D_{{p}}})({c_{\textit{in}}^{\prime}(\bar {t})}/{c_{\textit{in}}(\bar {t})})$
.
Figure 9. Validation of our analytical framework for an error-function-type slowly varying solute profile. (a) Final residual particle fraction, defined as the value when the average solute concentration in the dead-end pore reaches
$99\,\%$
, as a function of the effective solute-transition time. Results are obtained from continuous solutions for diffusive particles (diamond symbols) and from particle tracking under the non-diffusive assumption (black solid lines). The yellow solid line represents our analytical prediction that accounts for both diffusiophoresis and particle diffusion. (b) Breakthrough curves for
$\tau _{\textit{erf}}/\tau _{{s}} = 10^1, 10^2$
and
$10^3$
. Symbols represent numerical results, while solid lines show our analytical predictions, demonstrating good agreement over a wide range of solute-transition times.
The piecewise definition reflects the three distinct stages of the extraction process: early time pure diffusion, intermediate time coupled diffusiophoresis and diffusion and late time pure diffusion.
Figure 9(a) shows the final residual particle fraction (evaluated when the average solute concentration in the dead-end pore reaches 99 % of its final value) as a function of
$\tau _{\textit{erf}}/\tau _{{s}}$
. The black solid line corresponds to particle tracking results under the non-diffusive assumption, while the diamond symbols represent solutions of the advection–diffusion equation that include particle diffusion. Similar to the linear case, we identify a highest residual fraction
$\alpha _{{H}}=0.173$
and a lowest residual fraction
$\alpha _{{L}}=0.032$
. As
$\tau _{\textit{erf}}/\tau _{{s}}$
increases, the diffusive model diverges from the non-diffusive one, and our analytical prediction successfully captures the further reduction in the residual fraction introduced by particle diffusion.
Figure 9(b) shows BTCs for error-function-type solute profiles at various
$\tau _{\textit{erf}}/\tau _{{s}}$
values. Symbols represent numerical solutions of the advection–diffusion equations, and solid lines represent our theoretical predictions. The good agreement confirms the reliability of our three stage model for the error-function-type slowly varying solute profile.
A.3. Influence of salt type: LiCl versus NaCl
We repeated the dead-end pore experiments described in § 3 using NaCl as the solute. Because colloids in NaCl can adsorb onto channel walls at high salt concentrations, we used lower concentrations:
$c_0=0.1\,\textrm {mM}$
for the initial suspension and
$c_1=10\,\textrm {mM}$
for the introduced colloid-free solution. For direct comparison, we repeated the experiments with LiCl at the same concentrations (
$c_0=0.1\,\textrm {mM}$
,
$c_1=10\,\textrm {mM}$
).
We observe that, for both electrolytes, a diffuse solute front extracts more particles from the dead-end pore than a sharp front, as demonstrated in figure 10, where symbols denote experimental measurements, and solid lines are fits from the non-diffusive particle-trajectory model (3.4) with a constant diffusiophoretic mobility. The fitted mobilities are
$\varGamma _{{p}}^{\textit{NaCl}}\approx 0.4\times 10^{-9}\,\mathrm{m^2\,s^{-1}}$
and
$\varGamma _{{p}}^{\textit{LiCl}}\approx 0.7\times 10^{-9}\,\mathrm{m^2\,s^{-1}}$
, with NaCl exhibiting the lower value.
Figure 10. Dead-end pore experiments with (a) NaCl and (b) LiCl with
$c_0=0.1\,\mathrm{mM}$
and
$c_1=10\,\mathrm{mM}$
. Symbols show experimental measurements; solid lines are fits from the non-diffusive particle-trajectory model in (3.4) using a constant diffusiophoretic mobility. The fitted mobilities are
$\varGamma _{{p}}^{\textit{NaCl}}\approx 0.4\times 10^{-9}\,\mathrm{m^2\,s^{-1}}$
and
$\varGamma _{{p}}^{\textit{LiCl}}\approx 0.7\times 10^{-9}\,\mathrm{m^2\,s^{-1}}$
.
Figure 11. Comparison of simulations using the CSC mobility model versus a constant
$\varGamma _{{p}}$
. (a) Time evolution of the residual fraction
$N/N_0$
. Symbols: experiments; solid lines: 2-D simulations with constant
$\varGamma _{{p}}$
; dashed lines: 2-D simulations with the CSC model. (b) Longitudinal particle -density profiles along the dead-end pore for the sharp front, evaluated at
$t/\tau _{s}=0,\,0.1,\,0.2,\,0.5,\,1,\,2$
. (c) Longitudinal particle-density profiles for the diffuse front, evaluated at
$t/\tau _{s}=0,\,0.25,\,0.5,\,1,\,2,\,4$
. The CSC model captures the qualitative trends and matches the sharp-front case closely, while slightly underestimating the diffuse-front residuals relative to experiment.
A.4. Comparison between CSC and constant–
$\varGamma _{{p}}$
models
We compare 2-D dead-end pore with channel simulations (see § 3) using two expressions for the diffusiophoretic mobility: (i) a constant
$\varGamma _{{p}}$
, already described in § 3, and (ii) a concentration-dependent mobility given by the CSC model (§ 5). In the CSC runs, the slope
$a_1$
in the empirical
$\zeta (c)$
relation is treated as a single fitting parameter: it is calibrated to match the sharp-front experiment, yielding
$a_{1}\approx -42.0\,\mathrm{mV\,decade^{-1}}$
. The diffuse-front simulation then uses the same
$a_1$
and the inlet history inferred from the fluorescein signal with time rescaled as in § 3.
Figure 11(a) overlays the CSC results on the experimental data and the constant
$\varGamma _{{p}}$
simulations from figure 2(e). The CSC model reproduces the qualitative trends in both sharp and diffuse cases. Quantitatively, it closely matches the sharp-front residuals, while for the diffuse front it slightly underestimates the residual fraction relative to the measurements. Figure 11(b–c) shows the corresponding longitudinal density profiles along the pore. Compared with the experimental profiles in figure 2(f–g) and the constant
$\varGamma _{{p}}$
simulations in figure 2(i–j), the CSC model exhibits stronger inlet focusing in the sharp case, whereas for the diffuse case it likewise yields a nearly uniform profile at late times. Overall, for our dead-end pore experiments with LiCl concentration varying from
$c_0=1\,\textrm {mM}$
to
$c_1=100\,\textrm {mM}$
, both constant and variable mobility models capture the key qualitative behaviours observed in the experiments. It would certainly be interesting to explore further the impact of variable zeta potential and phoretic mobility over a wider range of salt concentrations and salt types as recent works offer evidence for their potential role (Akdeniz et al. Reference Akdeniz, Wood and Lammertink2023; Lee et al. Reference Lee, Lee and Ault2023).
Open access
