Transport of a passive scalar in wide channels with surface topography: An asymptotic theory

Transport processes at small scales play a central role in a number of fields, ranging from biology [1, 2], where solutes and blood cells are transported through vascular channels or bacteria spread through soil and tissues, to geophysics [3], where water flows through complex porous materials and eroded rock particulates sediment in rivers, to microfluidic applications [4–8], where various biological and synthetic samples are mixed, sorted, or focused while moving through microchannels. These systems often operate at low Reynolds number, where viscous forces dominate over inertia and strong stochastic fluctuations due to Brownian effects dictate the overall transport behavior [9]. In contrast to dilute, unconfined media, natural environments and microfluidic applications exhibit a variety of confining surfaces with different properties, including deformable and elastic materials [10], fluctuating boundaries [11], and structured topographies [12], which result in complex hydrodynamic flows. Understanding the interaction of Brownian suspensions with nearby boundaries and flows provides insight into microbiological processes and enables optimization and design of microfluidic devices [5].

The study of solute transport in microchannels dates back to the work of Taylor [13] and Aris [14] in the 1950s. In their seminal work, they predicted an enhancement compared to diffusion alone of the streamwise spreading of particles, as diffusion across the channel allows particles to sample different velocities of the hydrodynamic flows. Since then, numerous theoretical and experimental studies have extended their ideas, treating the transient profile of the solute concentration [15–17], the effect of a finite particle size [18], particle–particle interactions [19], channel corrugations [20], channel wall slip [21], and absorbing and pulsating channel walls [22]. The latter can lead to a reduction of particle dispersion due to an entropic slow-down resulting from channel wall constrictions. Such a slow-down has also been reported for particles diffusing through periodic, corrugated channel geometries [23]. On the other hand, channel wall topography may enhance dispersion due to the generation of additional hydrodynamic flows and velocity gradients [11].

Beyond dispersion through channel geometries, solute transport in porous materials has received substantial scientific attention dating back to Saffman in the 1960s [24]. Continuing the work of several decades, in the 1970s De Josselin De Jong noted that dispersion in porous materials is a tensorial quantity [25], which encodes the (anisotropic) properties of the medium. In 1980 Brenner developed a systematic theory to predict the second-order anisotropic dispersion tensor and drift of particles spreading through spatially periodic, porous media [26], which relies on the details of the flow field within one periodic cell. Brenner's theory has been extended to account for fluctuations in the environment that are of the order of the transport process itself [27]. At the same time, other researchers characterized the tensorial nature of diffusion in random porous media [28, 29]. Building on Brenner's pioneering work, the study of self-propelled agents in spatially periodic media has demonstrated that obstacles act as entropic barriers while producing additional velocity gradients, which can simultaneously reduce or enhance dispersion, depending on the local shear rate [30]. Reminiscent of the transport through fluctuating channels mentioned above, experiments revealed that fluctuations of the obstacles enhance transport of particles diffusing through the porous matrix [31]. Going beyond periodic structures, other experiments have demonstrated anomalous dispersal of solutes in random porous media, where complex flow patterns of Newtonian and non-Newtonian fluids emerge within the dead end pores [32, 33]. Other recent work has applied similar ideas to describe the evolution of volatile substances in droplets, which then relates to their shape through the generation of Marangoni forces [34].

While these studies provide immediate insights into biological and geophysical transport phenomena, they can also guide the design of new microfluidic devices, which rely on the precise control of transport of particulate suspensions [4–8]. For example, the method 'deterministic lateral displacement' relies on the use of arrays of posts within a microfluidic channel to efficiently separate biological and synthetic constituents of different sizes [35–37]. The effect of diffusion within this context has been analyzed both experimentally [38] and theoretically [39] in the realm of Brenner's theory. In particular, it has been shown that the interactions of finite-sized particles with the anisotropic obstacle field can generate long-time anisotropic dispersion [38, 39].

Another common approach is the use of tailored surfaces with particular surface topographies to achieve mixing [40], sorting [41–44], or focusing of particles in suspension [45, 46]. In addition, recently, it has been found that particles moving past herringbone structures have complex, three-dimensional trajectories due to the particle's interaction with the surface [47, 48], yet the effect of diffusion, which in many contexts is not negligible, on their long-time transport properties has not yet been addressed. While in several studies dispersion in narrow channels has been studied, in many of these microfluidic applications channels are significantly wider than they are tall, requiring a new, fully three-dimensional theoretical approach for the characterization of dispersion, which requires consideration of transport both along and perpendicular to the flow. It is important to note that in the long time limit, the effects of the channel's side walls have been shown to dominate dispersion and can lead to an increase in dispersion up to a factor of eight [49–52]. However, for any channel there will always be an intermediate time scale where dispersion is controlled solely by the influence of the channel's top and bottom surfaces with the side wall effects only mattering on much longer time scales, which are increased as the channel is widened [49–52].

Here, we revisit the classical Taylor dispersion theory and extend it to scalar transport in wide, structured channels. By 'structured' we refer to the shape of the channel walls. We present an asymptotic long-time, two-dimensional convection–diffusion equation, in contrast to Taylor's one-dimensional equation for dispersion through narrow channels. In particular, the three-dimensional nature of the surface structures no longer allows for a reduction to either a two-dimensional or axisymmetric description of the flow. Our asymptotic theory, which relies on a domain-perturbation method, valid for small surface amplitudes, for the hydrodynamic flows, provides an analytic prediction for the dispersion matrix and the overall drift as a function of the surface shape. We provide results for different surface structures, ranging from corrugated channel walls, as often used in microfluidic applications, to randomly structured topographies. Finally, we use stochastic simulations with asymptotic predictions for the flow fields as input to corroborate our coarse-grained theory.

Consider a Brownian particle at a position r at time t. The probability density $c(\mathbf{r},t)$ of the particle (or, equivalently, the solute concentration) is governed by the Fokker–Planck (or convection–diffusion) equation,

$$\begin{align} \partial_t c = - \boldsymbol{\mathbf{\nabla}} \cdot (\mathbf{u} c) + D_0 \nabla^2 c, \end{align} \tag{ 1 }$$

where D₀ is the short-time diffusivity of the particle and $\mathbf{u} = [u,v,w]^{\mathrm{T}}$ is a quasi-steady low-Reynolds-number flow field. The particle is suspended in a channel with a structured lower wall, S_w , and a planar upper wall, S_h , at z = h (measured from a reference surface S₀), see figure 1. The structured wall S_w is described by the profile $z = aH(\mathbf{r}_\parallel)$, where $\mathbf{r}_\parallel = [x,y,0]^{\mathrm{T}}$ are the in-plane coordinates, a denotes the characteristic surface amplitude, and $H(\mathbf{r}_\parallel)$ is the surface shape function. In addition, the upper wall is moving at a speed u₀ along the x-direction and the lower wall is stationary.

Zoom In Zoom Out Reset image size

We assume that the channel is not confined along the horizontal directions. At the top and bottom of the channel the probability density obeys a no-flux boundary condition,

$$\begin{align} \mathbf{n}\cdot\boldsymbol{\nabla} c(\mathbf{r},t) = 0 \qquad \mathrm{on}\; S_h \;\mathrm{and}\; S_w, \end{align} \tag{ 2 }$$

where n denotes the vector normal to the upper planar surface S_h and the lower corrugated surface S_w , respectively. We are interested in deriving the behavior of the height-averaged probability density

$$\begin{align} C(\mathbf{r}_\parallel,t) = \frac{1}{h-aH(\mathbf{r}_\parallel)} \int_{aH(\mathbf{r}_\parallel)}^h c(\mathbf{r},t)\,\mathrm{d} z \equiv\left\langle c(\mathbf{r},t)\right\rangle, \end{align} \tag{ 3 }$$

where $\langle \cdot \rangle$ represents the height-averaging. Unlike in traditional Taylor-dispersion theory, where averaging over the cross-section of the channel leads to an effective one dimensional transport equation, we retain two horizontal spatial dimensions. We decompose the full probability density as

$$\begin{align} c(\mathbf{r},t) = C(\mathbf{r}_\parallel,t) + \tilde{c}(\mathbf{r},t), \end{align} \tag{ 4 }$$

where $\tilde{c}$ is a perturbation term that accounts for the variations in the probability density along the z-direction.

Inserting (4) into (1) and defining the in-plane gradient $\nabla_\parallel = [\partial_x,\partial_y, 0]^{\mathrm{T}}$ we obtain

$$\begin{align} \partial_t (C + \tilde{c}) = - \boldsymbol{\mathbf{\nabla}}\cdot(\mathbf{u}(C+\tilde{c}))+D_0 (\nabla_\parallel^2 C + \nabla^2 \tilde{c}). \end{align} \tag{ 5 }$$

Height averaging (5) and recalling the no-flux boundary condition (2), we find

$$\begin{align} \partial_t C & = -\langle \boldsymbol{\mathbf{\nabla}} \cdot(\mathbf{u}(C+\tilde{c}))\rangle + D_0 \nabla_\parallel^2 C,\nonumber\\ & = -\boldsymbol{\mathbf{\nabla}}_\parallel \cdot(\langle \mathbf{u}\rangle C) - \langle \partial_z w\rangle C - \langle \boldsymbol{\mathbf{\nabla}}\cdot (\mathbf{u}\tilde{c})\rangle + D_0 \nabla^2_\parallel C, \end{align} \tag{ 6 }$$

where by definition $\langle \tilde{c}\rangle = 0$, $\langle\nabla^2\tilde{c}\rangle$ vanishes by the Leibniz rule, and we have:

$$\begin{align} \langle \partial_z w\rangle = \frac{w(z=h) - w(z=aH(\mathbf{r}_\parallel))}{h-aH(\mathbf{r}_\parallel)}. \end{align} \tag{ 7 }$$

To derive an equation for $\tilde{c}$ we subtract (6) from (5):

$$\begin{align} \partial_t \tilde{c}& = - \left[\boldsymbol{\mathbf{\nabla}} \cdot(\mathbf{u}C) - \left( \boldsymbol{\mathbf{\nabla}}_\parallel \cdot(\langle \mathbf{u}\rangle C) + \langle \partial_z w\rangle C\right) \right]\nonumber\\ & \quad - \left[\boldsymbol{\mathbf{\nabla}}\cdot (\mathbf{u}\tilde{c}) - \langle \boldsymbol{\mathbf{\nabla}} \cdot (\mathbf{u}\tilde{c}) \rangle \right] + D_0 \nabla^2 \tilde{c}. \end{align} \tag{ 8 }$$

As we are interested in the long-time behavior of the particle, we consider time scales greater than the time scale of diffusion in the vertical direction, $t \gtrsim h^2/D_0$ [53]. In this limit, diffusion across the channel averages out and we can consider the steady-state solution for $\tilde{c}$, $\partial_t \tilde{c} = 0$. Further, we expect that the perturbations in the probability density are smaller than the average, $\tilde{c} \lesssim C$, which leads us to neglect the second term of (8); this step is standard in this type of theoretical development. Finally, since $\tilde{c}$ encodes the entire fluctuations of the probability-density in the z-direction, we assume that the diffusion in z is dominant. These assumptions lead to a simplified form of (8):

$$\begin{align} D_0 \partial^2_z \tilde{c}= \boldsymbol{\mathbf{\nabla}} \cdot(\mathbf{u}C) - \left( \boldsymbol{\mathbf{\nabla}}_\parallel \cdot(\langle \mathbf{u}\rangle C) + \langle \partial_z w\rangle C\right). \end{align} \tag{ 9 }$$

Subsequently, we rescale times by $h^2/D_0$, lengths by h, and velocities by u₀. Introducing the Péclet number $\mathrm{Pe}\equiv u_0 h/D_0$, a dimensionless roughness $\epsilon = a/h$, and exploiting incompressibility, (6) and (9) become

$$\begin{align} \partial_t C &= -\mathrm{Pe}\left[\boldsymbol{\mathbf{\nabla}}_\parallel\cdot(\langle \mathbf{u}\rangle C) -\langle \partial_z w\rangle C - \langle \mathbf{u}\cdot \boldsymbol{\mathbf{\nabla}} \tilde{c}\rangle\right] + \nabla^2_\parallel C, \end{align} \tag{ 10a }$$

$$\begin{align} \partial^2_z\tilde{c} &= \mathrm{Pe}\left[\mathbf{u}\cdot\boldsymbol{\mathbf{\nabla}}_\parallel C - \boldsymbol{\mathbf{\nabla}}_\parallel \cdot\left(\langle \mathbf{u}\rangle C\right)- \langle \partial_z w\rangle C \right]. \end{align} \tag{ 10b }$$

We note here that the in-plane divergence $\nabla_\parallel\cdot\langle\mathbf{u}\rangle$ does not vanish in general. Specifying the surface shape $H(\mathbf{r}_\parallel)$ and the associated hydrodynamic flows u, allows us to recast (10a ) into a classical, two-dimensional convection–diffusion equation with effective transport parameters.

2.1. Perturbation flow due to surface roughness

The base shear flow $\mathbf{u}_0 = u_0 z/h \mathbf{e}_x$ that arises due to the motion of the upper boundary becomes perturbed by the presence of the surface topography. At low Reynolds numbers, the flow is governed by the continuity and Stokes equations. Non-dimensionalizing velocities u by u₀ and the pressure field p by $\mu u_0/h$, where µ denotes the viscosity of the fluid, the governing equations are

$$\begin{align} \nabla\cdot \mathbf{u} = 0, \end{align} \tag{ 11a }$$

$$\begin{align} -\nabla p + \nabla^2 \mathbf{u} = \mathbf{0}. \end{align} \tag{ 11b }$$

As we consider the case of a moving, smooth top boundary S_h and a static rough bottom boundary S_w , the boundary conditions are $\mathbf{u} = \mathbf{u}_0 = \mathbf{e}_x$ on S_h and $\mathbf{u} = \textbf{0}$ on S_w . For small surface roughness $\epsilon = a/h\ll1$, we can expand the flow field in terms of ε as

$$\begin{align} \mathbf{u} = \mathbf{u}^{(0)} + \epsilon \mathbf{u}^{(1)} + \epsilon^2 \mathbf{u}^{(2)} + \mathcal{O}(\epsilon^3), \end{align} \tag{ 12 }$$

and similarly the pressure field $p = p^{(0)} + \epsilon p^{(1)} + \epsilon^2 p^{(2)} + \mathcal{O}(\epsilon^3)$. We include terms up to ε², as we will show that these terms determine the leading-order effects of the surface roughness on dispersion. Next, we apply a domain perturbation to the boundary conditions [54, 55], which allows us to replace the structured wall S_w by a planar reference surface S₀, with the boundary conditions:

$$\begin{align} \mathbf{u}^{(0)}&=\mathbf{0}, \end{align} \tag{ 13a }$$

$$\begin{align} \mathbf{u}^{(1)}&=-H(\mathbf{r}_\parallel)\partial_z\mathbf{u}^{(0)}|_{z=0}, \end{align} \tag{ 13b }$$

$$\begin{align} \mathbf{u}^{(2)}& =-H(\mathbf{r}_\parallel)\partial_z\mathbf{u}^{(1)}|_{z=0}-\frac{H(\mathbf{r}_\parallel)^2}{2}\partial^2_z\mathbf{u}^{(0)}|_{z=0}. \end{align} \tag{ 13c }$$

From the boundary conditions we immediately find that the zeroth-order flow is simple shear $\mathbf{u}^{(0)} = z\mathbf{e}_x$ and therefore the boundary conditions reduce to: $\mathbf{u}^{(1)} = -H(\mathbf{r}_\parallel)\mathbf{e}_x$ and $\mathbf{u}^{(2)} = -H(\mathbf{r}_\parallel)\partial_z\mathbf{u}^{(1)}|_{z = 0}$ on S₀. Hence, we can model the perturbation velocity due to the presence of surface topography as an expansion in terms of ε, with an effective slip velocity on a flat surface S₀ entering into the boundary conditions for every order. This slip velocity effectively encodes the presence of the surface roughness while allowing us to consider flow between two flat parallel plates. The roughness-induced velocity field can be calculated analytically for surface shapes of the form

$$\begin{align} H(\mathbf{r}_\parallel) = \sum_{i=1}^N \alpha_{i} \cos(\mathbf{k}_{i}^\alpha\cdot\mathbf{r}_\parallel)+ \beta_{i}\sin(\mathbf{k}_{i}^\beta\cdot\mathbf{r}_\parallel), \end{align} \tag{ 14 }$$

with N modes, arbitrary coefficients α_i , β_i , and wave vectors $\mathbf{k}_{i}^\alpha$, $\mathbf{k}_{i}^\beta$ (see appendix A). This result implies that the flow fields can be obtained analytically for any shape function expandable in terms of a Fourier series.

We next consider a limiting case of (14) of a simple corrugated surface, which are commonly used in microfluidic settings [40] to help us develop the form of the dispersion tensor. Following this, to study tracer dispersion, we also consider more complex forms of (14), including periodic bumpy as well as randomly structured surfaces. See figure 2 for exemplary surface structures.

Zoom In Zoom Out Reset image size

2.2. Corrugated surface (single mode)

2.2.1. Details of the flow field.

We start our analysis by considering a surface shape consisting of a single Fourier mode

$$\begin{align} H(\mathbf{r}_\parallel) = \cos(\mathbf{k}\cdot\mathbf{r}_\parallel), \end{align} \tag{ 15 }$$

where $\mathbf{k} = [k_x,k_y,0]^{\mathrm{T}}$ denotes the in-plane wave vector (figure 2(a)). Due to the boundary condition $\mathbf{u}^{(1)}|_{z = 0} = -\cos(\mathbf{k}\cdot\mathbf{r}_\parallel)\mathbf{e}_x$ (see (13b )) the first-order flow fields must take the form

$$\begin{align} u^{(1)}(\mathbf{r}; \mathbf{k}) &=\hat{u}^{(1)}(z; \mathbf{k})\cos(\mathbf{k}\cdot\mathbf{r}_\parallel), \end{align} \tag{ 16a }$$

$$\begin{align} v^{(1)}(\mathbf{r}; \mathbf{k}) &=\hat{v}^{(1)}(z; \mathbf{k})\cos(\mathbf{k}\cdot\mathbf{r}_\parallel), \end{align} \tag{ 16b }$$

$$\begin{align} w^{(1)}(\mathbf{r}; \mathbf{k}) &= \hat{w}^{(1)}(z; \mathbf{k})\sin(\mathbf{k}\cdot\mathbf{r}_\parallel), \end{align} \tag{ 16c }$$

and consequently the pressure can be written as $p^{(1)}(\mathbf{r}; \mathbf{k}) = \hat{p}^{(1)}(z; \mathbf{k})\sin(\mathbf{k}\cdot\mathbf{r}_\parallel)$. Inserting this ansatz into the continuity equation (11a ), we can write the relationship between the z-component of the velocity field and the x- and y-components via

$$\begin{align} \frac{\mathrm{d} \hat{w}^{(1)}}{\mathrm{d} z} &= k_x\hat{u}^{(1)}+k_y\hat{v}^{(1)} = \mathbf{k}\cdot[\hat{u}^{(1)}, \hat{v}^{(1)},0]^{\mathrm{T}}. \end{align} \tag{ 17 }$$

Using this relation and the Stokes equations (11b ) we can obtain analytical expressions for the velocity field using Mathematica (see appendix A). It is convenient to introduce a vector field $\mathbf{Q}(z; \mathbf{k}) = [Q_1,Q_2,0]^{\mathrm{T}}$, which is related to the velocity field via $\mathrm{d}^2 Q_1/\mathrm{d} z^2 = \hat{u}^{(1)}$ and $\mathrm{d}^2 Q_2/\mathrm{d} z^2 = \hat{v}^{(1)}$. Then, using (17), the first-order roughness-induced velocities can be expressed in terms of Q via

$$\begin{align} \mathbf{u}^{(1)}(\mathbf{r}; \mathbf{k})&= \cos(\mathbf{k}\cdot\mathbf{r}_\parallel)\mathbf{Q}_{zz}+\sin(\mathbf{k}\cdot\mathbf{r}_\parallel)(\mathbf{k}\cdot\mathbf{Q}_z)\mathbf{e}_z, \end{align} \tag{ 18 }$$

where the subscript z denotes the derivative with respect to z. We note that the first-order flow fields are periodic functions of the in-plane coordinates x and y.

From the first-order velocities, the second-order boundary condition on z = 0 (13c ) becomes

$$\begin{align} \mathbf{u}^{(2)}(\mathbf{r}; \mathbf{k})\bigl|_{z=0} &= -\frac{1}{2}(1+\cos(2\mathbf{k}\cdot\mathbf{r}_\parallel))\mathbf{Q}_{zzz}\bigl|_{z=0}\nonumber\\ & \quad +\frac{1}{2}\sin(2\mathbf{k}\cdot\mathbf{r}_\parallel)(\mathbf{k}\cdot\mathbf{e}_x)\mathbf{e}_z, \nonumber\\ & \equiv \bar{\mathbf{u}}_0^{(2)}(\mathbf{k}) +\tilde{\mathbf{u}}^{(2)}(\mathbf{r}; \mathbf{k})|_{z=0}, \end{align} \tag{ 19 }$$

where we have written the second-order boundary condition as a combination of a non-periodic, $\bar{\mathbf{u}}^{(2)}_0(\mathbf{k})$, and a periodic, $\tilde{\mathbf{u}}^{(2)}(\mathbf{r}; \mathbf{k})$, contribution, which is periodic over the spatial coordinate $\mathbf{r}_\parallel$. Therefore, the second-order flow field $\mathbf{u}^{(2)}(\mathbf{r}; \mathbf{k})$ assumes the general, abbreviated form

$$\begin{align} \mathbf{u}^{(2)}(\mathbf{r}; \mathbf{k}) = \bar{\mathbf{u}}_0^{(2)}(\mathbf{k})(1-z)+\tilde{\mathbf{u}}^{(2)}(\mathbf{r};\mathbf{k}). \end{align} \tag{ 20 }$$

We note that the z-component of the velocity field is completely periodic with zero mean and hence $\bar{w}_0^{(2)}(\mathbf{k}) = 0$. The expressions for the velocity fields have been obtained analytically using Mathematica. The non-periodic contributions, $\bar{\mathbf{u}}^{(2)}_0(\mathbf{k})$, can be found in appendix A. Further, in appendix D we compare these asymptotic results to direct numerical simulations and discuss the validity of such approximations when applied to some of the geometries we consider in this work.

2.2.2. Derivation of the drift-diffusion equation.

In what follows we present a detailed derivation of the drift-diffusion equation for a 3D corrugated surface shape. We start by substituting the velocity field in (18) and (20) into (10a ) and (10b ):

$$\begin{align} \partial_t C = \mathrm{Pe}&\left[-\frac{1}{2}\left(\partial_x + \epsilon^2\bar{\mathbf{u}}_0^{(2)}\cdot \boldsymbol{\mathbf{\nabla}}_\parallel \right)C\right. \nonumber\\ & \;\;\! \left.- \left\langle \left(z \partial_x + (\epsilon \mathbf{u}^{(1)} + \epsilon^2\mathbf{u}^{(2)}) \cdot \boldsymbol{\mathbf{\nabla}} \right)\tilde{c}\right\rangle \right]\nonumber\\ & \;\; +\nabla^2_\parallel C+ \mathrm{p.t.} +\mathcal{O}(\epsilon^3), \end{align} \tag{ 21a }$$

$$\begin{align} \partial_z^2 \tilde{c} & = \mathrm{Pe}\Bigg[\left(z-\frac{1}{2}\right)\partial_x + \epsilon \left(\mathbf{u}^{(1)} \cdot \boldsymbol{\mathbf{\nabla}}_\parallel - \boldsymbol{\mathbf{\nabla}}_\parallel\cdot \left\langle \mathbf{u}^{(1)}\right\rangle\right)\nonumber\\ & \quad + \epsilon^2\left(\left(\frac{1}{2}-z\right) \bar{\mathbf{u}}^{(2)}_0 \cdot \boldsymbol{\mathbf{\nabla}}_\parallel \right)\Bigg]C + \mathrm{p.t.} +\mathcal{O}(\epsilon^3). \end{align} \tag{ 21b }$$

Here, 'p.t.' (periodic terms) includes all terms that vary periodically over the horizontal coordinate $\mathbf{r}_\parallel$. As we are concerned primarily with the long-time, average diffusivity of the particles we only consider the non-periodic contributions to the particle concentration. All of the periodic effects will average out when considering macroscopic motion, and so we collect all of these terms in 'p.t.' for further calculations. At this stage, this includes several terms of $\mathcal{O}(\epsilon^2)$ as well as the term $\epsilon \boldsymbol{\mathbf{\nabla}}_\parallel\cdot(\langle \mathbf{u}^{(1)} \rangle C)$ in (21a ). While some periodic terms included in p.t. in $\tilde{c}$ may lead to non-periodic contributions to C, such terms will be $\mathcal{O}(\epsilon^3)$ or smaller and so can be neglected.

We proceed by integrating (21b ) twice in z to obtain

$$\begin{align} \tilde{c} =&\,\, \mathrm{Pe}\Bigg[\left(\frac{z^3}{6}-\frac{z^2}{4}\right)\partial_x + \epsilon \Bigg(\cos(\mathbf{k}\cdot\mathbf{r}_\parallel) \left(\mathbf{Q}\cdot\boldsymbol{\mathbf{\nabla}}_\parallel - \frac{z^2}{2}\langle \mathbf{Q}_{zz}\rangle\cdot\boldsymbol{\mathbf{\nabla}}_\parallel\right)\nonumber \\ & \qquad + \sin(\mathbf{k}\cdot\mathbf{r}_\parallel)\frac{z^2}{2}\langle \mathbf{Q}_{zz}\rangle\cdot\mathbf{k} \Bigg) - \epsilon^2 \left(\frac{z^3}{6}-\frac{z^2}{4}\right) \bar{\mathbf{u}}_0^{(2)} \cdot \boldsymbol{\mathbf{\nabla}}_\parallel \Bigg]\nonumber \\ & \times C + \lambda_1 z + \lambda_2 + \mathrm{p.t.} + \mathcal{O}(\epsilon^3) \end{align} \tag{ 22 }$$

where λ₁ and λ₂ are constants of integration. By taking the height average of the continuity equation, we find that

$$\begin{align} \left\langle\mathbf{Q}_{zz}\right\rangle\cdot\mathbf{k} = \left\langle\partial_z \tilde{w}^{(1)}(z) \right\rangle= \tilde{w}^{(1)}(1)-\tilde{w}^{(1)}(0)=0, \end{align} \tag{ 23 }$$

which simplifies (22). The no-flux boundary conditions, $\partial_z \tilde{c} = 0$ on z = 0 and z = 1, provide the same condition for λ₁:

$$\begin{align} \lambda_1 = -\epsilon \cos(\mathbf{k}\cdot\mathbf{r}_\parallel) \mathbf{Q}_z(z=0) \cdot \boldsymbol{\mathbf{\nabla}}_\parallel C + \mathrm{p.t.}+ \mathcal{O}(\epsilon^3). \end{align} \tag{ 24 }$$

As by definition the vertical average vanishes, i.e. $\langle \tilde{c}\rangle = 0$, λ₂ is

$$\begin{align} \lambda_2 = \mathrm{Pe}&\Bigg[ \frac{1}{24}\partial_x\! -\! \epsilon \cos(\mathbf{k}\cdot\mathbf{r}_\parallel) \bigg(\langle \mathbf{Q}\rangle\!-\! \frac{1}{6}\langle \mathbf{Q}_{zz}\rangle -\frac{1}{2}\mathbf{Q}_z(0)\bigg)\cdot\boldsymbol{\mathbf{\nabla}}_\parallel\nonumber \\ &\quad - \frac{\epsilon^2}{24} \bar{\mathbf{u}}_0^{(2)} \cdot \nabla_\parallel\Bigg]C + \mathrm{p.t.} + \mathcal{O}(\epsilon^3). \end{align} \tag{ 25 }$$

Thus, the perturbed probability density $\tilde{c}(\mathbf{r},t; \mathbf{k})$ is given by

$$\begin{align} \tilde{c} &= \mathrm{Pe}\left[g(z) \partial_x + \epsilon \cos(\mathbf{k}\cdot\mathbf{r}_\parallel)\mathbf{f}\cdot \boldsymbol{\mathbf{\nabla}}_\parallel - \epsilon^2 g(z) \bar{\mathbf{u}}_0^{(2)} \cdot \boldsymbol{\mathbf{\nabla}}_\parallel \right]\nonumber \\ &\quad \times C + \mathrm{p.t.} + \mathcal{O}(\epsilon^3), \end{align} \tag{ 26 }$$

where we abbreviated the vector field

$$\begin{align} \mathbf{f}(z; \mathbf{k}) &= \mathbf{Q}(z; \mathbf{k}) - \left(\frac{z^2}{2}-\frac{1}{6}\right)\langle \mathbf{Q}_{zz}(z;\mathbf{k})\rangle -\langle \mathbf{Q}(z;\mathbf{k}) \rangle\nonumber \\ &\quad - \left(z-\frac{1}{2}\right)\mathbf{Q}_z(z=0;\mathbf{k}). \end{align} \tag{ 27 }$$

and introduced the function

$$\begin{align} g(z) = \frac{z^3}{6}-\frac{z^2}{4} + \frac{1}{24}. \end{align} \tag{ 28 }$$

To derive the effective drift-diffusion equation up to second order in the surface roughness $\mathcal{O}(\epsilon^2)$, we insert (26) into (10a ) and compute the contributions of the terms $\left\langle \left(z \partial_x + (\epsilon \mathbf{u}^{(1)} + \epsilon^2\mathbf{u}^{(2)})\cdot \boldsymbol{\mathbf{\nabla}}\right)\tilde{c}\right\rangle$ to $\mathcal{O}(\epsilon^2)$. The $\mathcal{O}(\epsilon^2)$ contribution to the first term is

$$\begin{align} \langle z \partial_x \tilde{c}\rangle = \frac{\mathrm{Pe}}{120}\left[-\partial_x^2 + \epsilon^2\bar{\mathbf{u}}_0^{(2)} \cdot \boldsymbol{\mathbf{\nabla}}_\parallel\partial_x \right]C + \mathrm{p.t.} + \mathcal{O}(\epsilon^3). \end{align} \tag{ 29 }$$

We further calculate non-vanishing terms to $\mathcal{O}(\epsilon)$ of

$$\begin{align} \langle \mathbf{u}^{(1)} \cdot \boldsymbol{\mathbf{\nabla}} \tilde{c}\rangle = \frac{1}{2}\epsilon\mathrm{Pe}\left(\boldsymbol{\mathbf{\nabla}}_\parallel \cdot \left\langle \mathbf{Q}_{zz} \mathbf{f}\right\rangle\cdot\boldsymbol{\mathbf{\nabla}}_\parallel\right) C + \mathrm{p.t.} + \mathcal{O}(\epsilon^2), \end{align} \tag{ 30 }$$

and terms to $\mathcal{O}(1)$ of

$$\begin{align} \langle \mathbf{u}^{(2)} \cdot \boldsymbol{\mathbf{\nabla}} \tilde{c} \rangle = \frac{\mathrm{Pe}}{120}\left(\bar{\mathbf{u}}^{(2)}_0 \cdot \boldsymbol{\mathbf{\nabla}}_\parallel\right) \partial_x C + \mathrm{p.t.} + \mathcal{O}(\epsilon). \end{align} \tag{ 31 }$$

By inserting these terms into (21a ), the convection–diffusion (Fokker–Planck) equation (1) reduces to a 2D problem when considering only the non-periodic effects. Therefore, we will drop the z component of the vectors $\boldsymbol{\mathbf{\nabla}}_\parallel$, $\mathbf{r}_\parallel$, and k, henceforth treating them as two-dimensional quantities.

Similar to previous work on dispersion through spatially periodic, porous media [26], we can divide the domain into subdomains of length L_x along the x-direction and L_y along the y-direction, where the length scales L_x and L_y are set by the underlying periodic surface structure. The drift-diffusion equation contains periodic terms (p.t.), which vary locally within one subdomain as a result of the periodic flows generated by the surface structure, and non-periodic terms, which are constant in the spatial coordinate $\mathbf{r}_\parallel$. We will consider (here in dimensional form) coarse-grained length scales $x,y\gg L_x,L_y$ and long times $t\gg \mathrm{max} (L_x,L_y)/(\epsilon u_0)$, so that the local effects of the periodic terms can be neglected, thus dropping all 'p.t.' from further consideration.

The 2D convection–diffusion equation for the concentration (probability density) $C\equiv C(\mathbf{r}_\parallel,t;\mathbf{k})$ of tracer particles considered on these coarse-grained scales is

$$\begin{align} \partial_t C = - \boldsymbol{\mathbf{\nabla}}_\parallel \cdot\left(\mathbf{u}_{\mathrm{eff}} C\right) + \boldsymbol{\mathbf{\nabla}}_\parallel \cdot \left(\mathbf{D}\cdot \boldsymbol{\mathbf{\nabla}}_\parallel C\right) +\mathcal{O}(\epsilon^3), \end{align} \tag{ 32 }$$

with effective drift velocity and dispersion tensor:

$$\begin{align} \mathbf{u}_{\mathrm{eff}} &= \mathbf{u}_{\mathrm{eff}}^{(0)} + \epsilon^2 \mathbf{u}_{\mathrm{eff}}^{(2)} +\mathcal{O}(\epsilon^3), \end{align} \tag{ 33a }$$

$$\begin{align} \mathbf{D} &= \mathbf{D}^{(0)} + \epsilon^2 \mathbf{D}^{(2)}+\mathcal{O}(\epsilon^3). \end{align} \tag{ 33b }$$

Note that we have defined $\mathbf{u}_{\mathrm{eff}}$ and D to include the dependence on Péclet number. In agreement with earlier work, the zeroth-order contributions to the effective drift velocity and dispersion are:

$$\begin{align} \mathbf{u}_{\mathrm{eff}}^{(0)}=\frac{1}{2}\mathrm{Pe}\left[\begin{array}{c}1 \\0\end{array}\right] \end{align} \tag{ 34a }$$

$$\begin{align} \mathbf{D}^{(0)} = \left[ \begin{array}{cc} 1 + \frac{\mathrm{Pe}^2}{120} & 0 \\ 0 & 1 \end{array} \right]. \end{align} \tag{ 34b }$$

The roughness-induced contributions, which represent a main finding of this work, evaluate to:

$$\begin{align} \mathbf{u}_{\mathrm{eff}}^{(2)}(\mathbf{k})=\frac{1}{2}\mathrm{Pe}\left[\begin{array}{c}\bar{u}_0^{(2)}(\mathbf{k})\\[3pt] \bar{v}_0^{(2)}(\mathbf{k})\end{array}\right], \end{align} \tag{ 35 }$$

and

$$\begin{align} \mathbf{D}^{(2)}(\mathbf{k}) &= - \frac{\mathrm{Pe}^2}{4} \left(\frac{1}{15} \left[ \begin{array}{cc} \bar{u}_0^{(2)}(\mathbf{k}) & \frac{1}{2}\bar{v}_0^{(2)}(\mathbf{k}) \\[3pt] \frac{1}{2}\bar{v}_0^{(2)}(\mathbf{k}) & 0 \end{array} \right] \right.\nonumber \\ &\quad \left. +\left(\left\langle \mathbf{Q}_{zz}(z;\mathbf{k}) \mathbf{f}(z;\mathbf{k}) \right\rangle +\left\langle \mathbf{Q}_{zz}(z;\mathbf{k}) \mathbf{f}(z;\mathbf{k}) \right\rangle^{\mathrm{T}}\right) \right). \end{align} \tag{ 36 }$$

The effective drift and dispersion tensor are functions of the wave vector k, the roughness ε, and the Péclet number (Pe). Note that only even powers in ε contribute to the effective drift and dispersion tensors. In particular, the zeroth-order contributions correspond to the classical Taylor dispersion in shear flow between two planar walls and the second-order contribution reflects the effects of the surface topography.

The form of the effective drift suggests that surface structures can induce particles to drift in directions that are different from the direction of applied forcing (flow), as surface corrugations can deflect the average flow direction by an angle β, see figure 3.

Zoom In Zoom Out Reset image size

Furthermore, we observe that the roughness-induced dispersion tensor is symmetric and contains non-vanishing off-diagonal terms, which represent an xy-coupling. This indicates that the hydrodynamic interactions of the fluid flow with the surface structure not only alters the dispersion of tracer particles but also changes the nature of the long-time anisotropic diffusion. In particular, we note that the solution C to the coarse-grained convection–diffusion equation (32) assumes a bivariate Gaussian distribution

$$\begin{align} C = \frac{1}{\mathcal{N}}\exp\left[-\frac{1}{2} \left(\mathbf{r}_\parallel-\boldsymbol{\mu}\right)^{\mathrm{T}}\boldsymbol{\cdot}\boldsymbol{\Sigma}^{-1}\boldsymbol{\cdot}\left(\mathbf{r}_\parallel-\boldsymbol{\mu}\right)\right], \end{align} \tag{ 37 }$$

with normalization $\mathcal{N} = 2\pi \sqrt{\mathrm{det}(\boldsymbol{\Sigma})}$, mean $\boldsymbol{\mu} = \mathbf{u}_{\mathrm{eff}} t$, and variance $\boldsymbol{\Sigma} = 2 t \mathbf{D}$ with $\mathbf{u}_{\mathrm{eff}}$ and D from (33a ) and (33b ). The level sets of the distribution are ellipsoids with major and minor axis along the eigenvectors of the dispersion matrix, which are proportional to the respective eigenvalue, Λ₁ and Λ₂, see figure 3. Hence, the particle diffuses along the eigenvectors of the dispersion matrix with diffusivities, Λ₁ and Λ₂ ($\Lambda_1\geqslant\Lambda_2$).

It is important to note that the eigenvectors of the dispersion matrix generally do not align with the effective drift, $\mathbf{u}_{\mathrm{eff}}^{(0)}+\epsilon^2\mathbf{u}_{\mathrm{eff}}^{(2)}$, as one may naively expect. The angle between the principle eigenvector and the drift, $\alpha + \beta$, where α is the angle between the applied forcing and the eigenvector corresponding to Λ₁ and β the deflection of the average flow from the forcing (see figure 3) depend on details of the surface structure.

We discuss results for a single surface mode in section 3.1.

2.3. Bumpy surface (two modes)

We can further generalize this approach by considering a bumpy surface

$$\begin{align} H(\mathbf{r}_\parallel) = H_1 \cos(\mathbf{k}_{1}\cdot \mathbf{r}_\parallel) + H_2 \sin(\mathbf{k}_{2}\cdot \mathbf{r}_\parallel), \end{align} \tag{ 38 }$$

where k₁, k₂ are parallel to the surface S₀ and H₁ and H₂ are two dimensionless parameters. Following our previous notation, the horizontal components of the first-order velocity $\mathbf{u}^{(1)}$ assume the form

$$\begin{align} \mathbf{u}^{(1)}(\mathbf{r};\{H_i,\mathbf{k}_i\}) &= H_1\cos(\mathbf{k}_1\cdot\mathbf{r}_\parallel)\mathbf{Q}_{zz}(z; \mathbf{k}_1) \nonumber\\ &\quad +H_2\sin(\mathbf{k}_{2}\cdot \mathbf{r}_\parallel)\mathbf{Q}_{zz}(z; \mathbf{k}_2)\nonumber\\ &\quad+\left( H_1\sin(\mathbf{k}_1\cdot\mathbf{r}_\parallel)\mathbf{k}_1\cdot \mathbf{Q}_z(z; \mathbf{k}_1)\right.\nonumber\\ &\quad \left.-H_2\cos(\mathbf{k}_{2}\cdot \mathbf{r}_\parallel)\mathbf{k}_2\cdot\mathbf{Q}_z(z;\mathbf{k}_2)\right)\mathbf{e}_z. \end{align} \tag{ 39 }$$

Consequently, the second-order boundary conditions on z = 0 reduce to

$$\begin{align} \mathbf{u}^{(2)}(\mathbf{r}; \{H_i,\mathbf{k}_i\})\bigl|_{z=0} = -\frac{1}{2}\sum_{i=1}^2(H_i)^2 \mathbf{Q}_{zzz}(z=0; \mathbf{k}_i) + \mathrm{p.t.} \end{align} \tag{ 40 }$$

where 'p.t.' contains all periodic terms. Note that the non-periodic terms do not include a coupling term $H_1H_2$. The second-order velocities then assume the form

$$\begin{align} \mathbf{u}^{(2)}(\mathbf{r}; \{H_i,\mathbf{k}_i\}) = \bar{\mathbf{u}}_0^{(2)}(\{H_i,\mathbf{k}_i\})(1-z) +\tilde{\mathbf{u}}^{(2)}(\mathbf{r};\{H_i,\mathbf{k}_i\}), \end{align} \tag{ 41 }$$

with $\bar{\mathbf{u}}_0^{(2)}(\{H_i,\mathbf{k}_i\}) = \sum_{i = 1}^2H_i^2\bar{\mathbf{u}}_0^{(2)}(\mathbf{k}_i)$. The derivation of the diffusion tensor proceeds as in the case of a single mode surface up to (21a ) and (21b ). Due to the presence of the additional terms in the second-order flow, the equation for $\tilde{c}$ becomes

$$\begin{align} \tilde{c} = \mathrm{Pe}&\Bigg[g(z)\partial_x + \epsilon \big(H_1 \cos(\mathbf{k}_1\cdot \mathbf{r}_\parallel)\mathbf{f}(z;\mathbf{k}_1)\nonumber \\ &\quad + H_2 \sin(\mathbf{k}_2\cdot \mathbf{r}_\parallel) \mathbf{f}(z;\mathbf{k}_2)\big)\cdot \boldsymbol{\mathbf{\nabla}}_\parallel\nonumber \\ &\quad- \epsilon^2 g(z) \left(\sum_{i=1}^2 H_i^2 \bar{\mathbf{u}}_0^{(2)}(\mathbf{k}_i)\right) \cdot \boldsymbol{\mathbf{\nabla}}_\parallel \Bigg]C + \mathrm{p.t.} + \mathcal{O}(\epsilon^3), \end{align} \tag{ 42 }$$

where the contributions of different modes to the velocity fields have been summed (compare with (22)) and the term $\langle \mathbf{Q}_{zz}(z;\mathbf{k}_{1,2})\rangle\cdot\mathbf{k}_{1,2} = 0$ vanishes for each mode. We further replace $\bar{\mathbf{u}}_0^{(2)}$ in (29) and (31) by $\bar{\mathbf{u}}_0^{(2)}(\{H_i,\mathbf{k}_i\})$ (from (41)) and obtain the non-vanishing terms up to $\mathcal{O}(\epsilon)$ of

$$\begin{align} \langle \mathbf{u}^{(1)} \cdot \boldsymbol{\mathbf{\nabla}} \tilde{c}\rangle &= \frac{1}{2}\epsilon\mathrm{Pe}\Bigg(\sum_{i=1}^2 H_i^2\boldsymbol{\mathbf{\nabla}}_\parallel \cdot \langle \mathbf{Q}_{zz}(z;\mathbf{k}_i) \mathbf{f}(z;\mathbf{k}_i)\rangle \cdot\boldsymbol{\mathbf{\nabla}}_\parallel\Bigg) \nonumber \\ &\quad \times C + \mathrm{p.t.} + \mathcal{O}(\epsilon^3), \end{align} \tag{ 43 }$$

Thus, the equation for the concentration $C\equiv C(\mathbf{r}_\parallel, t;\{H_i,\mathbf{k}_i\})$ in 2D reduces to (32), where the zeroth-order contributions are the same as before (34a ) and (34b ) and the roughness-induced drift velocity and dispersion tensor are:

$$\begin{align} \mathbf{u}_{\mathrm{eff}}^{(2)}=\frac{1}{2}\mathrm{Pe}\sum_{i=1}^2H_i^2 \left[\begin{array}{c} \bar{u}_0^{(2)}(\mathbf{k}_i)\\[3pt] \bar{v}_0^{(2)}(\mathbf{k}_i) \end{array}\right] = \sum_{i=1}^2 H_i^2\mathbf{u}_{\mathrm{eff}}^{(2)}(\mathbf{k}_i) \end{align} \tag{ 44a }$$

$$\begin{align} \mathbf{D}^{(2)} = \sum_{i=1}^2 H_i^2\mathbf{D}^{(2)}(\mathbf{k}_i). \end{align} \tag{ 44b }$$

We find that the second-order contributions are the sums over the effective drift (34a ) and dispersion (36) caused by a single wavy perturbation with wave vector k_i . Coupling of surface modes, $H_1H_2$, is only relevant when examining effects at higher orders.

2.4. General, randomly structured surface

More generally, we consider the surface shape of the form of a discrete Fourier series

$$\begin{align} H(\mathbf{r}_\parallel) = \frac{1}{N}\sum_{i,j=1}^N \alpha_{ij} \cos(\mathbf{k}_{ij}\cdot\mathbf{r}_\parallel)+ \beta_{ij}\sin(\mathbf{k}_{ij}\cdot\mathbf{r}_\parallel), \end{align} \tag{ 45 }$$

with $\mathbf{k}_{ij} = 2\pi [i,j,0]^{\mathrm{T}}/(\lambda N)$, where λ denotes the characteristic (smallest) wavelength of the surface and $\lambda_N = \lambda N$ the longest wavelength. The latter sets the periodicity of the surface along the x- and y-directions, respectively, and hence the coarse scales $L_x = L_y = \lambda_N$. The coefficients $\{\alpha_{ij}, \beta_{ij}\}$ are non-dimensional constants. We can extend our procedure for the case of two modes to the general case of a surface with arbitrarily many modes, and obtain the analytical predictions for the roughness-induced effective velocities and dispersion:

$$\begin{align} \mathbf{u}_{\mathrm{eff}}^{(2)} = \frac{1}{N^2}\sum_{i,j=1}^N\left((\alpha_{ij})^2 + (\beta_{ij})^2\right)\mathbf{u}_{\mathrm{eff}}^{(2)}(\mathbf{k}_{ij}), \end{align} \tag{ 46a }$$

$$\begin{align} \mathbf{D}^{(2)} = \frac{1}{N^2} \sum_{i,j=1}^N\left((\alpha_{ij})^2+(\beta_{ij})^2\right)\mathbf{D}^{(2)}(\mathbf{k}_{ij}). \end{align} \tag{ 46b }$$

Our findings demonstrate that the roughness-induced contributions to the dispersion at second order can be separated into linearly independent components for each surface mode.

We now use the theoretical structure developed in the previous section to examine the effects of surface topography on particle transport in a few particular cases of interest, including a corrugated surface consisting of a single mode, a bumpy surface composed of two balanced modes, and a randomly structured surface.

3.1. Corrugated surface (1 mode)

First, we consider surfaces that have a corrugation wavelength λ and are oriented relative to the applied forcing direction e_x by an angle θ. This corresponds to a mode with wave vector $\mathbf{k} = \frac{2\pi}{\lambda}[\cos\theta,\sin\theta]^{\mathrm{T}}$. To compare with our coarse-grained result for the effective drift and dispersion tensor, we perform Brownian dynamics simulations of tracer particles suspended in a corrugated channel, where we approximate the flow fields using our domain perturbation solutions (see appendix C for details of the simulation; see appendix D for a detailed comparison of the approximate flow field with results obtained using finite-element methods). Figure 4 depicts our simulation results and theoretical predictions for a planar surface and a corrugated surface (with $\lambda/h = 1/\sqrt{2}$, $\theta = \pi/4$, and ε = 0.1; see figure 3(a)) with $\mathrm{Pe}~= 100$. Compared to the planar surface, which represents classical Taylor dispersion, the presence of surface corrugations causes the particles to drift in the $-\mathbf{e}_y$ direction. Furthermore, the axes of the elliptical tracer distribution tilt relative to the lab coordinate system, which is consistent with our prediction of an anisotropic dispersion tensor with non-vanishing off-diagonal elements. The level curves of the coarse-grained analytical prediction for C (37) capture the simulation results well.

Zoom In Zoom Out Reset image size

Keeping $\theta = \pi/4$, a quantitative comparison of the effective transport parameters obtained from simulations across a range of λ shows good agreement with our theoretical predictions. In particular, the results in figures 5(a) and (b) confirm that our theory very closely captures the effective drift of the particles in both directions. This indicates that our extension of the coarse-graining procedure established by Taylor can be used to predict the effective drift and dispersion in the presence of more complex flow fields. Note that for a surface with a wavelength of $\lambda/h\sim 1$ we expect that the particle experiences an effective drift in the $-\mathbf{e}_y$ direction with a velocity of about 2% of the mean speed. However, the same surface can lead to reductions of almost 10% in the direction of applied forcing.

Zoom In Zoom Out Reset image size

In addition to the drift, the surface corrugations amplify the dispersion. Figure 5(c) shows the predicted values of Λ₁ and Λ₂ (corresponding to the eigenvalues of the dispersion matrix) normalized by the equivalent diffusion constant from classical Taylor dispersion theory (34b ). Note that we will maintain this normalization for the rest of this paper. We find that our simulation results match the theoretical predictions for Λ₁ well, with greater enhancement of diffusion at shorter wavelengths before settling to a constant value at longer wavelengths. Furthermore, the results show that $\Lambda_{1}\gtrsim D^{(0)}_{xx}$, which is indicative of the surface roughness continuing to have an effect on dispersion even when the wavelengths λ of the disturbances are very large. On the contrary, our theory predicts that Λ₂ is decreased at small λ and becomes slightly increased at large λ. It is important to note, however, that at the Péclet number and wavelength λ used here Λ₂ is comparable to molecular diffusion and hence the effect is so small that we could not reliably resolve it from our simulations.

Having verified our theoretical predictions for the coarse-grained effective drift and dispersion tensor with these simulations, in figure 6 we study the effects of surface roughness across a wide range of surface corrugations with different λ and θ. We note that our theory is expected to be valid as the ratio $\epsilon/\lambda\to 0$ [55], and therefore we expect that for small λ our perturbation approach may not be accurate. Further studies including finite-element solutions for more surfaces than addressed in appendix D will help to verify the range of validity of the results for other surface geometries.

Zoom In Zoom Out Reset image size

The surface structure tends to decrease the effective drift $|\mathbf{u}_{\mathrm{eff}}|$, although this effect gets smaller at higher wavelengths (figure 6(a)). The magnitude of drift is only weakly dependent on θ at low wavelengths but increases with θ at high wavelengths. Longer wavelength surfaces with a wave vector oriented perpendicular to the base flow direction (grooves aligned with the flow direction) tend to have the smallest effect on dispersion, as any given tracer particle effectively experiences a locally flat surface.

The orientation β of the drift velocity relative to the forcing direction e_x , is maximally reduced for θ = 0 and $\pi/2$ and for long wavelengths, with the highest deflections predicted at low wavelengths for $\theta\approx \pi/4$ (figure 6(b)). Note that this angle follows the sign convention depicted in figure 3, such that a positive surface orientation angle θ leads to drift in the negative e_y direction.

It is important to emphasize that the principal direction of dispersion does not align with the mean flow direction. In the absence of Brownian effects, as shown by [47] for sedimenting spheres and [48] for point particles in shear flow, particles moving near a wavy surface exhibit helical trajectories due to hydrodynamic couplings with the surface structure, which leads to transport across the grooves (see exemplary flow fields in figure D1). We expect that similar helical flow trajectories lead to the observed difference between the induced flow direction and the principal axes of diffusion. The variations in the angle between the principle eigenvector and the drift velocity, $\alpha + \beta$, is qualitatively similar to that of β (figure 6(c)). We note that the angle between the principle eigenvector and the base flow α is not depicted but behaves similarly, although with a smaller magnitude.

Turning now to the dispersion, small surface wavelength ($\lambda \sim 0.5$) generate more than 10% increased dispersion along the principal eigenvector (figure 6(d)). This effect is more pronounced when the corrugations are perpendicular to the applied forcing (θ = 0) and decreased when they are aligned with the applied forcing ($\theta = \pi/2$). At large λ the dispersion enhancement is reduced, as the perturbation flows driven by the surface topology weaken. Dispersion is observed to be enhanced only by a few percent at λ = 50. At shorter wavelengths, the dispersion is primarily controlled by changes in λ, while θ plays a secondary role. However, as λ gets large, θ plays the dominant role in determining the strength of dispersion, with less enhancement as the wave vector approaches an angle of $\theta = \pi/2$ (corresponding to surface grooves aligned along e_x ).

The behavior of Λ₂ for a range of surfaces (figure 6(e)) is markedly different from that of Λ₁. For shorter wavelengths, dispersion along the secondary eigenvector is less than classical Taylor dispersion (and hence molecular diffusion). This reduction is strongest near $\theta = \pi/4$ and falls off symmetrically. For $\lambda\sim1$ the difference between the surfaces with roughness and classical theory is small. At large $\lambda \gtrsim 1$ the dispersion in this direction is enhanced relative to classical theory, with maximal enhancement for $\theta = \pi/4$. However, this effect is small, contributing less than 1% to an increase in dispersion.

Finally, we consider the product of the two eigenvalues $\Lambda_1\Lambda_2$ (figure 6(f)), which is proportional to the area of the level curves of (37) describing the concentration of the tracer particles. For short wavelengths, $\lambda\lesssim 1$, the angle θ plays a large role in the dispersion of particles, with enhancements of over 20% or reductions in spreading relative to classical Taylor dispersion depending on the orientation angle θ. For longer wavelengths $\lambda \gtrsim 1$, the spreading is maximized for θ = 0 and roughness can provide enhancements to the overall spreading of around $5\%$ for $\lambda\gtrsim 5$.

3.2. Bumpy surface (two modes)

Having analyzed several important features of single-mode surfaces, we now turn to bumpy surfaces consisting of two modes. Here, we are particularly interested in symmetrical bumpy surfaces, such as that depicted in figure 2(b). In this case, the two wave vectors k₁ and k₂ are related by $k_{1,x} = k_{2,x}$ and $k_{1,y} = - k_{2,y}$ ($\mathbf{k}_1\cdot\mathbf{k}_2 = 0$). Furthermore, we take $H_1 = H_2 = 1/2$ in (38) so that the maximum height of the bumps is always 1. We also describe the surfaces in terms of an angle θ measured relative to k₁, such that $\theta = \arctan(k_{1,y}/k_{1,x})$, and a wavelength $\lambda = 2\pi/|\mathbf{k}_1|$.

Due to the reflectional symmetry of the wave vectors, the dispersion tensor for the resulting surface becomes diagonal when examined in the lab frame. Mathematically, the corrections to the dispersion tensor resulting from each surface mode are even in k_y for the diagonal terms but odd in k_y for the off-diagonal terms (see appendix B). As shown in (44b ), the dispersion for each mode is additive and hence the corresponding off-diagonal contributions cancel when choosing wave vectors of this form.

Combining wave vectors in this manner results in a surface that does not have a preferred direction. Therefore any drift along e_y is eliminated and we can analyze the diffusion tensor in the lab frame, as the two non-zero components are the eigenvalues of the tensor. While the diffusion is still anisotropic, as diffusion in the e_x direction is still much greater than that in the e_y direction, the result is an elliptical distribution of tracers whose axes are aligned with the lab frame. Examining D_xx , as displayed in figure 7(a), we find that the bumpy case is very similar to the single-mode case, as the first component of the tensor is even with respect to k_y . Therefore, we are just adding mirror-symmetric components, which results in the same profile, with a weak effect of θ at short wavelengths, while at long wavelengths changes in θ dominate the behavior of D_xx . The largest increase to diffusion is observed for θ = 0, which results in an increase larger than 5% above classical Taylor dispersion for λ = 1.

Zoom In Zoom Out Reset image size

We observe a different behavior for D_yy for the two-mode case than in the single-mode case, see figure 7(b). There exists a local maximum of D_yy for $\theta = \pi/4$ and $\lambda \sim 1$, which decreases for intermediate λ before increasing again at large λ. It is important to note that these effects are still less than a 1% increase relative to Taylor's classical result near a planar wall. However, the presence of this non-monotonic behavior allows one to potentially optimize cross-streamline dispersion of tracer particles by tuning the surface topography.

3.3. Random surface (N modes)

Finally, we address the case of a randomly structured surface consisting of N modes. We restrict ourselves to considering surfaces where the surface shape is described by

$$\begin{align} H(\mathbf{r}_\parallel) = \frac{1}{2N}\sum_{i,j=1}^N \gamma_{ij}\left( \cos(\mathbf{k}_{ij}\cdot\mathbf{r}_\parallel)+ \sin(\mathbf{k}^{^{\prime}}_{ij}\cdot\mathbf{r}_\parallel)\right), \end{align} \tag{ 47 }$$

with $\mathbf{k}_{ij} = 2\pi [i,j,0]^{T}/(\lambda N)$, $\mathbf{k}^{^{\prime}}_{ij} = 2\pi [i,-j,0]^{\mathrm{T}}/(\lambda N)$, and coefficients γ_ij . To model a random surface the latter are drawn from a Gaussian distribution with mean $\langle \gamma_{ij}\rangle = 0$ and unit variance $\langle \gamma_{ij}\gamma_{nm}\rangle = \delta_{in}\delta_{jm}$, where the brackets indicate an average over many surface realizations (see figure 2(c) for a realization of this surface). This description creates a surface consisting of a random distribution of symmetrical bumps similar to the two-mode case with no directional bias, which is what we would expect of a randomly chosen surface, such as from sanding lines or machine tool marks. From our analysis above, we find that the overall effect of each random mode is additive in the diffusion tensor (46b ) with coefficients $\alpha_{ij} = \beta_{ij} = \gamma_{ij}$.

Our theory allows us to characterize the average dispersion coefficients, $\langle D_{xx}\rangle$ and $\langle D_{yy}\rangle$, of tracers suspended near random surfaces, by taking the average over many surface realizations. As the number of modes in the surface increases, the diffusion constants appear to approach a limit larger than the diffusivities in classical Taylor dispersion, see figure 8(a). $\langle D_{xx}\rangle$ exhibits over 30% increase in diffusivity due to the random surface while the increase of $\langle D_{yy}\rangle$ is around 5%. These results clearly indicate that the introduction of a random surface topology, even with relatively small amplitude, can significantly increase the dispersion of tracer particles advected near that surface structure. The enhancement for large N is to leading order limited by the contribution of the smallest wave vector: $\langle \mathbf{D}\rangle \lesssim \mathbf{D}^{(0)}+\epsilon^2\mathbf{D}^{(2)}(\mathbf{k}_{NN})$.

Zoom In Zoom Out Reset image size

We can also use this random surface formulation to examine the effect of varying Péclet number. As shown in figure 8(b), which plots the average dispersion coefficients for random surfaces with N = 20, both components of dispersion are close to that of classical Taylor dispersion at small Pe, indicating that classical Taylor dispersion is dominant relative to the effects of surface roughness. As Pe increases the normalized $\langle D_{xx}\rangle/D_{xx}^{(0)}$ increases slightly but remains bounded. This is due to the normalization factor from classical Taylor dispersion, which includes terms of both $\mathcal{O}(1)$ (corresponding to the effect of molecular diffusion) and $\mathcal{O}(\mathrm{Pe}^2)$ (corresponding to the effect of advection). As the contribution from surface roughness scales quadratically with Pe, due to the fact that it also arises from advection, $\langle D_{xx}\rangle/D_{xx}^{(0)}$ remains bounded for all Pe and approaches a 30% increase due to roughness.

In contrast, $\langle D_{yy}\rangle/D_{yy}^{(0)}$ increases for $\mathrm{Pe}\gtrsim\mathcal{O}(1/\epsilon^2)$. Mathematically, this is due to the fact that $D_{yy}^{(0)}$ only depends on molecular diffusion while the advection contribution of $D_{yy}^{(2)}$ scales with $\mathrm{Pe}^2$. As the Péclet number increases advective contributions become more important than molecular diffusion and so the dispersion driven by this advection becomes more pronounced than the effects of molecular diffusion alone.

The impact of surface properties on the classical Taylor dispersion phenomenon [13, 14] in narrow tubes has been subject of several works [20–23]. Here, we consider wide rectangular channels, which are confined along the vertical direction only, under shear flow and study the impact of structured surface topographies on the transport behavior of diffusing tracer particles. By generalizing Taylor's dispersion theory [13], we derive analytical predictions for the probability distribution of tracer particles at coarse-grained length and time scales, in the limit of small surface roughness.

Our theory, valid for all surfaces expandable in terms of a Fourier series, shows that surface structures can induce asymmetric dispersion along a principal direction that is different from the main flow direction. This effect is prominent for simple surface corrugations. For a Péclet number of 100, these corrugations can lead to enhancements of 20% along the principle direction while reducing the dispersion in the perpendicular direction. This can generate enhanced particle spreading relative to classical Taylor dispersion theory. Furthermore, our results reveal that surface corrugations that are tilted with respect to the forcing direction induce drift of tracer particles along the corrugations and reduce the drift along the forcing direction. Tracer particle drift along corrugations has been observed in experiments of non-Brownian particles in shear flow [41, 42, 45] and under gravity [47]. Our study, however, provides the first quantitative treatment of Brownian effects. We further show that for surfaces that do not have a directional bias the dispersion tensor becomes diagonal within the lab frame.

In addition, we prove that the contributions of individual wave modes of the surface to the large-scale diffusivity and drift velocity can be decoupled to leading order in the surface roughness. This allows us to consider randomly structured surfaces by summing over the contributions of each wave number. Within the limits of our theory, we find that randomly structured surfaces lacking a directional bias always enhance dispersion in both directions due to resulting disturbance flows. Under some conditions, we have observed that relative enhancement can reach an asymptotic value around 30% in the direction of the forcing as the number of random modes grows; this effect varies slightly with the Péclet number as demonstrated in figure 8(b).

It is important to note that our asymptotic theory relies on approximate flow fields, obtained via a domain perturbation method, and an extension of Taylor's dispersion theory. The stochastic simulations performed in this work merely corroborate our coarse-graining procedure, while the approximate flow fields have been validated separately using numerical simulations of the full hydrodynamic flows. We have made some progress in validating these results using finite-element methods, as shown in appendix D; however, future work should include extending these simulations to more surface geometries and refining the simulation code to allow for computation at higher resolutions within a reasonable time.

Our analytical framework can be readily extended to study Taylor dispersion in pressure-driven flows. To compute the roughness-induced drift and dispersion tensor, our asymptotic solutions for the main and roughness-induced flows need to be adjusted. While the first-order perturbation flows due to the roughness remain unaffected, the second-order contributions change, as the second term in the boundary condition (13c ) cannot be dropped.

Our predictions, irrespective of the details of the surface, show that the impact of surface structure for motion along the main flow direction is always orders of magnitude larger than for motion perpendicular to it. The latter is always an effect coming from the (weak) perturbation flow generated by the surface roughness. Thus, this particular feature comes from our assumption of small surface roughness and is expected to become altered for surface roughness of the order of the channel width, as in porous media or membranes, where restrictions can reduce transport along the flow direction. It would be interesting to study the transition from flow-enhanced dispersion to a potential slow-down due to entropic constraints in this context. Our results further indicate that small amplitude surface roughness can slow down tracer transport across the entire width of a channel while maintaining the same flow conditions, which may be useful for applications, such as ensuring a chemical reaction to complete.

C K acknowledges support from the Austrian Science Fund (FWF) via the Erwin Schrödinger fellowship (Grant No. J4321-N27). H A S and J V R thank the NSF for support via CBET-2127563 and MCB-1853602. The authors thank Steffen Hardt for insightful comments on the manuscript as well as his student Sebastian Dillen for insights gained in his MSc thesis work. We thank Luc Deike and Jiarong Wu for their assistance with Basilisk and Juan Santiago for helpful feedback.

Most of the data can be reproduced by the formulas given in the appendix of the manuscript. The data that support the findings of this study are available upon reasonable request from the authors.

A.1. 3D wavy surface: one mode

Using $\mathbf{k} = k[\cos\theta, \sin\theta,0]^{\mathrm{T}}$ the roughness-induced velocities of first order assume the form:

$$\begin{align} \hat{u}^{(1)}(z; \mathbf{k})&\equiv \frac{\mathrm{d}^2Q_1}{\mathrm{d} z^2}\nonumber\\[3pt] &= -\frac{\mathrm{csch}(k)}{2\left(2 k^2-\cosh (2 k)+1\right)}\Big(k \big(\cos ^2(\theta ) (2 k \sinh (k (z+1))\nonumber\\[3pt] &\quad +z \cosh (k (z-3)) -2(z-1) \cosh (k (z-1)) \nonumber\\[3pt] &\quad+(z-2) \cosh (k (z+1)) -4 k z \sinh (k) \cosh (k z)) \nonumber\\[3pt] &\quad-k (\cos (2 \theta)-3) \sinh (k-k z)\big) \sinh (k (z-3))\nonumber\\[3pt] &\quad+\sinh (k (z+1))+2 \sinh (k-k z)\Big) \end{align} \tag{ 1.1a }$$

$$\begin{align} \hat{v}^{(1)}(z; \mathbf{k})& \equiv\frac{\mathrm{d}^2Q_2}{\mathrm{d} z^2}\nonumber\\ & =\frac{k \sin (\theta ) \cos (\theta )}{2\left(2 k^2-\cosh (2 k)+1\right)} \Big(-z (4 k+\coth (k)) \cosh (k z)\nonumber\\[1.6pt] &\quad+z \mathrm{csch}(k) \cosh(k (z-3))\nonumber\\ &\quad+\sinh (k z) (4 k \coth (k)+3 z-4)\Big) \end{align} \tag{ 1.1b }$$

$$\begin{align} \hat{w}^{(1)}(z;\mathbf{k}) = \frac{k \cos (\theta ) (2 k (z\!-\!1) \sinh (k z)\!+\!z \cosh (k (z\!-\!2))\!-\!z \cosh (k z))}{2 k^2\!-\!\cosh (2 k)\!+\!1} .\end{align} \tag{ 1.1c }$$

The non-periodic second-order contributions evaluate to

$$\begin{align} \bar{u}_0^{(2)}(\mathbf{k}) \!=\! \frac{k \left(\sinh (2 k) (\cos (2 \theta)\!+\!3)\!-\!4 k \left(k \coth (k) \sin ^2(\theta)\!+\!2 \cos ^2(\theta)\right)\right)}{8 k^2\!-\!4 \cosh (2 k)+4}, \end{align} \tag{ 1.2a }$$

$$\begin{align} \bar{v}_0^{(2)}(\mathbf{k}) = \frac{k \sin (2 \theta) (\sinh (2 k)+2 k (k \coth (k)-2))}{8 k^2-4 \cosh (2 k)+4}, \end{align} \tag{ 1.2b }$$

$$\begin{align} \bar{w}_0^{(2)}(\mathbf{k}) =0 . \end{align} \tag{ 1.2c }$$

A.2. General, randomly structured surface

By linearity the first-order velocities induced by a surface of shape as in (45) assume the form

$$\begin{align} &\mathbf{u}^{(1)}(\mathbf{r}; \lambda, \{\alpha_{ij},\beta_{ij}\}) \nonumber\\ &\quad= \frac{1}{N} \sum_{i,j=1}^N\left(\alpha_{ij}\cos(\mathbf{k}_{ij}\cdot\mathbf{r}_\parallel)+\beta_{ij}\sin(\mathbf{k}_{ij}\cdot \mathbf{r}_\parallel)\right)\mathbf{Q}_{zz}(z; \mathbf{k}_{ij})\nonumber\\ &\qquad+\frac{1}{N} \sum_{i,j=1}^N\left(\alpha_{ij}\sin(\mathbf{k}_{ij}\cdot\mathbf{r}_\parallel)-\beta_{ij}\cos(\mathbf{k}_{ij}\cdot \mathbf{r}_\parallel)\right)\mathbf{k}_{ij}\cdot\mathbf{Q}_{z}(z;\mathbf{k}_{ij})\mathbf{e}_z. \end{align} \tag{ 1.3 }$$

The boundary conditions at z = 0 for the second-order velocity field are

$$\begin{align} \mathbf{u}^{(2)}(\mathbf{r}; \lambda, \{\alpha_{ij},\beta_{ij}\})\Bigl|_{z=0} &= -\frac{1}{2N^2}\sum_{i,j=1}^N\left((\alpha_{ij})^2+ (\beta_{ij})^2\right)\nonumber\\ &\quad\times\mathbf{Q}_{zzz}(z=0;\mathbf{k}_{ij}) + \mathrm{p.t.}, \end{align} \tag{ 1.4 }$$

where we observe that the coupling of different surface modes only contributes to the periodic terms 'p.t.'. The second-order velocities are of the form as in (41).

Noting that the dispersion matrix is symmetrical ($D_{12} = D_{21}$), we may express the analytic form of leading order correction to the dispersion matrix for a given wave number $\mathbf{k} = [k\cos(\theta),k\sin(\theta)]^{\mathrm{T}}$ as:

$$\begin{align} D^{(2)}_{11} &= -\frac{1}{1920 k^4 \left(-2 k^2+\cosh (2 k)-1\right)^2}\Bigg[1440 k^6 \sin ^4(\theta ) \mathrm{csch}^2\left(\frac{k}{2}\right)+k \sinh (4 k) \left(4 \left(k^4+105\right) \cos (2 \theta )-3 \left(55 \cos (4 \theta )-4 k^4+45\right)\right)\nonumber \\ &\quad +960 \left(4 k^2-5\right) k \sin ^2(2 \theta ) \sinh (k)-160 \sin ^2(\theta ) \cosh (3 k) \big(\left(k^2-36\right) \cos (2 \theta ) -k^2-12\big)-80 \sin ^2(\theta ) \cosh (4 k) \nonumber \\ &\quad\times\left(\left(k^2+18\right) \cos (2 \theta )-k^2+6\right)+64 \left(2 k^2+1\right) k^6 \cos ^2(\theta ) +32 \left(k^4-90\right) k^5 \sin ^2(\theta ) \coth \left(\frac{k}{2}\right)\nonumber \\ &\quad-320 k^2\cos ^2(\theta ) \left(\left(2 k^4+22 k^2+15\right)\cos (2 \theta )-2\left(k^4+11 k^2+3\right)\right)+80 k^2 \left(k^2 \left(\left(2 k^2-7\right) \cos (4 \theta )-3 \left(2k^2+9\right)\right)\right.\nonumber \\ &\quad\left.+2 \left(-2 k^4+7 k^2\!+\!6\right) \cos (2 \theta )\right)\nonumber \! -\!\frac{320 \left(5 k^2\!+\!48\right) k^4 \sin ^4(\theta )}{\cosh (k)\!+\!1} \!+\!160 \sin ^2(\theta ) \cosh (k) \Big(-8k^4\!-\!k^2\!+\!\left(8 k^4\!-\!191 k^2\!-\!36\right) \cos (2 \theta )-12\Big) \nonumber \\ &\quad+40 \Big(8 k^6+37 k^4-123 k^2 +4 \left(-4 k^6-9 k^4+85 k^2+15\right) \cos (2 \theta )+\left(8 k^6+23 k^4-133k^2-45\right) \cos (4 \theta )-15\Big) \nonumber \\ &\quad-2 k \sinh (2 k) \big(32 k^6+72 k^4-1350 k^2+12 \left(7 k^4+90 k^2+35\right) \cos (2 \theta ) +5 \left(4 k^4-90 k^2+159\right) \cos (4 \theta)-1095\big)\nonumber \\ &\quad-4 \cosh (2 k) \big(8 k^6+120 k^4 +645 k^2+5 \left(16 k^4+59 k^2-72\right) \cos (4 \theta )+4 \Big(2 k^6-30 k^4-145 k^2+120\Big) \cos (2 \theta ) -120\big)\nonumber \\ &\quad+32 k^3\sin ^2(\theta ) \tanh \left(\frac{k}{2}\right) \Big(-480 \cos (2 \theta )+k^6 -90 k^2-480\Big) -960 k \sin ^2(2 \theta ) \sinh (3 k)\Bigg] \end{align} \tag{ 2.1 }$$

$$\begin{align} D^{(2)}_{12} &= \frac{1}{7680 k^4 \left(-2 k^2+\cosh (2 k)-1\right)^2}\Bigg[\mathrm{csch}^2(k) \Bigg(2 \sin (2 \theta ) \Big(44 k^5 \sinh (4 k)-k^5 \sinh (6 k)-40 \left(k^2-12\right) \left(8 k^2+3\right) \cosh (3 k)\nonumber \\ &\quad+40 \left(k^2-12\right) \cosh (5 k)+20\Big(k^2 +6\Big) \cosh (6 k)+\left(16 k^4-1525\right) k^5 \sinh (2 k)+4 \big(2 k^6-30 k^4-155 k^2+60\big) \cosh (4 k) \nonumber \\ &\quad+80 \Big(54 k^3 \sinh ^3(k)+\left(k^2-12\right) \left(8k^4+4 k^2+1\right) -42 k \sinh ^5(k)\Big) \cosh (k)+4 \big(8 k^8+446 k^6 +1850 k^4+1495 k^2+420\big)\nonumber \\ &\quad-4 \left(8 k^8\!-\!112 k^6\!-\!100 k^4\!+\!1345 k^2\!+\!510\right) \cosh (2 k)\Big) \!-\!5 \sin (4\theta ) \Big(k \Big(4 k \Big(\big(378 \!-\!16 k^2\big) \cosh (3 k)\!-\!61 \cosh (4 k)\!+\!2 \cosh (5 k) \nonumber \\ &\quad+\cosh (6 k)+k \Big(-480 \sinh (k) +294 \sinh (2 k)-96 \sinh (3 k)+45 \sinh (4 k) -16 k\Big(\cosh (4 k) +k \sinh ^3(k)\cosh(k)\Big)\Big)\Big)\nonumber \\ &\quad-768 \sinh ^5\left(\frac{k}{2}\right) \cosh ^3\left(\frac{k}{2}\right) (-10 \cosh (k)+11 \cosh (2 k)+43)\Big) +16\Big(8 k^6-92 k^4 -95 k^2-36\Big) \cosh (k) +4 \Big(32 k^6-108 k^4\nonumber \\ &\quad-533 k^2-306\Big) \cosh (2 k)+4 \big(80 k^6+508 k^4+593 k^2+216 \cosh (3 k) +36 \cosh (4 k)-72 \cosh (5 k)+18\cosh (6 k)+252\big)\Big)\Bigg)\Bigg] \end{align} \tag{ 2.2 }$$

$$\begin{align} D^{(2)}_{22} &= \frac{1}{192 k^4 \left(-2 k^2+\cosh (2k)-1\right)^2} \Bigg[\sin ^2(2 \theta ) \mathrm{csch}^2(k) \Big(-80 k^6-508 k^4 -593 k^2+k \Big(k \big(2 \left(8 k^2-189\right) \cosh (3 k) \nonumber \\ &\quad +61 \cosh (4 k)-2 \cosh (5 k)-\cosh (6 k)+k \big(480 \sinh(k)-294 \sinh (2 k)+96 \sinh (3 k)-45 \sinh (4 k) \nonumber \\ &\quad+16 k \left(\cosh (4 k)+k \sinh^3(k) \cosh(k)\right)\big)\big) +192 \sinh ^5\left(\frac{k}{2}\right) \cosh^3\left(\frac{k}{2}\right) (-10 \cosh (k)+11 \cosh (2 k)+43)\Big)\nonumber \\ &\quad+4 \left(-8 k^6+92 k^4+95 k^2+36\right) \cosh (k)+\Big(-32 k^6+108 k^4+533 k^2+306\Big) \cosh (2 k)-216 \cosh (3k)\nonumber \\ &\quad -36 \cosh (4 k)+72 \cosh (5 k) -18 \cosh (6 k)-252\Big)\Bigg] . \end{align} \tag{ 2.3 }$$

To validate our coarse-graining procedure to compute analytically the effective drift and diffusivities, we perform computer simulations of Brownian particles in shear flow near corrugated surfaces using the flow fields obtained via the domain perturbation method as input. Therefore, we assume that the particles are neutrally buoyant and small compared to the surface corrugations, so that they can be approximated by point particles, which, in addition to diffusion, follow the flow stream lines. The equation of motion reads

$$\begin{align} \frac{\mathrm{d} \mathbf{r}(t)}{\mathrm{d} t} &= \mathbf{u}(\mathbf{r}(t))+\sqrt{2D_0}\boldsymbol{\xi}(t), \end{align} \tag{ 3.1 }$$

where we have modeled Brownian motion in terms of an independent Gaussian white noise process $\boldsymbol{\xi}(t)$ of zero mean and delta correlated variance $\langle \xi_{i}(t)\xi_{j}(t^{^{\prime}})\rangle = \delta_{ij}\delta(t-t^{^{\prime}})$ with $i,j = 1,2,3$. We use the flow field u from our perturbation approach (as outlined earlier) and employ reflective boundary conditions at the walls. Further, we non-dimensionalize via

$$\begin{align} \mathbf{r}&= h\hat{\mathbf{r}} \qquad t= \frac{h^2}{D_0} \hat{t} \qquad \mathbf{u}= u_0\hat{\mathbf{u}}, \end{align} \tag{ 3.2 }$$

where $\hat{\mathbf{r}}$, $\hat{t}$, and $\hat{\mathbf{u}}$ are dimensionless quantities, and discretize the equation of motion

$$\begin{align} \hat{\mathbf{r}}(\hat{t}+\Delta\hat{t})&= \hat{\mathbf{r}}(\hat{t})+ \mathrm{Pe}\, \hat{\mathbf{u}}(\hat{\mathbf{r}}(\hat{t}))\Delta\hat{t} +\sqrt{2 \Delta \hat{t}}\ \mathbf{N}_\xi. \end{align} \tag{ 3.3 }$$

Here, $\Delta \hat{t}$ denotes the discretized time step and $\mathbf{N}_\xi$ is a vector with independent and normally distributed random variables with zero mean and unit variance. To corroborate our theoretical predictions with simulations, we use a Péclet number of Pe = 100, choose $H(\mathbf{r}_\parallel) = \cos(\mathbf{k}\cdot\mathbf{r}_\parallel)$ with $\mathbf{k} = (2\pi/\lambda)[\cos\theta,\sin\theta]^{\mathrm{T}}$ and $\theta = \pi/4$, and vary the wavelength $\lambda\in[0.5,12.5]/\sqrt{2}$. We evaluate 10⁴ trajectories using simulations performed with Python and using a time-step of $\Delta\hat{t} = 10^{-3}$.

We then compute the average displacements along $\langle \Delta x\rangle$ and $\langle \Delta y\rangle$ to extract the drift at long times via $\langle \Delta x\rangle \simeq u_{\mathrm{eff}}t$ and $\langle \Delta y\rangle \simeq v_{\mathrm{eff}}t$, see figure C1. Subtracting the drift from the displacements and realigning the displacements along the main direction of diffusion, we compute the mean-square displacements and extract the long-time diffusivities along the main direction of diffusion (and hence the dispersion matrix). The variances of displacements along the main directions of diffusion, $\mathrm{var}_1(t)$ and $\mathrm{var}_2(t)$, are shown in figure C2. They scale linearly in time ${\sim}t$ at long times, which, in principle, allows us to extract the eigenvalues of the dispersion matrix via $\mathrm{var}_1(t)\simeq 2\Lambda_1 t$ and $\mathrm{var}_2(t)\simeq 2\Lambda_2 t$. While we can reliably measure Λ₁, Λ₂ is smaller than the molecular diffusivity (see figure 5) and within our simulations is not resolvable. We note that we tested the number of realizations N and found that increasing it further does not significantly alter our results.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Our theory for the effective dispersion of solutes near rough surfaces relies on an asymptotic approximation for the velocity field. In order to provide some reasonable assurance that this approximation captures the important qualities of the flow, we numerically solved the Stokes equations in several channel geometries using the centered Navier–Stokes solver in the Basilisk package [56, 57], simulating the flow in one periodic unit of the surface.

We make use of Basilisk's adaptive tree grids [58] to limit the computational cost of our simulations. Our simulations employ a maximum grid resolution of $N = 2^9 = 512$ points in both horizontal axes. The vertical axis in z is constrained to give the appropriate ratio $\lambda/h$. We simulated the flow field for a single corrugated surface with θ fixed at $\pi/4$ but with $\lambda/h = 2$, 3, 4, 5, and 10 with ε varied between 0.01 and 0.2. We did not run simulations for $\lambda/h\lt 2$ due to issues with implementing the domain properly in the Basilisk environment. Our simulation iterated on the velocity field until the change in the velocity u_x was less than ε³ for each step. We then compared these numerical flow fields with the approximate flow field developed using our asymptotic approach. A qualitative comparison of these flow fields for $\lambda/h = 4$ and ε = 0.1 is shown in figure D1. Note that in this figure the velocity u_y has been enhanced for clarity. Both the simulation and analytic results agree qualitatively, with similar behavior as the flow encounters the peaks and valleys.

Zoom In Zoom Out Reset image size

In addition to a qualitative comparison, we can quantitatively compare our simulation results to our theory. As our theory predicts the effective motion on a coarse-grained scale, we evaluate the velocity averaged across one periodic cell. First, similar to our derivation using the approximate velocity field, we calculate the height averaged horizontal velocity field, $\mathbf{u}_\parallel(x,y) = \langle [u,v]^\mathrm{T}\rangle$. Next, we can average this velocity across one periodic cell C by computing

$$\begin{align} \mathbf{u}_{\mathrm{eff,sim}}=\frac{1}{C}\int_C \mathbf{u}_\parallel(x,y) \,\mathrm{d}\mathbf{r}_\parallel .\end{align} \tag{ 4.1 }$$

We compare this effective velocity field with that given by our asymptotic theory (33a ). We report the percent error of the asymptotic theory relative to the simulation results in figure D2. We observe in (a) that the error in u_x is relatively small (less than 2%) across all ε and λ. Generally the error increases as we increase the surface amplitude, consistent with the asymptotic nature of our theory, which is valid only for $\epsilon\ll 1$.

Zoom In Zoom Out Reset image size

The error in u_y is comparatively larger. However, note that it is never greater than 50% across the entire range of parameters tested. Given that u_y only enters the average velocity at $\mathcal{O}(\epsilon^2)$ we believe having less than a factor of 2 difference for all of our simulation results is quite good. Generally it appears that the error for longer wavelengths is lower, especially at higher ε. This is consistent with our theory, which is more approximately correct in the limit of long wavelengths. We generally observe the same trends as in (a), with the error mostly increasing as ε increases. At very low ε this general trend breaks down. However, this may be in part due to a lack of sufficient resolution to resolve these small velocities.

To test the resolution hypothesis, we performed an additional set of computations for $\lambda/h = 10$ at two different maximum resolutions, which are shown in figure D3. Except for an anomaly at ε = 0.15, where our simulation ended up very close to the analytic approximation, the higher resolution generally has lower error, as with more resolution we are better able to capture $\mathcal{O}(\epsilon^2)$ effects. At this higher resolution, our simulations and theory differ by only 3% even at ε = 0.2, which is quite large for an asymptotic approximation. We expect similar decreases in the overall error in u_y , given the computational resources to rerun the other wavelengths at this higher resolution.

Zoom In Zoom Out Reset image size

Based on these simulation results, we have reasonable confidence that our asymptotic theory can predict the coarse-grained behavior of the flow for $\lambda/h \geqslant 2$ and $\epsilon\ll 1$. Our theory captures both the qualitative behavior of the flow within one periodic cell as well as provides a reasonable approximation to the flow averaged over a periodic cell. While we have not conducted simulations for other surface geometries, including bumpy surfaces, given the linear nature of flows at low Re we are confident that our asymptotic theory can provide approximate descriptions of the drift and dispersion of particles over surfaces that conform our restrictions on the magnitude of λ and ε.