2.1. Microfluidic experiments

The microfluidic chip is fabricated with standard soft photolithography (Xia & Whitesides Reference Xia and Whitesides1998) using polydimethylsiloxane (PDMS) and plasma-bonded to a glass slide (channel height $\approx 50\, \rm \unicode{x03BC} m$ ). For each experiment, the porous-medium chip (200 × 60 posts or 42.6 mm × 10.7 mm) is filled with an aqueous solution of 0.01 mM LiCl or sodium dodecylsulphate (SDS) using a syringe pump (Harvard Apparatus). Blob creation is done by piercing and injecting the microfluidic chip in situ with particles suspended in 1 mM LiCl, or 0.01 mM LiCl, or 1 mM SDS for the attractive/control/repulsive cases, respectively (figure 1 a). Flow is begun ${\sim} 10\,\rm s$ after the introduction of the blob. Imaging is done using an inverted microscope (Nikon Eclipse Ti2 equipped with an ORCA-Fusion Digital CMOS camera with a spatial resolution of 2300 × 2300 pixels and dynamic range of 16 bit) and 200 nm fluorescent particles (Invitrogen), with excitation/emission wavelengths of 470/525 nm. We kept the particle concentration $n$ within the linear range of calibration correlating with pixel intensity. To image the entire length of the channel, we use large image stitching across six images at 0.2 frames per second.

Figure 1.

Solute gradients significantly modulate the macroscopic dispersion of colloidal blobs in porous media. (a) Set-up schematic showing the geometry of the porous medium consisting of a hexagonal lattice of circular posts with disorder strength $\chi =0.2$ and lattice spacing of $\lambda =213\,\rm \unicode{x03BC} m$ . (b–d) Two snapshots of the colloid blob at $\tilde {t}=0$ and $\tilde {t}=96$ , where $\tilde {t}=tU_m/\lambda$ . The colloidal blob in the attractive case splits into two blobs, while in the repulsive case the blob’s dispersion is suppressed. Colour scale is truncated for visual clarity. (e) Longitudinal dispersion $\tilde {\sigma }_{||}^{2}$ of the colloids as a function of time for three different cases: control, attractive and repulsive. Inset: the centre of mass of all blobs $\tilde {\bar {x}}$ moves with the mean flow velocity. (f) The transverse-averaged light intensity corresponding to panels (b–d).

The 2-D porous medium is designed using circular posts of diameter $158 \,\rm \unicode{x03BC} m$ arranged in a hexagonal lattice with porosity $\phi{\sim} 0.5$ . To create disorder, each obstacle is displaced from its lattice position in a random direction $0\lt \varTheta \lt 2 \pi$ by a perturbation drawn from a uniform random distribution between $(0,\chi r]$ , where $\chi$ is the disorder parameter and $r$ is the post radius. Any overlaps were removed by re-displacing the post with a new perturbation until there were no overlaps. The solute ambipolar diffusivities ( $2D_+D_-/(D_++D_-)$ , where $D_+$ and $D_-$ are cation and anion diffusivities) are $D_{LiCl}=1.37\times 10^{-9} \,\textrm {m}^2\,\textrm {s}^{-1}$ and $D_{{SDS}}=0.84\times 10^{-9} \textrm {m}^2\,\textrm {s}^{-1}$ and the particle diffusivity is $D_p=2\times 10^{-12}\,\textrm {m}^2\,\textrm {s}^{-1}$ . Since the Debye layer ( $\kappa ^{-1}\approx 14 \,\textrm {nm}$ at 298 K for 1 mM Z:Z = 1:1 electrolyte) is comparable to the particle radius, we use the finite Debye layer expression to calculate diffusiophoretic mobility (Lee, Lee & Ault Reference Lee, Lee and Ault2023): $\varGamma _p = ({\epsilon }/{\eta }) ( {k_BT}/{Ze} )^2 [( {Ze\zeta _p\beta }/{k_BT})g(\lambda ') + 4\ln [\cosh { ( {Ze\zeta _p}/{k_BT} )} ]h(\lambda ') ]$ , where $\epsilon$ is the permittivity of the medium, $\eta$ is the dynamic viscosity of the medium, $k_B$ is the Boltzmann constant, $T$ is the temperature, $Z$ is the valence of the solute, $e$ is the elementary charge, $\zeta _p$ is the zeta potential of the colloidal particle, $\beta$ is the diffusivity difference equal to $(D_+-D_-)/(D_++D_-)$ for a Z:Z electrolyte and $\lambda '=(\kappa a)^{-1}$ is the ratio of the Debye layer to the radius of the particle; $g(\lambda ')$ and $h(\lambda ')$ are functions given in Lee et al. (Reference Lee, Lee and Ault2023). Although $\zeta _p$ may depend on the solute concentration, pH and the species of counter-ion, we assume a constant $\zeta _p\approx -70\,\rm mV$ and therefore a constant $\varGamma _p$ in the simulations for simplicity (Gupta, Shim & Stone Reference Gupta, Shim and Stone2020). With this assumption, the corresponding constant diffusiophoretic mobilities are $\varGamma _{p, {LiCl}}=0.27\times 10^{-9}\, \textrm {m}^2\,\textrm {s}^{-1}$ and $\varGamma _{p,{SDS}}=-0.11 \times 10^{-9}\, \textrm {m}^2\,\textrm {s}^{-1}$ . For the control case, we set $\varGamma _p=0$ .

2.2. Experimental observations

In the absence of solute gradients, the colloidal blob is advected by the background flow and becomes dispersed due to velocity heterogeneities (figure 1 b). In the presence of solute gradients, colloids move up the gradient in the LiCl case, where $\varGamma _p\gt 0$ , and down the gradient in the SDS case, where $\varGamma _p\lt 0$ . The mean longitudinal flow velocity is $U_m$ $\approx 190\,\rm \unicode{x03BC} m\,s^{-1}$ (equivalent to the mean migration velocity of the blob), corresponding to $\textit{Pe}_c = U_m\lambda /D_c \approx 1.8\times 10^4$ and $\textit{Pe}_s= U_m\lambda /D_{\textrm {s}} \approx 30 \,(\text{LiCl, attractive}), 48\,(\text{SDS, repulsive})$ . In the attractive case, the colloidal blob suspended in $1$ mM LiCl solution is injected into the medium with background salt concentration $0.01$ mM LiCl. We expect the colloids to remain more coherent than in the control case as they are attracted towards the solute blob, suppressing the dispersion. Surprisingly, however, we observe the opposite: the colloidal blob becomes more dispersed in the longitudinal direction and even shows two peaks, indicating the splitting of the blob (figure 1 c). By contrast, in the repulsive case, where the colloidal blob suspended in $1$ mM SDS solution is injected into the medium filled with $0.01$ mM SDS, we expect the colloidal blob to become more dispersed, and yet we observe a suppressed dispersion (figure 1 d). We characterise the longitudinal dispersion of colloids as $\sigma ^2_{||} = {\iint (x-\bar {x})^2 n \,{\rm d}x {\rm d}y}/{\iint n \,{\rm d}x {\rm d}y}$ , where $\bar {x} = {\iint xn \, {\rm d}x {\rm d}y}/{\iint n \, {\rm d}x {\rm d}y}$ is the centre of mass of the blob, which demonstrates an enhanced/suppressed dispersion in the attractive/repulsive cases, respectively (figure 1 e). We non-dimensionalise these length scales using the lattice spacing $\lambda$ . The transverse-averaged light intensity further demonstrates the splitting in the attractive case and suppressed dispersion in the repulsive case (figure 1 f). To probe the mechanism behind these observations we resort to numerical simulations.

Figure 2.

Two-dimensional numerical simulations of the evolution of solute and colloid blobs. (a,b) A Gaussian blob of solute and colloid is introduced into the medium at $\tilde {t}=0$ . Flow direction is from left to right with $\textit{Pe}_c\approx 7200$ and $\chi =0.2$ . Colour scale is truncated for visual clarity. (c) The flow velocity distribution in the medium. (d–g) Snapshots of the solute, and control, attractive and repulsive cases of colloid fields at $\tilde {t}=tU_m/\lambda =97$ , respectively. (h) The corresponding cross-sectionally averaged concentration profiles show the bimodal splitting in the attractive case and inhibited dispersion in the repulsive case. (i, j) The corresponding longitudinal and transverse dispersion of the blob. Error envelopes depict averaging over different realisations of disorder. (k) The longitudinal dispersion coefficient.

2.3. Two-dimensional numerical simulations

The fluid flow is governed by the Stokes equations:

(2.1) \begin{align} -\boldsymbol{\nabla }\!p + \mu {\nabla} ^2 \boldsymbol {u} & = 0 , \\[-12pt]\nonumber \end{align}

(2.2) \begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol {u} & = 0, \end{align}

where $\boldsymbol{u}\equiv (u,v)$ , $p$ and $\mu$ are the fluid velocity, pressure and viscosity, respectively. We introduce a joint 2-D Gaussian blob of solute ( $c$ ) and colloids ( $n$ ) with characteristic length $l_b\ (\gt \lambda )$ , i.e. $c(x,y,t=0)=0.99\exp ({-(x^2+y^2)/2l_b^2})+0.01$ and $n(x,y,t=0){}=\exp ({-(x^2+y^2)/2l_b^2})$ , respectively (figure 2 a,b). Both solute and colloids are then advected by the background flow $\boldsymbol{u}(x,y)$ (figure 2 c). The evolutions of solute and colloid fields are governed by the advection–diffusion equations:

(2.3) \begin{align} \partial _t c + \boldsymbol{\nabla }\boldsymbol{\cdot }(c\boldsymbol {u}) &= D_{\textrm {s}} {\nabla} ^2 c, \\[-12pt]\nonumber \end{align}

(2.4) \begin{align} \partial _t n + \boldsymbol{\nabla }\boldsymbol{\cdot }\left [ n(\boldsymbol {u} + \boldsymbol {u}_{\textit{dp}}) \right ] &= D_c {\nabla} ^2 n, \end{align}

where the diffusiophoretic velocity $\boldsymbol{u}_{{{\rm d}p}} = \varGamma _p \boldsymbol{\nabla }\ln c$ . To lower the computational costs, values ofPe in the simulations are a factor of 2.5 times smaller than those in experiments, i.e. $\textit{Pe}_c \approx 7200$ , $\textit{Pe}_s \approx 12\,(\text{attractive}), 19\, (\text{repulsive})$ . The 2-D numerical simulations were performed using the finite volume method in OpenFOAM (Ault, Shin & Stone Reference Ault, Shin and Stone2018; Alipour et al. Reference Alipour, Li, Liu and Pahlavan2026). A triangular mesh of resolution ${\sim} 2.5\,\rm \unicode{x03BC} m$ is used. Post boundaries and walls are treated with no-slip boundary conditions. The diffusioosmotic mobility along PDMS surfaces is lower than the diffusiophoretic mobility of the colloidal particles. For glass surfaces, however, the two mobilities have been shown to be comparable in magnitude (Liu & Pahlavan Reference Liu and Pahlavan2025). Here, the effect of diffusioosmosis on the velocity profile is neglected in the simulations for simplicity and computational cost. The reasonable qualitative agreement between our 2-D simulations and experiments (as will be shown) suggests that diffusioosmotic flows may influence the observations quantitatively, but not qualitatively.

While the solute blob has a Gaussian profile at long times ( $t \approx 100 (\lambda / U_m) \approx 10 \lambda ^2/D_s$ ) (figure 2 d), the colloid blob shows a longitudinal bimodal splitting in the attractive case (figure 2 f) and inhibited dispersion in the repulsive case (figure 2 g). The cross-sectionally averaged colloid concentration profiles of the simulations (figure 2 h) and their longitudinal dispersion (figure 2 i) are also in agreement with those observed in the experiments. The transverse dispersion of the blob ( $\sigma^2_{\perp}$ ) is slightly enhanced in the repulsive case and slightly suppressed in the attractive case (figure 2 j), but marginal overall in comparison with the longitudinal dispersion. The longitudinal dispersion coefficient $D_{||} = {\rm d}\sigma^2_{||}/{\rm d}t$ evolves at early times as the blob samples the velocity heterogeneities of the medium and solute gradients, before reaching a plateau at longer times, showing that the attractive case disperses more than the control, while the repulsive case disperses less (figure 2 k). The increase of longitudinal dispersion in the attractive case is more pronounced in the simulations than in the experiments. We attribute this to: (i) the 2.5-D nature of the experiments, i.e. the presence of top/bottom walls with no-slip boundary conditions, and (ii) the asymmetry of the attractive and repulsive cases. In the repulsive case, colloids are pushed away from the low-velocity zones. In these zones, phoresis dominates over advection, and therefore has a strong footprint. In contrast, in the attractive case, colloids are being pulled from the high-velocity zones towards the low-velocity zones. As the colloids are in the high-velocity zones, their transverse motion due to diffusiophoresis will be much smaller than their longitudinal motion due to flow. Therefore, it will take a longer time until the difference between attractive and control case manifests itself.

Figure 3.

Phoretic migration between slow and fast flow zones leads to macroscopic changes in the longitudinal dispersion of the blob. (a) The magnitude of the advective velocity $|\boldsymbol{u}|$ depicting fast zones in the channels and slow zones between channels. (b) As the blob is displaced, higher-solute-concentration pockets emerge in the low-/high-velocity zones in the back/front of the blob. (c) In the attractive case, the transverse component of the phoretic velocity ( $v_{\textit{dp}}$ ) points towards the low-/high-velocity zones in the back/front of the blob, leading to the splitting of the blob in (d). (e) In the repulsive case, the transverse component of the phoretic velocity points towards the high-/low-velocity zones in the back/front of the blob, leading to a less dispersed blob in the longitudinal direction (f).

To gain a better understanding of the mechanism behind our observations, we probe the early-time evolution of the blob. The origin of small-scale transverse solute gradients is in the presence of shear, i.e. velocity gradient between slow and fast streamlines (figure 3 a). In the attractive case, near the front of the solute blob concentration gradients point towards the high-velocity pathways while pointing towards the low-velocity zones near the back of the blob (figure 3 b). These gradients lead to the focusing of the colloids in the high-velocity zones near the front, and in the low-velocity zones near the back of the blob, leading to enhanced dispersion and splitting of the colloidal blob (figure 3 c,d). The phoretic migration is reversed in the repulsive case, leading to an enhanced transverse dispersion near the front, and removing the colloids from the low-velocity zones near the back, and a longitudinally less dispersed colloidal (figure 3 e,f). Effectively, in the attractive case, diffusiophoresis acts against diffusion, lowering the transverse dispersion, while in the repulsive case, diffusiophoresis works with diffusion, enhancing the transverse dispersion. Phoretic transport that occurs at the scale of the pores in a porous medium is similar to the channel-wide phoretic transport between fast and slow streamlines in a channel, as studied in the configuration of a one-dimensional (1-D) Gaussian pulse dispersing in Poiseuille flow by Migacz & Ault (Reference Migacz and Ault2022). We may relate the transverse solute gradients with the longitudinal gradient as: ${\rm d}c/{\rm d}y \sim \textit{Pe}_s ({\rm d}c/{\rm d}x)$ ; therefore, as $\textit{Pe}_s$ decreases, the transverse solute gradients and the corresponding diffusiophoretic migration weaken. This explains the quantitative difference between experiments and the simulations at artificially small $\textit{Pe}_s$ (figures 1 e and 2 i). Further, the solute blob remains Gaussian with its peak concentration decaying as ${\sim} 1/(D_{||,s}t/l_b^2)^{1/2}$ , where $D_{||,s}$ is the longitudinal dispersion coefficient of the solute blob; therefore, the solute gradients decay slowly as ${\sim} t^{-1/2}$ . Even though at long times $t\gg l_b^2/D_{||,s}$ we expect diffusiophoretic drift to become negligible, before this time, significant redistribution of the colloid has occurred already. Since $D_c \ll D_s$ , the particles do not diffuse back into their ‘control case configuration’ even after solute gradients vanish, and therefore the disparity persists at long times.

Figure 4.

(a) We construct a two-layer model with $u_1\gt u_2$ to probe the evolution of a 1-D blob of solute and colloids. The two layers can communicate via diffusion and diffusiophoresis. (b) Full numerical model ((2.5) and (2.6)) depicting the evolution of colloid density field $n=(n_1+n_2)/2$ for $\textit{Pe}_c\sim 7200$ . The blob exhibits bimodal splitting in the attractive case and inhibited dispersion in the repulsive case. (c) The time evolution of second moment of the colloid blob ( $\tilde {\sigma }^2_{||} = \sigma ^2_{||}/l^2$ ) for all three cases and the corresponding dispersion coefficient (inset) show reasonable agreement with 2-D numerical simulations and experiments. (d,e) Corresponding colloid density field, second moment and dispersion coefficient (inset) obtained from the simplified model ((2.5) and (2.7)).

2.4. A minimal two-layer model

Using the insight provided by 2-D simulations, we develop a two-layer model with minimal ingredients to capture the main observations. While simplified models for dispersion in disordered media have been developed over the years (Bouchaud & Georges Reference Bouchaud and Georges1990; Guyon, Nadal & Pomeau Reference Guyon, Nadal and Pomeau2012), here we use a two-layer model, similar in spirit to generalised Taylor dispersion (Van den Broeck Reference Van den Broeck1990), to provide a mechanistic comparison with diffusiophoretic dispersion in Poiseuille flow. We consider two parallel layers with velocities $u_1$ and $u_2\lt u_1$ (figure 4 a), and solute and colloid fields $c_i$ , $n_i$ , where $i=1,2$ . The characteristic length scale between the two layers is $l$ , and is chosen to be the pore width in the porous medium. We introduce a 1-D Gaussian blob of solute and colloids, and observe how they evolve due to the background 1-D flow, longitudinal and transverse diffusion and diffusiophoresis. In the frame moving with mean flow velocity $\bar {u}$ , the governing equation for the solute field is

(2.5) \begin{align} \partial _t \underbrace {\begin{bmatrix} c_1 \\ c_2 \end{bmatrix}}_{\boldsymbol {c}} + \underbrace {\begin{bmatrix} u_1 -\bar {u} & 0 \\ 0 & u_2 -\bar {u} \end{bmatrix}}_{\boldsymbol {U}} \partial _x \begin{bmatrix} c_1 \\ c_2 \end{bmatrix} - D_{\textrm {s}} \partial _x^2 \begin{bmatrix} c_1 \\ c_2 \end{bmatrix} + & \underbrace {\frac {D_{\perp ,s}}{l^2}\overbrace {\begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}}^{\mathbb{A}}}_{\boldsymbol {D}_{\perp ,s}} \begin{bmatrix} c_1 \\ c_2 \end{bmatrix} = 0, \\[-12pt]\nonumber \end{align}

(2.6) \begin{align} \partial _t \underbrace {\begin{bmatrix} n_1 \\ n_2 \end{bmatrix}}_{\boldsymbol {n}} + \begin{bmatrix} u_1-\bar {u} & 0 \\ 0 & u_2-\bar {u} \end{bmatrix} \partial _x \begin{bmatrix} n_1 \\ n_2 \end{bmatrix} - D_c \partial _x^2 \begin{bmatrix} n_1 \\ n_2 \end{bmatrix} + & \underbrace {\frac {D_{\perp ,c}}{l^2}\begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}}_{\boldsymbol {D}_{\perp ,c}} \begin{bmatrix} n_1 \\ n_2 \end{bmatrix} \notag \\ - \underbrace {\frac {v_{\textit{dp}}}{l} \begin{bmatrix} H(-v_{\textit{dp}}) & H(v_{\textit{dp}}) \\[3pt] -H(-v_{\textit{dp}}) & -H(v_{\textit{dp}}) \end{bmatrix}}_{\boldsymbol {D}_{\perp ,{dp}}} \begin{bmatrix} n_1 \\ n_2 \end{bmatrix} = 0, \end{align}

where the tensors $\boldsymbol {D}_{\perp ,s}$ , $\boldsymbol {D}_{\perp ,c}$ and $\boldsymbol {D}_{\perp ,{dp}}$ represent the transverse dispersion of solute and colloids, and transverse phoretic exchange, respectively. Function $H(x)$ is the step function. Inter-layer transport by diffusiophoretic velocity is given by $v_{\textit{dp}}\partial _yn_i \sim v_{\textit{dp}}n_i/l$ , where $v_{\textit{dp}} = \varGamma _p( {\partial \log {\text{(solute conc.)}}}/{\partial y} )\approx \varGamma _p( {(\log {c_1}-\log {c_2})}/{l}) = ( {\varGamma _p}/{l} )\log ({ {c_1}/{c_2}})$ . Here we choose the length scale while approximating the spatial derivative to be $l$ . We assume $u_{\textit{dp}} \ll u_1, u_2$ , and is thus negligible. The initial conditions are: $c_{i}(x,0) = (c_0-c_\infty) {\rm e}^{-x^2/(2l_b^2)} + c_\infty$ , where $c_0=1, \, c_\infty =0.01$ and $n_{i}(x,0) = n_0 {\rm e}^{-x^2/(2l_b^2)}$ , where $n_0=1$ and $l_b$ is the characteristic blob size.

While it is clear that the diffusiophoretic term $\boldsymbol {D}_{\perp ,{dp}}$ in (2.6) is not a diffusive term, we can gain insight by approximating it as one (Raynal et al. Reference Raynal, Bourgoin, Cottin-Bizonne, Ybert and Volk2018). Our 2-D simulations showed that diffusiophoresis acts against the transverse dispersion of colloids in the attractive case, and acts together with it in the repulsive case. Thus, we can simplify $(\boldsymbol {D}_{\perp ,c}+\boldsymbol {D}_{\perp ,{dp}})\boldsymbol{\cdot }\boldsymbol {n} \approx \boldsymbol {D}^{\textit{eff}}_{\perp ,{\textit{att}}}\boldsymbol{\cdot }\boldsymbol {n}=(D^{\textit{eff}}_{\perp ,{\textit{att}}}/l^2)\mathbb{A}\boldsymbol{\cdot }\boldsymbol {n}=(D_{\perp ,c}/l^2+D_{\perp ,\textit{att}}(t)/l^2)\mathbb{A}\boldsymbol{\cdot }\boldsymbol {n}$ , with $D^{\textit{eff}}_{\perp ,{\textit{att}}}\lt D_{\perp ,{c}}$ for the attractive case, and $\approx \boldsymbol {D}^{\textit{eff}}_{\perp ,{\textit{rep}}}\boldsymbol{\cdot }\boldsymbol {n}=(D^{\textit{eff}}_{\perp ,{\textit{rep}}}/l^2)\mathbb{A}\boldsymbol{\cdot }\boldsymbol {n}=(D_{\perp ,c}/l^2+D_{\perp ,{\textit{rep}}}(t)/l^2)\mathbb{A}\boldsymbol{\cdot }\boldsymbol {n}$ , with $D^{\textit{eff}}_{\perp ,{\textit{rep}}}\gt D_{\perp ,{c}}$ for the repulsive case. Here, we take $D_{\perp , {att/rep}}(t)=0.5\varGamma _p\max [\log {c_1(t)/c_2(t)}]$ . Therefore, we expect the diffusiophoretic colloids to display behaviour similar to that of a solute with a modified transverse dispersion at early times: $D^{\textit{eff}}_{\perp ,{\textit{att}}}\lt D_{\perp ,{c}}\lt D^{\textit{eff}}_{\perp ,{\textit{rep}}}$ . We note that as $t\rightarrow \infty$ , $D^{\textit{eff}}_{\perp ,{\textit{att}}}\rightarrow D_{\perp ,c}$ and $D^{\textit{eff}}_{\perp ,{\textit{rep}}}\rightarrow D_{\perp ,c}$ .

We proceed to analyse the simplified two-layer transport equation with a transverse dispersion coefficient of $D_{\perp ,j}$ , dropping the superscript ‘eff’ for brevity:

(2.7) \begin{align} \partial _t \boldsymbol {n} + \boldsymbol {U}\boldsymbol{\cdot }\partial _x\boldsymbol {n} - D_c\partial _x^2\boldsymbol {n} + \big(D_{\perp ,j}/l^2\big)\mathbb{A}\boldsymbol{\cdot }\boldsymbol {n} = 0, \end{align}

where $j={\textit{att}}, {c}, {\textit{rep}}$ . Taking the Fourier transform of this equation, we then obtain the eigenmodes of the system, and analyse it in the diffusion- and advection-dominated regimes. On Fourier transformation, we obtain $ \partial _t \hat {\boldsymbol {n}} = ( -i\boldsymbol {U} k -k^2 D_c \boldsymbol {I} - (D_{\perp ,j}/l^2)\mathbb{A} ) \boldsymbol{\cdot }\hat {\boldsymbol {n}}$ . Seeking solutions of the form $\hat {\boldsymbol {n}} = e^{\lambda _\pm t} \hat {\boldsymbol {n}}_{\lambda _{{\pm }}}$ yields the following two eigenvalues:

(2.8) \begin{equation} \lambda _{\pm } = -D_ck^2 - ik\bar {u} - D_{\perp ,j}/l^2 \pm \beta (k), \end{equation}

where $\beta (k) = \sqrt {(\bar {\bar {u}}^2 - \bar {u}^2)k^2 + D_{\perp ,j}^2/l^4}$ , and $\bar {u}$ and $\bar {\bar {u}}$ are the arithmetic mean (AM) and geometric mean (GM) of $u_1$ and $u_2$ , respectively. The AM–GM inequality implies that the coefficient of the $k^2$ term in $\beta (k)$ is necessarily negative, and $\beta (k)$ can become complex. The difference between the AM and GM of the velocities is indicative of shear in the two-layer model, allowing us to define a ‘shear parameter’ $\gamma = \sqrt {\bar {u}^2 - \bar {\bar {u}}^2}$ . We then have $\beta (k)=\sqrt {-\gamma ^2k^2+D_{\perp ,j}^2/l^4}$ . In the limit when $\gamma k \ll D_{\perp ,j}/l^2$ , transverse dispersion dominates shear, leading to $\beta (k) \approx D_{\perp ,j}/l^2 - \gamma ^2l^2 k^2/2D_{\perp ,j}$ . The eigenvalues then are

(2.9) \begin{align} \lambda _+ & = -\big [ D_c + \gamma ^2l^2/2D_{\perp ,j} \big ] k^2 - i\bar {u}k , \\[-12pt]\nonumber\end{align}

(2.10) \begin{align} \lambda _- & = -\big [ D_c - \gamma ^2l^2/2D_{\perp ,j} \big ] k^2 - i\bar {u}k - 2D_{\perp ,j}/l^2, \end{align}

where the first term corresponds to the dispersion and the second term corresponds to advection. The constant term in (2.10) is a sink term. We find that both modes are advected with the same mean velocity $\bar {u}$ . In the $\lambda _+$ mode, $D_{||}=D_c+\gamma ^2l^2/2D_{\perp ,j}$ , i.e. an enhanced longitudinal dispersion. This mode persists as $t\rightarrow \infty$ , while the $\lambda _-$ mode decays due to the sink term. In the opposite limit, where $\gamma k \gg D_{\perp ,j}/l^2$ , we have $\beta (k) \approx i (\gamma k - ( {D_{\perp ,j}^2}/{2k\gamma l^4})) \approx i\gamma k$ . The eigenvalues are

(2.11) \begin{align} \lambda _{+} & = -D_ck^2 - i(\bar {u}-\gamma )k - D_{\perp ,j}/l^2 , \\[-12pt]\nonumber\end{align}

(2.12) \begin{align} \lambda _{-} & = -D_ck^2 - i(\bar {u}+\gamma )k - D_{\perp ,j}/l^2, \end{align}

where the existence of two modes migrating at different velocities ( $\bar {u}-\gamma$ and $\bar {u}+\gamma$ ) manifests as a ‘splitting’ of an initially unimodal distribution. Recalling that $D_{\perp ,{\textit{att}}}\lt D_{\perp ,{c}}\lt D_{\perp ,{\textit{rep}}}$ , we can gain some insight by assuming that the repulsive case corresponds to the diffusion-dominant regime $\gamma k \ll D_{\perp ,{\textit{rep}}}/l^2$ , leading to $D_{||,{\textit{rep}}}=D_c+\gamma ^2l^2/2D_{\perp ,{\textit{rep}}} \lt D_c+\gamma ^2l^2/2D_{\perp ,c}$ , i.e. suppressed longitudinal dispersion compared with the control case. In the opposite limit of shear-dominated transport, where $\gamma k \gg D_{\perp ,{\textit{att}}}/l^2$ , we obtain the two modes travelling with different velocities, leading to the splitting, and enhanced dispersion.

These results can also be interpreted in the broader context of chromatographic dispersion in porous media. In the simplified trapping model discussed by Bouchaud & Georges (Reference Bouchaud and Georges1990), a particle advected along a mobile backbone with velocity $V$ intermittently enters stagnant regions with probability $p$ , acquiring a residence time $\tau$ over a characteristic longitudinal distance $\xi$ . Using their notation, this gives an effective velocity $U^{-1}=V^{-1}+p\langle \tau \rangle /\xi$ and a longitudinal dispersion coefficient $D_{||}=pU^3(\langle \tau ^2\rangle -p\langle \tau \rangle ^2)/(2\xi )$ . Equivalently, in the small- $p$ limit, this can be written as $D_{||}=({1}/{2})fU^2\langle \tau ^2\rangle /\langle \tau \rangle$ , where $f$ is the fraction of time spent in stagnant regions, showing that the dispersion scale is set by the ratio $\langle \tau ^2\rangle /\langle \tau \rangle$ . Thus, longitudinal dispersion is controlled by fluctuations in the time spent in different velocity states. Our two-layer model retains this same physical ingredient, but replaces mobile–immobile exchange by exchange between two mobile layers with different velocities and transverse exchange rate $\kappa =D_{\perp ,j}/l^2$ (equivalently, an inter-layer residence time $\tau _\perp =\kappa ^{-1}$ ). In the long-wavelength limit, the model gives $D_{||}=D_c+\gamma ^2/(2\kappa )=D_c+\gamma ^2\tau _\perp /2$ , showing explicitly that reduced transverse exchange enhances longitudinal dispersion, whereas increased exchange homogenises the blob and suppresses dispersion. In the present system, diffusiophoresis dynamically and sign-dependently modifies this exchange: attraction lowers the effective transverse exchange and, in the shear-dominated regime $\gamma k\gg \kappa$ identified above, can lead to bimodal longitudinal splitting, whereas repulsion increases the exchange and suppresses dispersion.

To verify these predictions, we numerically solve the coupled solute and colloid fields ((2.5) and (2.6)), keeping the parameter values the same as in the 2-D simulations. The perpendicular dispersion coefficients are $D_{\perp ,s}\sim D_s$ and $D_{\perp ,c}\sim 10 D_c$ to reflect transverse dispersion due to geometric disorder. The layer thickness $l$ is chosen as the pore width for a hexagonal lattice given by $l=\lambda ( 1-(2\sqrt {3}(1-\phi )/\pi )^{1/2} )$ , which amounts to $l\approx \lambda /4$ for $\phi =0.5$ . Velocities $u_1$ and $u_2$ are such that $u_1+u_2=2U_m$ and $u_1-u_2 \sim 0.5 \langle |\partial u/\partial y|\rangle l$ . The mean shear rate is estimated as $\langle |\partial u/\partial y|\rangle \sim U_m/(l/2) \sim U_m/(\lambda /8)$ . The solution to the full numerical model ((2.5) and (2.6)) is plotted as the average concentration of both layers $n \equiv (n_1+n_2)/2$ (figure 4 b), along with that of the simplified model ((2.5) and (2.7)) (figure 4 d). Both the full numerical and simplified two-layer models capture the longitudinal split of the colloid concentration in the attractive case, and inhibited dispersion with a unimodal distribution in the repulsive case (figure 4 b,d). They also capture the observed qualitative trends for both longitudinal dispersion and dispersion coefficient in the presence of solute gradients (figure 4 c,e), showing that the minimal necessary ingredients are shear and the transverse diffusiophoretic exchange of colloids.

2.5. Impact of geometric disorder on colloid dispersion

To alter the velocity distribution in the porous medium, we introduce disorder into the ordered hexagonal lattice and study four different configurations of increasing disorder via simulations: $\chi = 0.0 \ \text{(ordered)},\ 0.1,\ 0.2,\ 0.36$ (figure 5 a), maintaining the mean flow velocity $\langle u \rangle$ constant (figure 5 b). As we increase the disorder, the mean transverse velocity magnitude, $\langle |v| \rangle$ , increases (figure 5 b). The mean transverse shear, however, decreases with disorder (figure 5 b inset). These trends are indicative of transition from shear-dominated to diffusion-dominated transport, associated with enhanced transverse dispersion and suppressed longitudinal dispersion. In the attractive case, this transition implies the suppression of splitting, and in the repulsive case, a more compact blob (figure 5 d,e). While dispersion is suppressed with disorder (figure 5 c), the relative impact of solute gradients on longitudinal dispersion remains constant (figure 5 c inset). Here, we have limited our study to weakly disordered media to avoid overlap of the posts. Many natural porous media, however, are strongly disordered, including regions of low permeability and dead-end zones, where diffusiophoretic transport is expected to dominate over advection (Park et al. Reference Park, Lee, Yoon and Shin2021; Jotkar et al. Reference Jotkar, de Anna, Dentz and Cueto-Felgueroso2024). These stagnant zones facilitate colloid trapping while preferential flow pathways advect colloids; together, this allows the observation of bimodal dispersion of colloids to persist for strongly disordered media. Beyond dead-end pores, recent work by Alipour et al. (Reference Alipour, Li, Liu and Pahlavan2026) showed that diffusiophoretic migration of colloids within preferential flow pathways of porous media leads to significant changes of macroscopic dispersion of colloids. Therefore, we expect the influence of phoretic migration on dispersion of colloidal blobs to persist for highly disordered porous media.

Figure 5.

The interplay of disorder and diffusiophoresis modulates the macroscopic dispersion of colloids. (a) Velocity distribution within the medium for different disorder strengths $\chi$ . (b) The mean longitudinal velocity $\langle \tilde {u} \rangle$ is held constant throughout the simulations. The mean absolute transverse velocity $\langle |\tilde {v}| \rangle$ increases with disorder, while the mean absolute transverse shear $\langle |\partial \tilde {u} / \partial \tilde {y}| \rangle$ decreases (inset). (c) Dispersion coefficient versus disorder at time $\tilde {t}=48$ for $\textit{Pe}_c \sim 7200$ . Inset: corresponding ratio of dispersion coefficient with respect to the control case. (d,e) Geometric disorder suppresses the influence of phoresis as the transverse velocity magnitude dominates phoretic exchange between fast and slow zones.