Structure induced laminar vortices control anomalous dispersion in porous media

Natural porous systems, such as soil, membranes, and biological tissues comprise disordered structures characterized by dead-end pores connected to a network of percolating channels. The release and dispersion of particles, solutes, and microorganisms from such features is key for a broad range of environmental and medical applications including soil remediation, filtration and drug delivery. Yet, owing to the stagnant and opaque nature of these disordered systems, the role of microscopic structure and flow on the dispersion of particles and solutes remains poorly understood. Here, we use a microfluidic model system that features a pore structure characterized by distributed dead-ends to determine how particles are transported, retained and dispersed. We observe strong tailing of arrival time distributions at the outlet of the medium characterized by power-law decay with an exponent of 2/3. Using numerical simulations and an analytical model, we link this behavior to particles initially located within dead-end pores, and explain the tailing exponent with a hopping across and rolling along the streamlines of vortices within dead-end pores. We quantify such anomalous dispersal by a stochastic model that predicts the full evolution of arrival times. Our results demonstrate how microscopic flow structures can impact macroscopic particle transport.

∂t for an instantaneous solute pulse c(x, t = 0) = δ(t) at the upper boundary of the roll. The solute advects and diffuses  To show this quantitatively, we first note that during a rapid initial phase, solute is spread uniformly along and 35 between streamlines due to shear dispersion [2,5]. In order to formalize this observation, we consider the streamline 36 coordinates (ψ, φ), where ψ(x) is the stream function, so that the superscript denotes the transpose and ∇ψ( Thus, the 38 coordinate ψ labels the streamline and φ denotes the distance along the streamline. We write c( and use for simplicity of notation the same letter for concentration in Cartesian and streamline coordinates. The 40 concentration is constant along streamlines, this means that c(x, t) can be represented by its average along a closed 41 streamline as .
The integral is along the closed contour ψ = ψ(x). Note that division by v in the integrand accounts for the variable 43 spacing between streamlines because the amount of solute between streamlines is constant. We now average the 44 advection-diffusion equation (2) along a streamline and use (4), which gives where we used the identity where we used the fact that v = |∇ψ|. Using the Gauss theorem, the integral term on the right side of Eq. (5) can be 47 written as Thus, Eq. (5) becomes where we defined 50 C(ψ) = dφv(ψ, φ).
In order to continue we need to specify T (ψ), which is the advection time along a closed contour ψ and C(ψ), boundary. There, we approximate the streamfunction by ψ(x) = σx 2 /2, where σ is the shear rate. This implies that 55 v(ψ) ∼ √ ψ and T (ψ) ∼ ψ −1/2 . The circulation C(ψ) is dominated by sections along the horizontal boundaries, close 56 to which the velocity field is approximately constant. This means C(ψ) ≈ C 0 . Equation (8) can then be written as where D 0 = DC 0 . We note that (10) is equivalent to the Langevin equation where ξ(t) is a Gaussian white noise. We define the operational time ds = dt/T (ψ) and write Thus, the time t(s) is given by where the angular brackets denote the noise average. As ψ(s) describes a Brownian motion, we obtain the scaling which implies that 62 t(s) ∼ s 3/4 , i.e., time increases slower than linearly with operational time.

63
A. Trapping times 64 We want to determine the distribution of trapping times, which implies the distribution of t(s) such that s is the operational time to escape the roll, this means to return to the streamline ψ = 0 after entering the role at s = 0.

66
This means we need to determine the return time PDF. As this problem is not well-defined in continous space, we This means, it scales as g ∼ s −3/2 . The trapping time distribution φ(t) is obtained from g(s) by variable transform This behavior is valid as long as the solute is distributed along the outer streamlines of the vortex. For times larger 74 than the diffusion time τ D = R 2 /D across the extension R of the roll, φ(t) is decaying exponentially fast.

B. Residence times 76
In order to obtain the residence time distribution, we consider a uniform distribution of particles inside the role.

77
Note that a uniform distribution p 0 (x) = 1 in equidistant space x implies in streamline coordinates that Thus, we obtain the residence time distribution by integration of (16) because it describes the residence time distri-79 bution of particles that originate from ψ. Thus, we obtain Thus, the residence time distribution behaves as The scaling of Φ is obtained substituting eq. (15) and eq. (20) into eq. (21), leading to This is the scaling observed for the residence time from the experimental and numerical data.

84
The breakthrough curves are fully determined by the particle retentention at the first step. Thus, we model the 85 breakthrough curves through the CTRW With the distribution p 0 (x, t) of initial particle positions and times (x 0 , t 0 ) given by where α is the proportion of particles initially located in the DEPs, and I(·) is an impulse function, which is 1 if 88 its argument is true and 0 else. This means that the initial positions x 0 are uniformly distributed in space, and the 89 residence time within a DEP at the initial position x 0 is given by a random time t 0 . Its distribution is given by 90 Φ(t), which is the residence time within a DEP of depth Λλ m and so it is the Γ distribution observed with our direct whose -2/3 scaling was we determined in the previous section. The time increment τ n is set constant and equal to 93 the advection time over the distance c as τ n = ∆t = c /u m with u m the average velocity in percolating pores. The 94 random displacement ξ n is a Gaussian random variable with 0 mean and unit variance, and D * is the fitted dispersion 95 coefficient of the transmitting pore network.

97
In order to determine the breakthrough curve at all times for this CTRW model, we first consider the BTC for 98 particles that originate at a distance x from the outlet. If the particle is not in a DEP, the time to reach the outlet 99 is distributed according to the inverse Gaussian distribution [6] because transport in this CTRW model is given by advection and dispersion. If the particle is initially trapped, the 101 time to arrive at the outlet is t 0 + t f , where t f is distributed according to the inverse Gaussian (26) because transport 102 is due to advection and dispersion once the particle has been released from the DEP. Thus, the BTC for particle 103 originating at a distance x from the outlet can be written as The full breakthrough curve F (t) is obtained by integration of f (t, x) over the full range L such that In order to extract the tailing behavior, we note that f (t, x) is peaked about t = x/v so that we can write Thus, we can write for the breakthrough curve For times t L/v, we obtain 108 The concentration of colloids after the saturation phase, that lasted 24 hours, resulted non-homogeneously dis-110 tributed across the medium. To quantify this distribution in terms of α, we consider a representative section (15 mm 111 × 4 mm) from the experiment and measure the number of suspended (mobile) colloids using the method described 112 in Methods e. Then, using the corresponding segregation map (see Method a), we compute the fraction of colloids 113 located in DEP areas (see Figure 9): the TPs have a slightly smaller colloids concentration with respect to the DEPs, 114 resulting in 22% of the total suspended colloids accumulated in the DEPs and 78% of the colloids in the TPs (while 115 DEPs represent the 9% of the total volume of the porous system and the TPs represent the 91% of it). This means 116 that the fraction of colloids initially located within a DEP is α = 0.22 and, as a consequence, 1 − α is the fraction 117 of colloids initially located within a TP. We took this α value into consideration to implement the CTRW model for 118 colloids transport and the BTC as discussed in section IV.  H 9 7 8 G z v L Q R M f F D z e q 6 K q X p Q K b s D z v p 2 l 5 Z X V t f X C R n F z a 3 t n 1 9 3 b b 5 g k 0 5 T V a S I S 3 Y q I Y Y I r V g c O g r V S z Y i M B G t G g + u x 3 7 x n 2 v B E 3 c E w Z a E k P c V j T g l Y q e M e + h c e D s 5 w I L M A 2 A N o m c t R x y 1 5 Z W 8 C v E j 8 G S m h G W o d 9 y v o J j S T T A E V x J i 2 7 6 U Q 5 k Q D p 4 K N i k F m W E r o g P R Y 2 1 J F J D N h P r l + h E + s 0 s V x o m 0 p w B P 1 9 0 R O p D F D G d l O S a B v 5 r 2 x + J / X z i C + C n O u 0 g y Y o t N F c S Y w J H g c B e 5 y z S i I o S W E a m 5 v x b R P N K F g A y v a E P z 5 l x d J 4 7 z s e 2 X / t l K q V m Z x F N A R O k a n y E e X q I p u U A 3 V E U W P 6 B m 9 o j f n y X l x 3 p 2 P a e u S M 5 s 5 Q H / g f P 4 A N 7 q U Z g = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " + 0 B W w O M V 1 z M z B b u N t O u i l q 7 h t U A = " > A A A B / X i c b V D J S g N B E O 1 x j X E b l 5 u X x i B 4 k D A j E T 0 G v H i M Y B b I D K G n 0 5 M 0 6 e 4 Z u m v E O A R / x Y s H R b z 6 H 9 7 8 G z v L Q R M f F D z e q 6 K q X p Q K b s D z v p 2 l 5 Z X V t f X C R n F z a 3 t n 1 9 3 b b 5 g k 0 5 T V a S I S 3 Y q I Y Y I r V g c O g r V S z Y i M B G t G g + u x 3 7 x n 2 v B E 3 c E w Z a E k P c V j T g l Y q e M e + h c e D s 5 w I L M A 2 A N o m c t R x y 1 5 Z W 8 C v E j 8 G S m h G W o d 9 y v o J j S T T A E V x J i 2 7 6 U Q 5 k Q D p 4 K N i k F m W E r o g P R Y 2 1 J F J D N h P r l + h E + s 0 s V x o m 0 p w B P 1 9 0 R O p D F D G d l O S a B v 5 r 2 x + J / X z i C + C n O u 0 g y Y o t N F c S Y w J H g c B e 5 y z S i I o S W E a m 5 v x b R P N K F g A y v a E P z 5 l x d J 4 7 z s e 2 X / t l K q V m Z x F N A R O k a n y E e X q I p u U A 3 V E U W P 6 B m 9 o j f n y X l x 3 p 2 P a e u S M 5 s 5 Q H / g f P 4 A N 7 q U Z g = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " + 0 B W w O M V 1 z M z B b u N t O u i l q 7 h t U A = " > A A A B / X i c b V D J S g N B E O 1 x j X E b l 5 u X x i B 4 k D A j E T 0 G v H i M Y B b I D K G n 0 5 M 0 6 e 4 Z u m v E O A R / x Y s H R b z 6 H 9 7 8 G z v L Q R M f F D z e q 6 K q X p Q K b s D z v p 2 l 5 Z X V t f X C R n F z a 3 t n 1 9 3 b b 5 g k 0 5 T V a S I S 3 Y q I Y Y I r V g c O g r V S z Y i M B G t G g + u x 3 7 x n 2 v B E 3 c E w Z a E k P c V j T g l Y q e M e + h c e D s 5 w I L M A 2 A N o m c t R x y 1 5 Z W 8 C v E j 8 G S m h G W o d 9 y v o J j S T T A E V x J i 2 7 6 U Q 5 k Q D p 4 K N i k F m W E r o g P R Y 2 1 J F J D N h P r l + h E + s 0 s V x o m 0 p w B P 1 9 0 R O p D F D G d l O S a B v 5 r 2 x + J / X z i C + C n O u 0 g y Y o t N F c S Y w J H g c B e 5 y z S i I o S W E a m 5 v x b R P N K F g A y v a E P z 5 l x d J 4 7 z s e 2 X / t l K q V m Z x F N A R O k a n y E e X q I p u U A 3 V E U W P 6 B m 9 o j f n y X l x 3 p 2 P a e u S M 5 s 5 Q H / g f P 4 A N 7 q U Z g = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " + 0 B W w O M V 1 z M z B b u N t O u i l q 7 h t U A = " > A A A B / X i c b V D J S g N B E O 1 x j X E b l 5 u X x i B 4 k D A j E T 0 G v H i M Y B b I D K G n 0 5 M 0 6 e 4 Z u m v E O A R / x Y s H R b z 6 H 9 7 8 G z v L Q R M f F D z e q 6 K q X p Q K b s D z v p 2 l 5 Z X V t f X C R n F z a 3 t n 1 9 3 b b 5 g k 0 5 T V a S I S 3 Y q I Y Y I r V g c O g r V S z Y i M B G t G g + u x 3 7 x n 2 v B E 3 c E w Z a E k P c V j T g l Y q e M e + h c e D s 5 w I L M A 2 A N o m c t R x y 1 5 Z W 8 C v E j 8 G S m h G W o d 9 y v o J j S T T A E V x J i 2 7 6 U Q 5 k Q D p 4 K N i k F m W E r o g P R Y 2 1 J F J D N h P r l + h E + s 0 s V x o m 0 p w B P 1 9 0 R O p D F D G d l O S a B v 5 r 2 x + J / X z i C + C n O u 0 g y Y o t N F c S Y w J H g c B e 5 y z S i I o S W E a m 5 v x b R P N K F g A y v a E P z 5 l x d J 4 7 z s e 2 X / t l K q V m Z x F N A R O k a n y E e X q I p u U A 3 V E U W P 6 B m 9 o j f n y X l x 3 p 2 P a e u S M 5 s 5 Q H / g f P 4 A N 7 q U Z g = = < / l a t e x i t > FIG. 9. Relative distribution of suspended particles between DEP and FC in experiment: Location of suspended particles inside the dead-end pores (yellow dots in magenta regions) and transmitting pores (blue dots in green regions) in a subsection of the porous medium. The fraction of total particle counts is α = 0.22 in DEP and 1 − α = 0.78 in TP.