Asymptotic transport and dispersion of active particles in periodic porous media

Asymptotic transport and dispersion of active particles in periodic porous media

Roberto Alonso-Matilla    Brato Chakrabarti    David Saintillan Department of Mechanical and Aerospace Engineering, University of California San Diego, 9500 Gilman Drive, La Jolla, CA 92093, USA
July 19, 2019

We elucidate the long-time transport of active Brownian particles flowing through a porous lattice using generalized Taylor dispersion theory and Langevin simulations. The asymptotic spreading of a dilute cloud of microswimmers is shown to obey an obstacle-free advection-diffusion equation, which we use to unravel the effects of motility, lattice geometry and fluid flow. Our model suggests that the interplay of particle self-propulsion, entropic trapping by obstacles and shear-induced dispersion can be exploited for controlling the transport of active swimmers through microstructured environments.


The transport of active self-propelled particles through complex microstructured media has important implications in microbial ecology as well as human health. Examples include the spreading of contaminants in soils and groundwater aquifers, bacterial filtering, biodegradation and bioremediation processes, and the transport of motile cells inside the body. Engineering areas such as medical diagnostics and biochemical analysis also often hinge on the manipulation and control of active particle motions through complex geometries C. Bechinger et al. (2016).

While the effective transport of passive tracers in porous media flows has been analyzed minutely in the past, the case of active particles remains largely unexplored. In free space, a microswimmer with constant speed and translational and rotational diffusivities and performs a spatial random walk with long-time dispersivity in dimensions Berg (1993). In a porous matrix and without flow, this active dispersion should intuitively be reduced, not only by geometric obstruction but also by the tendency of microswimmers to accumulate at boundaries Ezhilan and Saintillan (2015); Ezhilan et al. (2015); Yan and Brady (2015). Trapping near obstacles is indeed seen in experiments Takagi et al. (2014) and is thought to be enhanced by hydrodynamic interactions Spagnolie et al. (2015). When an external flow is applied, the coupling of particle orientations with the fluid shear further complicates transport. First, diffusion and swimming across streamlines alter long-time spreading by the classic phenomenon of shear-induced Taylor dispersion Taylor (1953); Aris (1956). Secondly, the rotation of particles in the strong shear near boundaries is expected to result in their alignment against the flow E. Altshuler et al. (2013); Creppy et al. (2018), possibly leading to a net reduction in mean transport by an effect similar to upstream swimming in pressure-driven flows Kaya and Koser (2012); Ezhilan and Saintillan (2015). The situation is yet more complex in the presence of asymmetric obstacles, where the gliding of microswimmers along curved boundaries can drive a net migration even in quiescent conditions Galajda et al. (2007); Yariv and Schnitzer (2014); M. S. Davies Wykes et al. (2017); Tong and Shelley (2017).

Some of these trends have been confirmed in experiments A. T. Brown et al. (2016); Sosa-Hernández et al. (2017); Morin et al. (2017); Creppy et al. (2018) as well as various computational models Chepizhko and Peruani (2013); Potiguar et al. (2014); Chamolly et al. (2017); Khalilian and Fazli (2016), yet a general theoretical framework for predicting active transport and dispersion in even simple porous matrices remains lacking. In this Letter, we present a macro-continuum transport model based on generalized Taylor dispersion theory Brenner (1980) that predicts the asymptotic transport velocity and dispersivity of active Brownian particles in periodic porous media. Starting from the Fokker-Planck equation for the configuration of a single particle in a periodic lattice, we derive theoretical expressions for its long-time mean velocity and dispersivity dyadic in an applied flow, where both quantities simply involve the solution of a boundary-value problem on the pore scale. We also propose a macroscale Eulerian interpretation for these coefficients, and use it to unravel the roles of particle motility, fluid flow, and lattice geometry on asymptotic transport properties.

Problem definition.—We analyze the long-time asymptotic transport of a dilute collection of self-propelled active particles dispersed in a viscous solvent and moving through the interstices of a doubly-periodic 2D porous material. The porous medium is modeled as an infinite square lattice comprised of rigid obstacles of characteristic dimension and generated through the discrete translation of a geometric unit cell of linear dimension  (Fig. 1). Each cell in the lattice is labeled by two integers , which identify its horizontal and vertical positions with respect to origin , chosen as the centroid of an arbitrary cell. The array is characterized by its porosity , or ratio of the interstitial fluid area of a cell over its total area.

In a dilute system, we can neglect interparticle interactions and need only consider the transport of a single swimmer. Its instantaneous configuration is described by its position with respect to point and by its swimming direction in the plane of motion, where . The global coordinate can be decomposed as , where is the position of the unit cell where the particle is located and is a local coordinate with respect to the center of that cell. The particle is released at at position with orientation . We introduce the following gradient operators: global , local and orientational .

Particle motion results from self-propulsion, from advection and rotation by the fluid flow, and from translational and rotational diffusion. We model its dynamics by two coupled Langevin equations


where and are Gaussian random vectors with zero mean and variance and , respectively. The velocity field is a 2D Stokes flow driven by a macroscopic pressure gradient applied across the array, with uniform speed and incoming angle upstream of the array; we solve for it numerically by the boundary integral method (see Supplemental Material 111Supplemental material can be found at XXXX.). The fluid vorticity also causes rotation of the particle in Eq. (2). Our aim is to describe the long-time statistics of , which we proceed to model within a continuum framework.

Figure 1: () Geometry of a representative lattice, unit cell and associated nomenclature. () Non-spherical obstacle shapes are generated by conformal mapping Note1 (); Avron et al. (2004), where parameters and control obstacle aspect ratio and fore-aft asymmetry, respectively.

Fokker-Planck description.—We introduce the conditional probability density of finding the particle at a position with orientation at time , given that it was initially released at . It satisfies the Fokker-Planck equation Chandrasekhar (1943); Doi and Edwards (1986)


where the initial condition is incorporated by the source term on the right-hand side. The spatial and orientational probability fluxes are easily inferred from the Langevin equations (1)–(2) as Note1 ()


For a given cell (), Eq. (3) can be recast solely in terms of local variables as


where . The probability density function is normalized as


which implies convergence of the previous sum and the subsequent convergence of the probability moments. A no-flux condition is enforced on the obstacle boundaries, as well as continuity at the edges of the unit cells Note1 (). We henceforth use dimensionless equations where we scale variables using length scale , time scale , and velocity scale , and we also normalize by . In addition to , non-dimensionalization of the governing equations yields three dimensionless groups: , , and Ezhilan and Saintillan (2015).

Figure 2: (a) Transport and spreading of a cloud of active particles around circular pillars in Brownian dynamics (BD) simulations at two different times (left two columns) compared to the Darcy-scale theoretical prediction obtained by solving Eq. (15) (third column). The two rows correspond to incoming flow angles and . See Note1 () for movies showing the temporal evolution. (b) Probability density function of the macroscale horizontal position for the two instantaneous times shown in (a) for , where analytical predictions are compared to BD simulation results. (c) Time evolution of the variances and of particle positions, showing convergence to the asymptotic diffusive regime. (d) Asymptotic solutions on the pore scale for two pillar shapes: streamlines and magnitude of the imposed flow (top), density field (middle), and polarization field (bottom). In all panels, , , , and .

Macrotransport model.—Following generalized Taylor dispersion theory Brenner (1980), we seek the mean asymptotic transport properties in terms of the th-order polyadic local moments of defined as


The mean asymptotic transport velocity can be obtained from the steady-state zeroth-order local moment , which describes the long-time probability of finding the particle at local position with orientation irrespective of which unit cell it is traversing. It satisfies the simplified steady Fokker-Planck equation , with fluxes


subject to normalization condition , to the no-flux condition on the surface of the pillar, and to jump conditions and on the unit cell edges and . We solve for numerically using a finite-volume method.  is then obtained by simple quadrature as Note1 ()


where is the outward unit normal along the cell edges.

Determining the long-time dispersivity tensor requires consideration of the first-order moment , whose asymptotic form can be derived as


The fluctuation field encapsulates information about dispersion and obeys the steady problem


subject to the no-flux condition at the obstacle walls, and to jump conditions and on and . We numerically solve for , and its knowledge provides as a second quadrature


where denotes the spatial and orientational average of with weighting factor Note1 ().

On long length and time scales, a dilute cloud of active particles is thus expected to be transported with mean velocity and to spread with dispersivity as provided by Eqs. (11) and (14). This suggests seeking a coarse-grained Eulerian interpretation of and at the so-called Darcy scale, where the fluid and solid obstacles that comprise the porous material become indistinguishable and distances shorter than the characteristic size of the unit cell are irrelevant. To this end, we introduce the discrete conditional probability density Brenner (1980) for a particle to be located in cell () at time . On large length scales, we can assimilate to a continuous variable and define a corresponding macro-scale gradient operator . The Darcy-scale probability density , which is the continuous counterpart of , is then expected to formally satisfy an obstacle-free advection-diffusion equation Note1 ()


As shown below, the analytical solution of Eq. (15) is in exquisite agreement with simulation results at long times, thereby solidifying the interpretation of and .

Results and discussion.—Results from Brownian Dynamics (BD) simulations of the spreading of a concentrated point source of active swimmers based on Eqs. (1)–(2) are compared to the analytical solution of Eq. (15) in Fig. 2. On short time scales, the cloud develops a complex shape that is continually distorted by the pillars and exhibits wakes, boundary layers, and other features typical of advective-diffusive transport around obstacles. At later times, the cloud shape becomes increasingly smooth and is very well captured by the theoretical model, which predicts spreading into a translating anisotropic Gaussian. This comparison asserts the strength of the Darcy-scale interpretation for describing asymptotic transport, but also highlights key differences. In particular, we observe tailed distributions at pre-asymptotic times in Fig. 2(b), a phenomenon also recently reported in experiments using E. coli Creppy et al. (2018). The distributions are found to be negatively skewed in the flow direction in strong flows of weakly swimming particles (high and low ), whereas they become positively skewed in the opposite case Note1 (). At longer times, the tails lose their asymmetry and the distributions converge to the asymptotic diffusive solution as shown in Fig. 2(,).

Asymptotic transport characteristics can be attributed to pore-scale features of the fluid flow and of the steady-state zeroth-order moment as illustrated in Fig. 2(). In the absence of flow, the interplay of self-propulsion and translational diffusion near obstacle boundaries creates a net swimmer polarization against the obstacles Ezhilan and Saintillan (2015); Ezhilan et al. (2015); Yan and Brady (2015), which vanishes at the cell edges by symmetry. This polarization is accompanied by a net increase in density , which is azimuthally symmetric near a circular obstacle and enhanced in regions of high curvature for non-circular pillars Yan and Brady (2018). When a flow is applied, the fore-aft symmetry is broken and particles now preferentially accumulate near stagnation points, with a maximum density reached in the rear of circular obstacles. This accumulation, which is observed in experiments G. Miño et al. (2018) and is akin to that occurring past funnel constrictions E. Altshuler et al. (2013), results from the combination of advection by the flow and of shear-rotation and swimming in the fast-flowing high-shear regions above and below the pillars. Advection and rotation also induce a density minimum ahead of the obstacle, which coincides with a stagnation point in the polarization field. As explained below, these features play a central role in determining the effective spreading behavior.

Figure 3: Dispersion in an imposed flow: ()–() Dependence of and on flow Péclet number for various values of in a lattice with circular pillars for an external flow in the -direction (). Inset in () highlights strong-flow scalings. () Maximum eigenvalue of the dispersivity dyadic and () corresponding direction of maximum dispersion as a function of incoming flow angle. Parameter values: , , . Symbols: BD simulations; lines: analytical model.

The dependence of dispersion coefficients on external flow strength and swimming activity is shown in Fig. 3(,) for a lattice comprised of circular pillars. In weak flows, and have similar magnitudes, with faster-swimming particles (high ) spreading more efficiently by active dispersion. Quite remarkably, increasing flow strength first causes a decrease in spreading. This decrease is unique to self-propelled particles, and indeed disappears at low values of . It is attributed to a reduction in active dispersion by the flow, which prevents particles from freely sampling all swimming directions. At higher values of , starts increasing again and displays an asymptotic scaling of , with an exponent of at intermediate flow strengths and of in very strong flows, which is the classical scaling for the porous dispersion of passive tracers Edwards et al. (1991). Interestingly, activity is seen to hinder in strong flows in Fig. 3(), as swimming enhances cross-stream migration, which reduces axial spreading within the framework of shear-induced dispersion. Unlike , the transverse dispersivity in Fig. 3() monotonically decays with as stronger flows cause faster rotations of the swimmers in the local shear and thus limit their ability to travel in the -direction.

Dispersion can be further controlled by tuning the incoming flow angle in Fig. 3(,). For arbitrary , the dispersion tensor is no longer diagonal and maximum dispersion occurs along the eigendirection associated with its maximum eigenvalue . The magnitude of is maximum for , which allows for easy gliding of the swimmers between rows of pillars with minimal collisions. Other flow directions result in more severe obstruction by the pillars and thus display much lower values of , which reaches its minimum for at high flow rates. While in strong flows the direction of maximum dispersion is close to that of the imposed flow, such is not the case in weak flows where and have opposite signs, indicating significant spreading transverse to the flow.

The effect of lattice porosity is analyzed in Fig. 4(,) and further underscores the subtle interplay of active vs shear-induced dispersion. In the absence of flow and at low porosity, dispersion is weak as the pillars act as entropic barriers that hinder active transport Laachi et al. (2007). Increasing porosity thus causes a monotonic increase in and as geometric obstruction plays a lesser role. In strong flows (), porosity instead causes a monotonic decrease in , as the obstacles act as regions of shear production, thus enhancing shear-induced dispersion. In weaker flows (), shows a non-monotonic behavior: as porosity is increased, axial dispersion first decreases due to a drop in local shear rate, before increasing again at low obstacle area fractions due to active swimming. Curiously, the limit of differs from the case of swimmers in free-space, as infinitesimal obstacles still induce fluid shear in an imposed flow while also modifying the swimmer distribution in their vicinity.

Figure 4: Effect of lattice porosity and obstacle geometry: ()–() Dependence of and on in the case of circular obstacles at and for various flow strengths . ()–() Effect of shape aspect ratio parameter and fore-aft asymmetry on at for different values of . Inset shows obstacle shapes. All plots are for and .

Pillar aspect ratio and fore-aft asymmetry both also affect transport in non-trivial ways. As shown in Fig. 4(), elliptical obstacles aligned with a weak flow cause higher streamwise dispersion than if aligned perpendicular to it, as the former allow swimmers to glide past while the latter act as entropic barriers that obstruct transport Laachi et al. (2007). These trends reverse at high flow rates, as perpendicular obstacles induce stronger velocity gradients thereby enhancing shear-induced dispersion. Fore-aft asymmetric obstacles also lead to complex trends in Fig. 4(). In the absence of flow, decreases as the shape parameter deviates from (circle) by the aforementioned entropic obstruction mechanism. When a flow is applied, fluid shear causes particles to swim along the walls and against the flow, resulting in an enhancement of particle accumulation around obstacles for and in a reduction for . Consequently, the effects of shear-alignment (in weak flows) and shear-induced dispersion (in strong flows) on transport are magnified for obstacle shapes that promote accumulation (), as seen in Fig. 4(). This asymmetry becomes negligible at very high values of , as activity becomes subdominant and the magnitude of the shear gradients that drive dispersion is independent of the sign of by reversibility of Stokes flow.

Concluding remarks.—We have analyzed the long-time transport properties of active particles in a porous lattice under both quiescent and flow conditions. The predictions of our continuum model, which agree perfectly with Langevin simulations at long times, highlight the complex interplay of particle motility, alignment under shear, cross-stream migration, lattice porosity and pillar geometry on asymptotic dispersion. Our theory provides a simple framework for analyzing microorganismal transport in natural structured environments and for designing engineered porous media in applications involving microswimmers. The fundamental premise of non-interacting active Brownian particles glosses over many details that may become relevant in specific systems. Future work should address the roles of hydrodynamic interactions with obstacles Spagnolie et al. (2015), swimmer-specific scattering dynamics Kantsler et al. (2013), rheological effects due to activity Saintillan (2018), and the potential emergence of spontaneous flows and active turbulence in denser systems Theillard et al. (2017).

D.S. gratefully acknowledges funding from NSF Grant DMS-1463965.


  • C. Bechinger et al. (2016) C. Bechinger et al., Rev. Mod. Phys. 88, 045006 (2016).
  • Berg (1993) H. C. Berg, Random Walks in Biology (Princeton University Press, 1993).
  • Ezhilan and Saintillan (2015) B. Ezhilan and D. Saintillan, J. Fluid Mech. 777, 482 (2015).
  • Ezhilan et al. (2015) B. Ezhilan, R. Alonso-Matilla,  and D. Saintillan, J. Fluid Mech. 781, R4 (2015).
  • Yan and Brady (2015) W. Yan and J. F. Brady, J. Fluid Mech. 785, R1 (2015).
  • Takagi et al. (2014) D. Takagi, J. Palacci, A. B. Braunschweig, M. J. Shelley,  and J. Zhang, Soft Matter 10, 1784 (2014).
  • Spagnolie et al. (2015) S. E. Spagnolie, G. R. Moreno-Flores, D. Bartolo,  and E. Lauga, Soft Matter 11, 3396 (2015).
  • Taylor (1953) G. I. Taylor, Proc. R. Soc. Lond. A 219, 186 (1953).
  • Aris (1956) R. Aris, Proc. R. Soc. Lond. A 235, 67 (1956).
  • E. Altshuler et al. (2013) E. Altshuler et al., Soft Matter 9, 1864 (2013).
  • Creppy et al. (2018) A. Creppy, E. Clément, C. Douarche, M. V. D’Angelo,  and H. Auradou, arXiv:1802.01879  (2018).
  • Kaya and Koser (2012) T. Kaya and H. Koser, Biophys. J. 102, 1514 (2012).
  • Galajda et al. (2007) P. Galajda, J. Keymer, P. Chaikin,  and R. Austin, J. Bacteriol. 189, 8704 (2007).
  • Yariv and Schnitzer (2014) E. Yariv and O. Schnitzer, Phys. Rev. E 90, 032115 (2014).
  • M. S. Davies Wykes et al. (2017) M. S. Davies Wykes et al., Soft Matter 13, 4681 (2017).
  • Tong and Shelley (2017) J. Tong and M. J. Shelley, to appear in SIAM J. Appl. Math.  (2017).
  • A. T. Brown et al. (2016) A. T. Brown et al., Soft Matter 12, 131 (2016).
  • Sosa-Hernández et al. (2017) J. E. Sosa-Hernández, M. Santillán,  and J. Santana-Solano, Phys. Rev. E 95, 032404 (2017).
  • Morin et al. (2017) A. Morin, D. L. Cardozo, V. Chikkadi,  and D. Bartolo, Phys. Rev. E 96, 042611 (2017).
  • Chepizhko and Peruani (2013) O. Chepizhko and F. Peruani, Phys. Rev. Lett. 111, 160604 (2013).
  • Potiguar et al. (2014) F. Q. Potiguar, G. Farias,  and W. Ferreira, Phys. Rev. E 90, 012307 (2014).
  • Chamolly et al. (2017) A. Chamolly, T. Ishikawa,  and E. Lauga, New J. Phys. 19, 115001 (2017).
  • Khalilian and Fazli (2016) H. Khalilian and H. Fazli, J. Chem. Phys. 145, 164909 (2016).
  • Brenner (1980) H. Brenner, Phil. Trans. R. Soc. Lond. A 297, 81 (1980).
  • (25) Supplemental material can be found at XXXX.
  • Avron et al. (2004) J. Avron, O. Gat,  and O. Kenneth, Phys. Rev. Lett. 93, 186001 (2004).
  • Chandrasekhar (1943) S. Chandrasekhar, Rev. Mod. Phys. 15, 1 (1943).
  • Doi and Edwards (1986) M. Doi and S. F. Edwards, The Theory of Polymer Dynamics (Oxford University Press, 1986).
  • Yan and Brady (2018) W. Yan and J. F. Brady, Soft Matter 14, 279 (2018).
  • G. Miño et al. (2018) G. Miño et al., Adv. Microbiol. 8, 451 (2018).
  • Edwards et al. (1991) D. Edwards, M. Shapiro, H. Brenner,  and M. Shapira, Trans. Porous Media 6, 337 (1991).
  • Laachi et al. (2007) N. Laachi, M. Kenward, E. Yariv,  and K. D. Dorfman, Europhys. Lett. 80, 50009 (2007).
  • Kantsler et al. (2013) V. Kantsler, J. Dunkel, M. Polin,  and R. E. Goldstein, Proc. Natl. Acad. Sci. USA 110, 1187 (2013).
  • Saintillan (2018) D. Saintillan, Annu. Rev. Fluid Mech. 50, 563 (2018).
  • Theillard et al. (2017) M. Theillard, R. Alonso-Matilla,  and D. Saintillan, Soft Matter 13, 363 (2017).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description