# Numerical Equivalence Between SPH and Probabilistic Mass Transfer Methods
for Lagrangian Simulation of Dispersion111This work was partially supported by
the US Army Research Office under Contract/Grant number W911NF-18-1-0338; the
National Science Foundation under awards EAR-1417145 and DMS-1614586; and the
Spanish Ministry of Economy and Competitiveness through project WE-NEED,
PCIN-2015-248.

Guillem Sole-Mari222Department of Civil and Environmental Engineering (DECA),
Universitat Politècnica de Catalunya, Barcelona, Spain333Hydrogeology Group
(GHS), UPC-CSIC, Barcelona, Spain<EMAIL_ADDRESS>Universitat
Politècnica de Catalunya
C/ Jordi Girona 1-3
08034 Barcelona, Spain  Michael J. Schmidt444Department of Applied Mathematics
and Statistics, Colorado School of Mines, Golden, CO, 80401, USA555Hydrologic
Science and Engineering Program, Department of Geology and Geological
Engineering, Colorado School of Mines, Golden, CO, 80401, USA
<EMAIL_ADDRESS>Stephen D. Pankavich666Department of Applied Mathematics
and Statistics, Colorado School of Mines, Golden, CO, 80401, USA
<EMAIL_ADDRESS>David A. Benson777Hydrologic Science and Engineering
Program, Department of Geology and Geological Engineering, Colorado School of
Mines, Golden, CO, 80401, USA<EMAIL_ADDRESS>Colorado School of Mines
1500 Illinois St.
Golden, CO 80401

###### Abstract

Several Lagrangian methodologies have been proposed in recent years to
simulate advection-dispersion of solutes in fluids as a mass exchange between
numerical particles carrying the fluid. In this paper, we unify these
methodologies, showing that mass transfer particle tracking (MTPT) algorithms
can be framed within the context of smoothed particle hydrodynamics (SPH),
provided the choice of a Gaussian smoothing kernel whose bandwidth depends on
the dispersion and the time discretization. Numerical simulations are
performed for a simple dispersion problem, and they are compared to an
analytical solution. Based on the results, we advocate for the use of a kernel
bandwidth of the size of the characteristic dispersion length
$\ell=\sqrt{2D\Delta t}$, at least given a “dense enough” distribution of
particles, for in this case the mass transfer operation is not just an
approximation, but in fact the exact solution, of the solute’s displacement by
dispersion in a time step.

###### keywords:

Lagrangian Modeling , Dispersion , Smoothed Particle Hydrodynamics , Mass
Transfer Particle Tracking , Kernel Bandwidth

††journal: Advances in Water Resources

##  1 Introduction

In recent years, a number of Lagrangian numerical schemes have been proposed
to simulate advection-dispersion processes in fluids. Some of these approaches
rely exclusively on traditional random walks to simulate dispersion [1, 2, 3,
4, 5, 6, 7, 8, 9, 10, 11, 12], whereas a second class represents dispersion
through mass transfer between particles that carry a given amount of fluid
[13, 14, 15, 16]. Other authors have hybridized random walks with mass
transfer [17, 18] in an approach that allows partitioning of total dispersion
between mixing (simulated by mass transfer) and non-mixed spreading (simulated
via random walks). Mass-transfer algorithms can be further subdivided into two
groups. The first group [13, 14] derives the mass exchange rates from the
well-established smoothed particle hydrodynamics (SPH) method [19], which,
besides solute transport, has been used in a variety of applications [20] such
as astrophysics, fluid dynamics, and solid mechanics. A second group of
approaches, often referred to as mass transfer particle tracking (MTPT)
algorithms [15, 16], derive the mass-exchange rate from stochastic rules
governing the co-location probability of particles moving via dispersion. To
date, a relationship between these two methodologies for mass transfer has not
been established. In this paper we analytically derive the connection between
the SPH and MTPT conventions and show that, for specific kernel choices and
provided that equivalent normalization and averaging conventions are used, the
SPH and MTPT approaches are numerically equivalent. Additionally, for the
fixed choice of a Gaussian kernel, we investigate the effect of differing
bandwidth choices on deviations from the analytical, well-mixed solution.

##  2 The link between SPH and MTPT

The SPH approach to approximating dispersion can be summarized by following
[13, 14]. Therein, the following equation describes the time evolution of the
concentration, $C_{i}(t)$, carried by a numerical particle labeled $i=1,..,N$,
assuming that all particles contain the same amount of fluid:

$\frac{dC_{i}}{dt}=2\sum^{N}_{j=1}{\frac{{\widehat{D}}_{ij}}{{\widehat{\rho}}_{ij}}\left(C_{i}-C_{j}\right)F\left({\boldsymbol{X}}_{i}-{\boldsymbol{X}}_{j};h\right)}.$
(1)

Here, $N$ is the number of particles, $\boldsymbol{X}_{i}$ is the position of
particle $i$, and $F\left(\boldsymbol{r};h\right)$ is a radial function
satisfying

$\boldsymbol{r}F\left(\boldsymbol{r};h\right)=\nabla
W\left(\boldsymbol{r};h\right),$ (2)

with $W$ representing a radially symmetric, translation-invariant kernel with
bandwidth $h$. Additionally, ${\widehat{D}}_{ij}$ is the effective dispersion
coefficient that, in the isotropic but spatially variable case, reduces to

${\widehat{D}}_{ij}\coloneqq
g\left(D({\boldsymbol{X}}_{i}),D({\boldsymbol{X}}_{j})\right),$ (3)

where $g$ is an averaging function (e.g., arithmetic or harmonic average). The
quantity ${\widehat{\rho}}_{ij}$, defined by

$\displaystyle\widehat{\rho}_{ij}$ $\displaystyle\coloneqq
g\left(\rho_{i},\rho_{j}\right),$ (4)
$\displaystyle\rho\left(\boldsymbol{X};h\right)$
$\displaystyle\coloneqq\sum_{k=1}^{N}W\left({\boldsymbol{X}}-{\boldsymbol{X}}_{k};h\right),\quad\rho_{q}\coloneqq\rho\left(\boldsymbol{X}_{q};h\right),\quad
q=i,j,$ (5)

is an average of the particle densities estimated at ${\boldsymbol{X}}_{i}$
and ${\boldsymbol{X}}_{j}$. A popular choice for $g$, in this case, is the
arithmetic average. Note that we make explicit the previously suppressed
dependence of $\widehat{\rho}_{ij}$ on the locations of the particles,
$\boldsymbol{X}_{i},\ i=1,\dots,N$, and the parameter $h$, which represents
the bandwidth of the kernel function $W$.

In the specific case that $W\left(\boldsymbol{r};h\right)$ is a Gaussian
kernel with the form

$W\left(\boldsymbol{r};h\right)={\left(2\pi h^{2}\right)}^{-d/2}\
\exp\left({-\frac{\left|{\boldsymbol{r}}\right|^{2}}{2h^{2}}}\right),$ (6)

where $d$ is the number of spatial dimensions, we have

$F\left(\boldsymbol{r};h\right)=-\frac{1}{h^{2}}W\left(\boldsymbol{r};h\right).$
(7)

Substituting (7) into (1), integrating the expression (first-order explicit),
and then rearranging, we arrive at

${C}_{i}\left(t+\Delta
t\right)={C}_{i}\left(t\right)+\sum^{N}_{j=1}{{\beta}_{ij}\mathcal{W}_{ij}\left({C}_{j}\left(t\right)-{C}_{i}\left(t\right)\right)},$
(8)

in which we define

$\displaystyle{\beta}_{ij}\left(h\right)$
$\displaystyle\coloneqq\frac{\ell_{ij}^{2}}{h^{2}},\quad\ell_{ij}\coloneqq\sqrt{2{\widehat{D}}_{ij}\Delta
t},$ (9) $\displaystyle\mathcal{W}_{ij}\left(h\right)$
$\displaystyle\coloneqq\frac{W\left({\boldsymbol{X}}_{i}-{\boldsymbol{X}}_{j};h\right)}{{\widehat{\rho}}_{ij}\left(h\right)}.$
(10)

Here, we once again denote the dependence of $\beta$ and $\mathcal{W}$ on $h$
because, for a different kernel bandwidth choice, these quantities will be
altered correspondingly. Note that $\ell_{ij}$ in (9) is equal to the
characteristic distance of the average dispersion of particles $i$ and $j$ in
a time step $\Delta t$.

We now consider, for the sake of comparison, the MTPT algorithm originally
formulated by _Benson and Bolster_ [15], further discussed in [16], and given
by

${C}_{i}\left(t+\Delta
t\right)={C}_{i}\left(t\right)+\frac{1}{2}\sum^{N}_{j=1}\mathcal{P}_{ij}\left({C}_{j}\left(t\right)-{C}_{i}\left(t\right)\right),$
(11)

where $\mathcal{P}_{ij}$ is the probabilistic weighting function for a mass
transfer from particle $j$ to particle $i$, with the form

$\mathcal{P}_{ij}=\frac{P(\boldsymbol{X}_{i}-\boldsymbol{X}_{j})}{\widetilde{\rho}_{ij}}.$
(12)

Here, the function $P$ is the probability density for the co-location of
particles $i$ and $j$, moving via dispersion,

$\displaystyle P(\boldsymbol{X}_{i}-\boldsymbol{X}_{j};D,\Delta t)$
$\displaystyle=\left(4\pi(D_{i}+D_{j})\Delta
t\right)^{-d/2}\exp\left[-\frac{|\boldsymbol{X}_{i}-\boldsymbol{X}_{j}|^{2}}{4(D_{i}+D_{j})\Delta
t}\right]$ (13) $\displaystyle\equiv
W\left(\boldsymbol{X}_{i}-\boldsymbol{X}_{j};\sqrt{2(D_{i}+D_{j})\Delta
t}\right)$
$\displaystyle=W\left(\boldsymbol{X}_{i}-\boldsymbol{X}_{j};\sqrt{4\widehat{D}_{ij}\Delta
t}\right)$
$\displaystyle=W\left(\boldsymbol{X}_{i}-\boldsymbol{X}_{j};\sqrt{2}\ell_{ij}\right),$

where $D_{k}\coloneqq D(\boldsymbol{X}_{k})$, and $\widetilde{\rho}_{ij}$ is a
normalizing factor that has classically been chosen to be $\rho_{j}$, as in
(5), with $h=\sqrt{2}\ell_{ij}$, in order to make the matrix
$\boldsymbol{\mathcal{P}}$ (with $i,j^{\text{th}}$ entry $\mathcal{P}_{ij}$) a
left stochastic matrix (i.e., a matrix where all columns sum to 1). However,
this does not guarantee that $\mathcal{P}_{ij}=\mathcal{P}_{ji}$. Hence, the
concentration increase (or decrease) at particle $i$ due to its interaction
with particle $j$ by (11) may not match the decrease (or increase) at particle
$j$ due to interaction with particle $i$. As a consequence of this asymmetry,
normalization by $\widetilde{\rho}_{ij}=\rho_{j}$ may not impose exact mass
conservation. Also, in the original paper [15], equation (11) is formulated in
terms of solute masses instead of concentrations, which are, in this case,
interchangeable since all particles carry an equal amount of fluid.

Comparing equations (8) and (11), it is evident that the co-location
probability-based mass exchange algorithm of _Benson and Bolster_ [15] is
numerically equivalent to the SPH formalism for $\beta_{ij}=1/2$ for all
$i,j=1,...,N$, with the standard deviation associated with particle co-
location by dispersion used for the bandwidth of $W$ in (6). Note that,
according to (9) and (10), imposing a constant value for ${\beta}_{ij}$
implies that the kernel bandwidth $h$ will change with the positions of
particles $i$ and $j$ for spatially-variable dispersion and will depend on
$\Delta t$, as can be seen from (13).

Expressions (8) and (11) can be written in a general matrix-vector form as

$\boldsymbol{{C}}\left(t+\Delta
t\right)=\boldsymbol{A}\left(t\right)\boldsymbol{{C}}\left(t\right),$ (14)

where $\boldsymbol{{C}}_{i}\coloneqq{C}_{i}$, and

$\boldsymbol{A}\coloneqq\boldsymbol{I}+\left[\boldsymbol{\beta}\circ\boldsymbol{\mathcal{W}}-\operatorname{diag}\left(\left[\boldsymbol{\beta}\circ\boldsymbol{\mathcal{W}}\right]\boldsymbol{1}\right)\right].$
(15)

Above, $\boldsymbol{I}$ is the $N\times N$ identity matrix, $\boldsymbol{1}$
is an $N\times 1$ vector of ones, $\circ$ denotes the entrywise, or Hadamard,
product, $\operatorname{diag}(\boldsymbol{x})$ is a square matrix with the
entries of $\boldsymbol{x}$ on its main diagonal, and the $i,j^{\text{th}}$
entries of the matrices $\boldsymbol{\beta}$ and $\boldsymbol{\mathcal{W}}$
are $\beta_{ij}$ and $\mathcal{W}_{ij}$, respectively. Note that, as mentioned
above and elsewhere (see [16]), choosing $\widetilde{\rho}_{ij}$ to be
$\rho_{j}$ in (12) ensures that $\boldsymbol{\mathcal{P}}$ (denoted as
$\boldsymbol{\mathcal{W}}$ in (15)) is a left stochastic matrix (but not
necessarily symmetric). On the other hand, we note from (15) that if
$\boldsymbol{\mathcal{W}}$ is symmetric, then $\boldsymbol{A}$ is also
symmetric with rows and columns that sum to 1, guaranteeing conservation of
mass. Thus, a better normalization approach is to choose
$\widetilde{\rho}_{ij}$ to be $\widehat{\rho}_{ij}$, as in (4), resulting in
symmetric $\boldsymbol{\mathcal{W}}$ and mass-conserving $\boldsymbol{A}$.

_Schmidt et al._ [16, 21] also present a discretized Green’s function approach
to simulating dispersion by mass transfer. For a time step $\Delta t$, this
algorithm is described as:

$\boldsymbol{{C}}\left(t+\Delta
t\right)=\boldsymbol{\mathcal{G}}\left(t\right)\boldsymbol{{C}}\left(t\right),$
(16)

with

$\displaystyle\mathcal{G}_{ij}$
$\displaystyle\coloneqq\frac{W\left({\boldsymbol{X}}_{i}-{\boldsymbol{X}}_{j};\ell_{ij}\right)}{\widetilde{\rho}_{ij}},$
(17)

where, once again, $\widetilde{\rho}_{ij}$ is traditionally defined to be
$\rho_{j}$ as in (5). We see that the matrix $\boldsymbol{\mathcal{G}}$ is
nearly identical to $\boldsymbol{\mathcal{W}}$ for $h=\ell_{ij}$ (and to
$\boldsymbol{\mathcal{P}}$ with twice the square bandwidth), the only
difference being the choice of non-symmetric normalization using
$\widetilde{\rho}_{ij}=\rho_{j}$.

We note that, for a sufficiently large $N$,
$\sum_{j=1}^{N}\mathcal{G}_{ij}\approx\int\rho\left(\boldsymbol{x}\right)\frac{W\left(\boldsymbol{X}_{i}-\boldsymbol{x}\right)}{\rho\left(\boldsymbol{x}\right)}\mathrm{d}\boldsymbol{x}=1$,
which implies
$\operatorname{diag}\left(\boldsymbol{\mathcal{G}}\boldsymbol{1}\right)\approx\boldsymbol{I}$.
Hence, knowing that, for $h=\ell_{ij}$,
$\boldsymbol{\mathcal{W}}\approx\boldsymbol{\mathcal{G}}$, we see that the
discretized Green’s function algorithm (16) is also nearly identical to the
SPH and particle co-location expression given in (14) and (15), under the
constraint that ${\beta}_{ij}=1$ for all $i,j$. Hereafter, for simplicity, we
refer to any matrix $\boldsymbol{\beta}$ with all-equal entries as a scalar
$\beta$.

Thus, we have unified the previously divergent approaches to simulating
dispersion that are employed by the SPH and MTPT algorithms. Namely, to frame
things in the SPH context, the MTPT algorithms hold the mass-transfer scaling
parameter $\beta$ constant ($1/2$ or $1$) and adapt the kernel, itself, to the
magnitude of dispersion over a time step. This is in contrast to the
traditional SPH approach, where the kernel bandwidth is independent from the
dispersion magnitude, and often set to contain a prescribed “number of
neighbors”, either locally or on average [22]. The kernel is then scaled in
amplitude by the parameter $\beta_{ij}$ to capture the magnitude of the
dispersion action.

Having established the link, through the parameter $\boldsymbol{\beta}$ in
(15) (alternatively viewed as the choice of kernel bandwidth $h$), between the
SPH and MTPT formalisms for simulating dispersion in a Lagrangian context, we
next consider the implications of varying this parameter. In the following
section, we conduct some numerical experiments to consider these effects.

##  3 Numerical investigations

Figure 1: Concentrations at $t=4$ for example realizations with  evenly-
spaced,  randomly-spaced, and  random-walking particles. The black line is the
analytical solution. For these simulations $N=255$, $\Delta t=0.01$ and
$h=h_{*}$ in each case. The initial condition is a Dirac delta positioned at
the center of the domain. Note that deviations from the analytical solution
are caused by irregular, and possibly wide, particle spacings, due to low
particle numbers in the non-equally-spaced cases. Figure 2: Numerical results
for evenly-spaced, stationary particles. $\mathrm{RMSE}$ (18), as a function
of the kernel bandwidth $h$, is given for different combinations of $N$ and
$\Delta t$. The dashed, semitransparent vertical lines indicate the values of
$h$ that correspond to $\beta=1$ ($h=\ell=\sqrt{2D\Delta t}$) for a given
value of $\Delta t$. Figure 3: Numerical results for evenly-spaced, stationary
particles. $\left(\mathrm{a}\right)$ Bandwidth $h_{*}$ associated with the
minimum $\mathrm{RMSE}$ plotted against the particle spacing $s=L/N$, for
different $\Delta t$ values (see color legend on Figure 2).
$\left(\mathrm{b}\right)$ Bandwidth $h_{*}$ associated with the minimum
$\mathrm{RMSE}$ plotted against the dispersion distance $\ell=\sqrt{2D\Delta
t}$ given different $N$ values (see marker legend on Figure 2).
$\left(\mathrm{c}\right)$ Lowest-error bandwidth $h_{*}$ against particle
spacing $s=L/N$, both normalized by the dispersion distance
$\ell=\sqrt{2D\Delta t}$. Figure 4: Numerical results for randomly-spaced,
stationary particles. Averaged $\mathrm{RMSE}$ (18), as a function of the
kernel bandwidth $h$, is given for different combinations of $N$ and $\Delta
t$. The dashed, semitransparent vertical lines indicate the values of $h$ that
correspond to $\beta=1$ ($h=\ell=\sqrt{2D\Delta t}$) for each value of $\Delta
t$. Figure 5: Numerical results for randomly-spaced, stationary particles.
$\left(\mathrm{a}\right)$ Bandwidth $h_{*}$ associated with the minimum
$\mathrm{RMSE}$ plotted against the average particle spacing $\langle
s\rangle=L/N$, for different $\Delta t$ values (see color legend on Figure 4).
$\left(\mathrm{b}\right)$ Bandwidth $h_{*}$ associated with the minimum
$\mathrm{RMSE}$ plotted against the dispersion distance $\ell=\sqrt{2D\Delta
t}$ given different $N$ values (see marker legend on Figure 4).
$\left(\mathrm{c}\right)$ Lowest-error bandwidth $h_{*}$ against average
particle spacing $\langle s\rangle=L/N$, both normalized by the dispersion
distance $\ell=\sqrt{2D\Delta t}$. Figure 6: Numerical results for random-
walking particles. Averaged $\mathrm{RMSE}$ (18), as a function of the kernel
bandwidth $h$, is given for different combinations of $N$ and $\Delta t$. The
dashed, semitransparent vertical lines indicate the values of $h$ that
correspond to $\beta=1$ for each value of $\Delta t$. Figure 7: Numerical
results for random-walking particles. $\left(\mathrm{a}\right)$ Bandwidth
$h_{*}$ associated with the minimum $\mathrm{RMSE}$ plotted against the
average particle spacing $\langle s\rangle=L/N$, for different $\Delta t$
values (see color legend on Figure 4). $\left(\mathrm{b}\right)$ Bandwidth
$h_{*}$ associated with the minimum $\mathrm{RMSE}$ plotted against the
dispersion distance $\ell=\sqrt{2D_{\text{MT}}\Delta t}$ given different $N$
values (see marker legend on Figure 4). $\left(\mathrm{c}\right)$ Lowest-error
bandwidth $h_{*}$ against average particle spacing $\langle s\rangle=L/N$,
both normalized by the dispersion distance $\ell=\sqrt{2D_{\text{MT}}\Delta
t}$.

To analyze the effect of the kernel bandwidth $h$ on SPH transport
simulations, we study a simple case of 1D constant dispersion, where the
initial condition is a Dirac delta pulse located at the center of the domain,
$x=0.5$ [L]. For simplicity, the model has no units. The dispersion
coefficient is fixed as $D=10^{-3}$ [L2T-1], and the total simulation time is
$T=4$ [T]. The analytical solution is then a Gaussian with variance
$\sigma^{2}=2DT$ (see Figure 1), or, to be more precise, the analytical
solution is a normalized, $N$-bin histogram populated with evaluations of the
density of a Normal distribution, $\mathcal{N}(0.5,\sigma^{2})$, at the
positions of the particles. We compare this analytical solution to the
numerical results for a range of values of $h$, $N$ and $\Delta t$, using
root-mean-squared error (RMSE) as the error metric, which is defined to be

$\mathop{\mathrm{RMSE}}\left(\boldsymbol{C}^{\mathrm{si}}\right)=\sqrt{\frac{1}{N}\sum_{i=1}^{N}{\left({C}_{i}^{\mathrm{si}}\left(T\right)-{C}_{i}^{\mathrm{an}}\left(T\right)\right)^{2}}},$
(18)

where $\boldsymbol{C}^{\mathrm{an}}\left(T\right)$ is the analytical solution
vector at time $T$, $\boldsymbol{C}^{\mathrm{si}}\left(T\right)$ is the
corresponding result from a given simulation.

For our numerical experiments, the $N$ particles are initially distributed
over a fixed interval $\left[0,L\right]$, with $L=1$ [L]. The Dirac delta
initial condition is represented in the numerical model by placing a particle
with concentration $N/L$ at the center of the domain. We compare three
different types of simulations: equally-spaced, stationary particles (Section
3.1), randomly-spaced, stationary particles (Section 3.2), and particles
moving by Brownian motion random walks (Section 3.3). For the latter two
cases, initial particle positions are assigned according to draws from a
uniform, $\mathcal{U}(0,1)$, distribution, and ensembles of $9\,520$ and
$1\,660$ realizations of each configuration, respectively, are performed in
order to obtain a smooth estimation of the expected error by averaging over
the ensemble. For fixed values of $N$ and $\Delta t$, we define $h_{*}$ as the
bandwidth for which the lowest average $\mathrm{RMSE}$ is obtained, i.e.,

$h_{*}=\mathop{\mathrm{argmin}}_{h>0}\left(\overline{\mathop{\mathrm{RMSE}}\left(\boldsymbol{C}^{\mathrm{si}};h\right)}\right),$
(19)

where
$\overline{\mathop{\mathrm{RMSE}}\left(\boldsymbol{C}^{\mathrm{si}};h\right)}$
is the average $\mathrm{RMSE}$ over all realizations.

###  3.1 Equally-spaced, stationary particles

Figure 2 shows $\mathrm{RMSE}$ (18) as a function of $h$ for different values
of $N$ and $\Delta t$, for simulations with evenly-spaced, stationary
particles. In this case we observe a high degree of overlap between the
curves, since marginal changes in $N$ and/or $\Delta t$ do not always have a
significant effect on the simulation results. The simple explanation for this
is that, for a fixed $\Delta t$ that implies a given dispersion distance,
$\ell=\sqrt{2D\Delta t}$, increasing $N$ beyond a certain point does nothing
to improve the “resolution” of the simulation, and the reverse also holds. We
see that, given a high enough density of particles ($N$ sufficiently large),
the closest possible representation of the dispersion equation (lowest
$\mathrm{RMSE}$) occurs for $\beta=1$. In other words, for evenly-spaced
particles, the smoothing kernel associated with $\beta=1$ is virtually free of
numerical error when used in the SPH algorithm, as it in fact matches the
analytical solution of the solute’s dispersion over a time step. It is worth
noting here that this value of $\beta=1$ does not correspond to the particle
co-location algorithm, given in (11) [see 15], but to the generalization of
the Green’s function algorithm instead [16], which is discussed in Section 2.
From a physical point of view, using a kernel bandwidth larger than $\ell$
($\beta<1$), could be seen as equivalent to assuming that the solute mass
carried by each particle is Gaussian-distributed in space over some support,
rather than a Dirac delta, prior to the start of the time step [10]. This is
consistent with the fact that, for low $N$, the $\mathrm{RMSE}$ can be reduced
(up to a certain point) by using a larger kernel; i.e., the assumption that
each particle is distributed over some support can mitigate the need for more
particles. Conversely, choosing a kernel bandwidth significantly smaller than
$\ell$ $(\beta>1$), in addition to not having a clear physical meaning,
generates numerical instabilities because the mass transfer between two
particles in one time step may be larger than the difference between their
masses (see (8)). As a result, these cases are excluded from the results shown
in Figure 2.

Some of the aforementioned relations can be better observed in Figure 3. Given
a coarse time discretization (Figure 3$\left(\mathrm{a}\right)$, green curves
and markers), $h_{*}$ does not depend on $s$, and $h_{*}=\ell$. Given a finer
time discretization and a low particle density, we have the relation
$h_{*}\propto s$ (see the linear trend, for large $s$, in the yellow curves of
Figure 3$\left(\mathrm{a}\right)$). This proportionality is consistent with
the known theoretical behavior for the truncation error of the SPH
interpolation, given evenly-spaced particles [23]. In examining the relation
of $h_{*}$ to the dispersion distance $\ell=\sqrt{2D\Delta t}$ in Figure
3$\left(\mathrm{b}\right)$, we observe that, for sufficiently high values of
$N$ and $\Delta t$, we have $h_{*}=\ell$ (corresponding to $\beta=1$, see the
clearly distinguished minima in Figure 2), and otherwise, $h_{*}\simeq
s/\sqrt{2}$ (the curves with less pronounced minima in Figure 2). All these
relations are summarized by the two distinguishable regimes that can be seen
in Figure 3$\left(\mathrm{c}\right)$, wherein $h_{*}$ and $s$ are non-
dimensionalized via scaling by the dispersion distance $\ell$.

###  3.2 Randomly-spaced, stationary particles

The numerical results for randomly-distributed particles show less distinct
trends, in terms of matching the analytical solution, than those seen for the
evenly-distributed particles of Section 3.1, and this can be seen in Figure 4.
In this case, the $\mathrm{RMSE}$ does not always have such a clearly
identifiable minimum in the vicinity of $h_{*}$, nor does $h_{*}$ reliably
correspond to $\beta=1$, as we saw in Section 3.1. Rather, its behavior
appears to roughly agree with the theoretical SPH truncation error for
randomly-spaced particles [22, 23], which can be expressed as the summation of
two terms: the smoothing error, which scales with $h$; and the quadrature
error, which scales with $\langle s\rangle/h$, where $\langle s\rangle$ is the
expected particle separation (here, $\langle s\rangle=L/N$). Balancing these
two terms results in $h_{*}\propto\sqrt{\langle s\rangle}$, and hence for that
choice of bandwidth the truncation error scales with $\sqrt{\langle
s\rangle}$. This is consistent with the results shown in Figure 4, where ,
given $h=h_{*}$ (i.e., considering only each curve’s minimum), the
$\mathrm{RMSE}$ scales with the particle number as $\mathrm{RMSE}\propto
N^{-1/2}$.

It is only when $\Delta t$ adopts large values that it appears to have a
noticeable influence on the $\mathrm{RMSE}$. This behavior is also evident in
the relative insensitivity of $h_{*}$ to $\ell$, as can be seen in Figure
5$\left(\mathrm{b}\right)$. In Figure 5$\left(\mathrm{a}\right)$ we see that
the relation of $h_{*}$ to the average particle spacing $\langle s\rangle$ is
not linear, not even for small $\Delta t$, unlike in the evenly-spaced
particle case. Instead, we observe a range of slopes in the log-log space
(about $1/2$ and lower), which can be related to the aforementioned truncation
error [23], which is minimized when $h\propto\sqrt{\langle s\rangle}$. Unlike
the equally-spaced case (Figure 3$\left(\mathrm{c}\right)$), we do not observe
a single linear trend in Figure 5$\left(\mathrm{c}\right)$ for the
relationship between $h_{*}/\ell$ and $\langle s\rangle/\ell$. Rather, we
observe the general tendency that $\langle s\rangle\to 0$ implies
$h_{*}\to\ell$. For the range of tested values, a relatively high particle
density, of $\langle s\rangle\lesssim 0.01\ell$, is required to observe the
relation $h_{*}\simeq\ell$.

###  3.3 Random-walking particles

The same set of simulations are also conducted for a hybrid model in which the
dispersion coefficient is partitioned as

$D=D_{\text{RW}}+D_{\text{MT}},$ (20)

where $D_{\text{MT}}$ is the dispersion coefficient used in the SPH/MTPT
algorithm described in the previous section, and particles move by Brownian
motion, according to the Langevin equation. For a time discretization
$\\{t_{1},t_{2},...,t_{n}\\}$, with $t_{k+1}=t_{k}+\Delta t$,

$X_{i}^{k+1}=X_{i}^{k}+\xi_{i}^{k}\sqrt{2D_{\text{RW}}\Delta t},$ (21)

where $X_{i}^{k}\coloneqq X_{i}(t_{k})$, and $\xi_{i}^{k}$ is a random number
drawn from a standard normal, $\mathcal{N}(0,1)$, distribution. With an
appropriate choice of $D_{\text{RW}}$ and $D_{\text{MT}}$, this type of
approach can be used to give a separate treatment to the non-mixed spreading
(RW) and the actual mixing (MT). Several authors [24, 25, 26, 27] have
suggested that these correspond to the anisotropic spreading (longitudinal
minus transverse hydrodynamic dispersion) and the isotropic mixing (molecular
diffusion plus transverse hydrodynamic dispersion) parts of the dispersion
tensor, respectively. Here we simply set $D_{\text{RW}}=D_{\text{MT}}=D/2$.
Note that, for this partitioning, random walks do not significantly perturb
spatial concentrations about their expected value. That is, the concentration
difference between two spatially coincident particles is negligible, meaning
that the concentrations at a given time vary “smoothly” with the particle
positions $X_{i}$ (see Figure 1, yellow markers). This is because particles
exchange mass at the same rate at which they diffuse by Brownian motion. For
this reason, we can study the influence of $h$ on the numerical results when
particles are random-walking and compare to the case where particles are
stationary (as in Sections 3.1 and 3.2), without introducing the concentration
variance that would be otherwise (purposefully) induced by setting
$D_{\text{RW}}\gg D_{\text{MT}}$. Since, at $t=0$, there is only one particle
with nonzero concentration, a strong variability in the results is introduced
by the random motion of that particle in the initial stages of the simulation,
when it is carrying nearly all the solute mass in the system. For this reason,
in order to favor faster convergence of the $\mathrm{RMSE}$ with the number of
simulations, we set that singular particle to be motionless and to use the
full dispersion coefficient in its mass-transfer calculations (i.e., for that
particle, $D_{\text{MT}}=D$ and $D_{\text{RW}}=0$). An alternative approach to
overcome the same issue could be to use more particles to represent the
initial Dirac delta condition.

The behavior of the $\mathrm{RMSE}$ in this case (Figure 6) can be seen as
occupying a middle ground between the equally-spaced (Figure 2) and the
randomly-spaced (Figure 4), stationary cases. The distribution of particle
spacings in the random-walking case at any given time is identical to the
stationary randomly-distributed case, but in the former, the expected, or
time-averaged, particle spacing distribution is much narrower, approximating
the stationary, evenly-spaced case in that sense. For that reason, we do
expect the value of $h_{*}$ for a random-walking model, in the context of this
specific example, to be bounded between the two extreme stationary cases,
which may be thought of as the most ordered and disordered systems,
respectively. Note, however, that the actual values of the $\mathrm{RMSE}$ in
Figure 6 are on the same order of magnitude as for the randomly-distributed,
stationary particles (Figure 4), and they can be even higher. This may be
attributed to the added natural variability of Brownian random walks used to
represent half of the dispersion, as opposed to the deterministic nature of
mass transfers. For high enough $N$ and $\Delta t$, we can see that
$\mathrm{RMSE}$ minima occur at $h_{*}=\ell$ and are strongly pronounced.
Otherwise, we see milder minima and $h_{*}>\ell$, similarly to what is
observed for equally-spaced particles (Figure 2). In these regions of milder
minima, we see the approximate scaling $\mathrm{RMSE}\propto N^{-1/2}$ given
$h=h_{*}$, which, in this behavior, is similar to the randomly-spaced,
stationary case (Figure 4).

We see that for a fine time discretization (blue line in Figure
7$\left(\mathrm{a}\right)$), we have $h_{*}\propto\sqrt{\langle s\rangle}$,
which, as mentioned in Section 3.2, indicates that $h_{*}$ in these regimes is
mainly controlled by the truncation error of the spatial interpolation. On the
other hand, we see a clear trend that $h_{*}=\ell$ for large enough $N$ and
$\Delta t$, as evidenced by the triangle symbols and green markers in Figure
7$\left(\mathrm{b}\right)$. As in the previous cases, $h_{*}$ departs from
$\ell$ at some threshold as the relative spacing $\langle s\rangle/\ell$
increases. Like in the stationary, randomly-spaced case, and unlike the
equally-spaced case, this threshold value for $\langle s\rangle/\ell$ appears
to depend on $\ell$ (i.e., no single linear trend is observed in Figure
7$\left(\mathrm{c}\right)$, unlike in Figure 3$\left(\mathrm{c}\right)$).
Nevertheless, for the range of tested values, $h_{*}\simeq\ell$ for $\langle
s\rangle\lesssim 0.1\ell$.

##  4 Summary and discussion

In this paper, we demonstrate an equivalence between the Lagrangian SPH
(smoothed particle hydrodynamics) and MTPT (mass transfer particle tracking)
methods for simulating dispersion, provided that the spatial kernel being
employed is Gaussian. These two methods originate from completely different
interpretations. The SPH community views their methods (classically speaking,
as recent work has included random walks in SPH simulations [18]) as solving
the dispersion equation by projecting the particles onto the continuum using
radial basis functions (kernels) and approximating the solution on that kernel
space. The random walk particle tracking community views the MTPT methods
considered in this paper in two ways: (i) a first-principles approach, wherein
mass-transfers between moving particles are scaled by the probability that
these particles co-locate via dispersion; (ii) a discretization of the Green’s
function for the dispersion equation, in which a particle’s solute mass is
spread in space via mass-transfers to its nearest neighbors. Previously, these
two MTPT methods were considered to be distinct approaches, and neither had
rigorous proofs associated with it. As a result of this work, however, both of
these MTPT methods now inherit a rigorous theoretical underpinning from the
SPH literature.

The numerical investigations we conduct yield compelling results regarding the
proper Gaussian kernel bandwidth for particle tracking simulations. We see
strong evidence that a kernel with bandwidth $h=\ell=\sqrt{2D\Delta t}$ (i.e.,
imposing $\beta=1$) is the ideal choice, provided there is a “dense enough”
spatial distribution of particles. This makes intuitive/physical sense
because, with bandwidth $\ell$, this Gaussian function is the fundamental
solution of the dispersion equation. In other words, aside from the error
introduced in the normalization step, using this kernel for mass transfer is
not an approximation, but rather a semi-analytical solution of the dispersion
in a time-step of length $\Delta t$. We also observe that, counter-
intuitively, a coarser time-discretization may be a better choice than a finer
one, if that allows one to use bandwidth $\ell$. However, there may be cases
in which the intent is to reproduce the dispersion equation without the
distortion associated with a low particle density (a subject that we discuss
below), but a high particle density cannot be afforded, computationally (as
may be likely to occur in multi-dimensional systems). If, in these cases, the
use of a long time-step would generate other forms of error (for instance, in
the chemical reactions), then a wider kernel bandwidth than $\ell$ (following
the traditional SPH bandwidth selection rules-of-thumb) may be a better choice
when seeking a compromise between accuracy and efficiency. One way to think of
this is to consider the wider-bandwidth particle to be a “macro-particle,” or
cluster of smaller particles, that is distributed in space over some support
volume.

Additional conclusions can be drawn from each of the individual cases tested
in Section 3. In the equally-spaced, stationary particle case, $h=\ell$ is
clearly the optimal bandwidth choice, provided that $N$ is sufficiently large,
as to capture the magnitude of dispersion, described by $\ell=\sqrt{2D\Delta
t}$ (i.e., particles must be close enough to “see” one another).

Considering the randomly-distributed, stationary particle case, we see a
different story, in that $\mathrm{RMSE}$ tends to be more related to average
inter-particle spacing, $\langle s\rangle=L/N$, than it is to the dispersion
distance, $\ell$. This is most likely because, for the range of $N$ and
$\Delta t$ values tested, the $\mathrm{RMSE}$ is dominated by the truncation
error of the SPH interpolation. Nevertheless, according to some authors in
particle methods [e.g., 3, 8], the distortion of the numerical solution caused
by heterogeneity in the inter-particle spacing and low particle densities can
represent incomplete mixing conditions, rather than being just a numerical
error. If we subscribe to this view, then the randomly-spaced case represents
areas in which particles are poorly-mixed and remain poorly mixed for the
duration of the simulation. From that perspective, using the $\ell$ bandwidth
would only be capturing the “average mixedness” of such a simulation, fully
simulating diffusive mixing in well-mixed areas and under-simulating mixing in
poorly mixed areas. In light of this, the increase in $\mathrm{RMSE}$ could be
thought of not as an error, but as desirable deviations from the well-mixed
solution, due to physically meaningful areas of poor mixing.

For the case of random-walking particles, we find that the qualitative
behavior of the $\mathrm{RMSE}$ with respect to the bandwidth $h$ can be
placed in a middle ground between the other two scenarios. In fact, the minima
($h_{*}$) are found to be bounded in this case between the two former cases.
It is clear from the results that, despite the particle disorder, the
dependence of the $\mathrm{RMSE}$ on $h$ should not be understood as a
function of the particle density alone. Instead, the error originated in
deviating from the dispersion kernel bandwidth $h=\ell$ should also be
considered. Again, if the effects of particle disorder on the numerical
solution are considered to be physically meaningful, it makes sense that
random-walking particles are closer to representing a well-mixed system
(distinguishable by $h_{*}=\ell$) than stationary randomly-distributed
particles, since in this case the poorly-mixed areas are not persistent in
time.

We believe the results of our numerical experiments are relevant in a general
sense, despite representing the specific simple case of a Dirac delta initial
condition in a one-dimensional setting. This particular dispersion problem,
where one initial concentration pulse spreads by dispersion, is no doubt the
simplest one; however, any more complex problem can be thought of as unions of
Dirac delta initial conditions, at least from a computational/discrete
standpoint. As long as the physics are being captured on a local, particle
level, as is demonstrated here, more complicated conditions will also be
properly simulated. Additionally, we expect the scaling with $s$ and $\ell$ to
be analogous for isotropic dispersion in higher dimensions because mass
transfers are merely a function of Euclidean distance between particles, and
hence not substantively different in higher spatial dimensions. However, the
scaling relations will likely need to be reformulated in terms of fill
distance, rather than the simple inter-particle spacing we see here in 1D.
Besides, the analysis performed in Section 3 would undoubtedly become more
complex in the case of anisotropic and spatially variable dispersion.

The traditional SPH extension to anisotropic dispersion entails a more
complicated expression for $\widehat{D}_{ij}$ in (1), while maintaining the
isotropy of the kernel $W$, and this approach may result in negative
concentrations [14]. This is in contrast to the more straightforward extension
of traditional MTPT to anisotropic dispersion, which would involve redefining
$W$ as an anisotropic multi-Gaussian with variance $2\Delta t\cdot
g\left(\boldsymbol{D}\left(\boldsymbol{X}_{i}\right),\boldsymbol{D}\left(\boldsymbol{X}_{j}\right)\right)/\beta$,
where $g$ is some averaging function. The subject of anisotropy is out of the
scope of this paper and should be addressed in future work. Nevertheless, as
mentioned in Section 3, another suitable approach to reproducing anisotropic
dispersion would be to split the dispersion tensor between an isotropic and an
anisotropic part, using the isotropic SPH/MTPT method addressed here to
simulate the former and reproducing the latter with random walks.

Open questions do remain in this area. For instance, we only consider the
Gaussian kernel in our analysis and results. Other kernels are commonly used
in the SPH literature, and compactly-supported kernels are known to result in
computational speedup. A standard choice is the compactly-supported Wendland
kernel that has been shown to approach a Gaussian in the infinitely-smooth,
limiting case [28]. How much error is introduced by this approximation, and
how does this compare to the common practice or imposing a cutoff distance of
$3h$ for mass transfers, as is commonly done in the particle tracking
literature?

The hybridization of SPH/MTPT with random walks is a very recent technique
that, to date, has not been studied in depth. In this work, we compare the
numerical results from one such model with an analytical solution in the
particular case wherein the simulation of the full dispersion tensor is
partitioned equally between random walks and mass transfers. If the purpose of
this hybridization is to simulate a two-scale system (as in [18]) in which the
random walk accounts for spreading and the mass transfer accounts for mixing,
it would be proper for the magnitude of mixing to be much smaller than that of
spreading, in order to generate states of local disequilibrium (as, for
instance, to simulate the effect of local heterogeneities in porous media).
Hence, further investigation is needed in this area, in order to: (i) analyze
the effect of using different spreading/mixing ratios, and (ii) evaluate the
capability of this kind of model to correctly reproduce the generation,
propagation, and decay of sub-scale concentration variance.

## References

  * [1] D. A. Benson, M. M. Meerschaert, Simulation of chemical reaction via particle tracking: Diffusion-limited versus thermodynamic rate-limited regimes, Water Resources Research 44 (12).  arXiv:https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2008WR007111, doi:10.1029/2008WR007111.
URL https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2008WR007111

  * [2] A. Paster, D. Bolster, D. A. Benson, Particle tracking and the diffusion-reaction equation, Water Resour. Res. 49 (2013) 1–6.  doi:10.1029/2012WR012444.
  * [3] A. Paster, D. Bolster, D. A. Benson, Connecting the dots: Semi-analytical and random walk numerical solutions of the diffusion–reaction equation with stochastic initial conditions, Journal of Computational Physics 263 (2014) 91 – 112.  doi:https://doi.org/10.1016/j.jcp.2014.01.020.
URL http://www.sciencedirect.com/science/article/pii/S0021999114000473

  * [4] D. A. Benson, T. Aquino, D. Bolster, N. Engdahl, C. V. Henri, D. Fernàndez-Garcia, A comparison of Eulerian and Lagrangian transport and non-linear reaction algorithms, Advances in Water Resources 99 (2017) 15 – 37.  doi:http://doi.org/10.1016/j.advwatres.2016.11.003.
  * [5] D. Bolster, A. Paster, D. A. Benson, A particle number conserving Lagrangian method for mixing-driven reactive transport, Water Resources Research 52 (2) (2016) 1518–1527.  doi:10.1002/2015WR018310.
  * [6] D. Ding, D. Benson, A. Paster, D. Bolster, Modeling bimolecular reactions and transport in porous media via particle tracking, Advances in Water Resources 53 (2012) 56–65.  doi:10.1016/j.advwatres.2012.11.001.
  * [7] D. Ding, D. A. Benson, Simulating biodegradation under mixing-limited conditions using Michaelis–Menten (Monod) kinetic expressions in a particle tracking model, Advances in Water Resources 76 (2015) 109 – 119.  doi:http://dx.doi.org/10.1016/j.advwatres.2014.12.007.
URL http://www.sciencedirect.com/science/article/pii/S0309170814002462

  * [8] D. Ding, D. A. Benson, D. Fernández-Garcia, C. V. Henri, D. W. Hyndman, M. S. Phanikumar, D. Bolster, Elimination of the reaction rate “scale effect”: Application of the Lagrangian reactive particle-tracking method to simulate mixing-limited, field-scale biodegradation at the Schoolcraft (MI, USA) site, Water Resources Research (doi:10.1002/2017WR021103).  doi:10.1002/2017WR021103.
URL http://dx.doi.org/10.1002/2017WR021103

  * [9] D. Bolster, A. Paster, D. A. Benson, A particle number conserving Lagrangian method for mixing-driven reactive transport, Water Resources Research 52 (2) (2016) 1518–1527.  doi:10.1002/2015WR018310.
URL http://dx.doi.org/10.1002/2015WR018310

  * [10] M. J. Schmidt, S. Pankavich, D. A. Benson, A kernel-based Lagrangian method for imperfectly-mixed chemical reactions, Journal of Computational Physics 336 (2017) 288 – 307.  doi:https://doi.org/10.1016/j.jcp.2017.02.012.
URL http://www.sciencedirect.com/science/article/pii/S0021999117301055

  * [11] G. Sole-Mari, D. Fernàndez-Garcia, P. Rodríguez-Escales, X. Sanchez-Vila, A KDE-based random walk method for modeling reactive transport with complex kinetics in porous media, Water Resources Research 53 (11) (2017) 9019–9039.
  * [12] G. Sole-Mari, D. Fernàndez-Garcia, Lagrangian modeling of reactive transport in heterogeneous porous media with an automatic locally adaptive particle support volume, Water Resources Research 54 (10) 8309–8331.  doi:10.1029/2018WR023033.
URL https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2018WR023033

  * [13] P. A. Herrera, M. Massabó, R. D. Beckie, A meshless method to simulate solute transport in heterogeneous porous media, Advances in Water Resources 32 (3) (2009) 413 – 429.  doi:https://doi.org/10.1016/j.advwatres.2008.12.005.
URL http://www.sciencedirect.com/science/article/pii/S0309170808002273

  * [14] P. A. Herrera, R. D. Beckie, An assessment of particle methods for approximating anisotropic dispersion, International Journal for Numerical Methods in Fluids 71 (5) 634–651.  arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1002/fld.3676, doi:10.1002/fld.3676.
URL https://onlinelibrary.wiley.com/doi/abs/10.1002/fld.3676

  * [15] D. A. Benson, D. Bolster, Arbitrarily complex chemical reactions on particles, Water Resources Research 52 (11) (2016) 9190–9200.  doi:10.1002/2016WR019368.
URL http://dx.doi.org/10.1002/2016WR019368

  * [16] M. J. Schmidt, S. D. Pankavich, D. A. Benson, On the accuracy of simulating mixing by random-walk particle-based mass-transfer algorithms, Advances in Water Resources 117 (2018) 115 – 119.  doi:https://doi.org/10.1016/j.advwatres.2018.05.003.
URL http://www.sciencedirect.com/science/article/pii/S0309170818301830

  * [17] N. B. Engdahl, D. A. Benson, D. Bolster, Lagrangian simulation of mixing and reactions in complex geochemical systems, Water Resources Research 53 (4) (2017) 3513–3522.  doi:10.1002/2017WR020362.
URL http://dx.doi.org/10.1002/2017WR020362

  * [18] P. A. Herrera, J. M. Cortínez, A. J. Valocchi, Lagrangian scheme to model subgrid-scale mixing and spreading in heterogeneous porous media, Water Resources Research 53 (4) 3302–3318.  arXiv:https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1002/2016WR019994, doi:10.1002/2016WR019994.
URL https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1002/2016WR019994

  * [19] R. A. Gingold, J. J. Monaghan, Smoothed particle hydrodynamics: theory and application to non-spherical stars, Monthly Notices of the Royal Astronomical Society 181 (3) (1977) 375–389.  doi:10.1093/mnras/181.3.375.
URL http://dx.doi.org/10.1093/mnras/181.3.375

  * [20] J. Monaghan, Smoothed particle hydrodynamics and its diverse applications, Annual Review of Fluid Mechanics 44 (1) (2012) 323–346.  doi:10.1146/annurev-fluid-120710-101220.
URL https://doi.org/10.1146/annurev-fluid-120710-101220

  * [21] M. J. Schmidt, S. D. Pankavich, A. Navarre-Sitchler, D. A. Benson, A Lagrangian method for reactive transport with solid/aqueous chemical phase interaction, Accepted, Journal of Computational Physics: X. .doi:https://arxiv.org/abs/1805.06072.
  * [22] A. M. Tartakovsky, N. Trask, K. Pan, B. Jones, W. Pan, J. R. Williams, Smoothed particle hydrodynamics and its applications for multiphase flow and reactive transport in porous media, Computational Geosciences 20 (4) (2016) 807–834.  doi:10.1007/s10596-015-9468-9.
URL https://doi.org/10.1007/s10596-015-9468-9

  * [23] N. J. Quinlan, M. Basa, M. Lastiwka, Truncation error in mesh-free particle methods, International Journal for Numerical Methods in Engineering 66 (13) (2006) 2064–2085.  arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1002/nme.1617, doi:10.1002/nme.1617.
URL https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.1617

  * [24] L. W. Gelhar, A. L. Gutjahr, R. L. Naff, Stochastic analysis of macrodispersion in a stratified aquifer, Water Resources Research 15 (6) (1979) 1387–1397.  doi:10.1029/WR015i006p01387.
URL http://dx.doi.org/10.1029/WR015i006p01387

  * [25] L. W. Gelhar, C. L. Axness, Three-dimensional stochastic analysis of macrodispersion in aquifers, Water Resources Research 19 (1) (1983) 161–180.  doi:10.1029/WR019i001p00161.
URL http://dx.doi.org/10.1029/WR019i001p00161

  * [26] O. A. Cirpka, E. O. Frind, R. Helmig, Numerical simulation of biodegradation controlled by transverse mixing, Journal of Contaminant Hydrology 40 (2) (1999) 159 – 182.  doi:https://doi.org/10.1016/S0169-7722(99)00044-3.
URL http://www.sciencedirect.com/science/article/pii/S0169772299000443

  * [27] C. J. Werth, O. A. Cirpka, P. Grathwohl, Enhanced mixing and reaction through flow focusing in heterogeneous porous media, Water Resources Research 42 (12) (2006) W12414, w12414.  doi:10.1029/2005WR004511.
URL http://dx.doi.org/10.1029/2005WR004511

  * [28] A. Chernih, I. H. Sloan, R. S. Womersley, Wendland functions with increasing smoothness converge to a Gaussian, Advances in Computational Mathematics 40 (1) (2014) 185–200.  doi:10.1007/s10444-013-9304-5.
URL https://doi.org/10.1007/s10444-013-9304-5