Roughness, Inertia, and Diffusion Effects on Anomalous Transport in Rough Channels

TL;DR

This paper investigates how roughness, inertia, and diffusion influence anomalous tracer transport in rough channels through simulations and stochastic modeling, revealing universal scaling laws and the importance of Lagrangian velocity statistics.

Contribution

It introduces a comprehensive study combining simulations and stochastic models to understand the combined effects of roughness, inertia, and diffusion on transport in rough channels, including universal scaling laws.

Findings

Pe controls first-passage time distributions

Roughness effects become prominent at higher Re

Universal power-law scaling at infinite Pe

Abstract

We study how the complex interplay between roughness, inertia, and diffusion controls the tracer transport in rough channels. We first simulate the flow and tracer transport over wide ranges of channel roughness, Reynolds number ($Re$), and P\'eclet number ($Pe$) observable in nature. $Pe$ exerts a first-order control on first-passage time distributions, and the effect of roughness on the tracer transport becomes evident with the increase in $Re$. The interplay between the roughness and $Re$ causes eddy flows, which intensify or suppress anomalous transport depending on $Pe$. At infinite $Pe$, the late-time scaling follows a universal power-law scaling, which is explained by conducting a scaling analysis. With extensive numerical simulations and stochastic modeling, we show that the roughness, inertia, and diffusion effects are encoded in Lagrangian velocity statistics represented by velocity distribution and correlation. We finally predict the anomalous transport using a stochastic model that considers the Lagrangian velocity statistics.

Equations

$$f_e(u)=-\frac{2}{b}\frac{dy}{du}=\frac{1}{6q\sqrt{\frac{1}{4b^2}-\frac{u}{6qb}}}.$$

$$f_t(t)\propto \frac{1}{t^3}\cdot \frac{1}{\sqrt{c_1-c_2/t}},$$

$$x^{(n+1)}=x^{(n)}+\Delta x,t^{(n+1)}=t^{(n)}+\frac{\Delta x}{v^{(n)}},$$

Full text

Fluid flow and mass transport in rough channels are ubiquitous phenomena occurring in numerous engineering applications and natural processes including microfluidics, biomedical devices, heat exchangers, and fractured geological media (Stroock et al., 2002; Squires and Quake, 2005; Sackmann et al., 2014; Morini, 2004; Berkowitz, 2002). Since Taylor’s seminal work on solute dispersion in shear flows (Taylor, 1953), the concept of effective dispersion has proven to be very powerful, and many studies have proposed various methods to quantify the effective dispersion in more complex flow fields (Aris, 1956; Koplik et al., 1988; Plumb and Whitaker, 1988; Frankel and Brenner, 1989; Sarracino et al., 2016; Detwiler et al., 2000; Gelhar and Axness, 1983; Drazer and Koplik, 2001; Bouquain et al., 2012; Måløy et al., 1988). Although these effective dispersion approaches are useful in many applications, the asymptotic regime is often not reached, and anomalous behaviors are widely observed  (Bijeljic and Blunt, 2006; Neuman and Tartakovsky, 2009; de Anna et al., 2013; Le Borgne et al., 2011; Kang et al., 2015; Nissan and Berkowitz, 2018; Cushman and Ginn, 2000; Benson et al., 2001; Singha et al., 2007; Berkowitz and Scher, 1997). The solute transport in finite-length channel flows at a relatively high Péclet number ($Pe$) is a representative example (Giona et al., 2009).

Recent studies on tracer transport through channel flow systems focused on the effects of eddies: eddies are shown to trap solute particles and cause heavy tailing of tracer breakthrough curves (Boutt et al., 2006; Cardenas et al., 2009; Bolster et al., 2014; Wang and Cardenas, 2014). However, the effects of eddies were investigated for a specific channel geometry and over narrow ranges of Reynolds number ($Re$) and $Pe$. In practice, the roughness, inertia, and diffusion effects can vary over wide ranges, leading to complex transport behaviors (Thompson and Brown, 1991; Lee et al., 2015; Nissan and Berkowitz, 2018, 2019). For example, the trapping mechanism induced by the diffusion of solute particles into eddies should strongly depend on $Pe$. Currently, we do not have a mechanistic understanding of how the complex interplay between channel roughness, $Re$, and $Pe$ controls the anomalous transport in rough channels.

In this study, we elucidate how this complex interplay governs the tracer transport in rough channels by conducting extensive numerical simulations and stochastic modeling. We consider self-affine rough walls, as rough surfaces in nature are often found to be self-affine (Mandelbrot, 1983; Kertesz et al., 1993; Ponson et al., 2006). Self-affine surfaces are scale invariant in that the height of the rough surface at a spatial location $x$ can be described using $z(x)=\lambda^{-H}z(\lambda x)$, where $H$ is the Hurst exponent that characterizes the roughness of the surface. We investigate the effects of roughness on the flow and transport by varying the Hurst exponent ($H$) in the range of $0.7$ – $0.9$, which is consistent with that observed in nature (Berkowitz, 2002; Bouchaud et al., 1990; Drazer et al., 2004). We use the successive random addition algorithm (Voss, 1988; Liu et al., 2004) to generate rough surfaces of length $L=10$ cm. The generated rough surfaces are duplicated and detached to have a constant aperture of $1$ mm.

We simulate a Newtonian fluid flow by solving steady-state incompressible Navier–Stokes equations for rough channels using the finite volume method (OPENFOAM, 2011). A constant flux boundary condition is imposed on the left boundary of the channel, and a zero-pressure gradient boundary condition is imposed on the right boundary. We discretize the channel with a resolution of $0.002$ mm, yielding $50,000\times 500$ grid cells within the channel domain. The inertial regimes are quantified with Reynolds number defined as $Re=\frac{ub}{\nu}$, where $b$ is the aperture, and $\nu$ is the kinematic viscosity of the fluid. We consider various laminar flow regimes at seven different Reynolds numbers: $Re=[1,10,20,40,60,80,100]$.

We simulate solute transport using a particle tracking method (Bijeljic et al., 2011). In this method, the advective transport is simulated using a streamline-based particle tracking algorithm that considers no-slip boundary conditions at solid–fluid interfaces (Mostaghimi et al., 2012). The diffusive displacement is modeled using a random walk method. The Lagrangian approach is free of numerical dispersion and can be used to accurately simulate particle transport at high $Pe$ regimes. We inject $10^{4}$ particles in each realization with a flux-weighted line injection, and the $x$-direnctional travel distance is $8$ cm. We investigate the effects of diffusion on the transport by varying the Péclet number defined as $Pe=\frac{ub}{2D}$, where $D$ is the molecular diffusivity. We consider five different $Pe$ regimes: $Pe=[10^{2},10^{3},10^{4},10^{5},\infty]$. We choose $Pe$ and $Re$ ranges such that they cover the observable Schmidt number, $Sc=\frac{Pe}{Re}$, in nature (Stroock et al., 2002; Delgado, 2006). In summary, we simulate the flow and transport by varying $H$, $Re$, and $Pe$ independently over wide ranges: $H=[0.7,0.8,0.9]$, $Re=[1,10,20,40,60,80,100]$, and $Pe=[10^{2},10^{3},10^{4},10^{5},\infty]$. The independent variation in $H$, $Re$, and $Pe$ helps discern the role of each factor. We simulate an ensemble of $10$ realizations for each combination of $H$, $Re$, and $Pe$.

Simulation Results.

To highlight the importance of understanding the complex interplay, we first present the effects of $H$ and $Pe$ on the tracer transport at $Re=100$. As shown in Figs. 1(a) and (b), an increase in the roughness increases the size and frequency of the eddies. The role of eddies on the tracer transport is sensitive to $Pe$. At $Pe=10^{3}$, tracers readily diffuse into eddies, and the eddies trap these tracers. Interestingly, the eddies play the opposite role at $Pe=10^{5}$. At $Pe=10^{5}$, the tracers can no longer easily enter the eddy zones, and the eddies instead act as slip boundaries that help in transporting the tracers near the channel walls (see supporting videos). This is clearly observed in the projected concentration profiles at a pore volume injection (PVI) of $0.3$ (Figs. 1(c) and (d)). The tailing in the profile becomes stronger with the increase in the roughness at $Pe=10^{3}$; however, the trend is opposite at $Pe=10^{5}$.

The normalized breakthrough curves or first-passage time distributions (FPTDs) for the combinations of $H=[0.7,0.8,0.9]$, $Re=[1,100]$, and $Pe=[10^{2},10^{3},10^{4},10^{5},\infty]$ are shown in Fig. 2. We show FPTDs up to $10$ PVI, though more than $90\%$ of the particles breakthrough within $1$ PVI. First, note that $Pe$ has a first-order control on the breakthrough curve shapes. The effects of roughness only become evident at higher $Re$. We focus on two key characteristics of the anomalous transport: early arrival and late-time tailing in the breakthrough curves. We observe an enhanced early arrival with the increase in the roughness and $Re$. A high roughness coupled with significant inertial effects (high $Re$) causes eddies to enlarge. The eddies narrow the mobile zone, or the main flow channel, and the decreased cross-sectional area of the main flow channel leads to an increase in the velocity along the main flow channel (Zhou et al., 2019). Fig. 1 shows this phenomenon. This explains the anomalously early arrival of the tracers, shown in Fig. 2, at high roughness and high $Re$.

The late-time behavior of the FPTDs is determined by the particles traveling through low-velocity regions. No-slip boundary conditions and eddies create low-velocity regions near the channel walls. To analyze the eddy trapping effects, we first delineate eddy zones based on the mass balance principle (Zhou et al., 2019) and quantify the effects of eddies on the first passage time by subtracting the total time spent in the eddies from the first passage time. The subtracted residence time distributions at $H=0.7$ and $Re=100$ clearly demonstrate eddy trapping effects (insets in Figs. 2(f)–(i)). However, with an increase in $Pe$, the trapping effects due to the eddies are observed much later. At $Pe=10^{5}$, the eddies appear to not delay but rather aid the transit of the particles until $\sim 3$ PVI. This is because the particles cannot easily enter the eddies due to the limited diffusion, and the eddies act as slip boundaries. Therefore, most of the tracers (approximately $99\%$) breakthrough before $3$ PVI. Consequently, the FPTDs in the smoother cases (e.g., $H=[0.8,0.9]$) show stronger tailing than that in the $H=0.7$ case until $t\sim 3$ PVI (Fig. 2(i)). We observe that although fewer particles are captured at high $Pe$ regimes, once captured, the trapped particles stay longer inside the eddies compared with that observed at lower $Pe$ regimes. This explains the change in scaling at PVI $\sim 3$ in Fig. 2(i).

Finally, at $Pe=\infty$, the late-time scaling of the FPTDs shows a universal power law scaling of $t^{-3}$ for all combinations of $Pe$, $Re$, and $H$ (Figs. 2(e) and (j)). We perform a scaling analysis to explain the observed universal power-law scaling. The low-velocity regions should determine the late-time scaling, and the tracers at $Pe=\infty$ cannot enter the eddies. Thus, the low-velocity distribution at $Pe=\infty$ should be determined by the no-slip boundary conditions. For a Poiseuille flow with an aperture $b$, the velocity profile across the channel follows the parabolic equation u(y)=6qb^{-3}\big{[}(b/2)^{2}-y^{2}\big{]}, where $q$ is the constant influx into the channel. The Eulerian velocity probability density function (PDF) corresponding to the parabolic profile is given as follows.

[TABLE]

The PDF of the Lagrangian velocities is related to the Eulerian velocity PDF through flux weighting as $f_{\mathcal{L}}(u)=\frac{uf_{e}(u)}{\int du\,uf_{e}(u)}$ (Dentz et al., 2016). The late-time scaling of the FPTD, $f_{t}(t)$, should be determined by the distribution of the low velocities in $f_{\mathcal{L}}(u)$. From $t\propto L/u$ and $f_{\mathcal{L}}(u)$, the late-time scaling of the FPTD follows,

[TABLE]

where $c_{1}=(4b^{2})^{-1}$ and $c_{2}=(6qb)^{-1}$. For large $t$, we obtain $f_{t}(t)\propto t^{-3}$. This confirms that the no-slip boundary condition indeed governs the late-time scaling at infinite $Pe$ regardless of the roughness and $Re$.

Lagrangian Velocity Statistics.

Recent studies reported that the pore structure governs the velocity PDF and that pore-scale velocities are strongly correlated (Bijeljic et al., 2011; de Anna et al., 2013; Datta et al., 2013; Kang et al., 2014; de Anna et al., 2017; Carrel et al., 2018; Crevacore et al., 2016). The transport in porous media has been successfully characterized by considering the interplay between the velocity PDF and the velocity correlation (Le Borgne et al., 2008; Dentz and Bolster, 2010; Kang et al., 2011; Bolster et al., 2014; Kang et al., 2015; Dentz et al., 2016). For channel flows, the roughness of the wall can lead to a significant difference in velocity between the main channel flow and the near-wall low-velocity zones (Méheust and Schmittbuhl, 2000; Cardenas et al., 2007; Lee et al., 2015). Moreover, the Lagrangian velocities sampled along a particle trajectory in the channel flows should retain the memory of the prior velocities because of the mass conservation constraint.

We characterize the Lagrangian velocity statistics using the Lagrangian velocity PDF (velocity distribution) and velocity correlation. We first quantify the motion of the solute particles in equidistance, $\Delta x$, in the mean flow direction. The velocity distribution is characterized using the PDF of the transition time $\tau=\frac{\Delta x}{v}$, where $v$ is the average Lagrangian velocity over $\Delta x$. The transition times are sampled at every $\Delta x=1$ mm from all the particle trajectories. Herein, we refer to the PDF of the transition times as the transition time distribution (TTD). We characterize the Lagrangian velocity correlation by quantifying the velocity correlation lengths conditional to initial velocity values. At preasymptotic regimes, velocity correlation can strongly depend on the initial velocity values (Le Borgne et al., 2007; Kang et al., 2019). We classify the initial Lagrangian velocities into $10$ classes, wherein each class is equidistantly spaced in a $\log$-scale. Based on the initial $\log$-velocity values, the particles are assigned into one of the $10$ classes, $i=[1,\dots,10]$, where $i=10$ is the class with the highest velocities. We estimate the characteristic correlation length for each class $i$ as ${\ell}_{i}=\int_{0}^{L}C(x|i)/C(0|i)dx$, where C(x|i)=\int_{-\infty}^{\infty}\big{|}P(\log v|i,x)-P(\log v,L)\big{|}d\log v (Le Borgne et al., 2007). Here, $P(\log v|i,x)$ is the conditional $\log$-velocity distribution for particles belonging to class $i$, and $P(\log v,L)$ is the marginal $\log$-velocity distribution at the outlet.

The underlying mechanisms of the early arrival and late-time tailing are accurately captured in the TTDs (Fig. 3) and velocity correlation (Fig. 4). The probability of having short transition times (high-velocity) increases with the increase in the roughness and $Re$ (Fig. 3). This implies that the flow channeling effect is encoded in the TTDs. The channeling effect is also encoded in the correlation length (Fig. 4). The correlation length of the velocity class $i=10$ controls the early arrival of the tracers, and we observe that, for $i=10$, the correlation length in the $H=0.7$ case is higher than those in the smoother cases as $Re$ increases (Fig. 4(c)). Eddies, which enlarge with the increase in $Re$ (inset of Fig. 4(c)), reduce the cross-sectional flow area and cause a strong preferential path (Fig. 1). The strong preferential flow reduces the velocity fluctuation and thereby increases the velocity correlation. The higher probability of having higher velocities combined with a strong velocity correlation explains the early arrival behavior at $H=0.7$.

Similar to the early arrival behavior, the late-time behavior is also encoded in the TTDs and velocity correlation. To quantify the tail of the TTDs, we fit the tail to a generalized Pareto distribution $G_{\xi,\sigma,\theta}(\tau)=\frac{1}{\sigma}(1+\frac{\xi(\tau-\theta)}{\sigma})^{-1-\frac{1}{\xi}}$ for $\tau>\theta$ when $\xi>0$ or for $\theta<\tau<\theta-\sigma/\xi$ when $\xi<0$. The three parameters $\xi$, $\sigma$, and $\theta$ are estimated using the maximum likelihood method. $\xi$, which is often called the tail index, gives an indication of the heaviness of a tail (the greater the value of $\xi$, the heavier the tail) (Coles, 2001). The threshold $\theta$ is set as the $0.9$ quantile of the cumulative distribution of the transition times, $\tau_{0.9}$. The insets in Fig. 3 show the estimates of the tail index. The method of maximum likelihood finds the parameter values that maximize the likelihood function $\mathcal{L}(\xi,\sigma,\theta\,|\mathbf{t})$ where $\mathbf{t}=\{\tau\,|\tau\geq\tau_{0.9}\}$.

When eddy flows are strong ($H=0.7$, $Re=100$), the tail index significantly increases with the increase in $Pe$ from $100$ to $10^{3}$ (inset of Fig. 3(d)). This is due to the stronger eddy trapping effect at $Pe=10^{3}$. At $Pe=100$, the particles easily enter the eddies but also easily exit. When $Pe$ increases beyond $10^{3}$, the tail index decreases significantly at $H=0.7$, confirming that the particles cannot efficiently sample velocities in the eddy zones. Moreover, at high $Pe$ regimes, the particles traveling near the walls will alternate between low (no-slip wall) and relatively high (eddy interface) velocities, thus decreasing the velocity correlation. For the velocity class $i=1$, the correlation length indeed decreases with the increase in $Re$ at $Pe=10^{5}$ and $H=0.7$ (Fig. 4(d)). This trend is well correlated with the increase in the relative eddy area as a function of $Re$ (inset of Fig. 4(c)). The decreased tailing in the TTDs and the loss of velocity correlation for the velocity class $i=1$ at high $Pe$, high roughness regimes explain how eddies suppress the anomalous transport. Finally, the role of $Pe$ in determining the overall shape of the FPTD is also evident from the TTDs and velocity correlation. $Pe$ has a strong effect on both TTD and velocity correlation, whereas the roughness at $Re=1$ has little effect on the TTD (inset of Fig. 3(b)) and velocity correlation (Fig. 4(a)).

Upscaled Stochastic Transport Model.

We hypothesize that the complex $H$–$Re$–$Pe$ interplay in rough channels is encoded in the velocity distribution and correlation, which in turn determine the effective transport. To test our hypothesis, we quantify the effective particle transport using an upscaled model that only takes the velocity distribution and correlation as input parameters. The continuous time random walk model with one-step correlation, often referred to as the spatial Markov model, has been successfully applied to predict anomalous transport across scales (Le Borgne et al., 2008; Kang et al., 2011; de Anna et al., 2013; Kang et al., 2014, 2015).

We test the hypothesis by running the spatial Markov model as an effective transport model. The effective particle transport can be characterized using the Langevin equations,

[TABLE]

where $\{\tau^{(n)}\}_{n=0}^{L/\Delta x}$ is a series of transition times. The transition time, $\tau^{(n)}$, in Eq. 3 is modeled as a Markov chain, whose transitions can be characterized using a transition matrix. We sample the Lagrangian velocity transitions at every $\Delta x=5$ mm, which is smaller than the estimated correlation length (Fig. 4). We then construct the $10\times 10$ transition matrix with the $10$ classes equidistantly spaced in a $\log$-scale. The spatial Markov model accurately predicts the projected concentrations and FPTDs for all combinations of $H$, $Re$, and $Pe$ (Figs. 1(c, d) and Fig. 2). This supports the hypothesis that the $H$–$Re$–$Pe$ interplay is effectively encoded in the velocity distribution and correlation and that they are sufficient for quantifying the effective transport.

In conclusion, we successfully established a mechanistic understanding of how the complex $H$–$Re$–$Pe$ interplay determines the effective tracer transport in rough channel flows. Wide ranges of $H$, $Re$, and $Pe$ observable in nature were investigated, and we found that eddy flows can either induce or nonintuitively suppress the anomalous transport depending on the $Pe$ value. Based on the improved understanding, we proposed a predictive model. For $Pe=\infty$ regime, we observed the universal power-law scaling, $t^{-3}$, in FPTDs and explained it with a scaling analysis.

.1 Acknowledgments

The authors gratefully acknowledge support from NSF via a grant EAR1813526. PKK also acknowledges the College of Science & Engineering at the University of Minnesota and the George and Orpha Gibson Endowment for its generous support of Hydrogeology. We thank the Minnesota Supercomputing Institute (MSI) at the University of Minnesota for computational resources and support.