Numerical Equivalence Between SPH and Probabilistic Mass Transfer Methods for Lagrangian Simulation of Dispersion

TL;DR

This paper demonstrates that mass transfer particle tracking methods can be mathematically unified with smoothed particle hydrodynamics by choosing an appropriate Gaussian kernel, leading to exact solutions for dispersion in certain conditions.

Contribution

It unifies MTPT and SPH methodologies, showing their equivalence with a specific kernel choice, and validates this with numerical simulations matching analytical solutions.

Findings

Kernel bandwidth of size $\ell=\sqrt{2D\Delta t}$ yields exact dispersion displacement.

Mass transfer operation becomes exact under dense particle distribution with the optimal kernel.

Numerical results agree with analytical solutions, confirming the theoretical unification.

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.

Figures

Equations

$$\frac{d{C}_i}{dt}=2\sum _{j=1}^N\frac{{\hat{D}}_{ij}}{{\hat{\rho }}_{ij}}\left(C_i-C_j\right)F\left({\mathbf{X}}_i-{\mathbf{X}}_j;h\right).$$

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

$${\hat{D}}_{ij}:=g\left(D({\mathbf{X}}_i),D({\mathbf{X}}_j)\right),$$

$${\hat{\rho }}_{ij}$$

$$\rho \left(\mathbf{X};h\right)$$

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

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

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

$${\beta }_{ij}\left(h\right)$$

$${\mathcal{W}}_{ij}\left(h\right)$$

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

$${\mathcal{P}}_{ij}=\frac{P({\mathbf{X}}_i-{\mathbf{X}}_j)}{{\tilde{\rho }}_{ij}}.$$

$$P({\mathbf{X}}_i-{\mathbf{X}}_j;D,\Delta t)$$

$$\equiv W\left({\mathbf{X}}_i-{\mathbf{X}}_j;\sqrt{2(D_i+D_j)\Delta t}\right)$$

$$=W\left({\mathbf{X}}_i-{\mathbf{X}}_j;\sqrt{4{\hat{D}}_{ij}\Delta t}\right)$$

$$=W\left({\mathbf{X}}_i-{\mathbf{X}}_j;\sqrt{2}{\ell }_{ij}\right),$$

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

$$\mathbf{A}:=\mathbf{I}+\left[\mathbf{\beta }\circ \mathcal{W}-\mathrm{diag}\left(\left[\mathbf{\beta }\circ \mathcal{W}\right]1\right)\right].$$

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

$${\mathcal{G}}_{ij}$$

$$\mathrm{RMSE}\left({\mathbf{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},$$

$$h_{\ast }={\mathrm{argmin}}_{h>0}\left(\overline{\mathrm{RMSE}\left({\mathbf{C}}^{\mathrm{si}};h\right)}\right),$$

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

$$X_i^{k+1}=X_i^k+{\xi }_i^k\sqrt{2D_{\text{RW}}\Delta t},$$

Full text

††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:

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

in which we define

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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:

[TABLE]

with

[TABLE]

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

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

[TABLE]

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.,

[TABLE]

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

[TABLE]

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$,

[TABLE]

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.