Diffusioosmotic dispersion

Jian Teng Center for Fluid Mechanics, Brown University, Providence, Rhode Island 02912, USA Bhargav Rallabandi Department of Mechanical Engineering, University of California, Riverside, CA 92521, USA Jesse T. Ault [email protected] Center for Fluid Mechanics, Brown University, Providence, Rhode Island 02912, USA

Abstract

Solute-surface interactions have garnered considerable interest in recent years as a novel control mechanism for driving unique fluid dynamics and particle transport with potential applications in fields such as biomedicine, the development of microfluidic devices, and enhanced oil recovery. In this study, we will discuss dispersion induced by the diffusioosmotic motion near a charged wall in the presence of a solute concentration gradient. Here, we introduce a plug of salt with a Gaussian distribution at the center of a channel with no background flow. As the solute diffuses, the concentration gradient drives a diffusioosmotic slip flow at the walls, which results in a recirculating flow in the channel; this, in turn, drives an advective flux of the solute concentration. This effect leads to cross-stream diffusion of the solute, altering the effective diffusivity of the solute as it diffuses along the channel. We derive theoretical predictions for the solute dynamics using a multiple-timescale analysis to quantify the dispersion driven by the solute-surface interactions. Furthermore, we derive a cross-sectionally averaged concentration equation with an effective diffusivity analogous to that from Taylor dispersion. In addition, we use numerical simulations to validate our theoretical predictions.

I Introduction

In fluid dynamics, dispersion is typically used to denote transport of a species from high to low concentrations due to non-uniform flow conditions. This is in contrast to diffusion, which denotes the similar transport of a species from high to low concentrations but due to Brownian motion. Since relatively few flow systems actually transport flow in a plug-like way, dispersion is a nearly ubiquitous transport mechanism in systems including microfluidic devices, filtration systems, chemical reactor systems, medical devices, lab-on-a-chip systems, and many others. A classical and simple example of dispersion is that which has come to be known as Taylor dispersion, which was first studied by Taylor [41] and later generalized by Aris [5]. Taylor dispersion describes the enhanced diffusivity that a diffusing species experiences in the presence of a shear flow, such as the diffusion of a solute concentration field in a pipe or channel flow. In both of these cases, the advection-diffusion equation governing the solute transport can be averaged over the cross-section to yield a 1D model for the depth-averaged concentration, which experiences an enhanced effective diffusivity that is a function of the Peclét number governing the transport. A simple physical understanding of this Taylor dispersion can be achieved by supposing we have a step initial condition in the concentration $c$ of some solute species. For example, suppose we have a Poiseuille flow in a pipe, and we suddenly add solute to the system such that $c(x<0)=c_{1}$ and $c(x>=0)=c_{2}$ , where $-\infty<x<\infty$ is the axial coordinate of our infinite pipe. Then, in the limit of no background flow, the interface between the two solute concentrations will smear out by diffusion in a purely 1D process. However, as the relative importance of fluid advection to solute diffusion (i.e., the Peclét number) increases, shear flows in the system distort the interface between the two solute concentrations as they diffuse, introducing cross-stream diffusion and enhancing the rate of axial diffusion of the cross-sectionally averaged concentration. For more rigorous background into the Taylor dispersion phenomenon, see (1) Barton [8] who extended the methods and results of Aris [5] to consider all times rather than just the asymptotic behavior, (2) Frankel and Brenner [18] who developed a generalized Taylor dispersion theory to greatly extend the ideas of Taylor and Aris to whole classes of other problems such as porous media flows or even sedimenting particles, and (3) Brenner and Edwards [12] who provide a comprehensive overview of the theory of macrotransport processes.

Many other studies and applications have built upon the theoretical work on Taylor dispersion by Taylor and Aris. Here we just briefly mention a few examples. Stone and Brenner [40] extended the theories of Taylor dispersion to consider laminar flows with velocity variations in the streamwise direction. Aminian et al. [3] studied the role of the channel boundary shape and aspect ratio on the dispersion as a means to control the delivery of chemicals in microfluidics. Salmon and Doumenc [31] studied the solute dispersion induced by buoyancy-driven flow and developed analytical methods analogous to Taylor dispersion. Chu et al. [14] added the ideas of oscillatory pressure-driven flow in a parallel-plate channel flow as well as patterned slip walls and investigated the role of both of these effects on the dispersion process. Finally, Chu et al. [15] coupled the Taylor dispersion in a pipe flow to the transport of a second species consisting of particles or bacteria via diffusiophoresis, which was also explored by Migacz and Ault [26]. This list is by no means exhaustive, but one unifying theme that is common to many of the studies related to dispersion is the role of imposed pressure gradients or moving boundaries to drive the shear flows that cause the dispersion. Typically, the transport of the solute is passive in the sense that it is not expected to couple to and alter the background fluid dynamics. However, there are many scenarios where the transport of solute is fully coupled to the fluid dynamics, which is the focus of this paper. In particular, we consider the case where there is no background pressure-driven flow and no moving boundaries, but where the solute interacts with the boundaries of the system to drive fluid flow via diffusioosmosis which is the spontaneous motion of fluid near a surface in the presence of a solute concentration gradient.

The physical origin of diffusioosmosis was discovered by Derjaguin et al. [16], where Derjaguin and his coworkers showed experimentally that a local solute concentration gradient near a boundary could induce a slip-flow boundary condition over the surface. Since then, significant theoretical progress has been made towards understanding and modeling this effect, and in the context of dispersion, a variety of studies have demonstrated the potential for diffusioosmosis to alter the transport of suspended species in confined geometries [30, 7, 35, 33, 23, 22, 13]. Diffusioosmosis is also closely related to the analogous phenomenon of diffusiophoresis. Whereas diffusioosmosis refers to the motion of fluid next to a surface in the presence of a chemical concentration gradient, diffusiophoresis refers to the reciprocal motion of suspended particles in a concentration gradient, which results due to the slip flow at the surface. Both diffusioosmosis and diffusiophoresis can have contributions due to chemi-osmosis/-phoresis that arise from the osmotic pressure gradient over the surface and electro-osmosis/-phoresis that arise in the case of electrolyte solutes with mismatched anion and cation diffusivities [4]. Coupled diffusioosmosis and diffusiophoresis have been the subject of a variety of recent studies, especially concerning the coupled transport of solutes and colloidal particles in confined geometries where convective transport is difficult to achieve [37, 22, 34, 38, 42, 1, 13].

Such coupled motions have also been the focus of studies in a variety of other natural and engineering settings, such as in underground porous media flows [29], water filtration systems [17, 35, 24, 10], microfluidic devices [28, 36], fabric cleaning systems [37], particle delivery methods [6], energy storage applications [21, 19], and many others. Typically, in such studies the role of diffusioosmosis has been a secondary effect, or neglected entirely, although a collection of studies on the motion of solutes and particles in narrow pores has considered both effects [34, 38, 6, 42, 1, 2]. In such systems, diffusioosmosis can play an essential role on the transport of both the solutes, the fluid, and any suspended particles [32]. In the present study, we investigate the dispersion of solute in a channel where the only fluid motion is driven by diffusioosmosis at the channel walls. That is, the diffusion of the initial solute concentration results in a local concentration gradient at the walls that in turn drives a fluid recirculation via diffusioosmosis. The resulting shear flows then in turn alter the transport of the solute concentration by a mechanism analogous to that of Taylor dispersion.

As mentioned, several studies have already considered the coupling between Taylor dispersion and diffusiophoresis. Specifically, Chu et al. [15] studied the diffusiophoretic dynamics of charged colloidal particles or bacteria in a solute concentration field that was experiencing Taylor dispersion. This work considered a one-way coupling where the fluid flow drives Taylor dispersion of the solute field and the particle concentration field, and the particles receive an additional contribution to their motion from diffusiophoresis via the solute field. They developed theoretical and numerical solutions in the long-time regime following an approach similar to that of Taylor and Aris. More recently, Migacz and Ault [26] built upon this work by developing solutions that are valid for both the early- and long-time regimes. However, neither of these studies considered the diffusioosmosis at the channel walls. One such study that did consider these effects was that of Alessio et al. [2], who considered the diffusioosmotic slip boundary condition and derived a multi-dimensional effective dispersion equation for solute transport into a dead-end pore similar to that of Taylor. In contrast to the previous studies, while the velocity profile in Alessio’s work is still approximately parabolic (as in classical Taylor dispersion), the magnitude is a function of position and time in the channel as the solute concentration evolves. We will find a similar behavior in our system.

To the best of our knowledge, no theoretical solutions have previously been developed to describe the diffusioosmosis-driven analog of Taylor dispersion. In this study, we take motivations from the works of Alessio et al. [2], Chu et al. [15], and Migacz and Ault [26] and develop analytical solutions for the diffusioosmosis-driven dispersion of a plug of solute in a channel. We consider a plug of solute that is initially normally distributed in a Gaussian peak at the center of the channel (see figure 1). As the solute diffuses, the local concentration gradient at the wall drives an effective slip boundary condition via diffusioosmosis that is dependent on the charge of the surface. This slip at the walls drives a recirculating flow in the channel. The recirculation contributes to the advection of the solute transport and introduces cross-stream diffusion, which alters the effective diffusivity of the solute along the channel. In the modeling process, we use a perturbation method to derive analytical solutions to the coupled fluid and solute dynamics. The theoretical analysis is performed for transport in a long, narrow 2D channel using a 2D Cartesian coordinate system and in a long, narrow cylindrical pipe using a 2D axisymmetric cylindrical coordinate system. In $\S$ II, we introduce the governing equations and boundary conditions for both systems. We then apply a perturbation method along with a multiple timescale analysis to theoretically solve for the fluid and solute dynamics. In $\S$ II.4, we derive the effective diffusivity of the cross-sectionally averaged solute concentration analogously to that of Taylor dispersion. In $\S$ III, we perform numerical simulations to solve for the fluid and solute dynamics and show good agreement with the theoretical predictions. In $\S$ IV, we analyze the dispersion behavior for various conditions and in different time regimes.

II Modeling diffusioosmotic dispersion in a long, narrow channel

In this section, we model the coupled fluid and solute transport for a diffusing plug of solute in a channel in the presence of diffusioosmosis-driven recirculation. We consider two configurations corresponding to planar and cylindrical channels and describe the flows in these configurations using Cartesian and cylindrical coordinates, respectively. The channel configurations for both systems are shown in figure 1. Initially, we introduce a plug of solute with a Gaussian distribution centered around the origin. In cases when the solute molecules/ions do not interact with the channel walls, the dynamics of the solute transport are governed by simple Brownian diffusion. Here, however, we consider the case where the channel walls have a non-zero surface charge and solute-surface interactions cannot be neglected. The local solute concentration gradient at the channel walls will induce a diffusioosmotic slip velocity boundary condition, which will in turn drive a recirculation in the channel as the solute diffuses. The magnitude of this slip velocity boundary condition is given by $u_{\textrm{wall}}=\Gamma_{w}\nabla_{\parallel}\ln c$ , where $u_{\textrm{wall}}$ is the velocity at the wall in the direction parallel to the wall, and the gradient is taken parallel to the surface. Here, $\Gamma_{w}$ is the diffusioosmotic mobility coefficient, which is a function of the surface charge of the channel walls, and $c$ is the solute ion concentration.

Refer to caption — Figure 1: Problem setup. We consider an initially Gaussian distributed plug of salt with characteristic length $\ell$ in both (a) 2D Cartesian coordinates and (b) axisymmetric cylindrical coordinates. Both channels are infinite in the axial direction, and $u_{\textrm{wall}}$ represents the slip velocity at the wall induced by diffusioosmosis.

II.1 Governing equations

The fluid and solute dynamics in the system are governed by the coupled continuity and incompressible Navier-Stokes equations and an advection-diffusion equation for the solute transport. Here, the fluid pressure, density, viscosity, and velocity are given respectively by $p$ , $\rho$ , $\mu$ , and $\boldsymbol{u}=(u,v)~{}\textrm{or}~{}(u_{r},u_{z})$ . The solute concentration and diffusivity are given by $c(x,y,t)$ and $D_{s}$ , respectively. Here, we consider $Re=\frac{\rho Uh}{\mu}\ll 1$ , where $U$ is some characteristic flow velocity in the axial direction, we neglect the influence of gravity, and we treat the fluid dynamics as quasi-steady. The dimensional form of the governing equations is given by


	$\displaystyle\nabla^{}\cdot\boldsymbol{u}^{}=0,$		(1a)
	$\displaystyle-\nabla^{}p^{}+\mu\nabla^{2}\boldsymbol{u}^{}=\boldsymbol{0},$		(1b)
	$\displaystyle\frac{\partial c^{}}{\partial t^{}}+\nabla^{}\cdot(\boldsymbol{u^{}}c^{})=D_{s}\nabla^{2}c^{*},$		(1c)

where asterisks denote dimensional quantities. The 2D axisymmetric cylindrical and the 2D Cartesian cases can be solved similarly and share similar boundary conditions. For the sake of simplicity, we only show the derivation for the 2D channel flow case in the main text and refer the reader to Appendix B for the derivation for the pipe flow case. The long, narrow channel has a height of $h$ in the Cartesian coordinate system. The channel is infinitely long, and the initial Gaussian solute distribution has a characteristic width of $\ell$ , and we will seek solutions in the limit that $h\ll\ell$ . To begin, we nondimensionalize the system of equations as follows,

\begin{split}&x=\frac{x^{*}}{\ell},\;y=\frac{y^{*}}{h},\;u=\frac{u^{*}}{U},\;v=\frac{v^{*}\ell}{Uh},\;p=\frac{p^{*}h^{2}}{\mu U\ell},\;\epsilon=\frac{h}{\ell},\;U=\frac{D_{s}}{\ell},\;t=\frac{t^{*}}{\ell^{2}/D_{s}},\end{split}

(2)

where $U=D_{s}/\ell$ is the characteristic speed of solute diffusion along the channel, and $\ell^{2}/D_{s}$ is the characteristic time of diffusion along the channel. With these scalings, we can rewrite the governing equations (1) in nondimensional form as,


	$\displaystyle 0=\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y},$		(3a)
	$\displaystyle 0=-\frac{\partial p}{\partial x}+\epsilon^{2}\frac{\partial^{2}u}{\partial x^{2}}+\frac{\partial^{2}u}{\partial y^{2}},$		(3b)
	$\displaystyle 0=-\frac{\partial p}{\partial y}+\epsilon^{4}\frac{\partial^{2}v}{\partial x^{2}}+\epsilon^{2}\frac{\partial^{2}v}{\partial y^{2}},$		(3c)
	$\displaystyle\epsilon^{2}\frac{\partial c}{\partial t}+\epsilon^{2}u\frac{\partial c}{\partial x}+\epsilon^{2}v\frac{\partial c}{\partial y}=\epsilon^{2}\frac{\partial^{2}c}{\partial x^{2}}+\frac{\partial^{2}c}{\partial y^{2}}.$		(3d)

The solution to the governing equations (3) is subject to boundary conditions on the fluid and solute. These boundary conditions can be summarized by:

\text{Quiescent far-field condition: }p=0,~{}~{}u=0,~{}~{}c=0~{}~{}\text{at}~{}~{}x=\pm\infty,

(4)

\text{No fluid penetration at the walls: }v=0~{}\text{at}~{}y=\pm\frac{1}{2},

(5)

\text{No-flux conditions at the channel walls: }\frac{\partial c}{\partial y}=0~{}\text{at}~{}y=\pm\frac{1}{2},

(6)

\text{Diffusioosmotic wall slip boundary condition: }u=\frac{\Gamma_{w}}{D_{s}}\frac{\partial\ln c}{\partial x}~{}\text{at}~{}y=\pm\frac{1}{2}.

(7)

The slip boundary condition given by (7) is induced by diffusioosmosis, which drives the flow inside the channel. In this study, we assume a constant zeta potential and diffusioosmotic mobility coefficient $\Gamma_{w}$ . This assumption is needed in order to achieve a final analytical solution, and is a reasonable approximation under many scenarios as discussed by several recent previous works (see, e.g., Ault et al. [6], Migacz and Ault [26], Lee et al. [25], Gupta et al. [20], and Alessio et al. [2]). We further focus our attention on cases where the initial condition is already spread out relative to the channel width, such that $\epsilon\ll 1$ , in which case the lubrication approximation can be used to simplify the governing equations.

II.2 Leading-order fluid dynamics

In our original nondimensionalization above, we chose a characteristic timescale that is the characteristic time for solute diffusion along the channel $\ell^{2}/D_{s}$ . This corresponds to the slow dynamics of the system. The other important timescale in the system is the characteristic time for solute diffusion across the channel, $h^{2}/D_{s}$ . This corresponds to the fast dynamics of the system, and the two timescales are separated by a factor of $\epsilon^{2}=h^{2}/\ell^{2}\ll 1$ . Here, in order to develop a solution that is uniformly valid across both the early and late dynamics, we use an approach similar to that of Migacz and Ault [26].

Following this approach, we introduce a multiple timescale analysis [9] in which we introduce a fast time variable $T=t/\epsilon^{2}$ . That is, $T=O(1)$ over dimensional times $\sim h^{2}/D_{s}$ corresponding to $t=O(\epsilon^{2})$ , whereas $t=O(1)$ on dimensional times $\sim\ell^{2}/D_{s}$ . Thus, we can map any time-dependent quantity as

f(t)\mapsto f(t,T)\Rightarrow\frac{\partial f}{\partial t}\mapsto\frac{\partial f}{\partial t}+\frac{\partial T}{\partial t}\frac{\partial f}{\partial T}=\frac{\partial f}{\partial t}+\epsilon^{-2}\frac{\partial f}{\partial T}.

(8)

With this mapping, the advection-diffusion equation (3 $d$ ), becomes

\epsilon^{2}\frac{\partial c}{\partial t}+\frac{\partial c}{\partial T}+\epsilon^{2}u\frac{\partial c}{\partial x}+\epsilon^{2}v\frac{\partial c}{\partial y}=\epsilon^{2}\frac{\partial^{2}c}{\partial x^{2}}+\frac{\partial^{2}c}{\partial y^{2}}.

(9)

We seek an analytical solution of the governing equations as perturbation expansions in the small parameter $\epsilon^{2}$ in the limit of $Re\ll 1$ . We seek solutions of the form


$\displaystyle c(x,y,t,T)$	$\displaystyle=c_{0}(x,y,t,T)+\epsilon^{2}c_{1}(x,y,t,T)+\epsilon^{4}c_{2}(x,y,t,T)+\ldots,$	(10a)
$\displaystyle p(x,y,t,T)$	$\displaystyle=p_{0}(x,y,t,T)+\epsilon^{2}p_{1}(x,y,t,T)+\epsilon^{4}p_{2}(x,y,t,T)+\ldots,$	(10b)
$\displaystyle u(x,y,t,T)$	$\displaystyle=u_{0}(x,y,t,T)+\epsilon^{2}u_{1}(x,y,t,T)+\epsilon^{4}u_{2}(x,y,t,T)+\ldots,$	(10c)
$\displaystyle v(x,y,t,T)$	$\displaystyle=v_{0}(x,y,t,T)+\epsilon^{2}v_{1}(x,y,t,T)+\epsilon^{4}v_{2}(x,y,t,T)+\ldots.$	(10d)

The initial condition of the solute concentration is $c(t=T=0)=\exp(-x^{2})$ . First, we need to obtain the leading-order velocity and pressure solutions. These can be obtained by substituting equations (10) into equations (3), which gives


	$\displaystyle 0=\frac{\partial^{2}u_{0}}{\partial y^{2}}+\frac{\partial p_{0}}{\partial x},$		(11a)
	$\displaystyle\frac{\partial p_{0}}{\partial y}=0,$		(11b)
	$\displaystyle\frac{\partial u_{0}}{\partial x}+\frac{\partial v_{0}}{\partial y}=0,$		(11c)
	$\displaystyle\frac{\partial c_{0}}{\partial T}=\frac{\partial^{2}c_{0}}{\partial y^{2}},$		(11d)

to leading order. Considering both the initial condition and the no-flux conditions at the channel walls, i.e., $\frac{\partial c_{0}}{\partial y}(y=\pm 1/2)=0$ , it must be true that $c_{0}(x,y,t,T)=c_{0}(x,t)$ , with $c_{0}(t=0)=\exp(-x^{2})$ . Furthermore, equation (11 $b$ ) indicates that $p_{0}$ is only a function of $x$ and $t$ .

To obtain the leading-order velocities, we first take equation (11 $a$ ) and solve for $u_{0}$ , which is subject to the slip boundary condition (7). This gives

u_{0}=\frac{\Gamma_{w}}{D_{s}c_{0}}\frac{\partial c_{0}}{\partial x}+\frac{1}{8}(-1+4y^{2})\frac{\partial p_{0}}{\partial x}.

(12)

Here, $u_{0}$ has a term that includes $p_{0}(x,t)$ , which can be found by considering the conservation of mass and integrating over the channel cross-section $A$ , i.e.,

\frac{1}{A}\iint_{A}u_{0}\,dA=0.

(13)

Following this approach, $p_{0}(x,t)$ , is found to be

p_{0}(x,t)=\frac{12\Gamma_{w}}{D_{s}}\ln c_{0}.

(14)

With this result, $u_{0}$ can be simplified and written as

u_{0}(x,y,t)=\frac{\Gamma_{w}}{2D_{s}}\frac{\partial\ln c_{0}}{\partial x}(-1+12y^{2}).

(15)

To solve for $v_{0}(x,y,t)$ , we integrate the continuity equation (11 $c$ ) and apply the boundary condition given by (5). This gives the leading-order $y$ -component of velocity to be

v_{0}(x,y,t)=\frac{\Gamma_{w}}{2D_{s}}\frac{\partial^{2}\ln c_{0}}{\partial x^{2}}(1-4y^{2})y.

(16)

As mentioned, the equivalent results for the flow in a cylindrical pipe can be found in Appendix B.

II.3 Higher-order solute transport

The higher-order solute concentration results from diffusioosmosis, which causes the deviation from pure diffusion. Using the leading-order velocity profiles found above, we seek a solution for the higher-order solute dynamics from (3 $d$ ). Substituting our asymptotic expansion into the advection-diffusion equation (9), we find, to leading-order, that

\frac{\partial c_{1}}{\partial T}+\frac{\partial c_{0}}{\partial t}+u_{0}\frac{\partial c_{0}}{\partial x}=\frac{\partial^{2}c_{0}}{\partial x^{2}}+\frac{\partial^{2}c_{1}}{\partial y^{2}}.

(17)

The term involving $v$ has disappeared since $\frac{\partial c_{0}}{\partial y}=0$ . To find a solution to this problem, we first consider long times such that $T\gg 1$ , but $t$ is finite. In this limit, $\frac{\partial c_{1}}{\partial T}$ is small. Then, averaging equation (17) over the channel cross-section gives

\frac{\partial c_{0}}{\partial t}=\frac{\partial^{2}c_{0}}{\partial x^{2}}.

(18)

The solution to this problem is given by

c_{0}(x,t)=\frac{1}{\sqrt{1+4t}}\exp{\left(-\frac{x^{2}}{1+4t}\right)},

(19)

which was previously developed by [20]. The advection-diffusion equation (17) can then be simplified to

\frac{\partial c_{1}}{\partial T}+\frac{\Gamma_{w}}{2D_{s}}\frac{\partial\ln c_{0}}{\partial x}(-1+12y^{2})\frac{\partial c_{0}}{\partial x}=\frac{\partial^{2}c_{1}}{\partial y^{2}},

(20)

subject to the initial condition $c_{1}(t=T=0)=0$ . To solve this, we note that at long times the fast-time dynamics should have all decayed such that the time derivative term can be ignored, and the equation can be integrated to yield

c_{1}(T\rightarrow\infty)\sim c_{1}^{\infty}(x,y,t)=\frac{y^{2}}{4c_{0}}\frac{\Gamma_{w}}{Ds}(-1+2y^{2})\left(\frac{\partial c_{0}}{\partial x}\right)^{2}+B(x,t)

(21)

where $B(x,t)$ is a yet unknown function that results from the integration. To solve for $B(x,t)$ , we substitute (21) into (9), and take the cross-sectional average of the equation, which gives

\frac{\partial^{2}B}{\partial x^{2}}=\frac{\Gamma_{w}}{D_{s}}\frac{e^{-\frac{-x^{2}}{1+4t}}}{420(1+4t)^{9/2}}\left[49(1+4t)^{2}+48\frac{\Gamma_{w}}{D_{s}}(1+4t)x^{2}-32\frac{\Gamma_{w}}{D_{s}}x^{4}\right]+\frac{\partial B}{\partial t}.

(22)

Note that, until now, we have left the results in terms of a general $c_{0}$ , but here and in the solution for $B$ below, we have substituted the specific solution for $c_{0}$ shown above since it is needed to solve the equation (22) using the Fourier transform approach. This procedure can be repeated for other arbitrary initial conditions as needed. Using a Fourier transform approach, $B(x,t)$ can be found to be

B=\frac{\Gamma_{w}}{D_{s}}\frac{e^{-\frac{-x^{2}}{\alpha}}}{840\alpha^{9/2}}\left(49(\alpha x)^{2}+16\frac{\Gamma_{w}}{D_{s}}t\left(3\alpha^{2}-12\alpha x^{2}+4x^{4}\right)-12\frac{\Gamma_{w}}{D_{s}}\alpha^{2}(\alpha-2x^{2})\ln(\alpha)\right),

(23)

where $\alpha(t)=1+4t$ . Note that $c_{1}^{\infty}$ only depends on the slow time $t$ and not the fast time $T$ . It does not satisfy the initial condition, so it is not yet the full solution for the higher-order solute dynamics. We seek a solution of the form $c_{1}(x,y,t,T)=c_{1}^{\infty}(x,y,t)+\hat{c}_{1}(x,y,t,T)$ . Substituting into equation (20), we find

\frac{\partial\hat{c}_{1}}{\partial T}=\frac{\partial^{2}\hat{c}_{1}}{\partial y^{2}},~{}~{}~{}~{}~{}\text{with the initial condition}~{}~{}~{}~{}\hat{c}_{1}(T=0)=-c_{1}^{\infty}(t=0),

(24)

the solution to which is

\hat{c}_{1}=\left.\frac{\Gamma_{w}}{4c_{0}D_{s}}\left(\frac{\partial c_{0}}{\partial x}\right)^{2}\right|_{t=0}\sum_{n=1}^{\infty}\frac{6(-1)^{n}}{n^{4}\pi^{4}}e^{-(2n\pi)^{2}T}\cos 2n\pi y.

(25)

We can then construct a composite solution $c_{1}=c_{1}^{\infty}+\hat{c}_{1}$ , which is valid for all $t$ . The solution of $c_{1}$ can be verified to satisfy the conservation of mass by considering

\int^{\infty}_{-\infty}\int^{1/2}_{-1/2}c_{1}(x,y,t,T)\mathrm{dy\,dx}=0.

(26)

Once again, analogous results for the coupled dynamics in a cylindrical pipe geometry can be found in Appendix B.

II.4 Effective diffusivities

As mentioned, the diffusioosmotic slip flow at the channel walls induces a recirculating flow that drives an advective transport of the solute, altering the effective diffusivity of its transport as it diffuses along the channel. Here, we seek to characterize the effective diffusivity of this transport by deriving a 1D transport equation for the cross-sectionally averaged solute transport. This approach is analogous to that of Taylor [41] and Aris [5] in their famous work on solute dispersion in the presence of pressure-driven shear flow (see also, e.g., Alessio et al. [2], Aminian et al. [3]).

To begin, we define $c^{\prime}(x,y,t)$ as the deviation of the solute concentration from its cross-sectionally averaged value $\overline{c(x,y,t)}$ , where overbars are used to denote the average over a cross-section,

c^{\prime}(x,y,t)=c(x,y,t)-\overline{c(x,t)}.

(27)

Substituting this definition into the solute advection-diffusion equation gives

\frac{\partial\overline{c}}{\partial t}+\frac{\partial c^{\prime}}{\partial t}+u\frac{\partial c^{\prime}}{\partial x}+u\frac{\partial\overline{c}}{\partial x}+v\frac{\partial c^{\prime}}{\partial y}=\frac{\partial^{2}\overline{c}}{\partial x^{2}}+\frac{\partial^{2}c^{\prime}}{\partial x^{2}}+\frac{1}{\epsilon^{2}}\frac{\partial^{2}c^{\prime}}{\partial y^{2}}.

(28)

Next, averaging this equation over the cross-section gives

\frac{\partial\overline{c}}{\partial t}+\overline{u\frac{\partial c^{\prime}}{\partial x}}+\overline{v\frac{\partial c^{\prime}}{\partial y}}=\frac{\partial^{2}\overline{c}}{\partial x^{2}}.

(29)

Here, substituting in our solutions for $u$ and $v$ from above, this can be rewritten into a 1D diffusion equation with a known forcing term given by

\frac{\partial\overline{c}}{\partial t}+\frac{\partial}{\partial x}\left[\left(\frac{\Gamma_{w}}{D_{s}}\right)^{2}\frac{\epsilon^{2}}{210}\left(\frac{\partial}{\partial x}\ln c_{0}\right)^{2}\frac{\partial c_{0}}{\partial x}\right]=\frac{\partial^{2}\overline{c}}{\partial x^{2}}.

(30)

This equation can easily be solved numerically. For the purposes of making an analogy to the classic Taylor dispersion problem, this can be rewritten into a 1D pure diffusion problem by recognizing that $\frac{\partial c_{0}}{\partial x}-\frac{\partial\overline{c}}{\partial x}=\mathcal{O}(\epsilon^{2})$ , which gives

\frac{\partial\overline{c}}{\partial t}=\frac{\partial}{\partial x}\left(D_{\textrm{eff}}\frac{\partial\overline{c}}{\partial x}\right),

(31)

where the effective diffusivity $D_{\textrm{eff}}$ is given by

D_{\textrm{eff}}=1+\left(\frac{\Gamma_{w}}{D_{s}}\right)^{2}\frac{\epsilon^{2}}{210}\left(\frac{\partial}{\partial x}\ln c_{0}\right)^{2}+\mathcal{O}(\epsilon^{4}).

(32)

From equation (32) we see that to leading-order the effective nondimensional diffusivity is $D_{\text{eff}}=1$ , and the effects of diffusioosmosis on the dispersion are $\mathcal{O}(\epsilon^{2})$ . That is, as $\epsilon\rightarrow 0$ , the initial condition of the solute plug becomes more spread out, the concentration gradients get weaker, and the diffusioosmosis becomes negligible. The same behavior occurs as $t\rightarrow\infty$ as the solute spreads out over long times. Curiously, the contribution from diffusioosmosis is also $\mathcal{O}((\Gamma_{w}/D_{s})^{2})$ . Thus, despite the non-linearity of the diffusioosmotic boundary condition, flipping the sign of $\Gamma_{w}$ (and thus the direction of the recirculation) results in the same effective diffusivity. This will be made slightly more surprising when we visualize the results below and see that the deviations in the solute concentration from $c_{0}$ are not mirror images of each other when the sign of the mobility is flipped.

Following a similar approach and using the results from Appendix B, the analogous averaged transport equation for the cylindrical pipe case is given by

\frac{\partial\overline{c}}{\partial t}+\frac{\partial}{\partial z}\left[\left(\frac{\Gamma_{w}}{D_{s}}\right)^{2}\frac{\epsilon^{2}}{48}\left(\frac{\partial}{\partial z}\ln c_{0}\right)^{2}\frac{\partial c_{0}}{\partial z}\right]=\frac{\partial^{2}\overline{c}}{\partial z^{2}}.

(33)

This can also be written into a form analogous to Taylor dispersion as

\frac{\partial\overline{c}}{\partial t}=\frac{\partial}{\partial z}\left(D_{\textrm{eff}}\frac{\partial\overline{c}}{\partial z}\right),

(34)

where the effective diffusivity $D_{\textrm{eff}}$ to leading-order is given by

D_{\textrm{eff}}=1+\left(\frac{\Gamma_{w}}{D_{s}}\right)^{2}\frac{\epsilon^{2}}{48}\left(\frac{\partial}{\partial z}\ln c_{0}\right)^{2}+\mathcal{O}(\epsilon^{4}).

(35)

Note that the contribution of diffusioosmosis to the effective diffusivity in a 2D channel flow is a factor of 48/210 weaker than in a cylindrical pipe system. This is a consequence of the fact that for a given slip velocity at the walls, the centerline velocity in a cylindrical pipe must be greater than in a 2D channel flow, leading to greater velocity gradients, greater distortion of the solute profile, and a greater contribution to the effective diffusivity enhancement. Using these effective diffusivities along with the 1D diffusion equation, the evolution of the cross-sectionally averaged solute concentration is quite efficient to compute numerically.

III Numerical methods

To validate the theoretical results above, we performed numerical simulations for the coupled transport in both geometries. In this section, we describe the numerical methods used and show that the theoretical results above accurately match the numerical results. The velocity profiles are solved for the 2D channel flow case using a Fourier transform approach and for the cylindrical pipe flow case using a multigrid relaxation approach. Because we are dealing with the low Reynolds number regime and the fluid dynamics can be treated as quasi-steady, the numerical approach can be outlined as follows. First, we initialize the solute concentration to the previously-described Gaussian distribution. We then solve for the quasi-steady velocity profile using either a numerical multigrid relaxation approach or the theoretical Fourier transform approach. Next, using the determined velocity profiles, we update the solute concentration profile using a numerical solution of the governing advection-diffusion equation. This procedure is then repeated until the final time.

III.1 Velocity solver

For the case of incompressible Newtonian Stokes flow, the governing equations can be simplified to a biharmonic equation for the streamfunction, $\psi$ , which is given by

\nabla^{4}\psi=0.

(36)

This is a convenient formulation, as it allows for the solution of the velocity profiles without needing to solve for the pressure. In the Cartesian coordinate system, a Fourier transform approach can be used to greatly accelerate the numerical solution of this equation. In particular, the biharmonic equation can be Fourier transformed in the $x$ direction as

k^{4}\hat{\psi}(k,y)-2k^{2}\frac{\partial^{2}\hat{\psi}(k,y)}{\partial y^{2}}+\frac{\partial^{4}\hat{\psi}(k,y)}{\partial y^{4}}=0.

(37)

Here, we ignore the functional dependence of $\psi$ on time because the velocity can be treated as quasi-steady. The boundary conditions for (37) are $\hat{\psi}(k,0)=0$ , $\frac{\partial^{2}\hat{\psi}}{\partial y^{2}}\bigr{|}_{y=0}=0$ , $\hat{\psi}(k,\pm 1/2)=0$ , and $\frac{\partial\hat{\psi}}{\partial y}\bigr{|}_{y=\pm 1/2}=\hat{u}_{\textrm{wall}}$ . Here, $\hat{u}_{\textrm{wall}}$ can be found by using the Fast Fourier Transform fft() in MATLAB, and $\psi$ can be found by using the inverse Fast Fourier Transform of $\hat{\psi}$ using ifft(). Finally, the velocity components can be obtained directly from $\psi$ using

u=\frac{\partial\psi}{\partial y},~{}~{}~{}~{}~{}~{}v=-\frac{\partial\psi}{\partial x}.

(38)

For the case in cylindrical coordinates, we have not found a similar approach to use a Fourier transform to rapidly solve the biharmonic equation, so we instead do this directly using a multigrid solution to a set of coupled Poisson equations, and we then use the direct finite difference method to solve for the velocity. In particular, the biharmonic equation can be written as follows


$\displaystyle D^{2}\psi$	$\displaystyle=r\frac{\partial}{\partial r}\left(\frac{1}{r}\frac{\partial\psi}{\partial r}\right)+\frac{\partial^{2}\psi}{\partial z^{2}}=\phi,$	(39a)
$\displaystyle D^{2}\phi$	$\displaystyle=r\frac{\partial}{\partial r}\left(\frac{1}{r}\frac{\partial\phi}{\partial r}\right)+\frac{\partial^{2}\phi}{\partial z^{2}}=0,$	(39b)

where $D^{2}$ is a special operator in cylindrical coordinates that is useful for the biharmonic equation, which is given by $D^{2}=r\frac{\partial}{\partial r}\left(\frac{1}{r}\frac{\partial}{\partial r}\right)+\frac{\partial^{2}}{\partial z^{2}}$ [39, 11]. The wall has boundary conditions of $\frac{1}{r}\frac{\partial\psi}{\partial r}=u_{\text{wall}}$ and $\phi=-\frac{1}{r}\frac{\partial\psi}{\partial r}+\frac{\partial^{2}\psi}{\partial r^{2}}$ . All of the other boundaries have $\psi=0$ and $\phi=0$ because of symmetry conditions and no-flux conditions at the end of the channel.

This set of coupled Poisson equations is solved by using the Gauss-Seidel relaxation method with second-order accuracy in space and time for both governing equations and boundary conditions. As mentioned, the multigrid approach was used to rapidly accelerate the convergence of this iterative approach.

III.2 Concentration solver

The solute concentration profiles can be numerically solved using the advection-diffusion equation. We use finite difference with the approximate factorization method to solve the advection-diffusion equation for both geometries [27]. The semi-discretized governing equations applicable to the approximate factorization method for the Cartesian and cylindrical coordinate systems are given, respectively, by

\begin{split}\left[I-\frac{D_{s}\Delta t}{2}\delta_{xx}\right]\left[I-\frac{D_{s}\Delta t}{2}\delta_{yy}\right]c^{n+1}=-\frac{\Delta t}{2}\left(3u^{n}\frac{\partial c^{n}}{\partial x}+u^{n-1}\frac{\partial c^{n-1}}{\partial x}\right)\\ -\frac{\Delta t}{2}\left(3v^{n}\frac{\partial c^{n}}{\partial y}+v^{n-1}\frac{\partial c^{n-1}}{\partial y}\right)+\left[I+\frac{D_{s}\Delta t}{2}\delta_{xx}\right]\left[I+\frac{D_{s}\Delta t}{2}\delta_{yy}\right]c^{n},\end{split}

(40)

and

\begin{split}\left[I-\frac{D_{s}\Delta t}{2}\left(\delta_{rr}+\frac{1}{r}\delta_{r}\right)\right]\left[I-\frac{D_{s}\Delta t}{2}\delta_{zz}\right]c^{n+1}=-\frac{\Delta t}{2}\left(3u_{r}^{n}\frac{\partial c^{n}}{\partial r}-u_{r}^{n-1}\frac{\partial c^{n-1}}{\partial r}\right)\\ -\frac{\Delta t}{2}\left(3u_{z}^{n}\frac{\partial c^{n}}{\partial z}-u_{z}^{n-1}\frac{\partial c^{n-1}}{\partial z}\right)+\left[I+\frac{D_{s}\Delta t}{2}\left(\delta_{rr}+\frac{1}{r}\delta_{r}\right)\right]\left[I+\frac{D_{s}\Delta t}{2}\delta_{zz}\right]c^{n},\end{split}

(41)

where $\Delta t$ is the timestep, $n$ indicates the current time index, $I$ represents the identity matrix, and the $\delta_{xx}$ , etc. are the spatial derivative matrices. Here, implicit second-order accurate time-stepping is used for the diffusive terms, and second-order Adams-Bashforth is used for the advective terms. For the cylindrical case, second-order finite differencing was used for the spatial derivatives, whereas fourth-order finite differencing was uses for the Cartesian case. To start the simulation, we first take a series of very small timesteps using a first-order Euler method to calculate the advective terms. For each time step, the solute concentration profiles are updated in time by solving equations (40) or (41), and then the fluid velocities are recalculated as described above for use in the next timestep. Note that in all cases we add a small offset background concentration of $10^{-7}$ to prevent the so-called ‘ballistic motion’ described by [20] which represents the background ion concentration typically present in solution due to dissolved CO₂ or other factors. In this section, we have developed and presented the numerical methods used for both systems. In the following section, we explore the range of results and the physical evolution of such systems with numerical validation, and we show that the cross-sectionally averaged approach closely approximates the results of fully 2D simulations and can greatly simplify the analysis as in the case of Taylor dispersion.

IV Results

To begin, we first use the numerical simulations to validate the theoretical predictions. An example comparison between the numerical and theoretical results is shown in figure 2.

The analytical predictions of velocities in the parallel-plate channel are calculated by using equations (15) and (16), and those in the cylindrical channel are calculated from (54) and (55). Results are compared for $\Gamma_{w}/D_{s}=-1$ and $\epsilon=0.1$ at $t=0.2$ . With a negative diffusioosmotic mobility, the wall slip velocity is away from the peak solute concentration, driving flow away from the centerplane ( $x=0$ or $z=0$ ) at the walls and toward the centerplane along the channel centerline ( $y=0$ or $r=0$ ). The theoretical and numerical results agree well. One feature of the results to notice is that the velocity along the centerline in the cylindrical case is enhanced relative to the Cartesian case. We will see later how this alters the effective disperion in such configurations.

Figure 3 shows the comparison between (a,c) numerical simulations and (b,d) theoretical predictions of the higher-order solute concentration in the parallel-plate (a,b) and the cylindrical (c,d) channels, respectively. Here, the results are plotted over one quarter of the domain due to symmetry. The comparison is again calculated with $\Gamma_{w}/D_{s}=-1$ and $\epsilon=0.1$ , and the results demonstrate that the numerical results closely match the theoretical predictions.

Having validated the theoretical solutions using numerical simulations, we now provide a detailed examination of the diffusioosmotic dispersion process. We first consider the early-time dynamics of the system over which the initial deviation of the solute profile from the purely 1D dynamics forms. Figure 4 illustrates these early time dynamics by presenting visualizations of both $\hat{c}_{1}$ , $c_{1}^{\infty}$ , and $c_{1}$ in the early-time regime for times up to $t=1\times 10^{-3}$ with $\epsilon=0.1$ and $\Gamma_{w}/D_{s}=-1$ .

Recall that the fast timescale corresponds to the characteristic time for diffusion to occur across the channel and is represented by $\hat{c}_{1}$ from equation (25). In contrast, the slow timescale corresponds to the diffusion along the channel and is represented by $c_{1}^{\infty}$ from equation (21). Here, $c_{1}$ is the total deviation of the solute concentration profile from the purely 1D dynamics and is formed by the sum of both $c_{1}^{\infty}$ and $\hat{c_{1}}$ . As can be seen, the purpose of $\hat{c}_{1}$ is to cancel out the initial condition of $c_{1}^{\infty}$ such that the initial condition of $c_{1}$ can be zero. Then, in this example, $\hat{c}_{1}$ has almost entirely decayed by $t=10^{-3}$ after which the solution is dominated by the slow-time dynamics.

The early-time dynamics are expected to decay over the timescale $\epsilon^{2}$ . This can be verified by plotting the peak value of $\hat{c}_{1}$ over a range of mobilities and $\epsilon$ values, which is shown in figure 5. Solid dots correspond to the results of 2D simulations which are calculated as $\hat{c}_{\textrm{num}}=(c-c_{0})/\epsilon^{2}-c_{1,\textrm{theory}}^{\infty}$ .

Specifically, figure 5 $(a)$ shows the peak of $\hat{c_{1}}$ with $\epsilon=0.1$ at various $\Gamma_{w}/D_{s}$ values. As can be seen, higher $\Gamma_{w}/D_{s}$ values correspond to larger peak $\hat{c}_{1}$ values, reflecting the enhanced dispersion in those cases. For all $\Gamma_{w}/D_{s}$ , the peak $\hat{c}_{1}$ values have all apparently decayed before $t$ reaches $\epsilon^{2}=10^{-2}$ . Figure 5(b) extends these results by considering cases with different $\epsilon$ values at a constant $\Gamma_{w}/D_{s}=-1$ . Here, the dashed lines are the locations where $t=\epsilon^{2}$ . As expected, in every case the peak $\hat{c}_{1}$ value vanishes over the timescale $t=O(\epsilon^{2})$ as predicted by the theory.

In §II, we developed a theoretical model for the higher-order correction to the solute concentration profile due to diffusioosmosis. As discussed, the diffusioosmotic slip flow drives a recirculating flow in the channel that alters the transport of the solute. The vorticity and recirculating fluid flow in the channels are shown in figure 6 at $t=1$ for both coordinate systems.

Here, (a) and (c) show the nondimensional vorticity on the cross-section for the parallel-plate and cylindrical channels, respectively, with $\Gamma_{w}/D_{s}=-1$ . Streamlines illustrating the recirculation on the cross-section are superposed on top of the color map. Panels (b) and (d) show the same nondimensional vorticity data but with scaled velocity vector maps superposed instead. With, $\Gamma_{w}/D_{s}=-1$ , the diffusioosmosis at the channel walls drives a slip flow away from the centerplane, driving a recirculating flow that is towards the centerplane along the channel centerline. The flow directions and the signs of the vorticity will be reveresed in cases with positive mobilities.

Next, we visualize the higher-order solute dynamics on the cross-section to better understand the role of the diffusioosmotic dispersion on the solute transport. Figure 7 shows the theoretical values of $c_{1}$ for both parallel-plate (a) and cylindrical (b) channels as functions of $\Gamma_{w}/D_{s}$ at a fixed time of $t=1$ .

In order to interpret the figure, recall that negative $\Gamma_{w}$ values correspond to slip flow away from the centerplane along the walls and toward the centerplane along the centerline of the channel, while positive $\Gamma_{w}$ values correspond to flows that recirculate in the opposite direction. Regions of positive $c_{1}$ (red) indicate locations that have increased solute concentration relative to the 1D pure diffusion case ( $c_{0}$ ), and regions of negative $c_{1}$ (blue) have relatively less concentration relative to $c_{0}$ . We can interpret the formation of these regions as follows. For example, consider the case with $\Gamma_{w}/D_{s}=-1$ in the parallel-plate channel. Along the channel walls, the flow is away from the centerplane $x=0$ , pulling the relatively higher concentration fluid away from $x=0$ and enhancing the solute concentration somewhat away from the centerplane. This is enhanced by the recirculating nature of the flow that also pulls flow away from the centerline of the channel and towards the walls. Along the centerline, the flow is towards the centerplane, pulling relatively lower concentration fluid toward the centerplane and resulting in a depletion region. These effects are flipped in the case of positive $\Gamma_{w}/D_{s}$ . One last point to note is that the cylindrical case has relatively greater solute depletion and enhancement along the channel centerline due to the relatively greater centerline velocity in a cylindrical pipe compared to a 2D channel flow for the same wall slip velocity.

Next, we visualize the long-time behavior of the solute concentration profile as modeled by $c_{1}^{\infty}$ . Recall that for $t>\mathcal{O}(\epsilon^{2})$ , the higher-order solute profile is $c_{1}\approx c_{1}^{\infty}$ , as $\hat{c}_{1}$ has already decayed. The long-time results are shown in figure 8 for times up to $t=1000$ with $\Gamma_{w}/D_{s}=-1$ . Here, the horizontal axis is scaled by $\sqrt{1+4t}$ , since this represents the rate of spread of $c_{0}$ along the channel.

Recall that with our nondimensionalization, this corresponds to 1000 times the characteristic time for diffusion along the channel, such that the initial pulse of solute has well decayed by this time. As shown in the figure, over long times the higher-order solute profile also spreads significantly in the axial direction, retaining qualitatively the same shape, and ultimately decaying. As the initial pulse of solute decays, the concentration gradient at the walls likewise decays such that the diffusioosmosis and recirculating flow also decay with time leading to the ultimate decay of $c_{1}$ . Here, panel (a) corresponds to the 2D channel flow case, and panel (b) corresponds to the axisymmetric pipe flow case. The only significant notable difference between the two cases is that the magnitude of the higher-order solute concentration is a factor of 4-5 higher for the cylindrical case. As mentioned, this is due to the relatively greater centerline velocity in the axisymmetric geometry, which yields greater velocity gradients and enhanced dispersion.

Before proceeding to investigate the cross-sectionally averaged dynamics, we consider the effect of varying $\epsilon$ and the breakdown of the theoretical solution at large $\epsilon$ . Figure 9 shows the theoretical and numerical predictions of the higher-order solute concentration profiles in the parallel-plate channel with $\Gamma_{w}/D_{s}=-1$ at $t=1$ .

As shown above, the theoretical predictions show good agreement with the numerical results in the limit of small $\epsilon$ . Here, we see that reasonable quantitative agreement between the theory and simulations exists up to about $\epsilon=2$ , above which the theoretical predictions break down. In particular, as can be seen for $\epsilon=10$ in Figure 9, the numerical results show a much more uniform depletion along the channel centerline near $x=0$ compared to the theoretical results.

Finally, following the strategy of Taylor and Aris, we consider the dynamics of the cross-sectionally averaged concentration profile. In $\S$ II.4, we derived a cross-sectionally averaged concentration equation that can be used to model the net effects of diffusioosmotic dispersion. These equations are relatively simpler 1D diffusion equations with an effective diffusion coefficient that depends on the channel geometry and captures the effects of the dispersion. Unlike in the relatively simpler case of pure Taylor dispersion, here, the effective diffusivity is no longer a constant but rather a function of both axial position in the channel and time, as the net effects of the dispersion evolve with time and space. Thus, several methods are available for characterizing the cross-sectionally averaged solute dynamics. First, this 1D variable diffusivity model can be easily numerically integrated in time to yield the dynamics. Second, the results of the 2D numerical simulations can be averaged over the cross-section. Finally, the theoretical results for $c=c_{0}+\epsilon^{2}c_{1}$ can be averaged over the cross-section to yield the theoretical leading-order dynamics. All of these approaches should be expected to match up to $\mathcal{O}(\epsilon^{2})$ . A visualization of the cross-sectionally averaged solute dynamics is shown in figure 10 for the parallel-plate channel. Here, the results are the deviation of the cross-sectionally averaged solute concentration from the 1D results i.e., $\overline{c}-c_{0}$ . Figure 10(a) shows results as a function of $|\Gamma_{w}/D_{s}|$ with $\epsilon=0.1$ and $t=1$ . The solid lines indicate the theoretical predictions. The square markers correspond to the results of numerically solving the 1D forced diffusion equation given by (30). The star markers represent the results of averaging the 2D numerical simulation results over the cross-section. All three methods show close agreement. Recall that in the effective diffusivity coefficients derived in $\S$ II.4, the contribution from diffusioosmosis is $\mathcal{O}(\Gamma_{w}/D_{s})^{2}$ , such that the sign of the mobility does not affect the cross-sectionally averaged evolution.

However, increasing the magnitude of $\Gamma_{w}/D_{s}$ enhances the diffusioosmosis relative to the solute diffusion, leading to greater dispersion and greater deviations from $c_{0}$ . Figure 10(b) shows the cross-sectionally averaged solute concentration as a function of $\epsilon$ with $|\Gamma_{w}/D_{s}|=1$ and $t=1$ . The solid lines indicate the numerical simulations, and the square markers correspond to the theoretical predictions. As the value of $\epsilon$ increases, the magnitude of diffusioosmotic dispersion increases, enhancing the deviations from the pure diffusion case. The theoretical predictions show a good agreement with the simulations up to $\epsilon=2$ . Figures 10(c,d) present the cross-sectionally averaged concentration profile as a function of time. Figure 10(c) shows the theoretical (solid lines) and numerical (star symbols) predictions of the cross-sectionally averaged solute concentration for constant $\Gamma_{w}/D_{s}=-1$ and $\epsilon=0.1$ and (d) shows the cross-sectionally averaged solute concentration from numerical simulations for constant $\Gamma_{w}/D_{s}=-1$ and $\epsilon=10$ . As time progresses, a relative depletion zone forms near the channel center that is balanced by accumulation regions to the left and right of the solute peak. Note that the depletion at $x=0$ does not form quite as quickly as that near $x=0$ , resulting in a small bump that decays with time. This is due to the fact that the axial concentration gradient is zero at $x=0$ due to symmetry, such that the wall slip velocity and any recirculating flow is negligible very near the centerplane. Thus, some time is still required for the solute to diffuse and adjust.

In addition, we study the width of the spread of the cross-sectionally averaged concentration profiles. Here, we define the width of the solute distribution as $\mathcal{L}_{99}$ , which corresponds to the axial position ( $x$ or $z$ ), where the concentration has decayed by 99% of its current peak value. Numerical results for a range of $\epsilon$ values are shown in figure 11 and compared to the theoretical predictions. Panel (a) shows the rescaled width of the concentration distributions as functions of time with constant $\Gamma_{w}/D_{s}=-1$ for different values of $\epsilon$ . The colored markers represent the numerical results of solute concentration, and the solid and light blue dashed lines correspond to theoretical predictions of the distribution width of $c0+\epsilon^{2}c_{1}^{\infty}$ and $c_{0}$ , respectively. Here, the distribution widths are rescaled by $\sqrt{1+4t}$ , which represents the expected spreading rate of $c_{0}$ . As can be seen, the results tend toward constant values, indicating that the spread of the higher-order solute profile is set by the spread of $c_{0}$ . This is due to the fact that while $c_{0}$ becomes small at the edge of the distribution, the gradient of the logarithm of $c_{0}$ may remain large, allowing the diffusioosmotic dispersion to act over a relatively wider distribution. As $\epsilon$ increases, the width of the higher-order solute distribution increases. Figure 11(b) shows $\mathcal{L}_{99}$ as a function of $\epsilon$ with constant $\Gamma_{w}/D_{s}=-1$ at a fixed time of $t=1$ . Here, the theoretical predictions of $c_{0}+\epsilon^{2}c_{1}^{\infty}$ and $c_{0}$ are shown as a solid line and a dashed line, respectively. The symbol markers correspond to the numerical results. At a given time, the width of the distribution increases with a larger $\epsilon$ value. The theoretical predictions show a good agreement with numerical simulation for $\epsilon<1$ , and the error between simulations and theory for this metric remains less than 10% even at $\epsilon=10$ .

V Conclusion

In this study, we investigated the diffusioosmotic dispersion in a long, narrow channel with an initial Gaussian plug of solute at the center of the channel. The concentration gradients associated with the diffusion of this solute induce diffusioosmotic slip flow at the channel walls. The slip flow, in turn, drives recirculation inside the channel so that the transport of the solute is not purely diffusive, but also advective. The recirculating fluid flow in the channel introduces a shear flow that distorts the solute concentration profile and alters the effective transport of it analogously to the process of Taylor dispersion. We derived theoretical solutions for the coupled fluid and solute dynamics in such systems for both 2D channel flows and axisymmetric pipe flows. The theoretical derivation utilized a multiple-timescale analysis that incorporates both the fast and slow time dynamics. By averaging the governing equations over the cross section and incorporating the leading-order solutions, the system was reduced to a 1D diffusion equation with a variable effective diffusivity coefficient that depends on the geometry of the system and incorporates the effects of diffusioosmotic dispersion in an averaged sense. This approach is analogous to that of Taylor and Aris in the study of Taylor dispersion. Numerical solutions of this effective 1D model were compared with cross-sectionally averaged results of the 2D numerical simulations as well as the theoretical averaged results, all of which show good agreement. As much as possible, we have kept this analysis independent of the specific initial condition, except for the limitation that it be uniform over the cross section and have a characteristic distribution that is reasonably spread out relative to the channel width. The specific Gaussian initial condition chosen here was only needed for the calculation of $B(x,t)$ due to the Fourier transform approach needed to solve for it. Thus, while the full derivation of the effective diffusivity coefficients is not completely general, it can easily be adapted to other systems and initial conditions by modifying this one section of the calculation. It may also be possible to seek an explicit solution to $B(x,t)$ that remains general of the initial solute concentration, but this is a problem that we leave for future work. The results and analysis presented here have documented the role of diffusioosmosis-driven recirculation on the transport of solute in a straight channel or pipe flow and have shown how the cross-sectionally averaged effective dynamics can be reduced to a 1D effective diffusion problem analogously to Taylor dispersion.

Appendix A Numerical details

In this section, we will discuss the numerical convergence tests. In the numerical simulations for the Cartesian case, we implemented a numerical scheme that is second-order accurate in time and fourth-order accurate in space. For the cylindrical case, we implemented a second-order accurate scheme in both time and space. Figure 12 shows the numerical convergence study results with $\Gamma_{w}/D_{s}=-1$ . The solid lines are the best-fit power law curves, and the matching color equations are the corresponding fitted power law functions. Figure 12 $a$ shows the relative norm error between the theoretical and numerical prediction of $c_{0}+\epsilon^{2}c_{1}$ at $t=0.01$ as we increase the value of $\epsilon$ . The theoretical prediction of solute concentration $c_{0}+\epsilon^{2}c_{1}$ has correction terms due to diffusioosmosis, which are order $\epsilon^{4}$ . The theoretical predictions converge to $c_{0}$ as $\epsilon$ goes to zero, which can also be observed in figure 12 $a$ . Figure 12 $b$ shows the order of convergence in the time step, $dt$ , for both Cartesian and cylindrical cases with $\epsilon=0.1$ . Here, we fixed the cylindrical case grid size as $2049\times 1025$ and the Cartesian case grid size as $1600\times 200$ . In both cases, the simulation ran until $t=0.1$ . The relative norm error is calculated in reference to the smallest $dt$ case in the simulation. Both Cartesian and cylindrical cases are second-order accurate in time as shown in the best-fit curve in figure 12 $b$ , where the absolute norm error is shown to decrease as $O(\textrm{dt}^{2})$ , as expected. We chose $\textrm{dt}=10^{-4}$ for both cases as the relative error is less than $10^{-4}$ . Figure 12 $c$ and 12 $d$ are the spatial convergence studies with $\epsilon=0.1$ for the Cartesian and cylindrical cases, respectively. Here, we used $\textrm{dt}=1\times 10^{-5}$ and $t_{\textrm{final}}=0.1$ for both Cartesian and cylindrical spatial convergence studies. The relative norm error is calculated in reference to the finest grid in the simulation. As shown in figures 12 $c$ and 12 $d$ , the errors are $O(\textrm{dx}^{4})$ in the Cartesian case and $O(\textrm{dz}^{2})$ in the cylindrical case.

Appendix B Derivation of cylindrical channel

In this section, we demonstrate the theoretical solution to the fluid and solute dynamics in the axisymmetric cylindrical coordinate system. The general procedure is the same as that for the 2D Cartesian coordinate system. The governing equations for the system are provided by (1). Here, we consider an axisymmetric configuration. We introduce the following nondimensionalizations:

z=\frac{z^{*}}{\ell},\;r=\frac{r^{*}}{a},\;u_{z}=\frac{u_{z}^{*}}{U},\;u_{r}=\frac{u_{r}^{*}\ell}{Ua},\;P=\frac{P^{*}a^{2}}{\mu U\ell},\;\epsilon=\frac{a}{\ell},\;\;U=\frac{D_{s}}{\ell},\;t=\frac{t^{*}}{\ell^{2}/D_{s}}.

(42)

With these scalings the nondimensional form of the governing equations in cylindrical coordinates are given by


	$\displaystyle\frac{1}{r}\frac{\partial(ru_{r})}{\partial r}+\frac{\partial u_{z}}{\partial z}=0,$		(43a)
	$\displaystyle-\frac{\partial p}{\partial r}+\epsilon^{2}\frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial u_{r}}{\partial r}\right)-\epsilon^{2}\frac{u_{r}}{r^{2}}+\epsilon^{4}\frac{\partial^{2}u_{r}}{\partial z^{2}}=0,$		(43b)
	$\displaystyle-\frac{\partial p}{\partial z}+\frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial u_{z}}{\partial r}\right)+\epsilon^{2}\frac{\partial^{2}u_{z}}{\partial z^{2}}=0,$		(43c)
	$\displaystyle\epsilon^{2}\frac{\partial c}{\partial t}+\epsilon^{2}u_{r}\frac{\partial c}{\partial r}+\epsilon^{2}u_{z}\frac{\partial c}{\partial z}=\frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial c}{\partial r}\right)+\epsilon^{2}\frac{\partial^{2}c}{\partial z^{2}}.$		(43d)

The solution to the governing equation (43) is subject to boundary conditions on the fluid and solute. These boundary conditions can be summarized by:

\text{Quiescent far-field conditions:}~{}p=0~{}~{}\text{and}~{}~{}~{}u_{z}=0~{}~{}~{}\text{and}~{}~{}~{}c=0~{}~{}~{}\text{at}~{}z=\pm\infty,

(44)

\text{No fluid penetration at the walls: }u_{r}=0~{}\text{at}~{}r=1,

(45)

\text{No-flux conditions at the channel walls: }\frac{\partial c}{\partial r}=0~{}\text{at}~{}r=1,

(46)

\text{Diffusioosmotic wall slip boundary condition: }u_{z}=\frac{\Gamma_{w}}{D_{s}}\frac{\partial\ln c}{\partial z}~{}\text{at}~{}r=1.

(47)

Similar to the Cartesian case, we introduce a multiple timescale approach [9] in which we introduce a fast time variable $T=t/\epsilon^{2}$ . That is, $T=O(1)$ over dimensional times $\sim a^{2}/D_{s}$ corresponding to $t=O(\epsilon^{2})$ , whereas $t=O(1)$ on dimensional times $~{}\ell^{2}/D_{s}$ . Thus, we rewrite the advection-diffusion equation as,

\epsilon^{2}\frac{\partial c}{\partial t}+\frac{\partial c}{\partial T}+\epsilon^{2}u_{r}\frac{\partial c}{\partial r}+\epsilon^{2}u_{z}\frac{\partial c}{\partial z}=\frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial c}{\partial r}\right)+\epsilon^{2}\frac{\partial^{2}c}{\partial z^{2}}.

(48)

Again, similar to the approach used in Cartesian coordinates, we seek the analytical solution of the governing equations as perturbation expansions with the small parameter $\epsilon=a/\ell$ using an expansion that has the form


$\displaystyle c(r,z,t,T)$	$\displaystyle=c_{0}(r,z,t,T)+\epsilon^{2}c_{1}(r,z,t,T)+\epsilon^{4}c_{2}(r,z,t,T)+\ldots,$	(49a)
$\displaystyle p(r,z,t,T)$	$\displaystyle=p_{0}(r,z,t,T)+\epsilon^{2}p_{1}(r,z,t,T)+\epsilon^{4}p_{2}(r,z,t,T)+\ldots,$	(49b)
$\displaystyle u_{z}(r,z,t,T)$	$\displaystyle=u_{z0}(r,z,t,T)+\epsilon^{2}u_{z1}(r,z,t,T)+\epsilon^{4}u_{z2}(r,z,t,T)+\ldots,$	(49c)
$\displaystyle u_{r}(r,z,t,T)$	$\displaystyle=u_{r0}(r,z,t,T)+\epsilon^{2}u_{r1}(r,z,t,T)+\epsilon^{4}u_{r2}(r,z,t,T)+\ldots.$	(49d)

We first need to obtain the leading-order velocity and pressure solutions. These can be obtained by substituting equations (49) into the governing equations, which gives


	$\displaystyle\frac{\partial u_{z0}}{\partial z}+\frac{u_{r0}}{r}+\frac{\partial u_{r0}}{\partial r}=0,$		(50a)
	$\displaystyle\frac{\partial p_{0}}{\partial r}=0,$		(50b)
	$\displaystyle-\frac{\partial p_{0}}{\partial z}+\frac{1}{r}\frac{\partial u_{z0}}{\partial r}+\frac{\partial^{2}u_{z0}}{\partial r^{2}}=0,$		(50c)
	$\displaystyle\frac{\partial c_{0}}{\partial T}=\frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial c_{0}}{\partial r}\right),$		(50d)

to leading order. Considering the fact that the initial condition is given by $c_{0}(T=0,t=0)=\exp(-z^{2})$ and that there are no-flux conditions at the channel walls, i.e., $\frac{\partial c_{0}}{\partial r}(r=\pm 1)=0$ , it must be true that $c_{0}(r,z,t,T)=c_{0}(z,t)$ , with $c_{0}(t=0)=\exp(-x^{2})$ .

To solve for the leading-order velocities, we first take equation (50 $a$ ) and solve for $u_{z0}$ , which is subject to the slip boundary condition (47). The leading-order $u_{z0}$ becomes

u_{z0}=\frac{\Gamma_{w}}{D_{s}c_{0}}\frac{\partial c_{0}}{\partial z}+\frac{1}{4}(-1+r^{2})\frac{\partial p_{0}}{\partial z}.

(51)

Here, $u_{z0}$ has a term that includes $p_{0}(z,t)$ , which can be found by considering the conservation of mass and integrating over the channel cross-section,

\int_{0}^{1}u_{z0}2\pi r\,dr=0.

(52)

Following this approach, $p_{0}(z,t)$ , is found to be

p_{0}=8\frac{\Gamma_{w}}{D_{s}}\ln(c_{0}).

(53)

With equation (53), $u_{z0}(r,z,t)$ can be simplified and written as

u_{z0}=\frac{\Gamma_{w}}{D_{s}c_{0}}(-1+2r^{2})\frac{\partial c_{0}}{\partial z}.

(54)

To solve for $u_{r0}(r,z,t)$ , we integrate the continuity equation (50a) and apply the boundary condition given by (5). This gives the leading-order $u_{r0}$ to be

u_{r0}=\frac{\Gamma_{w}r(-1+r^{2})}{D_{s}2c_{0}^{2}}\left(\left(\frac{\partial c_{0}}{\partial z}\right)^{2}-c_{0}\frac{\partial^{2}c_{0}}{\partial z^{2}}\right).

(55)

Substituting the asymptotic expansion into the advection-diffusion equation, we find, to leading-order, that

\frac{\partial c_{0}}{\partial t}+\frac{\partial c_{1}}{\partial T}+u_{z0}\frac{\partial c_{0}}{\partial z}-\frac{\partial^{2}c_{0}}{\partial z^{2}}-\frac{1}{r}\frac{\partial c_{1}}{\partial r}-\frac{\partial^{2}c_{1}}{\partial r^{2}}=0,

(56)

subject to the initial condition $c_{0}(t=0)=\exp(-z^{2})$ and $c_{1}(t=T=0)=0$ . To solve this, we note that at long times the fast-time dynamics should have all decayed such that the time derivative term can be ignored, and the equation can be averaged across the channel to obtain

\frac{\partial c_{0}}{\partial t}=\frac{\partial^{2}c_{0}}{\partial z^{2}}.

(57)

As in the Cartesian case, the solution to this is

c_{0}(z,t)=\frac{1}{\sqrt{1+4t}}\exp{\left(-\frac{z^{2}}{1+4t}\right)}.\\

(58)

Substituting this into (56) gives

\frac{\partial c_{0}}{\partial t}+\frac{\Gamma_{w}}{D_{s}c_{0}}(-1+2r^{2})\left(\frac{\partial c_{0}}{\partial z}\right)^{2}+\frac{\partial c_{1}}{\partial T}-\frac{\partial^{2}c_{0}}{\partial z^{2}}-\frac{1}{r}\frac{\partial c_{1}}{\partial r}-\frac{\partial^{2}c_{1}}{\partial r^{2}}=0.

(59)

At long times, the fast time dynamics have decayed such that derivatives with respect to $T$ can be neglected, and the equation can be integrated to yield

c_{1}(T\rightarrow\infty)\sim c_{1}^{\infty}(r,z,t)=\frac{1}{8}r^{2}\left[2\frac{\partial c_{0}}{\partial t}+\frac{\Gamma_{w}}{D_{s}c_{0}}(-2+r^{2})\left(\frac{\partial c_{0}}{\partial z}\right)^{2}-2\frac{\partial^{2}c_{0}}{\partial z^{2}}\right]+B(z,t),

(60)

where $B(z,t)$ is a yet unknown function obtained from integration. This can be determined by applying conservation of mass and taking the cross-sectional average of the advection-diffusion equation, which gives

\frac{\partial B}{\partial t}=\frac{\partial^{2}B}{\partial z^{2}}-\frac{e^{\frac{-z^{2}}{1+4t}}\Gamma_{w}/D_{s}(2+32t^{2}+4t(4+\Gamma_{w}/D_{s}z^{2})+\Gamma_{w}/D_{s}(z^{2}-z^{4}))}{3(1+4t)^{9/2}}.

(61)

Using a Fourier transform approach, the solution can be found to be

B=\frac{2\Gamma_{w}/D_{s}\left[16z\alpha^{2}+4\Gamma_{w}/D_{s}t(3\alpha^{2}-12\alpha z^{2}+4z^{4})-3\Gamma_{w}/D_{s}\alpha^{2}(\alpha-2z^{2})\ln\alpha\right]}{e^{\frac{z^{2}}{\alpha}}(\alpha^{7/2}(96+384t))}.

(62)

This equation does not satisfy the initial condition, so it is not yet the full solution for the higher-order solute dynamics. Now, we look for a solution $c_{1}(r,z,t,T)=c_{1}^{\infty}(r,z,t)+\hat{c}_{1}(r,z,t,T)$ . Substituting into equation (48), we find

\frac{\partial\hat{c_{1}}}{\partial T}=\frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial\hat{c}_{1}}{\partial r}\right),

(63)

with the initial condition given by

\hat{c}_{1}(T=0)=-c_{1}^{\infty}(t=0)=-\frac{1}{6}e^{-z^{2}}\frac{\Gamma_{w}}{D_{s}}(2-6r^{2}+3r^{4})z^{2}.

(64)

The solution to this is in the form of a Bessel series $\hat{c}_{1}(r,T)=\sum_{n=0}^{\infty}a_{n}e^{-\lambda_{n}^{2}T}J_{0}(\lambda_{n}r)$ , with initial condition of $f(r)=-\frac{1}{6}e^{-z^{2}}\frac{\Gamma_{w}}{D_{s}}(2-6r^{2}+3r^{4})z^{2}$ and boundary condition of $\frac{\partial\hat{c_{1}}}{\partial r}=0$ at $r=1$ . Using the boundary condition, the $\lambda_{n}$ can be found to be the roots of $J_{1}$ , which we denote as $j_{1,n}$ . The coefficients $a_{n}$ can be determined as follows

\begin{split}a_{n}&=\frac{2}{J_{0}(j_{1,n})^{2}}\int_{0}^{1}rf(r)J_{0}(j_{1,n}r)\,dr\\ &=-\frac{e^{-z^{2}}\frac{\Gamma_{w}}{D_{s}}z^{2}(96J_{2}(j_{1,n})-J_{1}(j_{1,n})j_{1,n}(24+j_{1,n}^{2}))}{3J_{0}(j_{1,n})^{2}j_{1,n}^{4}}.\end{split}

(65)

Thus, the solution of $\hat{c}_{1}$ is as follows

\hat{c}_{1}(r,z,t,T)=2\int_{0}^{1}rf(r)\,dr+\sum_{n=1}^{\infty}a_{n}e^{-j_{1,n}^{2}T}J_{0}(j_{1,n}r).

(66)

We can then construct a composite solution $c_{1}=c_{1}^{\infty}+\hat{c}_{1}$ that is valid for all $t$ using the fact that $T=t/\epsilon^{2}$ .

References

Alessio et al. [2021] B. M. Alessio, S. Shim, E. Mintah, A. Gupta, and H. A. Stone. Diffusiophoresis and diffusioosmosis in tandem: Two-dimensional particle motion in the presence of multiple electrolytes. Phys. Rev. Fluids, 6(5):054201, 2021.
Alessio et al. [2022] B. M. Alessio, S. Shim, A. Gupta, and H. A. Stone. Diffusioosmosis-driven dispersion of colloids: a Taylor dispersion analysis with experimental validation. J. Fluid Mech., 942, 2022.
Aminian et al. [2016] M. Aminian, F. Bernardi, R. Camassa, D. M. Harris, and R. M. McLaughlin. How boundaries shape chemical delivery in microfluidics. Science, 354(6317):1252–1256, 2016.
Anderson and Prieve [1984] J. L. Anderson and D. C. Prieve. Diffusiophoresis: Migration of colloidal particles in gradients of solute concentration. Separ. Purif. Method., 13(1):67–103, 1984.
Aris [1956] R. Aris. On the dispersion of a solute in a fluid flowing through a tube. Proc. R. Soc. Lond., 235(1200):67–77, 1956.
Ault et al. [2017] J. T. Ault, P. B. Warren, S. Shin, and H. A. Stone. Diffusiophoresis in one-dimensional solute gradients. Soft Matt., 13(47):9015–9023, 2017.
Ault et al. [2019] J. T. Ault, S. Shin, and H. A. Stone. Characterization of surface–solute interactions by diffusioosmosis. Soft Matt., 15(7):1582–1596, 2019.
Barton [1983] N.G. Barton. On the method of moments for solute dispersion. J. Fluid Mech., 126:205–218, 1983.
Bender and Orszag [1999] C. M. Bender and S Orszag. Advanced mathematical methods for scientists and engineers I: Asymptotic methods and perturbation theory, volume 1. Springer Science & Business Media, 1999.
Bone et al. [2020] S. E. Bone, H. G. Steinrück, and M. F. Toney. Advanced characterization in clean water technologies. Joule, 4(8):1637–1659, 2020.
Brenner [1961] H. Brenner. The slow motion of a sphere through a viscous fluid towards a plane surface. Chem. Eng. Sci., 16(3-4):242–251, 1961.
Brenner and Edwards [1993] H. Brenner and D. A. Edwards. Macrotransport Processes. Butterworth-Heinemann, 1993.
Chakra et al. [2023] Adnan Chakra, Naval Singh, Goran T Vladisavljević, François Nadal, Cécile Cottin-Bizonne, Christophe Pirat, and Guido Bolognesi. Continuous manipulation and characterization of colloidal beads and liposomes via diffusiophoresis in single-and double-junction microchannels. arXiv preprint arXiv:2302.05800, 2023.
Chu et al. [2019] H. C. W. Chu, S. Garoff, T. M. Przybycien, R. D. Tilton, and A. S. Khair. Dispersion in steady and time-oscillatory two-dimensional flows through a parallel-plate channel. Phys. Fluids, 31(2):022007, 2019.
Chu et al. [2021] H. C.W. Chu, S. Garoff, R. D. Tilton, and A. S. Khair. Macrotransport theory for diffusiophoretic colloids and chemotactic microorganisms. J. Fluid Mech., 917, 2021.
Derjaguin et al. [1961] B. V. Derjaguin, S. S. Dukhin, and A. A. Korotkova. Diffusiophoresis in electrolyte solutions and its role in mechanism of film formation from rubber latexes by method of ionic deposition. Kolloid. Zhur., 23(1):53, 1961.
Florea et al. [2014] D. Florea, S. Musa, J. Huyghe, and H. M. Wyss. Long-range repulsion of colloids driven by ion exchange and diffusiophoresis. Proc. Natl Acad. Sci., 111(18):6554–6559, 2014.
Frankel and Brenner [1989] I. Frankel and H. Brenner. On the foundations of generalized taylor dispersion theory. J. Fluid Mech., 204:97–119, 1989.
Gupta et al. [2020a] A. Gupta, A. G. Rajan, E. A. Carter, and H. A. Stone. Ionic layering and overcharging in electrical double layers in a Poisson-Boltzmann model. Phys. Rev. Lett., 125(18):188004, 2020a.
Gupta et al. [2020b] A. Gupta, S. Shim, and H. A. Stone. Diffusiophoresis: From dilute to concentrated electrolytes. Soft Matt., 16(30):6975–6984, 2020b.
Gupta et al. [2020c] A. Gupta, P. J. Zuk, and H. A. Stone. Charging dynamics of overlapping double layers in a cylindrical nanopore. Phys. Rev. Lett., 125(7):076001, 2020c.
Kar et al. [2015] A. Kar, T. Chiang, I. Ortiz Rivera, A. Sen, and D. Velegol. Enhanced transport into and out of dead-end pores. ACS Nano, 9(1):746–753, 2015.
Keh and Ma [2005] H. J. Keh and H. C. Ma. Diffusioosmosis of electrolyte solutions along a charged plane wall. Langmuir, 21(12):5461–5467, 2005.
Lee et al. [2018] H. Lee, J. Kim, J. Yang, S. W. Seo, and S. J. Kim. Diffusiophoretic exclusion of colloidal particles for continuous water purification. Lab Chip, 18(12):1713–1724, 2018.
Lee et al. [2022] S. Lee, J. Lee, and J. T. Ault. The role of variable zeta potential on diffusiophoretic and diffusioosmotic transport. Colloids Surfaces A, page 130775, 2022. ISSN 0927-7757.
Migacz and Ault [2022] R. E. Migacz and J. T. Ault. Diffusiophoresis in a Taylor-dispersing solute. Phys. Rev. Fluids, 7(3):034202, 2022.
Moin [2010] P. Moin. Fundamentals of Engineering Numerical Analysis. Cambridge University Press, 2010.
Palacci et al. [2010] J. Palacci, B. Abécassis, C. Cottin-Bizonne, C. Ybert, and L. Bocquet. Colloidal motility and pattern formation under rectified diffusiophoresis. Phys. Rev. Lett., 104(13):138302, 2010.
Park et al. [2021] S. W. Park, J. Lee, H. Yoon, and S. Shin. Microfluidic investigation of salinity-induced oil recovery in porous media during chemical flooding. Energy & Fuels, 35(6):4885–4892, 2021.
Rasmussen et al. [2020] M. K. Rasmussen, J. N. Pedersen, and R. Marie. Size and surface charge characterization of nanoparticles with a salt gradient. Nat. Commun., 11(1):1–8, 2020.
Salmon and Doumenc [2020] J.B. Salmon and F. Doumenc. Buoyancy-driven dispersion in confined drying of liquid binary mixtures. Phys. Rev. Fluids, 5(2):024201, 2020.
Shim [2022] S. Shim. Diffusiophoresis, diffusioosmosis, and microfluidics: Surface-flow-driven phenomena in the presence of flow. Chem. Rev., 122(7):6986–7009, 2022.
Shin et al. [2016] S. Shin, E. Um, B. Sabass, J. T. Ault, M. Rahimi, P. B. Warren, and H. A. Stone. Size-dependent control of colloid transport via solute gradients in dead-end channels. Proc. Natl Acad. Sci., 113(2):257–261, 2016.
Shin et al. [2017a] S. Shin, J. T. Ault, J. Feng, P. B. Warren, and H. A. Stone. Low-cost zeta potentiometry using solute gradients. Adv. Mat., 29(30):1701516, 2017a.
Shin et al. [2017b] S. Shin, J. T. Ault, P. B. Warren, and H. A. Stone. Accumulation of colloidal particles in flow junctions induced by fluid flow and diffusiophoresis. Phys. Rev. X, 7(4):041038, 2017b.
Shin et al. [2017c] S. Shin, O. Shardt, P. B. Warren, and H. A. Stone. Membraneless water filtration using CO₂. Nat. Commun., 8(1):1–6, 2017c.
Shin et al. [2018] S. Shin, P. B. Warren, and H. A. Stone. Cleaning by surfactant gradients: Particulate removal from porous materials and the significance of rinsing in laundry detergency. Phy. Rev. Appl., 9(3):034012, 2018.
Singh et al. [2020] N. Singh, G. T. Vladisavljević, F. Nadal, C. Cottin-Bizonne, C. Pirat, and G. Bolognesi. Reversible trapping of colloids in microgrooved channels via diffusiophoresis under steady-state solute gradients. Phys. Rev. Lett., 125(24):248002, 2020.
Stimson and Jeffery [1926] M. Stimson and G. B. Jeffery. The motion of two spheres in a viscous fluid. Proc. R. Soc. Lond., 111(757):110–116, 1926.
Stone and Brenner [1999] H. A. Stone and H. Brenner. Dispersion in flows with streamwise variations of mean velocity: Radial flow. Ind. Eng. Chem. Res., 38(3):851–854, 1999.
Taylor [1953] G. I. Taylor. Dispersion of soluble matter in solvent flowing slowly through a tube. Proc. R. Soc. Lond., 219(1137):186–203, 1953.
Williams et al. [2020] I. Williams, S. Lee, A. Apriceno, R. P. Sear, and G. Battaglia. Diffusioosmotic and convective flows induced by a nonelectrolyte concentration gradient. Proc. Natl Acad. Sci., 117(41):25263–25271, 2020.