# The influence of finite size particles on fluid velocity and transport
though porous media

M. Residori  S. Praetorius  Institute of Scientific Computing, TU Dresden,
01062 Dresden, Germany  P. de Anna  Institute of Earth Sciences, University of
Lausanne, Lausanne 1015, Switzerland<EMAIL_ADDRESS>A. Voigt
Institute of Scientific Computing, TU Dresden, 01062 Dresden, Germany and
Dresden Center for Computational Materials Science (DCMS), TU Dresden, 01062
Dresden, Germany<EMAIL_ADDRESS>
###### Abstract

Understanding the coupling between flow, hydrodynamic transport and dispersion
of colloids with finite-size in porous media is a long-standing challenge.
This problem is relevant for a broad range of natural and engineered
subsurface processes, including contaminant and colloidal transport, mixing of
bio-chemical compounds, kinetics of reactions and groundwater bio-remediation,
but also transport phenomena related to different systems like membranes, or
blood flow. While classical models for colloidal transport rely on macro-
dispersion theory and do not take into consideration the complex and
heterogeneous structure of the porous host medium, recent studies take into
consideration the detailed structure of the porous system and its impact on
the fluid velocity. However, the impact of confinement condition, represented
by the ratio of particles radius $a$ and pore throat size $\lambda$, has been
overlooked. Here, we use numerical simulations of fluid particle dynamics in
resolved porous media to demonstrate that particles confinement affects the
fluid macroscopic velocity field $\mathbf{u}$ which in turn affects the
particles transport itself. Our results show that even for small confinement
conditions ($a/\lambda\sim 2$ %), fluid and transported particles are
dynamically re-routed towards more permeable paths. This leads to the
emergence of ephemeral laminar vortexes at pore throat entrances and affects
the variance and mean fluid velocity.

##  I Introduction

Most engineered and natural systems characterized by a porous structure can
host fluids that, driven by a macroscopic pressure drop, move through their
network of pores beardynamics1988 . In most scenarios, the concerned fluid is
water that, being an excellent solvent, is carrying dissolved substances or
suspended particles, whose transport rate controls several phenomena,
including mixing, reaction kinetics and filtration Molnar2005 . The
characteristics of these phenomena are a direct consequence of the underlying
pore structure. A common issue of such pore structures is their inherent
structural heterogeneity (Dentz2011, ), or spatial variability, that has been
shown to emerge from sub-micron to meter scales (freezecherry, ). A wide range
of variability appears in the distribution of the individual pores size (i.e.
the space among solid grains, available for fluid motion).

Since within such porous systems fluid velocities are typically low (on the
order of $1\,m/d$, equivalent to about $20\,\mu m/s$) the flow of water can
assumed to be laminar, controlled by viscous forces rather than inertia, and
stationary (beardynamics1988, ). Under such conditions, an analytical
relationship between the pore throat size distribution
$f_{\lambda}\sim\lambda^{-\beta}$ and the distribution of the low fluid
velocities $f_{u}\sim u^{-\beta/2}$ has been suggested in (deAnnaPRF2017, ).
Under the assumption of Stokes flow it has been shown that each pore-throat
hosting a net transfer of fluid has the same statistical distirbution of
velocities as the one of a pipe, called the porelet, so that the overall
probability density function (PDF) of fluid velocities within the
heterogeneous structure can be determined by the superposition of all porelets
in two-dimensional (deAnnaPRF2017, ) and three-dimensional (DentzJFM2018, )
systems. This relationship, which allows to make predictions within a
continuous-time random-walk framework for the asymptotic statistics of the
spreading of fluid particles along their own trajectories, is based on the
linearity of the Stokes equation and assumes the dissolved substances and
suspended particles to have no impact on the stationary flow field.

The spatial variability of fluid velocity, which measures the overall velocity
contrast between pores or between the channels of high velocity and zones of
stagnation, has a major control on the transport of substances, dissolved or
suspended (Dentz2011, ). While classical macro dispersion theories have been
shown to adequately represent relatively homogeneous media (e.g. the ones
characterized by a well defined scale for average pore size) beardynamics1988
, it is known that such model predictions diverge from the observation of
transport through complex porous systems that results in anomalous transport
properties as early arrival times and long tailing Berkowitz2006 ; Dentz2011 .
In such media the observed transport also impacts mixing kinetics
deAnnaEST2014 ; Heyman2020 , and mixing-driven processes as filtration
Nishiyama2012 ; miele2019 , or microbial dispersal (Scheidweiler2020, ;
deAnnaNaturePhys2021, ). In several scenarios, flowing water carries suspended
colloidal particles and microbial cells or aggregates. The size of such
suspended particles and cell aggregates can vary over several orders of
magnitude, affecting their sedimentation (morrisbook, ; bergbook, ) and
overall transport (Sirivithayapakorn2003, ; Auset2004, ). In particular,
enhanced suspended particle velocity in confined environments has been
observed and associated to their finite-size, also known as size exclusion
effect (Prieve1978, ; Babakhani2019, ). This phenomenon forces larger
particles to remain within channels of high flow while preventing them from
access pores of low velocity whose size is comparable to the particle
diameter. The macroscopic effect is to rise the average suspended particle
velocity modifying the macroscopic breakthrough curves (Keller2004, ) and
dispersion coefficient (James2003, ). These phenomena clearly breakdown the
assumption that transported particles do not affect the stationary velocity
field, as in (deAnnaPRF2017, ; miele2019, ; Scheidweiler2020, ) that relate
the medium physical structure to the fluid and suspended particles velocity to
predict macroscopic transport.

Here, we investigate the impact of transported finite-size particles on the
fluid velocity distribution and dynamics, by means of a novel fluid particle
numerical scheme for simulating flow and transport through complex porous
structures that takes particle-particle, particles-fluid interactions and
inertia effects into account. We consider the local feedback of particles with
finite radius $a$, smaller but comparable to the local pore throat size
$\lambda$, on the fluid velocity distribution. Strictly, finite-size particles
when facing a constriction (pore opening) of comparable size experience the
viscous drag exerted by the no-slip boundary conditions at grain walls. This
effect is enhanced by particle agglomerates of larger size. Thus, a dynamical
interaction between the passage of finite-size particles through such pores
and the flow field is expected. We numerically investigate the dynamical
change of the fluid and particles velocity in laminar conditions. We compare
our simulation results with Stokes flow and point-particles. We show that for
even small confinement conditions ($a/\lambda\sim 2$ %), fluid and transported
particles are dynamically re-routed towards more permeable paths. This leads
to the emergence of ephemeral laminar vortexes at pore throat entrances and
affects the variance and mean fluid velocity.

##  II Methods

To study the impact of suspended and finite-size particles on flow and
transport we consider a two-dimensional porous domain $\Omega$ modeled by a
rectangle $[x_{0},x_{1}]\times[y_{0},y_{1}]$ with cutout non-overlapping disks
of random position and diameter, see Figure 1(a). This disordered arrangement
of disks is characterized by a Delaunay triangulation of the disk centers to
identify the nearest neighbors: each triangle defines a pore and each edge
defines a throat, see Figure 1(b). The statistical properties of the structure
are characterized in terms of pore throat size $\lambda=d-r_{1}-r_{2}$
distribution: $d$ being the distance between the two neighboring disk centers
and $r_{1}$ and $r_{2}$ are the respective radii. To consider a heterogeneous
medium the domain structure is generated such that the PDF of the pore throat
size is roughly following the power law distribution
$p_{\lambda}(\lambda)\sim\lambda^{-\beta}$, with $\beta=0.17$ as in
deAnnaPRF2017 .

Figure 1: _(a)_ Geometry $\Omega$ and frames of the sub-geometries
$\Omega_{i}$ considered for Stokes flow and fluid particle dynamics,
respectively. The frames are color coded with respect to the velocity
magnitude averaged over the sub-geometry, see Figure 3. _(b)_ Delaunay
triangulation of $\Omega$, where edges are marked in red. Each edge
corresponds to $d$. _(c)_ Rescaled velocity magnitude
$u_{r}=\log_{10}(u/\langle u\rangle);$ $u=\lVert\mathbf{u}\rVert$, from eqs.
(1)-(2). _(d)_ Streamlines of point-particles initialized at the left side of
$\Omega$. The different colors are for visualization purposes. They correspond
to three equal size subdomains of the left boundary.

###  II.1 Stokes dynamics

Within this geometry we consider the Stokes equations

$\displaystyle-\nu\nabla^{2}\mathbf{u}+\nabla p$ $\displaystyle=\mathbf{0},$
(1) $\displaystyle\nabla\cdot\mathbf{u}$ $\displaystyle=0,$ (2)

where $\mathbf{u}=\mathbf{u}(\mathbf{x})$ is the fluid velocity,
$p=p(\mathbf{x})$ the fluid pressure, $\nu$ the fluid viscosity and
$\mathbf{x}\in\Omega\subset\mathbb{R}^{2}$. Boundary conditions at
$\partial\Omega$ are specified at the rectangle edges and the grain walls
(disks perimeter). We impose $\mathbf{u}=\mathbf{u}_{\text{inflow}}$ at the
left side of the rectangle $\partial\Omega_{L}$,
$\nabla\mathbf{u}\cdot\mathbf{n}=\mathbf{0}$ at the right side of the
rectangle $\partial\Omega_{R}$, $\mathbf{u}\cdot\mathbf{n}=\mathbf{0}$ at top
$\partial\Omega_{T}$ and bottom $\partial\Omega_{B}$ rectangle boundaries and
$\mathbf{u}=\mathbf{0}$ at the remaining boundaries
$\partial\Omega_{\text{disk}}$. Thereby, $\mathbf{n}$ denotes the outward
normal to $\partial\Omega$. We consider

$\mathbf{u}_{\text{inflow}}(x_{0},y)=\left[\frac{1}{2}\left(1+\tanh\left(\alpha-\left|y-\frac{y_{1}+y_{0}}{2}\right|\right)\right),0\right]^{T},$

where $\alpha$ is chosen such that
$\lVert\mathbf{u}_{\text{inflow}}(x_{0},y_{0})\rVert=\lVert\mathbf{u}_{\text{inflow}}(x_{0},y_{1})\rVert\leq
10^{-16}$ in order to numerically satisfy the compatibility condition with the
top and bottom boundary conditions.

The problem is solved by finite elements. To implement no-slip boundary
conditions at solid interfaces and inflow conditions, we consider weak
solutions of eqs. (1) and (2) with $\mathbf{u}\in\mathbf{V}:=\\{\mathbf{v}\in
H^{1}(\Omega)^{2}:\mathbf{v}=\mathbf{0}\text{ at
}\partial\Omega_{\text{disk}},\,\mathbf{v}=\mathbf{u}_{\text{inflow}}\text{ at
}\partial\Omega_{L},\,\nabla\mathbf{u}\cdot\mathbf{n}=\mathbf{0}\text{ at
}\partial\Omega_{R}\\}$ and $p\in L^{2}(\Omega)$, such that

$\displaystyle\int_{\Omega}\nabla\mathbf{v}:\nu(\nabla\mathbf{u}+(\nabla\mathbf{u})^{T})-p\nabla\cdot\mathbf{v}\,\mathrm{d}\mathbf{x}$
$\displaystyle=\mathbf{0},$ (3)
$\displaystyle\int_{\Omega}q\nabla\cdot\mathbf{u}\,\mathrm{d}\mathbf{x}$
$\displaystyle=0,$ (4)

for every $\mathbf{v}\in H^{1}_{\mathbf{0}}(\Omega)^{2}$ and $q\in
L^{2}(\Omega)$. We partition the domain $\Omega$ by a conforming triangulation
$\mathcal{T}_{h}$. Then, the continuous spaces $\mathbf{V}\times L^{2}$ are
approximated by the Taylor-Hood space $\mathbb{T}_{h}$ defined as

$\mathbb{T}_{h}=(\mathbb{V}_{2,g_{0}}\times\mathbb{V}_{2,g_{1}})\times\mathbb{V}_{1},\quad\\!\\!\\!\mathbb{V}_{m,g}=\\{v\in\mathbb{V}_{m}:\operatorname{tr}_{\partial\Omega^{\prime}}v=g\\},\quad\\!\\!\\!\mathbb{V}_{m}=\\{p\in
C^{0}(\Omega):p|_{T}\in\mathbb{P}_{m}(T),\forall T\in\mathcal{T}_{h}\\},$

where $\mathbb{P}_{m}$ is the space of polynomials of order at most $m$ and
$\partial\Omega^{\prime}=\partial\Omega\setminus\partial\Omega_{R}$. The
velocity $\mathbf{u}$ and pressure $p$ are approximated by functions from
$\mathbb{T}_{h}$, i.e. $m=2$ piece-wise quadratic for $\mathbf{u}$ and $m=1$
piece-wise linear for $p$. In the simulations we consider $\nu=10^{2}$, see
Section II.2.4.

The resulting velocity field $\mathbf{u}$ is heterogeneous and organized in
channels of high and fluctuating velocity and zones of stagnation. Figure 1(c)
shows the logarithm of $u_{r}=u/\langle u\rangle$ with
$u=\lVert\mathbf{u}\rVert$ and $\langle\cdot\rangle$ the average. The color-
map is such that blue (red) regions are associated to velocity magnitude below
(above) its average value (gray). In order to visualize the transport
properties within the medium we track the displacement of point-particles
along the streamlines, as shown in Figure 1(d) where for each trajectory the
color denotes its own initial vertical location, separated in three regions.

Figure 2: _Left_. Double logarithmic plot of the probability density function
(PDF) of the pore throat sizes distribution $f_{\lambda}$. The black line is a
fit of the data. The theoretical distribution is $\lambda^{-\beta}$,
$\beta=0.166$, which match the fitted exponent. _Right_. Double logarithmic
plot of the probability density function of the velocities distribution
$f_{u_{r}}$. The black line is a fit of the data for small velocities, i.e.
$u_{r}<10^{-2}$. The fitted line has slope $-0.100\pm 0.022$ in agreement with
the prediction $f_{u_{r}}\sim u_{r}^{-\beta/2}$.

Figure 2 shows the pore throat size distribution (left) and the distribution
of the low velocities (right). The results indicate a related power law
distribution in the low range of pore throat size
$f_{\lambda}\sim\lambda^{-\beta}$ and low velocities $f_{u_{r}}\sim
u_{r}^{-\beta/2}$, with $\beta=0.166$. The exponents are within the 95%
confidence interval of the analytic theory of (deAnnaPRF2017, ).

###  II.2 Fluid particle dynamics

####  II.2.1 Mathematical model

To develop a model that takes into account the fluid-particle dynamics we
follow the method propose by tanaka00 . We consider the Navier-Stokes
equations

$\displaystyle\rho(\mathbf{u}_{t}+(\mathbf{u}\cdot\nabla)\mathbf{u})-\nabla\cdot(\nu(\nabla\mathbf{u}+\nabla\mathbf{u}^{T}))+\nabla
p$ $\displaystyle=\mathbf{f},$ (5) $\displaystyle\nabla\cdot\mathbf{u}$
$\displaystyle=0,$ (6)

where $\mathbf{u}=\mathbf{u}(t,\mathbf{x})$ is the fluid velocity,
$p=p(t,\mathbf{x})$ the fluid pressure, $\rho=const$ the fluid density,
$\nu=\nu(t,\mathbf{x})$ the viscosity, $\mathbf{f}=\mathbf{f}(t,\mathbf{x})$
the external volume forces and $(t,\mathbf{x})\in(0,T]\times\Omega$. The
external boundary conditions are the same as in the previous Section II.1 with
$\Omega$ replaced by sampled subdomains $\Omega_{i}$ to reduce the
computational cost, see Figure 1(a) and Figure 3.

Figure 3: Different geometries $\Omega_{i}$, $i=1,\dots,12$ obtained from
$\Omega$ in Figure 1(a) ordered with respect to the average value of the
Stokes flow $\langle u\rangle$ over the domain. The geometries consider
different throat sizes but are too small to define a reasonable throat size
distribution.

The fluid-particle dynamics approach describes particles as high viscosity
regions within the fluid. This implies that the momentum transfer within the
particle is much slower than through the fluid outside and it thus behaves
like a solid particle. In order to deal with the viscosity contrast the
particle is described by a diffuse domain approach. Let
$\mathbf{x}_{j}=\mathbf{x}_{j}(t)$ be the position of the center of mass of
particle $j$ of radius $a$. Then, the particle is represented through the
concentration field

$\phi_{j}(\mathbf{x})=\frac{1}{2}\left(1+\tanh\left(\frac{a-|\mathbf{x}-\mathbf{x}_{j}|}{\xi}\right)\right),\quad\mathbf{x}\in\Omega_{i},$
(7)

where $\xi$ is the width of the diffuse interface. For $\xi\to 0$ the
concentration field $\phi_{j}(\mathbf{x})$ approaches a characteristic
function to represent the particle. Let $\nu_{F}$ and $\nu_{P}$ with
$\nu_{F}\ll\nu_{P}$ be the fluid and particle viscosity, respectively. Then,
the viscosity field $\nu$ is defined as

$\nu(\mathbf{x})=\nu_{F}+\sum_{j=1}^{N}(\nu_{P}-\nu_{F})\phi_{j}(\mathbf{x}),\quad\mathbf{x}\in\Omega_{i},$
(8)

where $N$ is the number of particles. In the limit case
$\nu_{P}/\nu_{F}\to\infty$ and $\xi\to 0$ we recover solid particle dynamics,
see tanaka00 . In this way, interactions between fluid and particles are
encoded in the viscosity field. Particle-particle interactions are encoded in
the force term $\mathbf{f}$. We consider the repulsive part of the Lennard-
Jones potential, i.e.

$V(\ell)=\varepsilon\left(\frac{\sigma}{\ell}\right)^{12},$ (9)

where $\varepsilon$ is the strength of the potential, $\sigma$ the interaction
range and $\ell$ the distance. The force acting on particle $j$ is given by

$\mathbf{f}_{j}=-\frac{\partial\sum_{k\neq
j}V(|\mathbf{x}_{j}-\mathbf{x}_{k}|)}{\partial\mathbf{x}_{j}},\quad j=1,\dots
N.$ (10)

Finally, the continuous force field, which enters the Navier-Stokes equations,
is given by

$\mathbf{f}(\mathbf{x})=\sum_{j=1}^{N}\mathbf{f}_{j}\phi_{j}(\mathbf{x}).$
(11)

The position of particle $j$ is determined by

$\frac{\mathrm{d}\mathbf{x}_{j}(t)}{\mathrm{d}t}=\mathbf{u}_{j}(t)\quad\mbox{with}\quad\mathbf{u}_{j}(t)=\frac{\int_{\Omega}\mathbf{u}(t,\mathbf{x})\phi_{j}(\mathbf{x})\,\mathrm{d}\mathbf{x}}{\int_{\Omega}\phi_{j}(\mathbf{x})\,\mathrm{d}\mathbf{x}}.$
(12)

Eqs. (5) - (12) together with the described boundary conditions and
appropriate initial conditions define the model to solve. The model contains
several length scales (diffuse interface width $\xi$, particle radius $a$,
average throat size $\langle\lambda\rangle$ and domain size $|\Omega_{i}|$),
for which $\xi<a<\langle\lambda\rangle<|\Omega_{i}|$, and several time scales
which range from fast particle-particle interactions to slow penetration
times. All these different scales need to be resolved by the numerical method,
which makes it computational expensive.

In the limit of point-particles (particle radius $a\to 0$) also the
interaction range $\sigma$, which is related to $a$, goes to zero and the
fluid-particle and particle-particle interactions vanish. In this limit eqs.
(5) - (12) can be approximated by the Stokes equations (1) and (2).

####  II.2.2 Time discretization

As the finite-size of the particles dynamically modify the local flow field,
the governing equation (5) explicitly depends on time $t$. To solve the flow
evolution through time, we discretize time into $K$ steps of duration
$\tau=T/K$ so that $t_{k}=k\tau\,\in[0,T]$, $k=0,\dots,K$. We consider a semi-
implicit Euler method and an operator splitting approach to systematically
solve first for the fluid flow and then for the particle positions integrating
their motion equations (12). Given the value of
$\mathbf{u}^{k}=\mathbf{u}(t_{k},\mathbf{x})$, $\nu^{k}=\nu(t_{k},\mathbf{x})$
and $\mathbf{f}^{k}=\mathbf{f}(t_{k},\mathbf{x})$, the velocity field and
pressure at the next time step, $\mathbf{u}^{k+1}$ and $p^{k+1}$,
respectively, are approximated by

$\displaystyle\rho\left(\frac{1}{\tau}\mathbf{u}^{k+1}+(\mathbf{u}^{k}\cdot\nabla)\mathbf{u}^{k+1}+(\mathbf{u}^{k+1}\cdot\nabla)\mathbf{u}^{k}\right)-\nabla\cdot(\nu^{k}(\nabla\mathbf{u}^{k+1}+$
$\displaystyle(\nabla\mathbf{u}^{k+1})^{T}))+\nabla p^{k+1}$
$\displaystyle=\mathbf{f}^{k}+\rho\left(\frac{1}{\tau}\mathbf{u}^{k}+(\mathbf{u}^{k}\cdot\nabla)\mathbf{u}^{k}\right),$
(13) $\displaystyle\nabla\cdot\mathbf{u}^{k+1}$ $\displaystyle=0.$ (14)

The velocity field $\mathbf{u}^{k+1}$ is used to update the particles position
at the next time step by solving eq. (12) with the explicit Euler method

$\mathbf{x}_{j}^{k+1}=\mathbf{x}_{j}^{k}+\tau\mathbf{u}_{j}^{k+1}.$ (15)

####  II.2.3 Space discretization

The spatial discretization is accomplished by finite elements. We consider
weak solutions of eqs. (II.2.2) and (14) with $\mathbf{u}\in\mathbf{V}$ and
$p\in L^{2}$, and $\Omega$ replaced by $\Omega_{i}$, such that

$\displaystyle\int_{\Omega_{i}}\rho\left(\frac{1}{\tau}\mathbf{u}^{k+1}+(\mathbf{u}^{k}\cdot\nabla)\mathbf{u}^{k+1}+(\mathbf{u}^{k+1}\cdot\nabla)\mathbf{u}^{k}\right)\mathbf{v}+$
$\displaystyle\nabla\mathbf{v}:\nu^{k}(\nabla\mathbf{u}^{k+1}+(\nabla\mathbf{u}^{k+1})^{T})-p^{k+1}\nabla\cdot\mathbf{v}\,\mathrm{d}\mathbf{x}$
$\displaystyle=\int_{\Omega_{i}}\left(\mathbf{f}^{k}+\rho\left(\frac{1}{\tau}\mathbf{u}^{k}+(\mathbf{u}^{k}\cdot\nabla)\mathbf{u}^{k}\right)\right)\mathbf{v}\,\mathrm{d}\mathbf{x},$
(16)
$\displaystyle\int_{\Omega_{i}}q\nabla\cdot\mathbf{u}^{k+1}\mathrm{d}\mathbf{x}$
$\displaystyle=0,$ (17)

for every $\mathbf{v}\in H^{1}_{0}(\Omega_{i})$ and $q\in L^{2}(\Omega_{i})$.
As for the Stokes equations we partition the domain $\Omega_{i}$ by a
conforming triangulation $\mathcal{T}_{h}$. Then, the continuous spaces
$\mathbf{V}\times L^{2}_{0}$ are approximated by the Taylor-Hood space
$\mathbb{T}_{h}$ with $\Omega$ replaced by $\Omega_{i}$. The velocity
$\mathbf{u}^{k}$ and pressure $p^{k}$ are approximated by functions from
$\mathbb{T}_{h}$, i.e. $m=2$ piece-wise quadratic for $\mathbf{u}^{k}$ and
$m=1$ piece-wise linear for $p^{k}$.

####  II.2.4 Particle constraints and physical parameters

Expressing time in $s$, space in $cm$ and mass in $g$, we used for the density
the value of $\rho=1$, as water, and for the fluid viscosity $\nu_{F}=10^{2}$,
which corresponds to a very viscous fluid. Thus, since the average pore size
$\langle\lambda\rangle$ is about 50 and the average fluid velocity magnitude
$\langle u\rangle$ is about 1, the characteristic Reynolds number of the
system is $Re=\frac{\langle\lambda\rangle\langle u\rangle\rho}{\nu_{F}}=0.5$.
Particles are represented as disks of high viscosity $\nu_{P}=10^{4}$ with
radius $a$ and center of mass $\mathbf{x}_{j}$. The values we have chosen for
$\nu_{F}$ and $\nu_{P}$ are a compromise between computational feasibility and
representing the system physical properties. The relation between $\nu_{F}$
and $\nu_{P}$ follows previous studies tanaka00 .

Particle positions $\mathbf{x}_{j}$ are updated through eq. (15). The new
position results by averaging the velocity field $\mathbf{u}^{k+1}$ over the
high viscosity region associated to the particle. As a consequence the
distance between $\mathbf{x}_{j}$ and $\partial\Omega_{i}$ could be less than
$a$ at the new time instance. To avoid this we impose a constraint on eq.
(15). As soon as the distance between $\partial\Omega_{i}$ and
$\mathbf{x}_{j}$ is less or equal to $a$, we adjust the velocity
$\mathbf{u}^{k}_{j}$ by taking only the tangential component with respect to
the boundary $\partial\Omega_{i}$.

The fluid and particles densities have been matched and set to $\rho=1$ to
neglect any gravitational effect. We consider the potential strength
$\varepsilon=\nu_{P}$ and the interaction range $\sigma=2a$ in order to ensure
effective repulsion. In order to deal with the singularity of Lennard-Jones
potential $V(x)$ at $x=0$ we consider $V(x)=V(2a)$ for $x<2a$. We further
truncate the potential such that $V(x)=0$ for $x>(5/2a)^{2}$.

The domain partition $\mathcal{T}_{h}$ is chosen in such a way that its
elements remain as coarse as possible but are small enough in order to resolve
the flow field. For the fluid particle dynamics simulations we employ adaptive
refinement at the diffused particle interfaces to ensure the mesh size to be
significantly smaller than $\xi$. We chose $\xi=a/4$, where $a$ is the
particle radius. The time step size $\tau$ has been chosen small enough so
that $\tau\mathbf{u}_{j}^{k+1}\sim a$: this prevents particles to move further
then their radius $a$. Considering that at every time step the whole velocity
field is recomputed and particles are advected, we have chosen $\tau=0.5$. A
similar finite element implementation of the fluid particle method, with
similar parameters but simpler geometries can be found in PV_JCP_2015 .

###  II.3 Implementation of the numerical scheme

The numerical simulations are performed by using AMDiS Vey_CVS_2007 ;
witkowski15 . In its current version amdis2  it is based on the DUNE sander20
; Bastian2021  framework https://www.dune-project.org/. Particular features
are adaptive refinement and wrappers to linear system solvers as PETSc. Both
features are extensively used. We choose to discretize the domains $\Omega$
and $\Omega_{i}$ with the grid manager ALUGrid, see alugrid , and use the
solver library MUMPS mumps1 ; mumps2  imported by PETSc as linear solver, see
petsc-web-page . However, also with these advanced tools the fluid particle
dynamics model cannot be solved with reasonable effort on the full domain
$\Omega$ and is therefore considered on sampled subdomains $\Omega_{i}$, shown
in Figure 3.

##  III Results

Our goal is to show how particle confinement, defined by the particle radius
to averaged pore throat size ratio $a/\langle\lambda\rangle$, affects the
velocity field $\mathbf{u}$ and the particles transport itself. As the
particles radius is comparable to the constriction the particles have to pass
through, the pore-throat, the fluid velocity is expected to dynamically change
and to re-route fluid through other porous paths free of solid particles. This
pore-scale phenomenology is also expected to breakdown the assumptions
classically used that the laminar flow, dominated by viscous forces,
transporting suspended particles is stationary (Molnar2005, ; Dentz2011, ;
deAnnaPRF2017, ). We perform several simulations by systematically varying the
particle radius $a\in\\{0,1,2,4\\}$. In Table 1 we report the average pore
throat size for all considered geometries $\Omega_{i}$. We thus impose
confinement conditions $a/\langle\lambda\rangle$ of $\sim 2\%$ ($a=1$), $\sim
5\%$ ($a=2$) and $\sim 10\%$ ($a=4$), corresponding to realistic scenarios for
suspensions moving through soil systems (e.g. a bacterial aggregate can easily
reach a diameter of $10\,\mu$m moving through pores of about $100\,\mu$m).

Geometry | $\Omega$ | $\Omega_{1}$ | $\Omega_{2}$ | $\Omega_{3}$ | $\Omega_{4}$ | $\Omega_{5}$ | $\Omega_{6}$ | $\Omega_{7}$ | $\Omega_{8}$ | $\Omega_{9}$ | $\Omega_{10}$ | $\Omega_{11}$ | $\Omega_{12}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---
$\langle\lambda\rangle$ | 56.791 | 42.539 | 42.673 | 47.410 | 47.924 | 53.106 | 53.729 | 53.633 | 42.748 | 48.904 | 52.591 | 52.127 | 52.996
Table 1: Average throat size $\langle\lambda\rangle$ for $\Omega$ and the
sampled sub-geometries.

In order to compare our simulations with each other we introduce the rescaled
time $t/t^{*}$, where $t^{*}=\langle\lambda_{i}\rangle/\langle u_{S,i}\rangle$
with $\lambda_{i}$ the throat size and

$\langle
u_{S,i}\rangle=\frac{1}{\Omega_{i}}\int_{\Omega_{i}}\lVert\mathbf{u}_{S,i}\rVert\,\mathrm{d}\mathbf{x}$
(18)

the average velocity magnitude of the Stokes flow over the geometry
$\Omega_{i}$, respectively. Then, we study the evolution of the system for
$t\in[0,1000]$, corresponding to a final time $1000/t^{\ast}\in[10,29]$.

Figure 4: _Left_. Snapshots of $u_{r}$ of the fluid particle dynamics model at
different times $t/t^{*}=6.42,19.26,28.90$ (t = 200, 600, 900, from top to
bottom). The velocity field is visualized with the line integral contour
technique (LIC) and color coded by the rescaled velocity using eq. (19).
_Right_. Magnification of the region identified by the black box. In
particular, particles are visualized as dark dots with radius $a$. On top we
observe typical fluid flow through the channel. In the middle magnified
particle-particle interactions are highlighted and the effect on the fluid
flow is visible, e.g. by vortexes on the left side of the channel. On bottom a
particle almost blocks the flow through the pore throat. These results
correspond to $\Omega_{9}$ and confinement condition
$a/\langle\lambda\rangle\sim 10\%$.

Figure 4 shows snapshots of the rescaled velocity magnitude

$u^{k}_{r}:=\log_{10}(\lVert\mathbf{u}^{k}\rVert/\langle u_{S,i}\rangle),$
(19)

for $\Omega_{9}$ at three time instances and a confinement
$a/\langle\lambda\rangle\sim 10\%$. The color-code is such that blue (red)
zones are regions where the velocity magnitude is below (above) the average
velocity magnitude (gray). As for the Stokes solution under no confinement
($a/\langle\lambda\rangle=0$), the velocity field exhibits high velocity
channels and stagnation zones. However, this heterogeneity now changes over
time. In particular, the magnified pictures show that fluid-particle and
particle-particle interactions significantly influences the velocity field.
Laminar vortex structures appear for $t/t^{*}=19.26$ ($t=600$) left of the
pore-throat. They result from the fluid re-routing towards other, more
permeable, paths. The fluid-particle and particle-particle interactions can
also lead to a dramatic flow reduction almost stopping fluid motion in some
pore-throats, as seen at $t/t^{*}=28.90$ ($t=900$), where this situation
results from a particle almost blocking the flow through the pore-throat.

###  III.1 Velocity field

To analyze the temporal variability of the velocity field $\mathbf{u}$, we
consider the mean $u^{k}=\lVert\mathbf{u}^{k}\rVert$ and plot it against the
rescaled time $t/t^{*}$. The PDFs are constructed by considering a fixed bin
range $[b_{l},b_{r}]$ for all $\Omega_{i}$, where $b_{r}=\max_{i}u_{S,i}$ and
$b_{l}=b_{r}/10^{4}$. The bins are logarithmically distributed in
$[b_{l},b_{r}]$. The mean and variance of the PDFs are shown in Figure 5.
These results suggest that the mean velocity only oscillates around a constant
value. However, the variance strongly depends on the confinement and increases
significantly for increasing confinement $a/\langle\lambda\rangle$.

Figure 5: _Top_. From left to right the mean of the PDF of $u^{k}$ for the
different geometries $\Omega_{i}$ for the confinement
$a/\langle\lambda\rangle$ $\sim 2\%$ (_left_), $\sim 5\%$ (_mid_), $\sim 10\%$
(_right_). On the x-axis the rescaled time $t/t^{*}$,
$t^{*}=\langle\lambda\rangle/\langle u_{S,i}\rangle$, where $\langle
u_{S,i}\rangle$ is the average velocity field of the static Stokesian flow
corresponding to the domain $\Omega_{i}$. _Bottom_. The variance of the PDF of
$u^{k}$. It is evident that for small confinement ($\sim 2\%$) the mean and
variance are almost constant, while for larger confinements ($\sim 5\%$ and
$\sim 10\%$) the values oscillate with increasing intensity.

To confirm this we consider the time-average velocity field $\bar{\mathbf{u}}$
and its variance $\sigma_{\mathbf{u}}$, as

$\bar{\mathbf{u}}=\frac{1}{K}\sum_{k=1}^{K}\mathbf{u}^{k},\qquad\sigma_{\mathbf{u}}=\sqrt{\frac{1}{K}\sum_{k=1}^{K}\lVert\mathbf{u}^{k}-\bar{\mathbf{u}}\rVert^{2}}.$
(20)

Figure 6 shows both quantities as functions of the confinement
$a/\langle\lambda\rangle$ for different subdomains $\Omega_{i}$. The space-
time averaged velocities are constant with respect to the confinement, but the
velocity field changes over time and these changes increase with confinement.
While the actual values depend on the considered subdomain $\Omega_{i}$ and
with it the actual pore structure, the qualitative behavior remains the same
for all $\Omega_{i}$.

Figure 6: For the geometries $\Omega_{i}$ in Figure 3 the quantities
$\langle\bar{\mathbf{u}}\rangle$ (_left_) and
$\langle\sigma_{\mathbf{u}}\rangle$ (_right_) are displayed. The numbering and
color coding of the subdomains is ordered with respect the the averaged
velocity.

Both quantities are also shown in Figure 7 for the geometry $\Omega_{9}$ and
particle confinement $a/\langle\lambda\rangle\sim 10\%$. For comparison we
consider also the stationary Stokes equations with point-like particles for
the same geometry and parameters setting. The time average velocity field
corresponds to the solution of the Stokes equations, while its variance
strongly differs in the high velocity channels, which reflects the fact that
particles have been flux-weighted injected.

Figure 7: _Top._ Color plot in logarithmic scale of the rescaled velocity
$u_{r}$ for the Stokes flow with a magnification of the region identified by
the black box. _Middle._ Average flow $u_{r}$ obtained by the fluid particle
dynamics approach. The resulting velocity fields are similar. This is
justified by the fact that high velocity channels of the dynamic flow flip
repetitively, resulting, on average, in similar flow profiles as the Stokes
flow. _Bottom._ Color plot of $\sigma_{\mathbf{u}}$. Notice that the velocity
field has mostly changed within the high velocity channels, see also the
magnified region.

###  III.2 Transport properties

To illustrate the impact of the suspended particles confinement condition not
only on the velocity field, but also on their own transport properties, we
track them over time. Figure 8 visualizes particle streamlines at $1000/t^{*}$
for the geometry $\Omega_{9}$: for point-like particles, and particles of
radius $a=1,2,4$. The starting locations of particles on the left boundary are
the same for all simulations and are flux-weighted (deAnnaPRF2017, ). However,
the trajectories develop in a very different manner depending on the pore
structure. Due to the finite-size of particles and their interaction with the
fluid, trajectories tend to spread over the computational domain. The color
coding is as in Figure 1(d). While for point-like particles the divergence
free velocity field implies that trajectories never cross, for finite size
particles this is no longer true, see the mixing colors in Figure 8. This
effect increases with the suspended particle radius. The trajectory spreading
clearly invades low velocity regions which where poorly invaded by point-
particles.

Figure 8: Particle streamlines for the geometry $\Omega_{9}$ and different
particle radius $a$: point-particles ($a=0$) (top-left), $a=1$ (top-right),
(bottom-left) $a=2$, (bottom-right) $a=4$. The different colors are for
visualization purposes. They correspond to three equal parts of the left
boundary.

We next conduct the same analysis as for the velocity field but now restricted
to the particles. In particular, for each geometry $\Omega_{i}$ we consider
the PDF of the particles velocities magnitude
$u_{j}^{k}=\lVert\mathbf{u}_{j}(t^{k})\rVert$ (see 15) and compute the mean
and variance at every time step $t^{k}$, shown in Figure 9. The mean particle
velocity magnitude, Figure 9 (_top_), is larger than the corresponding mean of
the velocity field in Figure 5 (_top_). This can be explained as particles
most likely are found in high velocity channels. The numerical simulations
suggest that the mean particle velocity only oscillates around a constant
value during the time evolution for $t/t^{*}$ sufficiently large. Differently
from the variance of the velocity field magnitude in Figure 5 (_bottom_), the
variance of the particles velocity magnitude (Figure 9 (_bottom_)) does not
increase systematically with the particle confinement
$a/\langle\lambda\rangle$. This indicates that particles remain in high
velocity channels even if the flow changes. In order to give a cleaner view to
these effects we further average over time and plot the results in Figure 10.
The space-time averaged velocity and its variance remains essentially constant
with respect to the particle confinement.

Figure 9: At each time step the PDF of the particle velocity magnitude is
computed. On the x-axis the rescaled time $t/t^{*}$, where
$t^{*}=\langle\lambda\rangle/\langle u\rangle$ with $u$ magnitude velocity of
the Stokesian flow. Notice that we have different time scales for different
geometries since the throat size average varies with $\Omega_{i}$. On the
y-axis the mean (_top_) and variance (_bottom_) of the PDF for different
particle confinements $\sim 1\%$ (_left_), $\sim 5\%$ (_middle_) and $\sim
10\%$ (_right_) are shown. Figure 10: The space-time average of the particle
velocities (_left_) and their variance (_right_) is shown for all sub-
geometries $\Omega_{i}$ as a function of the particle confinement
$a/\langle\lambda\rangle$.

To give a quantitative measure of the particles displacement, we compute the
average propagation particles velocity along the longitudinal, $x$, direction.
Let $t^{0}_{j}$ and $t^{*}_{j}$ be the spawn and stop time of the individual
particle $j$. The spawn time corresponds to the injection time of particle $j$
at the left side of the domain, the stop time $t^{*}_{j}$ corresponds either
to end of the simulation time or the time at which the particle $j$ reaches
the right side of the computational domain. Let $[x_{j}(t),y_{j}(t)]$ be the
position of the particle $j$ at time $t$. Then, the average propagation
velocity along the $x$-direction is given by

$\bar{v}=\frac{1}{J}\sum_{j=1}^{J}v_{j},\quad
v_{j}=\frac{x_{j}(t^{*}_{j})-x_{j}(t^{0}_{j})}{t^{*}_{j}-t^{0}_{j}}.$

Table 2 gives $\bar{v}$ for the different geometries and radii. Notice that as
the particle confinement increases, the propagation velocity decreases for
most of the considered geometries.

Radius | $\Omega_{1}$ |  | $\Omega_{2}$ |  | $\Omega_{3}$ |  | $\Omega_{4}$ |  | $\Omega_{5}$ |  | $\Omega_{6}$ |
---|---|---|---|---|---|---|---|---|---|---|---|---
$0$ | $0.790\,607\,159\,698\,867\,5$ |  | $1.542\,765\,536\,383\,29$ |  | $1.265\,744\,386\,972\,096\,8$ |  | $1.697\,531\,690\,790\,071\,5$ |  | $1.905\,573\,392\,489\,862$ |  | $1.469\,853\,205\,658\,250\,3$ |
$1$ | $0.721\,349\,562\,610\,800\,3$ | ($0.087\,600\,52$) | $1.500\,168\,091\,515\,123$ | ($0.027\,611\,09$) | $1.175\,621\,818\,440\,004\,2$ | ($0.071\,201\,24$) | $1.687\,223\,507\,247\,289\,9$ | ($0.006\,072\,45$) | $1.765\,176\,159\,764\,037\,4$ | ($0.073\,677\,16$) | $1.460\,523\,196\,709\,400\,3$ | ($0.006\,347\,58$)
$2$ | $0.707\,961\,188\,062\,446\,5$ | ($0.104\,534\,81$) | $1.517\,819\,191\,873\,465$ | ($0.016\,169\,89$) | $1.191\,876\,900\,004\,521\,7$ | ($0.058\,358\,93$) | $1.734\,134\,035\,747\,872\,3$ | ($-0.021\,562\,1$) | $1.789\,152\,244\,063\,490\,2$ | ($0.061\,095\,07$) | $1.445\,044\,607\,286\,734\,7$ | ($0.016\,878\,28$)
$4$ | $0.606\,451\,564\,921\,749\,2$ | ($0.232\,929\,33$) | $1.383\,742\,333\,544\,798\,3$ | ($0.103\,076\,71$) | $1.106\,712\,562\,096\,999\,8$ | ($0.125\,642\,92$) | $1.688\,467\,546\,921\,600\,2$ | ($0.005\,339\,6$) | $1.887\,946\,101\,198\,029$ | ($0.009\,250\,39$) | $1.303\,872\,851\,411\,751$ | ($0.112\,923\,08$)

Radius | $\Omega_{7}$ |  | $\Omega_{8}$ |  | $\Omega_{9}$ |  | $\Omega_{10}$ |  | $\Omega_{11}$ |  | $\Omega_{12}$ |
---|---|---|---|---|---|---|---|---|---|---|---|---
$0$ | $1.832\,855\,034\,108\,713$ |  | $1.871\,619\,794\,646\,815\,6$ |  | $2.247\,769\,531\,474\,388\,4$ |  | $2.362\,920\,712\,143\,502\,3$ |  | $1.956\,479\,289\,056\,703\,5$ |  | $1.805\,464\,197\,911\,871\,3$ |
$1$ | $1.846\,836\,151\,799\,182\,5$ | ($-0.007\,628\,05$) | $1.961\,274\,951\,056\,896$ | ($-0.047\,902\,44$) | $2.220\,330\,848\,844\,713\,8$ | ($0.012\,207\,07$) | $2.283\,732\,214\,992\,107$ | ($0.033\,512\,97$) | $1.943\,131\,649\,125\,236\,6$ | ($0.006\,822\,28$) | $1.776\,941\,447\,054\,544$ | ($0.015\,798\,02$)
$2$ | $1.737\,394\,822\,791\,188\,4$ | ($0.052\,082\,79$) | $1.854\,612\,679\,764\,501$ | ($0.009\,086\,84$) | $2.020\,958\,351\,258\,437\,3$ | ($0.100\,905$) | $2.234\,303\,321\,560\,587\,6$ | ($0.054\,431\,53$) | $1.965\,344\,041\,150\,803\,9$ | ($-0.004\,530\,97$) | $1.775\,582\,651\,977\,670\,5$ | ($0.016\,550\,62$)
$4$ | $1.494\,538\,644\,317\,411$ | ($0.184\,584\,37$) | $1.588\,451\,043\,014\,832\,4$ | ($0.151\,296\,09$) | $1.494\,998\,111\,065\,284\,4$ | ($0.334\,897\,07$) | $2.046\,028\,142\,185\,807$ | ($0.134\,110\,54$) | $1.813\,442\,877\,216\,675\,3$ | ($0.073\,109\,09$) | $1.745\,707\,663\,168\,257\,3$ | ($0.033\,097\,6$)

Table 2: Average propagation velocity $\bar{v}$ along the $x$-direction of
particles for different geometries $\Omega_{i}$ and radii $a$. In brackets the
ratio between the average propagation velocity of the particle of radius $a$
and the corresponding average velocity of the point-particle:
$1-\bar{v}_{a}/\bar{v}_{0}$, where $\bar{v}$ is the average propagation
velocity along the $x$-direction. A negative ratio means that the particles
with radius $a$ propagated faster than the corresponding point-particle.

##  IV Conclusions

We have simulated the transport of soft particles of finite size through 2D
porous geometries systematically changing the confinement conditions,
represented by the ratio between the particle radius and the average pore-
throat size of the pore structure. The employed numerical approach is the
fluid particle dynamics approach introduced in tanaka00 . The numerical scheme
allows to simulate the Navier-Stokes equations with particle interactions
without explicitly treating particle boundaries. A limiting factor of this
approach is that the diffused particle interfaces need to be resolved. Our
results show that the particles confinement affects the fluid velocity field
$\mathbf{u}$ which in turn affects the particles transport. We show that these
fluid-particles interactions impact the macroscopic transport. As shown in
Figure 8, even for small confinement conditions, fluid and transported
particles are dynamically re-routed towards more permeable paths. This leads
to the emergence of ephemeral laminar vortexes at pore throat entrances and
affects the variance and mean fluid velocity. While these results have been
proven for a fluid more viscous than water, we expect these confinement
effects to be even more pronounced when solving the problem for water.

The complexity of the fluid-suspended particles interaction, flow
hydrodynamics, the coupling between particles size, local fluid velocity and
associated transport have been widely overlooked. Our results estimate the
effect of these phenomena on macroscopic flow kinematics and transport which
are relevant for mixing, reaction kinetics and filtration. We expect these
results to be also of relevance in more complex scenarios where also
morphological variability, in terms of grains shape (Alhashmi2016, ;
Xiong2016, ; Wu2019, ) plays critical roles, e.g. in groundwater contamination
and remediation Kahler2019 , enhanced hydrocarbon recovery Kar2015 , transport
through river sediments that create a closely packed pore network Lei2022 ,
water filtration systems Kosvintsev2002  and extra-cellular transport in brain
tissue Nicholson2017 . Also, for pores surrounded by a single grain that
cannot host a net transfer of fluid, the so called dead-end pores, which have
been recently shown to host laminar flow vortexes and are able to trap fluid
for long times Bordoloi2022 , it would be interesting to estimate the impact
of fluid-suspended particles interaction and flow hydrodynamics.

Acknowledgments: This work was funded by the EU H2020 program within FET-OPEN
project NARCISO (Grant Agreement no. 828890). We further acknowledge computing
resources at FZ Jülich under grant PFAMDIS and at ZIH under grant WIR.

## References

  * [1] J. Bear.  Dynamics of Fluids in Porous Media.  Elsevier, 1972.
  * [2] I. L. Molnar, W. P. Johnson, J. I. Gerhard, C. S. Willson, and D. M. O’Carroll.  Predicting colloid transport through saturated porous media: A critical review.  Water Resour. Res., 51:6804–6845, 2005.
  * [3] M. Dentz, T. Le Borgne, A. Englert, and B. Bijeljic.  Mixing, spreading and reaction in heterogeneous media: A brief review.  J. Contam. Hydrol., 120–121, 2011.
  * [4] R. A. Freeze and J. A. Cherry.  Groundwater.  Prentice Hall, 1979.
  * [5] P. de Anna, B. Quaife, G. Biros, and R. Juanes.  Prediction of velocity distribution from pore structure in simple porous media.  Phys. Rev. Fluids, 2:124103, 2017.
  * [6] M. Dentz, M. Icardi, and J. J. Hidalgo.  Mechanism of dispersion in a porous medium.  J. Fluid Mech., 841:851–882, 2018.
  * [7] B. Berkowitz, A. Cortis, M. Dentz, and H. Scher.  Modeling non-Fickian transport in geological formations as a continuous time random walk.  Rev. Geophys., 44:1–49, 2006.
  * [8] P. de Anna, J. Jimenez-Martinez, H. Tabuteau, R. Turuban, T. Le Borgne, M. Derrien, and Y. Méheust.  Mixing and reaction kinetics in poroous media: an experimental pore scale quantification.  Environ. Sci. Technol., 48:508–516, 2014.
  * [9] J. Heyman, D. R. Lester, R. Turuban, Y. M’eheust, and T. Le Borgne.  Stretching and folding sustain microscope chemical gradients in porous media.  Proc. Nat. Acad. Sci. (USA), 117:13359–13365, 2020.
  * [10] N. Nishiyama, T. Yokoyama, and S. Takeuchi.  Size distributions of pore water and entrapped air duringdrying-infiltration processes of sandstone characterized bywater-expulsion porosimetry.  Water Resour. Res., 48:W09556, 2012.
  * [11] F. Miele, P. de Anna, and M. Dentz.  Stochastic model for filtration by porous materials.  Phys. Rev. Fluids, 4:094101, 2019.
  * [12] D. Scheidweiler, F. Miele, H. Peter, T. J. Battin, and P. de Anna.  Trait-specific dispersal of bacteria in heterogeneous porous environments: from pore to porous medium scale.  J. Roy. Soc. Interface., 17:20200046, 2020.
  * [13] P. de Anna, A. A. Pahlavan, Y. Yawata, R. Stocker, and R. Juanes.  Chemotaxis under flow disorder shapes microbial dispersion in porous media.  Nat. Phys., 17:68–73, 2021.
  * [14] E. Guazzelli and J. F. Morris.  A physical introduction to suspension dynamics.  Cambridge University Press, 2012.
  * [15] H. C. Berg.  Random Walks in Biology.  Princeton University Press, 1993.
  * [16] Sirivithayapakorn S. and A. Keller.  Transport of colloids in saturated porous media: A pore-scale observation of the size exclusion effect and colloid acceleration.  Water Resour. Res., 9:1109, 2003.
  * [17] M. Auset and A. A. Keller.  Pore‐scale processes that control dispersion of colloids in saturated porous media.  Water Resour. Res., 40:W03503, 2004.
  * [18] D. C. Prieve and P. M. Hoysan.  Role of colloidal forces in hydrodynamic chromatography.  J. Colloid Interface Sci., 64:201–213, 1978.
  * [19] P. Babakhani.  The impact of nanoparticle aggregation on their size exclusion during transport in porous media: one- and three-dimensional modelling investigations.  Sci. Rep., 9:14071, 2019.
  * [20] A. A. Keller, S. Sirivithayapakorn, and C. V. Chrysikopoulos.  Early breakthrough of colloids and bacteriophage MS2 in a water‐saturated sand column.  Water Resour. Res., 40:W08304, 2004.
  * [21] S. C. James and C. V. Chrysikopoulos.  Effective velocity and effective dispersion coefficient for finite-sized particles flowing in a uniform fracture.  J. Colloid Interface Sci., 263:288–295, 2003.
  * [22] H. Tanaka and Ta. Araki.  Simulation method of colloidal suspensions with hydrodynamic interactions: Fluid particle dynamics.  Phys. Rev. Lett., 85:1338–1341, 2000.
  * [23] S. Praetorius and A. Voigt.  A Navier-Stokes phase-field crystal model for colloidal suspensions.  J. Chem. Phys., 142:154904, 2015.
  * [24] S. Vey and A. Voigt.  AMDiS: Adaptive multidimensional simulations.  Comput. Vis. Sci., 10:57–67, 2007.
  * [25] T. Witkowski, S. Ling, S. Praetorius, and A. Voigt.  Software concepts and numerical algorithms for a scalable adaptive parallel finite element method.  Adv. Comput. Math., 41:1145–1177, 2015.
  * [26] S. Praetorius.  The Adaptive Multi-Dimensional simulation toolbox (AMDiS), a discretization module on top of the Dune framework.  https://gitlab.com/amdis/amdis.
  * [27] O. Sander.  Dune — The Distributed and Unified Numerics Environment.  Springer, 2020.
  * [28] The Dune framework: Basic concepts and recent developments.  Comput. Math. with Appl., 81:75–112, 2021.
  * [29] M. Alkämper, A. Dedner, R. Klöfkorn, and M. Nolte.  The Dune ALUGrid module.  Arch. Numerical Software, 4:1–28, 2016.
  * [30] P. R. Amestoy, I. S. Duff, J. Koster, and J.-Y. L’Excellent.  A fully asynchronous multifrontal solver using distributed dynamic scheduling.  SIAM J. Matrix Anal. Appl., 23:15–41, 2001.
  * [31] P. R. Amestoy, A. Guermouche, J.-Y. L’Excellent, and S. Pralet.  Hybrid scheduling for the parallel solution of linear systems.  Parallel Comput., 32:136–156, 2006.
  * [32] S. Balay, S. Abhyankar, M. F. Adams, S. Benson, J. Brown, P. Brune, K. Buschelman, E. M. Constantinescu, L. Dalcin, A. Dener, V. Eijkhout, W. D. Gropp, V. Hapla, T. Isaac, P. Jolivet, D. Karpeev, D. Kaushik, M. G. Knepley, F. Kong, S. Kruger, D. A. May, L. C. McInnes, R. T. Mills, L. Mitchell, T. Munson, J. E. Roman, K. Rupp, P. Sanan, J. Sarich, B. F. Smith, S. Zampini, H. Zhang, H. Zhang, and J. Zhang.  PETSc Web page.  https://petsc.org/, 2022.
  * [33] Z. Alhashmi, M. J. Blunt, and B. Bijeljic.  The impact of pore structure heterogeneity, transport, and reaction conditions on fluid–fluid reaction rate studied on images of pore space.  Transport in Porous Media, 115:215–237, 2016.
  * [34] Q. Xiong, T. G. Baychev, and A. P. Jivkov.  Review of pore network modelling of porous media: Experimental characterisations, network constructions and applications to reactive transport.  J. Contam. Hydrol., 192:101–117, 2016.
  * [35] H. Wu, W.-Z. Fang, Q. Kang, W.-Q. Tao, and R. Qiao.  Predicting effective diffusivity of porous media from images by deep learning.  Sci. Rep., 9:20387, 2019.
  * [36] D. M. Kahler and Z. J. Kabala.  Acceleration of groundwater remediation by rapidly pulsed pumping: Laboratory column tests.  J. Environmental Eng., 145:06018009, 2019.
  * [37] A. Kar, T. Chiang, I. Ortiz Rivera, A. Sen, and D. Velegol.  Enhanced transport into and out of dead-end pores.  ACS Nano, 9:746–753, 2015.
  * [38] L. Lei, T. Park, K. Jarvis, L. Pan, I. Tepecik, Y. Zhao, Z. Ge, J.-H. Choi, X. Gai, S. A. Galindo-Torres, R. Boswell, S. Dai, and Y. Seol.  Pore-scale observations of natural hydrate-bearing sediments via pressure core sub-coring and micro-CT scanning.  Sci. Rep., 12:3471, 2022.
  * [39] S. Kosvintsev, R.G. Holdich, I.W. Cumming, and V.M. Starov.  Modelling of dead-end microfiltration with pore blocking and cake formation.  J. Membr. Sci., 208:181–192, 2002.
  * [40] C. Nicholson and S. Hrabetova.  Brain extracellular space: The final frontier of neuroscience.  Biophys. J., 113:2133–2142, 2017.
  * [41] A. D. Bordoloi, D. Scheidweiler, M. Dentz, M. Bouabdellaoui, Abbarchi M., and P. de Anna.  Structure induced laminar vortexes control anomalous dispersion in porous media.  Nat. Commun., 13:3820, 2022.