Particle Density Estimation with Grid-Projected Adaptive Kernels

TL;DR

This paper introduces a grid-projected adaptive kernel density estimation method that combines binning and KDE for efficient, boundary-aware particle density estimation in fluid transport simulations.

Contribution

It develops a hybrid KDE approach with local bandwidth optimization, improving efficiency and boundary handling over traditional KDE methods.

Findings

Enhanced computational efficiency demonstrated in examples

Effective boundary condition handling achieved

Adaptive bandwidth improves density estimation accuracy

Abstract

The reconstruction of smooth density fields from scattered data points is a procedure that has multiple applications in a variety of disciplines, including Lagrangian (particle-based) models of solute transport in fluids. In random walk particle tracking (RWPT) simulations, particle density is directly linked to solute concentrations, which is normally the main variable of interest, not just for visualization and post-processing of the results, but also for the computation of non-linear processes, such as chemical reactions. Previous works have shown the superiority of kernel density estimation (KDE) over other methods such as binning, in terms of its ability to accurately estimate the "true" particle density relying on a limited amount of information. Here, we develop a grid-projected KDE methodology to determine particle densities by applying kernel smoothing on a pilot binning; this may be seen as a "hybrid" approach between binning and KDE. The kernel bandwidth is optimized locally. Through simple implementation examples, we elucidate several appealing aspects of the proposed approach, including its computational efficiency and the possibility to account for typical boundary conditions, which would otherwise be cumbersome in conventional KDE.

Equations

$$c\left(\mathbf{x}\right)=\frac{m\rho \left(\mathbf{x}\right)}{\varphi \left(\mathbf{x}\right)},$$

$$\rho \left(\mathbf{x}\right):=\sum _{p=1}^{N}W\left(\mathbf{x}-{\mathbf{X}}_p;{\mathbf{h}}_p\right),$$

$$W\left(\mathbf{r};\mathbf{h}\right)=\prod _{i=1}^d\frac{1}{\sqrt{2\pi }{h}_i}\exp \left(-\frac{r_i^2}{2h_i^2}\right),$$

$${\rho }_u:=\frac{1}{\Lambda }{\int }_{{\mathbf{x}}_u-\text{uplambda}/2}^{{\mathbf{x}}_u+\text{uplambda}/2}\rho \left(\mathbf{x}\right)\mathrm{d}\mathbf{x}\approx \frac{1}{\Lambda }\sum _{\omega =1}^{\nu }{\mu }_{\omega }\overline{W}\left({\mathbf{x}}_{\omega }-{\mathbf{x}}_u;{\mathbf{h}}_{\omega },\text{uplambda}\right),$$

$$\overline{W}\left(\mathbf{r};\mathbf{h},\text{uplambda}\right):={\int }_{\mathbf{r}-\text{uplambda}/2}^{\mathbf{r}+\text{uplambda}/2}W\left({\mathbf{r}}^{\prime };\mathbf{h}\right)\mathrm{d}{\mathbf{r}}^{\prime }.$$

$$\overline{W}\left(\mathbf{r};\mathbf{h},\text{uplambda}\right)=\frac{1}{2^d}\prod _{i=1}^d\left[\mathrm{erf}\left(\frac{r_i+{\lambda }_i/2}{\sqrt{2}{h}_i}\right)-\mathrm{erf}\left(\frac{r_i-{\lambda }_i/2}{\sqrt{2}{h}_i}\right)\right].$$

$$\mathbf{r}=\text{uplambda}\odot \mathbf{z},$$

$$\overline{W}\left(\text{uplambda}\odot \mathbf{z};\mathbf{h},\text{uplambda}\right)=\frac{1}{2^d}\prod _{i=1}^d\left[\mathrm{erf}\left({\lambda }_i\frac{z_i+1/2}{\sqrt{2}{h}_i}\right)-\mathrm{erf}\left({\lambda }_i\frac{z_i-1/2}{\sqrt{2}{h}_i}\right)\right].$$

$$c_u:=\frac{1}{\Lambda }{\int }_{{\mathbf{x}}_u-\text{uplambda}/2}^{{\mathbf{x}}_u+\text{uplambda}/2}c\left(\mathbf{x}\right)\mathrm{d}\mathbf{x}\approx \frac{m{\rho }_u}{{\varphi }_u},$$

$$\hat{h}:={\left(\prod _{i=1}^d{h}_i\right)}^{\frac{1}{d}},s_i:=h_i/\hat{h},$$

$${\hat{h}}_u={\left[\frac{d{n}_u}{{\left(4\pi \right)}^{\frac{d}{2}}{T}_u}\right]}^{\frac{1}{d+4}},$$

$$n_u:=\int \rho \left(\mathbf{x}\right)W\left(\mathbf{x}-{\mathbf{x}}_u;{\text{upsigma}}_u\right)\mathrm{d}\mathbf{x}\approx \sum _{\omega =1}^{\nu }{\rho }_{\omega }\overline{W}\left({\mathbf{x}}_{\omega }-{\mathbf{x}}_u;{\text{upsigma}}_u,\text{uplambda}\right),$$

$$T_u=\{\begin{array}{ll}{\Psi }_{11,u},& \mathrm{for}d=1,\\ 2{\left({\Psi }_{11,u}{\Psi }_{22,u}\right)}^{\frac{1}{2}}+2{\Psi }_{12,u},& \mathrm{for}d=2,\\ 3{\left({\Psi }_{11,u}{\Psi }_{22,u}{\Psi }_{33,u}\right)}^{\frac{1}{3}}+{\sum }_{i\ne j\ne k}{\Psi }_{ij,u}{\left(\frac{{\Psi }_{kk,u}^2}{{\Psi }_{ii,u}{\Psi }_{jj,u}}\right)}^{\frac{1}{6}},& \mathrm{for}d=3.\end{array}$$

$$\begin{array}{rl}{\Psi }_{ij,u}:=\int {\kappa }^{(i)}\left(\mathbf{x}\right){\kappa }^{(j)}\left(\mathbf{x}\right)W& \left(\mathbf{x}-{\mathbf{x}}_u;{\text{upsigma}}_u\right)\mathrm{d}\mathbf{x}\\ & \approx \sum _{\omega =1}^{\nu }{\kappa }_{\omega }^{(i)}{\kappa }_{\omega }^{(j)}\overline{W}\left({\mathbf{x}}_{\omega }-{\mathbf{x}}_u;{\text{upsigma}}_u,\text{uplambda}\right),\end{array}$$

$${\kappa }^{(i)}\left(\mathbf{x}\right):=\frac{{\partial }^2\rho }{\partial x_i^2}\left(\mathbf{x}\right)=\sum _{p=1}^N{V}^{(i)}\left(\mathbf{x}-{\mathbf{X}}_p;{\mathbf{g}}_p^{\left(i\right)}\right),$$

$$V^{(i)}\left(\mathbf{r};\mathbf{g}\right):=\frac{{\partial }^{2}W}{\partial r_i^2}=\left(\frac{r_i^2}{g_i^4}-\frac{1}{g_i^2}\right)W\left(\mathbf{r};\mathbf{g}\right).$$

$${\kappa }_u^{(i)}:=\frac{1}{\Lambda }{\int }_{{\mathbf{x}}_u-\text{uplambda}/2}^{{\mathbf{x}}_u+\text{uplambda}/2}{\kappa }^{(i)}\left(\mathbf{x}\right)\mathrm{d}\mathbf{x}\approx \frac{1}{\Lambda }\sum _{\omega =1}^{\nu }{\mu }_{\omega }{\overline{V}}^{(i)}\left({\mathbf{x}}_{\omega }-{\mathbf{x}}_u;{\mathbf{g}}_{\omega }^{(i)},\text{uplambda}\right),$$

$${\overline{V}}^{(i)}\left(\mathbf{r};\mathbf{g},\text{uplambda}\right):={\int }_{\mathbf{r}-\text{uplambda}/2}^{\mathbf{r}+\text{uplambda}/2}{V}^{(i)}\left({\mathbf{r}}^{\prime };\mathbf{g}\right)\mathrm{d}{\mathbf{r}}^{\prime }.$$

$${\overline{V}}^{(i)}\left(\mathbf{r};\mathbf{g},\text{uplambda}\right)=-\frac{1}{2^{(d-\frac{1}{2})}\sqrt{\pi }{g}_i^3}\times$$

$$\left[\left(r_i+\frac{{\lambda }_i}{2}\right)\exp \left(-\frac{{\left(r_i+{\lambda }_i/2\right)}^2}{2g_i^2}\right)-\left(r_i-\frac{{\lambda }_i}{2}\right)\exp \left(-\frac{{\left(r_i-{\lambda }_i/2\right)}^2}{2g_i^2}\right)\right]\times$$

$$\prod _{j\ne i}\left[\mathrm{erf}\left(\frac{r_j+{\lambda }_j/2}{\sqrt{2}{g}_j}\right)-\mathrm{erf}\left(\frac{r_j-{\lambda }_j/2}{\sqrt{2}{g}_j}\right)\right].$$

$$s_{i,u}={\left(\frac{{\hat{\Psi }}_u}{{\Psi }_{ii,u}}\right)}^{\frac{1}{4}},{\hat{\Psi }}_u:={\left(\prod _{j=1}^d{\Psi }_{jj,u}\right)}^{\frac{1}{d}}.$$

$${\hat{h}}_u\equiv {\hat{h}}_u^{\sigma }={\left[\frac{4}{\left(d+2\right)N_u^{\sigma }}\right]}^{\frac{1}{d+4}}{\hat{\sigma }}_u.$$

$$N_u^{\sigma }=\frac{{\left(\sqrt{8\pi }{\hat{\sigma }}_u\right)}^d{n}_u^2}{{\rho }_u}.$$

$${\hat{\sigma }}_u={\left[\frac{\left(d+2\right){\left(8\pi \right)}^{\frac{d}{2}}{n}_u^2{\hat{h}}_u^{d+4}}{4{\rho }_u}\right]}^{\frac{1}{4}}.$$

$$g_j^{(i)}\equiv {\hat{g}}^{(i)},\forall j=1,\dots ,d.$$

$${\hat{g}}_u^{(i)}\equiv {\hat{g}}_u^{(i),\sigma }={\left[\frac{4+2^{\frac{d}{2}+4}}{3\left(d+4\right)N_u^{\sigma }}\right]}^{\frac{1}{d+6}}{\hat{\sigma }}_u{\vartheta }_i\left({\mathbf{s}}_u\right),$$

$${\vartheta }_i\left(\mathbf{s}\right)={\left[\frac{1}{d+4}\sum _{j=1}^d\frac{1+4{\delta }_{ij}}{s_i^4{s}_j^2}\right]}^{-\frac{1}{d+6}},$$

$${\gamma }_u^{(i)}:=\frac{{\hat{g}}_u^{(i)}}{{\hat{h}}_u}=\alpha \cdot {\left(N_u^{\sigma }\right)}^{\beta }\cdot {\vartheta }_i\left({\mathbf{s}}_u\right),$$

$$\alpha :={\left[\frac{1+2^{\frac{d+4}{2}}}{3\cdot 2^{\frac{4}{d+4}}}\right]}^{\frac{1}{d+6}}\cdot \frac{{\left(d+2\right)}^{\frac{1}{d+4}}}{{\left(d+4\right)}^{\frac{1}{d+6}}},\beta :=\frac{2}{\left(d+4\right)\left(d+6\right)},$$

$${\mathbf{n}}^{\mathrm{T}}\left[{\left(\varphi \mathbf{D}\nabla c\right)}^{+}-{\left(\varphi \mathbf{D}\nabla c\right)}^{-}\right]={\mathbf{n}}^{\mathrm{T}}\left[{\left(\mathbf{q}c\right)}^{+}-{\left(\mathbf{q}c\right)}^{-}\right],$$

Full text

††journal: Advances in Water Resources

1 Introduction

Random Walk Particle Tracking (RWPT) methods are a family of methods commonly used in the hydrologic sciences to simulate transport. They are appealing as they can accurately emulate many different physical processes that occur in porous media such as diffusion, hydrodynamic dispersion, mass transfer across multiple porosity systems and linear sorption [1, 2]. They are also conducive to simulating anomalous non-Fickian transport that arises due to medium heterogeneities below the scale of resolution [e.g. 3]. With RWPTs, the solute mass is discretized into a large number of discrete particles that move across the porous medium following deterministic and probabilistic rules, which account for the processes of advection, dispersion, matrix diffusion, etc. Lagrangian methods for simulating scalar transport, among which RWPTs are some of the most common, have been shown to be particularly useful when modeling transport in advection-dominated systems, where Eulerian methods can suffer from numerical dispersion and instabilities [1, 4].

However, the main shortcoming of RWPT methods is that, without modification, they can result in very noisy concentration fields due to subsampling effects associated with the finite number of particles in the system. This can be particularly troublesome when simulating solute transport in systems where processes are governed by nonlinearities or tight coupling and interactions between different solute concentrations, of which nonlinear chemical reactions are a prime example. Linear processes such as simple degradation, slow sorption or chain reactions can efficiently be incorporated to the RWPT algorithm by means of additional probabilistic rules with little additional computational cost. On the other hand, nonlinear reactions involve interactions between neighboring particles, which adds complexity to the problem as particles do not just need to know where they are, but where all others are also, which in naive implementations would mean an $\mathcal{O}(N^{2})$ numerical cost. Even with more optimized approaches that use better search algorithms [e.g. 5], the additional numerical cost can become significant for high particle numbers.

A problem that clearly highlights these issues and has received considerable recent attention in the literature is the simulation of bimolecular reactions of the type $\mathrm{A}+\mathrm{B}\rightarrow\mathrm{C}$ via RWPT [6, 7, 8, 9]. In many such cases it has been shown that the noise associated with the particles can fundamentally change the large scale behavior of the system; in some instances, it may reflect a true noise in the system [10, 11], but in others it may be a numerical artifact that leads to incorrect predictions [12]. Most attempts to simulate other, more complex reactions from a Lagrangian perspective have chosen to attribute concentrations to particles instead of fixed masses [13, 14], and to represent non-advective processes by means of mass transfer, thus allowing arbitrarily complex reactive processes to be simulated on-particle. This approach, however, hinders some of the intrinsic advantages of RWPT. For instance, species-dependent transport properties cannot be readily implemented, at least not in a quick and straightforward manner. Moreover, since solute mass tends to occupy new fluid as the simulation advances, empty particles need to be included in all those areas where it is anticipated that the solute may reach by dispersion [15].

A first attempt at simulating arbitrarily complex kinetic chemical reactions with traditional RWPT methods was recently presented by [16]. The method smooths concentration fields using an optimal kernel density estimator. In the paper, the authors derive an expression for the probability of reaction of a particle for any kinetic rate expression, using the optimal kernel density estimator as the particle support volume. This was extended in [17], where using a locally adaptive optimal kernel to determine the solute concentrations and ultimately the particle reaction probabilities were demonstrated to be superior over methods based on simple binning, where concentrations are computed by defining fixed representative volumes or bins over which the contained particle mass is assumed to be uniformly distributed. Despite potential shortcomings, the definition of a spatial discretization as a set of fixed bins is very convenient; it allows establishing links between the Lagrangian representation of the solute mass, and other properties that may be defined in space, thus readily enabling nonlinear and coupled interactions. Moreover, particle counting in a regular mesh can be very efficient from a computational viewpoint, typically much more so than a kernel density approach. Thus, one has to weigh disadvantages and benefits when deciding which approach to take. Ideally, one would have a method designed to obtain the best of all worlds.

Another challenge that often arises in the use of RWPT is the application of nontrivial boundary conditions, best defined in an Eulerian framework. Several methods have been proposed in the recent literature to incorporate different kinds of boundary conditions to RWPT, allowing us to simulate impermeable, Dirichlet [18], or Robin [19, 20, 21] boundary conditions. However, current kernel methods for the reconstruction of the concentrations [e.g., 22] do not address the subject of boundary conditions.

In this paper we propose a robust methodology to obtain solute concentrations from particle positions in RWPT simulations. These concentration fields can then be used to visualize the results, and perhaps more importantly, to incorporate nonlinear fate and transport processes. Although we focus on advection-dispersion in relation to a porous medium, this methodology is broadly applicable to other similar transport processes and systems. In fact, any particle-based approach that relies on some form of particle density estimation could benefit from it. Our proposed approach can be seen as a hybrid technique, where the nonreactive part of the transport (advection-dispersion) is simulated following classical principles of RWPT, and the reactive part is decoupled from the former and assisted by a grid on which the concentrations are estimated following an improved form of the local optimal kernel density estimation technique introduced in [17]. This new on-grid kernel-based method combines the practical efficiency of binning techniques with the accuracy gains of kernel methods, while also allowing us to impose Neumann, Dirichlet and Robin boundary conditions to the density estimation.

The paper is structured as follows. First, in $\S$2, we describe the methodology by which we can adapt the kernel density estimation procedure described in [17] to the case of a gridded domain. Then, in $\S$3, we explain how to correct the resulting concentration estimations to account for a variety of boundary conditions typically considered in porous media. Then, in $\S$4 we perform several computational experiments to test and also illustrate the proposed methodology, and in $\S$5 we finish with a summary and conclusions.

2 On-Grid Concentration Estimation From Particle Positions With Local Optimal Kernels

Let us consider a numerical particle cloud, made up of $N$ particles of identical mass $m$, that represent the spatial distribution of the total mass of a given chemical compound (solute), over a $d$-dimensional continuum. Particle positions are given by $\left\{\mathbf{X}_{1},\dots,\mathbf{X}_{N}\right\}$. The concentration of the compound at position $\mathbf{x}\equiv\left[x_{1},\dots,x_{d}\right]^{\mathrm{T}}$ would then be

[TABLE]

where $\rho\left(\mathbf{x}\right)$ is the particle density (particles per unit volume of medium), and $\phi\left(\mathbf{x}\right)$ is the volumetric content of the fluid per unit volume of medium. The particle density $\rho\left(\mathbf{x}\right)$ needs to be estimated from particle positions, and one method to do so is via kernel density estimation (KDE) such that

[TABLE]

where $W$ is the kernel, chosen here as a “product” multi-Gaussian:

[TABLE]

with $\mathbf{h}\equiv\left[h_{1},\dots,h_{d}\right]^{\mathrm{T}}$ being the vector of directional kernel bandwidths. Note that in (2), every particle $p$ can have a different bandwidth $\mathbf{h}_{p}$. Expression (2) has been used in previous works [22, 23, 12, 16, 17] to link particle positions to solute concentrations. Recently, we [17] proposed a technique to determine the optimal bandwidth $\mathbf{h}_{p}$ based on the minimization of the root mean squared error (RMSE) on a local environment. Here, we adapt the density estimation and the bandwidth optimization methods to apply them within a regular grid.

2.1 Concentration in a Bin

Let us discretize our entire domain of interest into $\nu$ regular bins of size $\boldsymbol{\uplambda}\equiv\left[\lambda_{1},\dots,\lambda_{d}\right]^{\mathrm{T}}$, labeled as $u=1,\dots,\nu$. Let us also group particle positions into the centers of the containing bins; i.e., if particle $p$ falls into bin $u$, then $\mathbf{X}_{p}\approx\mathbf{x}_{u}$, where $\mathbf{x}_{u}$ is the position of the center of bin $u$. Then we can define a discrete (mean) value of the particle density $\rho$ in the $u$th bin as

[TABLE]

where $\Lambda=\prod_{i=1}^{d}\lambda_{i}$ is the size of the bin, $\mu_{\omega}$ is the particle count (number of particles) in bin $\omega$, $\mathbf{h}_{\omega}$ is its associated kernel bandwidth, and

[TABLE]

The kernel $\overline{W}$ is a projected form of $W$. By combining (3) and (5) we obtain the closed form:

[TABLE]

Note from (4) that, since all $\mathbf{x}_{u}$, $\mathbf{x}_{\omega}$ belong to a regular grid of size $\boldsymbol{\uplambda}$, then $\mathbf{r}$ in (6) can be written as

[TABLE]

where $\mathbf{z}=\left[z_{1},\dots,z_{d}\right]^{\mathrm{T}}$ is a vector of integers ($z_{i}\in\mathbb{Z}$, $\forall i=1,\dots,d$), and $\odot$ is the Hadamard (element-wise) product. That is, $\mathbf{z}$ is a cell index, with $\mathbf{z}=\mathbf{0}$ corresponding to the cell where the kernel is centered. Then expression (6) can be rewritten as

[TABLE]

By examining equation (8), we see that for any given ratio $\mathbf{h}\oslash\boldsymbol{\uplambda}$ (where $\oslash$ is the Hadamard division) the set of required evaluations of $\overline{W}$ in (4) (setting a cutoff distance) can be fully defined as a $d$-dimensional matrix (see leftmost illustrations in Figure 1). Then, by limiting the possible values that $\mathbf{h}$ can adopt to a discrete set, it is possible to greatly speed up evaluation of (4) by storing these matrices after their first computation, hence never having to re-evaluate (6) for similar values of $\mathbf{h}_{\omega}$. More details on the properties, generation, cut-off correction, storage and use of matrix kernels can be found in A.

Now, the kernel-based evaluation of particle densities on bins consists of two steps. First, a simple binning is performed to obtain the particle count $\mu_{\omega}$ in every bin, which can be seen as an initial, perhaps noisy, “pilot” density estimation. Then, smoothing (4) is performed using the matrix kernels corresponding to bandwidths $\mathbf{h}_{\omega}$ to obtain densities $\rho_{u}$.

Finally, the concentration in bin $u$ can be calculated as

[TABLE]

where $\phi_{u}\coloneqq\phi\left(\mathbf{x}_{u}\right)$. Details on the local selection of parameter $\mathbf{h}_{\omega}$ to use in (4) are given in the following section.

2.2 The Optimal Kernel Bandwidth $\mathbf{h}$

Let us decompose the bandwidth $\mathbf{h}$ by defining the bandwidth scale $\widehat{h}$ and the vector of “shape” parameters $\mathbf{s}\equiv\left[s_{1},\dots,s_{d}\right]^{\mathrm{T}}$:

[TABLE]

i.e., $\mathbf{h}\equiv\widehat{h}\mathbf{s}$, with $\prod_{i=1}^{d}s_{i}=1$. Then, following [17], the local optimal bandwidth scale is

[TABLE]

where

[TABLE]

and

[TABLE]

In (13), the functionals $\varPsi_{ij,u}$ are defined as

[TABLE]

where

[TABLE]

with the kernel function $V^{(i)}$ being defined as

[TABLE]

In (14), $\kappa^{(i)}_{u}$ is defined as the mean value of $\kappa^{(i)}$ in bin $u$:

[TABLE]

with

[TABLE]

The notation $\mathbf{g}_{\omega}^{\left(i\right)}$ in lieu of $\mathbf{h}_{\omega}$ indicates that this bandwidth can adopt a different value than that of kernel $W$, and also different values for different directions of derivative. By combining (16) and (18) we obtain the following closed form:

[TABLE]

Similar to $\overline{W}$, the quantity $\lambda_{i}^{2}\overline{V}^{(i)}$ can be stored as a matrix (Figure 1) with values that only depend on the ratio $\mathbf{g}\oslash\boldsymbol{\uplambda}$. Additional corrections are performed on kernel $\overline{V}^{(i)}$ to ensure that it keeps the main properties of the original kernel $V^{(i)}$ despite the projection. Details on these corrections, as well as the generation, storage and use of matrix kernels can be found in A.

The shape vector is also determined by the “roughness” functionals $\varPsi_{ij,u}$:

[TABLE]

Detailed information on the theory behind these expressions can be found in [17]. Nevertheless, here we give an intuitive explanation: In (11), $n_{u}$ is a smooth average of the particle density over a local environment, whereas $T_{u}$, termed as “effective roughness”, is a measure of the square of the second spatial derivatives of the densities, also averaged over a local environment. Thus, typically, those areas where particle densities are higher, or where they form “peaks” and “valleys” that are more pronounced, will yield smaller kernel bandwidths (lower values of $\widehat{h}_{u}$). In (20), we see that the elongation of the kernel bandwidth in a direction $i$ will be inversely correlated to the (region-averaged) squared derivatives of the density in that direction, normalized by its geometric average over all dimensions; that is, the kernel will tend to stretch along the direction of minimum curvature. In sum, the local kernel size and shape adapts to “mimic” the features of the surrounding particle cloud.

At this juncture, the computation of the optimal bandwidth $\mathbf{h}_{u}$ to use in (4) requires the input of two additional kernel bandwidths: $\boldsymbol{\upsigma}_{u}$ in (12) and (14) ; and $\mathbf{g}_{\omega}^{\left(i\right)}$ in (17). These are addressed in the following two sections, respectively.

2.3 The Integration Support $\boldsymbol{\upsigma}$

Computation of (12) and (14) requires the definition of a Gaussian integration support, represented here by its vector of directional bandwidths $\boldsymbol{\upsigma}_{u}$. In the original development of the local optimization methodology [17], an isotropic, spatially constant support $\sigma_{i}\equiv\widehat{\sigma}$, $\forall i=1,\dots,d$ was proposed, such that $\widehat{\sigma}=3\widehat{h}^{\mathcal{G}}$, with $\widehat{h}^{\mathcal{G}}$ defined as the global AMISE-optimal kernel bandwidth scale. Not only is this approach completely heuristic in nature, but it also renders the local kernel indirectly dependent on a global feature, compromising local benefits. Here, we overcome this problem by introducing the concept of an equivalent normal particle distribution. We assume that the local optimal kernel scale $\widehat{h}_{u}$ is also the global optimal kernel bandwidth $\widehat{h}_{u}^{\sigma}$ associated with a virtual Gaussian distribution of variance $\widehat{\sigma}_{u}^{2}$, composed of a virtual number of particles $N_{u}^{\sigma}$; the classic expression for the AMISE-optimal kernel bandwidth in this case is [24]:

[TABLE]

We then impose that this virtual distribution has locally matching values with the actual particle distribution for $\rho_{u}$ and $n_{u}$, the latter being defined as in (12), using $\boldsymbol{\upsigma}_{u}$ as the integration support. Then, it can be shown (see B) that the following relation holds,

[TABLE]

Combining (21) with (22), and after some algebraic manipulation, we obtain:

[TABLE]

The integration support used in (12) and (14) is then $\boldsymbol{\upsigma}_{u}=\widehat{\sigma}_{u}\mathbf{1}$, where $\mathbf{1}$ is a $d\times 1$ vector of ones. Note that expression (23) is recursive, in the sense that $h_{u}$, $n_{u}$ and $\rho_{u}$ are all affected by some previous choice of $\widehat{\sigma}_{u}$. Nonetheless, it can be implemented iteratively. More details on the iterative implementation are given in $\S$2.5.

2.4 The Curvature Kernel Bandwidth $\mathbf{g}^{\left(i\right)}$

Computation of the optimal bandwidth as presented in $\S$2.2 requires a bandwidth for the estimation of the particle density curvatures ($\mathbf{g}^{(i)}_{\omega}$ in (15)). Previously [17], the Improved Sheather-Jones plug-in method by Botev et al. [25] was used to determine this bandwidth. Here, we rely on the equivalent Gaussian particle distribution described in $\S$2.3 to determine a local value for $\mathbf{g}^{(i)}_{\omega}$ recursively based on $\mathbf{h}_{\omega}$.

Above, we have defined both these bandwidths as diagonal. Here we further assume that $\mathbf{g}^{(i)}$ is isotropic,

[TABLE]

We also assume that the anisotropy of $\boldsymbol{\upsigma}_{u}$ is identical to that of bandwidth $\mathbf{h}_{u}$ (represented by vector $\mathbf{s}_{u}$). Then, in a similar fashion to what was done in $\S$2.3 for $\widehat{h}_{u}$, the value of bandwidth scale $\widehat{g}_{u}^{(i)}$ is assumed to match its $\mathrm{AMISE}$-optimal magnitude within the virtual Gaussian distribution of particles defined in $\S$2.3 through $\boldsymbol{\upsigma}_{u}$ and $N_{u}^{\sigma}$. This magnitude (see derivation in C) is

[TABLE]

where $\vartheta_{i}$ is a function of the anisotropy of $\mathbf{s}_{u}$, with a value of 1 in the isotropic case (i.e., $s_{i}=1,\ \forall i=1,\dots,d$):

[TABLE]

where $\delta_{ij}$ is the Kronecker delta. By combining (25) and (21), we can determine the ratios between optimal bandwidth scales:

[TABLE]

[TABLE]

with $N_{u}^{\sigma}$ obtained from (22). It is worth noting that $\alpha$ in (27) is a constant value close to 1 (ranging from $1.04$ for $d=1$ to $1.12$ for $d=3$), and $\beta\ll 1$ (meaning that the ratio $\gamma^{(i)}_{u}$ is very rigid with respect to $N_{u}^{\sigma}$). Hence $\widehat{g}_{u}^{(i)}$ is typically just somewhat larger than $\widehat{h}_{u}$ for a wide range of values of $N_{u}^{\sigma}$. For instance, in 2D, assuming isotropy ($\vartheta_{i}=1$), $\gamma^{(i)}_{u}\approx 1.31$ for $N_{u}^{\sigma}=10^{2}$, and $\gamma^{(i)}_{u}\approx 1.75$ for $N_{u}^{\sigma}=10^{5}$. Using definition (27), we can recursively obtain $\widehat{g}_{u}^{(i)}=\gamma^{(i)}_{u}\widehat{h}_{u}$, and then, in equation (17), we use $\mathbf{g}_{u}=\widehat{g}_{u}\mathbf{1}$.

2.5 Optimization Algorithm

The approach that we propose to determine the locally optimal bandwidth can be seen as a fixed-point iteration method: A set of local kernel bandwidths is given as an input, and a new set is obtained as an output, until convergence. Below, for any variable “$a$”, the notation $\{a\}$ indicates the set of all $a_{u}$. Before starting the iteration process, the particle counts $\{\mu\}$ are computed, and, if unavailable, an initial pilot $\{\widehat{\sigma}\}$ is defined such that $\widehat{\sigma}_{u}=3\widehat{h}_{u}$. Then, the structure of one iteration can be summarized as follows:

Compute $\{\rho\}$ by (4) using the input $\{\mathbf{h}\}$. 2. 2.

To obtain $\{\widehat{\sigma}\}$ and $\{n\}$:

Compute pilot $\{n\}$ by (12) using $\{\rho\}$ and $\{\widehat{\sigma}\}$.

Compute $\{\widehat{\sigma}\}$ by (23) using $\{\widehat{h}\}$, $\{\rho\}$ and pilot $\{n\}$.

Compute $\{n\}$ by (12) using $\{\rho\}$ and $\{\widehat{\sigma}\}$. 3. 3.

For $i=1,\dots,d$:

Compute $\{{\widehat{g}}^{(i)}\}$ by (27) using $\{\rho\}$, $\{n\}$, $\{\widehat{\sigma}\}$ and $\{\mathbf{s}\}$.

Compute $\{\kappa^{(i)}\}$ by (17) using $\{g^{(i)}\}$. 4. 4.

For $i=1,\dots,d$, for $j=i,\dots,d$:

Compute $\{\varPsi_{ij}\}$ by (14) using $\{\widehat{\sigma}\}$, $\{\kappa^{(i)}\}$ and $\{\kappa^{(j)}\}$. 5. 5.

Compute the new $\{\mathbf{h}\}$ by (11) and (20) using $\{n\}$ and all $\{\varPsi_{ij}\}$.

The iteration process may be exited when relative changes in $\{\mathbf{h}\}$ are below some tolerance level. In the context of a full RWPT simulation, since changes in the optimal kernel bandwidth typically occur at a much slower pace than solute transport [17], the optimization does not need to be performed at every time step of the transport time discretization, and if performed often enough, a single iteration may suffice. In between optimizations, particles should “carry” the identifier of the bin where they were located at the time of the latest optimization; then, the bandwidth at a bin is the average of the bandwidths associated to the contained particles.

3 Boundary Condition Corrections

Kernel methods typically suffer from boundary bias problems [e.g., 26]. Near no flux boundaries, the standard KDE as written in (2) may for instance produce an underestimation bias; as its support can cross the boundary in question, it may unphysically project nonzero solute concentrations on the other side of such boundaries, where no such mass can actually exist. Similar problems would occur for constant concentration or fixed flux boundaries, where without some correction, kernels will project incorrect and unphysical concentrations at and across the boundaries. Such problems also arise with the proposed discrete estimator in (4). However, for the discrete case, the position of the bins with respect to the boundaries is known and does not change over the course of a simulation. This allows us to efficiently introduce corrections to account for the presence and influence of boundaries, that would otherwise be unfeasible (or at least extremely cumbersome) directly from (2).

Here we propose an approach based on the assumption that the Gaussian kernel that represents a particle’s support volume behaves like a purely diffusive process, whose interactions with boundaries are well understood. In other words, we treat the kernel associated with a particle as if it were the result of a very fast diffusive process that has resulted in the spatial spreading of the particle’s mass, with an initial condition of a Dirac delta located at the particle’s physical position. If the particle is far from the boundary, this results in the standard Gaussian kernel (3). Previous works [e.g., 25] acknowledge this link between KDE and diffusion. The correction on the kernel that we propose is independent of the implementation of the boundary condition itself on the RWPT algorithm [see 18]; that is, the particles undergoing the random walks in the system are unaware of the kernel and so their motion must still account for the presence of a boundary (e.g. reflection on a no flux boundary). Note also that, in principle, our proposed approach could be extended to any kernel-based Lagrangian method [e.g., 27] in a bounded domain with physical boundary conditions. For most commonly used boundary conditions, we derive the simple reflection principles that can be used to modify the kernel near a regular boundary that is aligned with the principal directions of the kernel. Then, we also propose an approach to extend this procedure to the case of irregularly-shaped boundaries, which present unique challenges.

3.1 Regular Boundaries

All boundary conditions most typically used in transport simulations can be written in terms of a balance of mass fluxes between the inner ($+$) and the outer ($-$) side of the boundary; i.e.:

[TABLE]

where $\mathbf{n}$ is the unit vector normal to the boundary. Depending on the assumptions made for the outer side of the boundary, this balance can result in a variety of common boundary conditions.

3.1.1 Impermeable/Outlet

If the boundary is assumed to be a no flux boundary (impermeable), we set $\mathbf{n}^{\mathrm{T}}\mathbf{q}=0$, and $\left(\phi\mathbf{D}\boldsymbol{\nabla}c\right)^{-}=0$. In this case, equation (29) becomes the homogeneous Neumann boundary condition, which, for boundary $\mathfrak{N}$, can be written as

[TABLE]

Note that (30) would also be obtained at an outlet, by assuming $\mathbf{q}c^{-}=\mathbf{q}c^{+}$, and $\left(\phi\mathbf{D}\boldsymbol{\nabla}c\right)^{-}=0$.

Let us consider a particle $p$ located at $\mathbf{X}_{p}$. The density kernel associated with this particle is $W\left(\mathbf{x}-\mathbf{X}_{p};\mathbf{h}_{p}\right)$, if we neglect the effect of the boundary. As noted above, let us assume that this density distribution can be seen as the Green’s function of a virtual, fast diffusive process. Then, if the boundary is regular and aligned with the principal directions of the kernel, the kernel is altered by the proximity to $\mathfrak{N}$ by a reflection principle such that

[TABLE]

where $\widetilde{\mathbf{X}}_{p}$ is the mirror of $\mathbf{X}_{p}$ through $\mathfrak{N}$:

[TABLE]

Because of the regularity assumption,

[TABLE]

where

[TABLE]

which allows the corrected density estimation (i.e. replacing $W$ with $W_{\mathfrak{N}}$ in (2)) to be expressed as a function of the uncorrected density estimation:

[TABLE]

Hence, densities can be computed conventionally by (2), then the corrected density at a point is obtained by adding the density for the mirror symmetry point. It follows directly, that for the on-grid methodology presented in this work,

[TABLE]

where $\rho_{\widetilde{u}}$ is the uncorrected density (4) computed at the mirror bin $\widetilde{u}$, whose center is located at

[TABLE]

3.1.2 Dirichlet

Assuming in (29) that the dispersive flux, $\phi\mathbf{D}\boldsymbol{\nabla}c$, is the same on the outer and inner sides of the boundary, and that the outer concentration is prescribed, we obtain a Dirichlet boundary condition,

[TABLE]

where $c_{\mathrm{o}}$ is the prescribed concentration. In this case, the boundary can be permeable to a diffusive process. This is true in both directions, so part of a particle’s mass may fall outside the domain, and mass may enter the domain as well. Based again on the assumption that the kernel represents a fast diffusive process, we can separate this process into two parts following the superposition principle, first accounting for the initial condition with a homogeneous boundary condition (a perfectly absorbing boundary), and second accounting for the inhomogeneous part:

[TABLE]

The solution to the homogenous part, $W_{\mathfrak{D}\mathrm{H}}$, given a regular boundary aligned with the principal directions of the kernel, is [28]:

[TABLE]

The inhomogeneous part can also be solved under the same assumptions. Constant diffusion from a regular boundary follows the 1D analytical solution in the direction normal to the boundary [29]:

[TABLE]

with

[TABLE]

assuming a constant $\phi\approx\phi\left(\mathbf{x}_{\mathfrak{D}}\right)$ near the boundary. Expression (41) is also the Green’s function of diffusion for an initial condition of uniform density $2\rho_{\mathrm{o}}$ at the outer side of the boundary [18]. Thus, we note that equation (41) is approximately equivalent to a sum of mirror kernels (as in (39)) provided that

[TABLE]

where $\rho_{*}$ is the unknown true density. We replace it by a pilot estimate at the containing bin $\omega$, i.e.,

[TABLE]

with $\mu_{\mathrm{o}}\coloneqq\Lambda\rho_{\mathrm{o}}$. Then, substituting (40) and (44) into (39), and integrating as in (4), the on-grid corrected density estimate becomes:

[TABLE]

where

[TABLE]

where mirror bin $\widetilde{u}$ is defined in analogy with (37). Expression (45) is somewhat similar to (36); however now the simple reflection $\rho_{\widetilde{u}}$ is replaced by a new term, $\rho_{\widetilde{u}}^{\mathrm{o}}$. Examining expression (46) we see that the only modification in the smoothing process is that, when transferring density from any non-empty bin $\omega$ to an external bin $\widetilde{u}$, it must be done as if the bin $\omega$ contained a virtual number of particles equal to $2\mu_{\mathrm{o}}-\mu_{\mathrm{\omega}}$, instead of the actual value of $\mu_{\omega}$.

3.1.3 Robin

If we assume in (29) that the outer side of the boundary is a reservoir with a prescribed concentration, i.e., $\left(\phi\mathbf{D}\boldsymbol{\nabla}c\right)^{-}=0$ and $c^{-}=c_{\mathrm{o}}$, what we obtain is a Robin (or third-type) boundary condition, which can be written for boundary $\mathfrak{R}$ as

[TABLE]

where $\mathbf{v}=\mathbf{q}/\phi$. Like in the previous cases, we treat the kernel as a very fast “diffusive” process, which interacts with the boundary analogously to the Green’s function of the simulated dispersion. Since, for very large $\mathbf{D}$, condition (47) converges to (30), a Robin boundary affects the kernel identically to a homogeneous Neumann boundary (see $\S$3.1.1):

[TABLE]

This would also apply to similar Robin boundary conditions originating from other assumptions, such as reactive walls [e.g., 30].

3.2 Extension to Irregular Boundaries

The reflection principles given in $\S$3.1 are derived based on two assumptions: that the boundaries are regular and that they are aligned with the principal directions of the kernel. The latter can always be fulfilled by using an isotropic kernel (restricting the degrees of freedom of the bandwidth, see [17]). However, the assumption of regular boundaries may just not be fulfilled, and then expressions (36), (45) and (48) are not valid as they rely on substituting the mirroring of the kernel for that of the evaluation point through (33). Moreover, in that case, the true Green’s function of the diffusion process associated to bandwidth $\mathbf{h}_{\omega}$ is much more complex than the solutions given by (31), (40) and (43). Nevertheless, as an approximation, one can use these solutions to modify the kernel directly, defining the mirror points by reflection through the closest boundary point, and with a weighting parameter to compensate for the irregularity. For impermeable or Robin, we propose the following kernel:

[TABLE]

where $\eta_{\omega}$ is a weighting parameter such that the mass of the “reflection” kernel $\overline{W}\left(\widetilde{\mathbf{x}}_{\omega}-{\mathbf{x}}_{u};\mathbf{h}_{\omega},\boldsymbol{\uplambda}\right)$ entering the domain matches the mass of the “regular” kernel $\overline{W}\left(\mathbf{x}_{\omega}-\mathbf{x}_{u};\mathbf{h}_{\omega},\boldsymbol{\uplambda}\right)$ falling outside the domain. Thus, $\eta_{\omega}$ will be higher than one for a convex boundary and lower than one for a concave boundary. Essentially, the problem of non-bijection of the mirror points is fixed through this parameter. For Dirichlet boundary conditions:

[TABLE]

with the same definition for $\eta_{\omega}$. In the event that, because of the irregularity, the use of kernel (50) may generate negative density estimations in some bins, they should be set to zero.

See $\S$4.3.4 for an example implementation of the proposed approach.

4 Computational Investigations and Discussion

In this section we investigate the performance of the proposed methodology on hypothetical RWPT models, with a focus on the features that are new with respect to the original approach introduced in [17]. This includes the Eulerian grid projection, the iterative nature of the local kernel optimization, and the accounting for boundary conditions, both regular and irregular. In D we also discuss the use of local auxiliary parameters in the optimization process. We study each of these aspects separately, and then we focus on fixed-time concentration estimations rather than full reactive simulations, as it should be clear that a better concentration estimation will result in improved reactive transport modeling [17].

4.1 Grid-Projected KDE

The use of kernel functions in particle-based methodologies is common for non-linear processes that involve interaction between individual particles, such as chemical reactions. Traditionally, the representation of the support volume of a particle as a smoothing kernel results in a loop through all pairs of potentially interacting numerical particles [16]. Let us consider the simple example of an irreversible bimolecular reaction of type,

[TABLE]

where $\mathrm{A}$, $\mathrm{B}$ and $\mathrm{C}$ are chemical compounds, and the “batch” reaction kinetics (that is, neglecting transport) can be described as:

[TABLE]

where $k$ is the kinetic reaction constant and $c_{\mathrm{A}}$, $c_{\mathrm{B}}$ and $c_{\mathrm{C}}$ are the concentrations of compounds $\mathrm{A}$, $\mathrm{B}$ and $\mathrm{C}$, respectively. In this case, for a particle $a$ of compound $\mathrm{A}$, we could estimate its probability of reaction in a time step $\Delta t$ as [17]:

[TABLE]

For the sake of simplicity, $\phi$ in (53) is assumed constant. Here, the density of B-particles $\rho_{\mathrm{B}}$ is estimated from equations (1) and (2) at position $\mathbf{X}_{a}$ of the A-particle $a$. Doing this for all potentially reactive $\mathrm{A}$-particles would involve a double loop to see interaction between all particle pairs, and therefore scales in number of calculations as $N_{\mathrm{A}}N_{\mathrm{B}}$. In addition, it also requires a search algorithm, which would scale at best as $N_{\mathrm{A}}\log N_{\mathrm{B}}$. Aside from incorporating several kinds of boundary conditions, as addressed in $\S$3, a powerful argument for performing a pilot binning before the kernel smoothing is to increase computational efficiency, by pre-grouping the particles and avoiding the need for search algorithms; additionally one can pre-compute and store the matrix kernels, thus not having to evaluate the kernel function continuously. In this way, we simultaneously benefit from the low computational demands of binning as well as the accuracy of KDE.

In order to exemplify this, consider a simple hypothetical 2D problem where, at a given time, compounds $\mathrm{A}$ and $\mathrm{B}$, each represented by $N_{\mathrm{A}}=N_{\mathrm{B}}=N$ particles, are distributed in space as partially overlapping multi-Gaussian distributions. Both these distributions have isotropic, unit variance $\sigma^{2}=1$, and their centers $\langle\mathbf{X}_{\mathrm{A}}\rangle$, $\langle\mathbf{X}_{\mathrm{B}}\rangle$, are separated by a distance $0.8\sigma$.

To perform a reactive time step given the described conditions, we consider three possible alternatives to estimate $\rho_{\mathrm{B}}(\mathbf{X}_{a})$ in (53), for all $a=1,\dots,N$: (i) Binning, (ii) Traditional KDE (eq. (2)), and (iii) The KDE method proposed herein (eq. (4)).

To independently evaluate and compare these techniques, we use a constant (global), isotropic kernel bandwidth of size $\widehat{h}=\sigma N^{-\frac{1}{6}}$ (see equation (10)). In the limit of $N\rightarrow\infty$ and $\boldsymbol{\uplambda}\rightarrow\mathbf{0}$, the estimated density at a point should converge to the true solution

[TABLE]

We measure the difference between $\rho_{\mathrm{B}}$ and $\rho_{\mathrm{B}}^{*}$ through the normalized root mean squared error (NRMSE):

[TABLE]

Additionally, we compute the time spent to perform the density estimation. Figure 2 compares $\mathrm{NRMSE}$, $N$ and CPU time, given different values of isotropic bin size $\lambda$.

For fixed $N$ and $\lambda$, the proposed density estimation technique is always more accurate than binning and more efficient than the traditional KDE. As a result, given a desired degree of accuracy, or equivalently, a fixed spatial scale of interest $\lambda$, we have strong evidence that the proposed technique is the optimal choice, in terms of attainable computational efficiency, for estimating concentrations and reaction rates. For a very high value of $N$, as the kernel bandwidth $\widehat{h}$ becomes small in comparison to $\lambda$, the proposed method and binning converge in terms of performance and accuracy, because there is in fact no effective smoothing. Prior to reaching that point (and mostly, in areas of lower particle density), the kernels are able to successfully and efficiently make up for the lack of particles.

4.2 Local Bandwidth Optimization

Next, the local kernel bandwidth optimization and density estimation algorithm is tested in a synthetic example of advective-dispersive transport in a 2D heterogeneous porous medium (see Figure 3(a)). The spatial distribution of log hydraulic conductivity $K$, $Y=\log K$, in this domain of size ($80\times 50\ \mathrm{m}$) is built in 3 steps: first, (i) a random multi-Gaussian field is generated with mean $\langle Y\rangle=0$, variance $\langle Y^{2}\rangle=1$, and an exponential isotropic variogram with integral scale $I_{Y}=3\ \mathrm{m}$; (ii) it is then multiplied with a field that evolves linearly in the vertical direction from 0 at the top and bottom boundaries to 1 at the center; (iii) finally, some low-conductivity ($Y=-2.5$) inclusions measuring $2.5\times 13.5\ \mathrm{m}$ are added as shown in Figure 3(a). The objective of this multistep procedure is to generate multiple different local features, that could be used as a test for the local bandwidth selection algorithm.

Water flows through the porous medium following:

[TABLE]

where $\mathbf{q}$ is the Darcy velocity and $H$ is the hydraulic head. The top and bottom boundaries are impermeable, and the head is prescribed at the left and right boundaries to force a mean hydraulic gradient of $2.5\%$, generating a mean flow from left to right. The domain is discretized in square cells of $0.5\times 0.5\ \mathrm{m}$, and $\mathbf{q}$ at the cell interfaces is obtained via MODFLOW 2005 [31].

As an initial condition for transport, a uniform, rectangular distribution of $1.8\cdot 10^{5}$ particles (representing a solute) is injected at $t=0$ near the left boundary, as shown in Figure 3. For time intervals $\left[t,t+\Delta t\right]$, particles move by random walk particle tracking (RWPT) [1]:

[TABLE]

which is equivalent to solving the advection-dispersion equation (ADE). In (57), $\mathbf{A}\coloneqq\frac{1}{\phi}\left[\mathbf{q}+\nabla\cdot\left(\phi\mathbf{D}\right)\right]$, with $\mathbf{D}$ being the dispersion tensor; $\mathbf{B}$ is a $d\times d$ matrix such that $\mathbf{B}\mathbf{B}^{\mathrm{T}}=2\mathbf{D}$; and $\boldsymbol{\upxi}$ is a $d\times 1$ vector of standard-normally distributed random numbers, uncorrelated in time. The spatially variable, anisotropic dispersion tensor $\mathbf{D}$ is determined following [32], with a longitudinal dispersivity of $\alpha_{\ell}=2\cdot 10^{-3}\ \mathrm{m}$, a transverse dispersivity of $\alpha_{t}=3\cdot 10^{-4}\ \mathrm{m}$, and a molecular diffusion of $D_{m}=2\cdot 10^{-4}\ \mathrm{m}^{2}/\mathrm{day}$. The spatial interpolation of velocities and the dispersion tensor is done using the hybrid linear-bilinear method proposed by [33].

Figure 3(b) shows the spatial distribution of the particle density $\rho$ after 600 simulated days, estimated from the particle position information using the methodology presented in $\S$2. In Figure 4, we show the iterative process of bandwidth differentiation and convergence (see $\S$2.5), starting from an arbitrary, uniform isotropic bandwidth with $\widehat{h}=0.45\ \mathrm{m}$, within a zoomed-in region of the domain (delimited by the dashed lines shown in Figure 3). The estimation of the curvatures $\kappa^{(i)}$, included in Figure 4, is crucial for the correct determination of the locally optimal scale ($\widehat{h}$) and elongation ($s_{1}$) of the smoothing kernel bandwidth. We see for this particular example that, after 4 iterations, the solution has nearly stabilized. The smoothing of the “pilot” (binning) concentrations through this optimal local kernel is able to visibly reduce the noise, without generating excessive over-smoothing (see the two plots on the left of Figure 4). As a result, the NRMSE was reduced from a 93.7% to a 7.0% after the smoothing, with the NRMSE in this case being defined as in (55), but for $\rho$ on all $\mathbf{x}_{u}$, and with the “true” solution in this case being approximated as the binning solution obtained for $N=8.64\cdot 10^{7}$ ($\sim 2.5$ orders of magnitude more than the test case).

Another relevant property of the iterative optimization algorithm is its robustness, i.e. its ability to reach the same solution given different initial values. This is particularly important in the context of particle tracking simulations, where the particles can “carry” the support volume for some number of steps and then use it as an input to update the next optimal support volume, conferring it an evolutionary nature. For this to make sense, it is critical that for any given current distribution of particles, the solution always converges to a unique set of values, regardless of the history of the local kernel bandwidth.

To test this, we repeat the previous numerical experiment, but this time, the initial bandwidth is random and uncorrelated in space, with $\widehat{h}$ drawn from a uniform distribution with lower and upper bounds of $0.1$ and $0.8$; and with $s_{1}$ drawn from a log-uniform distribution (the logarithm is uniformly distributed) with lower and upper limits of $0.5$ and $2.0$. These are deliberately chosen to be challenging for the convergence of the algorithm. As can be seen from Figure 5, in the first iteration this results in a curvature estimation ($\kappa^{(1)},\kappa^{(2)}$) that is far from the correct solution, and characterized by high absolute values and strong variations with no apparent spatial correlation. As a consequence, the first estimation of $\widehat{h}$ and $s_{1}$ gives a small bandwidth without a well-defined elongation direction. Yet, after a few iterations, we see that the identification of the curvatures substantially improves, and concurrently the bandwidth scale and elongation start differentiating spatially distinct regions. For this highly adverse case of a locally random and uncorrelated input value, after only 7 iterations we have nearly reached the same stable solution as in Figure 4.

4.3 Boundary Conditions

Here we implement and evaluate the boundary condition correction techniques described in $\S$3. First, we perform a concentration estimation near the boundary in three simple 1D RWPT settings. In all cases, the solute moves by dispersion with $D=0.1\ \mathrm{m}^{2}/\mathrm{day}$ for a total time of $\tau=1000\ \mathrm{days}$; the medium porosity is $\phi=0.25$. Each particle has a mass (or more precisely, an amount of substance) of $m=10^{-4}\ \mathrm{mol}$. The bin size is $\lambda=0.5\ \mathrm{m}$. We compare the three estimation methods (the pilot binning, the uncorrected KDE, and the corrected KDE) to the analytical solution for each case.

4.3.1 Diffusion near impermeable wall

In the first example, the initial condition is a Dirac delta pulse of solute (with a total mass $M=1\ \mathrm{mol}$) located at $x=10\ \mathrm{m}$, near an impermeable boundary located at $x_{\mathfrak{N}}=0$. The results are shown in Figure 6(a). Since here diffusion is the only physical process, and the corrected kernel emulates diffusion, the concentration estimation has an excellent agreement with the analytical solution.

4.3.2 Diffusion through boundary

In the second example, all initial concentrations are zero inside the domain. There is a Dirichlet boundary condition such that $c_{\mathrm{o}}=0.18\ \mathrm{mol}/\mathrm{m}$ at $x_{\mathfrak{D}}=0$. The results are shown in Figure 6(b). The Dirichlet “reflection” technique (45) is able to correctly reconstruct the concentration field and gradient near the boundary. It is worth noting that a simple reflection as in (36) would have resulted in a zero-gradient instead, as shown by the curve labeled as “Simple Reflection”. Nevertheless, the error involved in using (36) instead of (45), at least in this specific case, was relatively small (compared, for instance, to the error associated with not performing a boundary correction at all).

4.3.3 Column with inlet and outlet reservoirs

In this third example we consider a constant rightwards advection ($q=0.055\ \mathrm{m}/\mathrm{day}$) along with a linear degradation:

[TABLE]

The reaction is simulated stochastically by a particle reaction probability $P=k\Delta t$ (i.e., this is the probability a particle dies in any given time step). At the inlet we have a Robin boundary condition as described in $\S$3.1.3, and at the outlet we set a zero-gradient (homogeneous Neumann) condition. This problem has a stationary solution. From Figure 6(c), we see that the corrected density estimation has a zero-gradient near the boundaries, and that the gradient is not zero in the true solution. This is because the kernel reflection method is based on pure diffusion and unable to account for the gradient generated by the combined action of advection and reaction. Nonetheless, one can appreciate in Figure 6(c) a substantial improvement for the corrected KDE method with respect to the uncorrected one.

4.3.4 Irregular Boundaries

In order to illustrate the boundary correction technique for the case of irregular boundaries explained in $\S$3.2, we designed a problem consisting of a domain with the shape of an arched tube (Figure 7) with impermeable boundaries. The region shown has dimensions $115\times 115$, and the spatial discretization is square with $\lambda=1$. The particle cloud (Figure 7(a)) is the result of a random injection of $100\,000$ particles. In polar coordinates (the origin being the center of curvature of the tube), the square of the radial coordinate of the particle positions is uniformly-distributed, whereas the angular coordinate is Gaussian-distributed. Such a distribution could be thought of as a hypothetical result of advection-dispersion through the tube.

Despite the high number of particles used, we observe in Figure 7(b) a considerable amount of noise in the binning estimation. This issue is fixed by our locally adaptive kernel smoothing technique, as shown in Figure 7(c), but at the expense (prior to any boundary correction) of an artificial density loss near the boundaries (see zoomed-in region). As illustrated by the representation of the kernel $\overline{W}$ near the boundary, this is caused by the lack of a reflection of the mass lost through the boundary by the smoothing. In Figure 7(d), we see that the modified kernel (see zoomed-in representation) successfully fixes the loss of mass through the boundary, yielding a satisfactory reconstruction of the concentrations.

5 Summary and Conclusions

We have presented a novel technique to estimate particle densities using the limited amount of information provided by a finite sample of particle positions. Although the spectrum of possible applications is wide, the focus of this work is the reconstruction of solute concentrations in random walk particle tracking (RWPT) simulations, which is relevant for the visualization of results and the incorporation of chemical reactions. Our technique relies on the accuracy of locally adaptive kernel density estimation (KDE), which is implemented in combination with a spatial discretization, resulting in benefits including computational efficiency and accurate implementation of boundary conditions. The method is valid in 1, 2 or 3 spatial dimensions. In principle, it can be applied to Lagrangian modeling techniques other than RWPT, or even to other applications of density estimation. An open-source MATLAB code, “Bounded Adaptive Kernel Smoothing” (BAKS), has been developed and made available to the community as a result of this work [34].

Our computational investigations provide strong evidence that our proposed methodology deals well with the dilemma between the accuracy of kernel methods and the low computational requirements of binning. From our tests, we see that the desired degree of accuracy (characterized by the choice of bin size $\boldsymbol{\uplambda}$) is always achieved faster (in terms of invested computational effort) with our proposed methodology compared with binning. The two methods converge, both in accuracy and performance, for high particle numbers, but prior to reaching that point, the kernel approach reaches an acceptable level of error earlier. Likewise, direct Lagrangian utilization of the kernel functions involves a worse scaling of CPU time with particle number than our method. As a result, the proposed method can achieve the same level of accuracy with lower computational effort.

Through a sample simulation of conservative transport in a heterogeneous porous medium, we show the convergence of the local kernel bandwidth iterative optimization method towards a stable result. Using that optimal local bandwidth to smooth the pilot binning density estimation resulted in a reduction of normalized error from $93.7\%$ to $7.0\%$. We also demonstrate the robustness of the method, in terms of being able to reach this same solution regardless of the initial input values provided. The fully local nature of the optimization process, which is a novel aspect with respect to the original methodology [17], allows the kernel to achieve a higher degree of local differentiation in terms of size and shape.

Near boundaries, the kernel is assumed to emulate the Green’s function of a fast diffusion process. Hence, the kernel is affected by the boundary following the analytical solution of pure diffusion corresponding to the relevant boundary condition, which results in simple and efficient reflection rules. The case of irregular boundaries is slightly more complicated, as it requires individual modification of each kernel, with a weighting of the reflection to ensure proper mass conservation. The simple implementation examples illustrate the importance of including these boundary corrections in the density estimation.

Appendix A Details on treatment of matrix kernels

The projected kernels $\overline{W}$ and $\overline{V}$, introduced in $\S$2.1 and $\S$2.2, respectively, require some corrections to ensure that they keep the main properties of the original kernels $W$, $V$, after the on-grid projection. In the case of $\overline{W}$, the original Gaussian kernel $W$ integrates exactly to 1 only over $\mathbb{R}^{d}$. If a cutoff distance is imposed, a normalization is needed to impose mass conservation:

[TABLE]

where $\overline{W}\left(\boldsymbol{\uplambda}\odot\mathbf{z};\mathbf{h},\boldsymbol{\uplambda}\right)$ is the unmodified kernel at cell index $\mathbf{z}$ (equation (8)), $\mathbf{z}=\mathbf{0}$ being the cell where the center of the kernel is located, and $\overline{W}^{\prime}$ is the modified kernel. In (59), the denominator is a sum over all cells within the cutoff limits.

In the case of $\overline{V}$, there are two corrections that need to be performed in order to conserve the original purpose of $V$. On one hand, the positive values must be weighted so that the kernel integrates to zero within the cutoff limits:

[TABLE]

[TABLE]

A simple intuitive example of the importance of this correction is that a constant particle density must necessarily yield a zero-curvature estimation, and this will only occur if the curvature kernel has zero-mean.

Besides, the final purpose of $\overline{V}^{(i)}$ is the estimation of the squared second spatial derivatives of the particle density ($\kappa^{(i)}_{\omega}\kappa^{(j)}_{\omega}$ in (14)). The spatial averaging involved in the grid projection could lead to a systematic under-prediction of these values. To prevent that, we scale the kernel so that it keeps the $\mathrm{L}^{2}$-norm of the original kernel $V$ after the projection:

[TABLE]

where $\left\lVert\ \right\rVert^{2}$ is the $\mathrm{L}^{2}$-norm operator, i.e.,

[TABLE]

$\overline{V}^{\prime\prime(i)}$ are the final values of the projected curvature kernel.

As mentioned in the main body of this work, the values of $\overline{W}$ and $\lambda_{i}^{2}\overline{V}^{(i)}$ only depend on the directional ratios between the bandwidth ($\mathbf{h}$ or $\mathbf{g}$) and the grid size ($\boldsymbol{\uplambda}$). Discretizing the values that these ratios are allowed to adopt, and imposing a cutoff distance (as in Figure 1), $\overline{W}$ and $\lambda_{i}^{2}\overline{V}^{(i)}$ become finite sets of matrices with a finite number of entries. Repeated use of the same (or very similar) kernel bandwidth results in the exact same matrix kernel, which can be stored in the memory after its first generation and correction. Then, computation of (4), (17), (12) and (14) only requires accessing the pre-computed matrix entries in the memory and performing the relevant weighted summation, avoiding redundant computational efforts.

Appendix B Derivation of expression (22)

Within a virtual Gaussian distribution of $N_{u}^{\sigma}$ particles, centered at $\boldsymbol{\upmu}_{u}\equiv[\mu_{u,1},\dots,\mu_{u,d}]^{\mathrm{T}}$ with the vector of directional standard deviations $\boldsymbol{\upsigma}_{u}\equiv[\sigma_{u,1},\dots,\sigma_{u,d}]^{\mathrm{T}}$,

[TABLE]

If $\boldsymbol{\upsigma}_{u}$ is also the bandwidth of the integration kernel, then following (12),

[TABLE]

Here, the approximation $\rho(\mathbf{x}_{u})\approx\rho_{u}$ is equivalent to the one in (12). Taking the square on both sides of expression (65) and rearranging we obtain

[TABLE]

which is the expression given in (22).

Appendix C Curvature kernel bandwidth for Gaussian distribution

Following the approach of [25, Appendix E], here generalized for $d$ dimensions, the optimal isotropic $\widehat{g}^{(i)}$, given a distribution of particles $\rho(\mathbf{x})$, will be

[TABLE]

with $N$ being the total number of particles. Taking $\rho(\mathbf{x})$ to be a Gaussian distribution of $N$ particles with a vector of directional standard deviations such that $\boldsymbol{\upsigma}=\widehat{\sigma}\mathbf{s}$, with $\prod_{i=1}^{d}s_{i}=1$, we have

[TABLE]

and then substituting (68) in (67) and rearranging we obtain

[TABLE]

which is equivalent to the combination of equations (25) and (26) given in $\S$2.4.

Appendix D Local vs Global selection of $\boldsymbol{\upsigma}$ and $\mathbf{g}$

Within the procedure presented in $\S$2, the selection of the curvature kernel bandwidth ($\mathbf{g}_{u}^{(i)}$) at a bin $u$ is completely local and independent of all particles located outside the range determined by the integration support $\boldsymbol{\upsigma}_{u}$. In the original methodology presented in [17], a global $\mathbf{g}^{(i)}$ was used over the whole domain, as noted in $\S$2.4. As a result, the local kernel bandwidth would be indirectly conditioned by global features. For instance, if the plume was dominantly formed by strongly elongated shapes in one specific direction, the curvature estimation anywhere would be biased towards detecting those kinds of features. This would limit the ability of the kernel bandwidth to optimally adapt strictly to the nearby particle distribution. To illustrate this, we compare the distribution of local bandwidths $\mathbf{h}$ obtained with the method presented in $\S$2 to the one that we obtain when using global (instead of local) values for $\boldsymbol{\upsigma}$ and $\mathbf{g}$, following what is described in [17]. We use the same example particle distribution of $\S$4.2.

Figure 8 shows the local bandwidth scale ($\widehat{h}$) and elongation ($s_{1}$) values obtained with the two described approaches. We clearly observe a higher degree of differentiation, that is, an increased ability of the bandwidth scale and elongation to reach extreme values when using the novel, fully local methodology presented here. This is particularly true for large bandwidths, which can only exist in the absence of noise in the curvature estimation (see Figure 5 as an example). The distribution of bandwidth scales and elongations also appears to be significantly smoother (less affected by “spurious” fluctuations) with the novel methodology, which suggests a more accurate identification of the dominant shapes of the local particle distribution.

Quantitatively, the NRMSE of $\rho$ changes from 7.4% for the global parameter choice to the aforementioned 7.0% when using the novel fully local methodology. Although this may seem like a small reduction, the difference would probably become larger in the case of an even stronger spatial differentiation of the local particle distributions.