# Smoothed Particle Hydrodynamics for anisotropic dispersion in heterogeneous porous media

Rodrigo Pérez-Illanes, Guillem Sole-Mari, Daniel Fernàndez-Garcia

<figure><img src="image_1.png" /></figure>

## To cite this version:

Rodrigo Pérez-Illanes, Guillem Sole-Mari, Daniel Fernàndez-Garcia. Smoothed Particle Hydrodynamics for anisotropic dispersion in heterogeneous porous media. Advances in Water Resources, 2024, 183, pp.104601. 10.1016/j.advwatres.2023.104601. insu-04342000v2

HAL Id: insu-04342000

https://insu.hal.science/insu-04342000v2

Submitted on 11 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Smoothed Particle Hydrodynamics for anisotropic dispersion in heterogeneous porous media

Rodrigo Pérez-Illanes, Guillem Sole-Mari, Daniel Fernàndez-Garcia

► To cite this version:

Rodrigo Pérez-Illanes, Guillem Sole-Mari, Daniel Fernàndez-Garcia. Smoothed Particle Hydrodynamics for anisotropic dispersion in heterogeneous porous media. Advances in Water Resources, 2024, 183, pp.104601. 10.1016/j.advwatres.2023.104601. insu-04342000v2

HAL Id: insu-04342000

https://insu.hal.science/insu-04342000v2

Submitted on 11 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

# Smoothed Particle Hydrodynamics for anisotropic dispersion in heterogeneous porous media

Rodrigo Pérez-Illanes<sup>a,b,*</sup>, Guillem Sole-Mari<sup>c</sup>, Daniel Fernández-Garcia<sup>a,b</sup>

## ARTICLE INFO

**Keywords:**

SPH

Anisotropy

Dispersion

Heterogeneous

Particles

## ABSTRACT

In the context of modeling solute transport through heterogeneous porous media, particle methods possess inherent advantages with respect to mesh-based (Eulerian) methods. In Smoothed Particle Hydrodynamics (SPH), particles represent fluid volumes exchanging concentrations with their neighbors to emulate hydrodynamic dispersion, and advection is simulated by the particles' displacement. This crucially prevents problems that are otherwise typically associated with Eulerian advection schemes (especially at high grid-Péclet numbers), such as numerical diffusion. Despite the advantages of SPH, modeling dispersion with anisotropic coefficients remains a challenge for the approach, with studies reporting unphysical negative concentrations in conservative problems. This has likely hindered its practical use because dispersion is intrinsically anisotropic in porous media. This article provides a review and numerical evaluation of SPH for simulating dispersion, focusing on three formulations compatible with anisotropic dispersion coefficients. The analysis includes a scheme for which negative concentrations have been formerly reported, plus two more recently developed methods which are applied here for the first time to the problem of anisotropic dispersion in heterogeneous porous media. The SPH schemes are tested under different degrees of dispersion anisotropy for both homogeneous and heterogeneous velocity fields. The results indicate that the newer SPH schemes can produce accurate results without negative concentrations while considering anisotropic dispersion, providing a valid alternative to simulate solute transport through heterogeneous domains.

## 1. Introduction

Smoothed Particle Hydrodynamics (SPH) is a formulation for solving continuous conservation equations based on particles. The method was originally developed in the context of astrophysics (Lucy, 1977; Gingold and Monaghan, 1977) and remains a popular alternative for simulating problems with complex dynamics on a wide spectrum of scientific and engineering applications (e.g., Monaghan and Gingold, 1983; Cleary and Monaghan, 1999; Springel, 2010; Monaghan, 2012; Vacondio et al., 2012; Tartakovsky et al., 2015; Ye et al., 2019; Vacondio et al., 2021). Due to its Lagrangian nature, SPH presents inherent advantages for the simulation of solute transport through heterogeneous porous media with spatially variable flow velocities. These can be challenging conditions for classical grid solvers when advection is the locally dominant transport mechanism, yielding results that are strongly affected by numerical diffusion (Cirpka et al., 1999; Herrera et al., 2009; Boso et al., 2013; Benson et al., 2017). For reactive transport, an accurate estimation of concentrations and gradients is

needed because these quantities determine metrics of mixing and reaction rates (e.g., Kitanidis, 1994; Cirpka and Valocchi, 2007; De Simoni et al., 2007; Le Borgne et al., 2010; Fernández-Garcia and Sanchez-Vila, 2011). The SPH method is tightly connected to Kernel Density Estimation (KDE), naturally leading to continuous concentrations for each particle, which is attractive for the modeling of chemical reactions (e.g., Tartakovsky et al., 2007, 2009).

In the context of porous media, different applications of SPH have illustrated the versatility of the method for a variety of reactive and conservative transport processes. From a pore-scale perspective, the method has been applied to the simulation of multi-phase flow (Bardara et al., 2013; Kunz et al., 2015; Tartakovsky et al., 2015), the study of hydrodynamic dispersion through granular media (e.g., Zhu and Fox, 2001, 2002), reactive transport with mineral precipitation (e.g., Tartakovsky et al., 2007, 2008) and with kinetic models for biomass growth (e.g., Tartakovsky et al., 2009). From a macroscopic point of view, the method has been discussed for modeling variable-density

<footnote><sup>a</sup> Department of Civil and Environmental Engineering (DECA), Universitat Politècnica de Catalunya, Barcelona, Spain</footnote>

<footnote><sup>b</sup> Hydrogeology Group (UPC-CSIC), Barcelona, Spain</footnote>

<footnote><sup>c</sup> Géosciences Rennes, Université de Rennes, Rennes, France</footnote>

<footnote>* Corresponding author at: Department of Civil and Environmental Engineering (DECA), Universitat Politècnica de Catalunya, Barcelona, Spain. E-mail address: rodrigo.alfonso.perez@upc.edu (R. Pérez-Illanes).</footnote>
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

# Smoothed Particle Hydrodynamics for anisotropic dispersion in heterogeneous porous media

Rodrigo Pérez-Illanes<sup>a,b,*</sup>, Guillem Sole-Mari<sup>c</sup>, Daniel Fernàndez-Garcia<sup>a,b</sup>

## ARTICLE INFO

**Keywords:**

SPH

Anisotropy

Dispersion

Heterogeneous

Particles

## ABSTRACT

In the context of modeling solute transport through heterogeneous porous media, particle methods possess inherent advantages with respect to mesh-based (Eulerian) methods. In Smoothed Particle Hydrodynamics (SPH), particles represent fluid volumes exchanging concentrations with their neighbors to emulate hydrodynamic dispersion, and advection is simulated by the particles' displacement. This crucially prevents problems that are otherwise typically associated with Eulerian advection schemes (especially at high grid-Péclet numbers), such as numerical diffusion. Despite the advantages of SPH, modeling dispersion with anisotropic coefficients remains a challenge for the approach, with studies reporting unphysical negative concentrations in conservative problems. This has likely hindered its practical use because dispersion is intrinsically anisotropic in porous media. This article provides a review and numerical evaluation of SPH for simulating dispersion, focusing on three formulations compatible with anisotropic dispersion coefficients. The analysis includes a scheme for which negative concentrations have been formerly reported, plus two more recently developed methods which are applied here for the first time to the problem of anisotropic dispersion in heterogeneous porous media. The SPH schemes are tested under different degrees of dispersion anisotropy for both homogeneous and heterogeneous velocity fields. The results indicate that the newer SPH schemes can produce accurate results without negative concentrations while considering anisotropic dispersion, providing a valid alternative to simulate solute transport through heterogeneous domains.

## 1. Introduction

Smoothed Particle Hydrodynamics (SPH) is a formulation for solving continuous conservation equations based on particles. The method was originally developed in the context of astrophysics (Lucy, 1977; Gingold and Monaghan, 1977) and remains a popular alternative for simulating problems with complex dynamics on a wide spectrum of scientific and engineering applications (e.g., Monaghan and Gingold, 1983; Cleary and Monaghan, 1999; Springel, 2010; Monaghan, 2012; Vacondio et al., 2012; Tartakovsky et al., 2015; Ye et al., 2019; Vacondio et al., 2021). Due to its Lagrangian nature, SPH presents inherent advantages for the simulation of solute transport through heterogeneous porous media with spatially variable flow velocities. These can be challenging conditions for classical grid solvers when advection is the locally dominant transport mechanism, yielding results that are strongly affected by numerical diffusion (Cirpka et al., 1999; Herrera et al., 2009; Boso et al., 2013; Benson et al., 2017). For reactive transport, an accurate estimation of concentrations and gradients is

needed because these quantities determine metrics of mixing and reaction rates (e.g., Kitanidis, 1994; Cirpka and Valocchi, 2007; De Simoni et al., 2007; Le Borgne et al., 2010; Fernández-Garcia and Sanchez-Vila, 2011). The SPH method is tightly connected to Kernel Density Estimation (KDE), naturally leading to continuous concentrations for each particle, which is attractive for the modeling of chemical reactions (e.g., Tartakovsky et al., 2007, 2009).

In the context of porous media, different applications of SPH have illustrated the versatility of the method for a variety of reactive and conservative transport processes. From a pore-scale perspective, the method has been applied to the simulation of multi-phase flow (Bardara et al., 2013; Kunz et al., 2015; Tartakovsky et al., 2015), the study of hydrodynamic dispersion through granular media (e.g., Zhu and Fox, 2001, 2002), reactive transport with mineral precipitation (e.g., Tartakovsky et al., 2007, 2008) and with kinetic models for biomass growth (e.g., Tartakovsky et al., 2009). From a macroscopic point of view, the method has been discussed for modeling variable-density

<footnote><sup>a</sup> Department of Civil and Environmental Engineering (DECA), Universitat Politècnica de Catalunya, Barcelona, Spain</footnote>

<footnote><sup>b</sup> Hydrogeology Group (UPC-CSIC), Barcelona, Spain</footnote>

<footnote><sup>c</sup> Géosciences Rennes, Université de Rennes, Rennes, France</footnote>

<footnote>* Corresponding author at: Department of Civil and Environmental Engineering (DECA), Universitat Politècnica de Catalunya, Barcelona, Spain. E-mail address: rodrigo.alfonso.perez@upc.edu (R. Pérez-Illanes).</footnote>
flow (e.g., Basser et al., 2017, 2019), the transport of bacterial concentrations by chemotaxis (e.g., Avesani et al., 2016), and applied to the simulation of solute transport through heterogeneous domains (e.g., Herrera et al., 2009; Boso et al., 2013). In particular, these last examples considered isotropic solute dispersion coefficients, and remarked the ability of SPH to provide accurate results for complex transport conditions, without the influence of numerical artifacts. However, macroscopic hydrodynamic dispersion is intrinsically anisotropic in porous media, with transverse components (both in the horizontal and vertical directions) much smaller than the longitudinal (Gelhar et al., 1992; Lichtner et al., 2002). The degree of dispersion anisotropy should be expected to cover a wide range of magnitudes in complex geological formations. Studies have discussed the scale- and time-dependence of dispersion parameters, and the influence on transport of hydraulic aspects of the system such as the spatial variability of the medium properties or the relative alignment of the flow direction with the preferential structures of the solid matrix (e.g., Pickens and Grisak, 1981; Silliman and Simpson, 1987; Lichtner et al., 2002; Fernández-Garcia et al., 2005). Furthermore, transverse dispersion plays a critical role in mixing, which ultimately determines the outcome and extent of reactive transport processes (e.g., Cirpka et al., 1999; Klenk and Grathwohl, 2002; Fernández-Garcia et al., 2008; Rolle et al., 2009; Hochstetler et al., 2013), hence, an accurate characterization of dispersion anisotropy is needed. This context underlines that if the SPH method is to be considered for practical field-scale studies, it should be verified that formulations of dispersion are able to properly handle scenarios with anisotropic coefficients for a wide range of anisotropy ratios. It is at this point where the method has encountered some challenges that may have hindered its adoption. Particularly, simple transport simulations in homogeneous porous media have reported until recently the occurrence of unphysical negative concentrations when considering anisotropic dispersion coefficients (e.g., Herrera and Beckie, 2013; Avesani et al., 2015; Alvarado-Rodríguez et al., 2019; Klapp et al., 2022). Therefore, the discussion has remained constrained to these results and no applications have been reported in the literature where the method is employed to simulate solute transport with more realistic porous media configurations, simultaneously considering an explicit representation of the medium heterogeneity and of dispersion anisotropy.

This context motivates a review of the SPH method in order to evaluate its potential for modeling anisotropic dispersion in heterogeneous domains. In this regard, the following major objectives have been established for this work: (i) to understand the formulation and numerical properties of SPH methods for dispersive transport, (ii) to address the origin of negative concentrations observed in scenarios of anisotropic dispersion, and (iii) to apply the SPH method to a problem of solute transport through heterogeneous media considering anisotropic transport coefficients. After an initial evaluation, three different SPH schemes have been selected for discussion: the *Español and Revenga* method (ER; *Español and Revenga*, 2003), Anisotropic SPH for Anisotropic Diffusion (ASPHAD; Tran-Duc et al., 2016) and Two First Derivatives (TFD; Biriukov and Price, 2018). In subsurface applications, the ER scheme has been the common choice for simulations with anisotropic dispersion and it is known to yield negative concentrations (e.g., Herrera and Beckie, 2013; Alvarado-Rodríguez et al., 2019; Klapp et al., 2022). Briefly, the scheme can reverse the physical direction of dispersive fluxes in some sections within the kernel, and the solute is then transported unphysically from low to high concentration particles (see Tran-Duc et al., 2016; Biriukov and Price, 2018). This effect occurs for ratios of anisotropy exceeding a relatively low threshold determined by the number of spatial dimensions, hence it is very likely that the scheme will lead to instabilities for most practical cases. Two novel alternative formulations, compatible with anisotropic transport coefficients, have been recently proposed in the literature after identifying the unphysical transport effect, reportedly leading to results without negative concentrations. These formulations have not

been discussed yet for the transport of solutes through heterogeneous porous media, and this work evaluates them for the first time for said purpose. By means of a change of the coordinates system, the ASPHAD scheme (Tran-Duc et al., 2016) transforms the anisotropic dispersion equation into an equivalent isotropic form, which is then integrated following a corrected version of the isotropic SPH method (Brookshaw, 1985). A different approach is taken in the TFD method (Biriukov and Price, 2018), where a double SPH integration is performed: first to compute the spatial concentration gradients and then to obtain the divergence of dispersive fluxes. Furthermore, these authors provide a comprehensive physical explanation for the origin of the instabilities displayed by the ER scheme.

Two numerical experiments are conducted in order to evaluate the performance of the selected SPH schemes, both of them exploring different scenarios of dispersion anisotropy. Firstly, the anisotropic dispersion of a Gaussian plume in a homogeneous medium is revisited (as in Herrera and Beckie, 2013; Avesani et al., 2015; Alvarado-Rodríguez et al., 2019; Klapp et al., 2022) and compared with the available analytical result. In the second numerical experiment, a solute is transported through a randomly heterogeneous aquifer. The focus is placed on the SPH integration of dispersive transport and particle velocities are simply interpolated from a pre-calculated steady, divergence-free groundwater flow (as in Herrera et al., 2009; Boso et al., 2013). SPH results from this test are compared with those from a high-resolution Random Walk Particle Tracking (RWPT) model of the same problem, which is assumed to represent the exact solution. The RWPT method has a well-established trajectory in hydrogeology as a particle formulation capable of providing results not influenced by numerical dispersion in heterogeneous media (e.g., LaBolle et al., 1996; Lichtner et al., 2002; Salamon et al., 2006; Fernández-Garcia et al., 2008; Le Borgne et al., 2010; Boso et al., 2013).

The structure of this article is as follows. In *Methods*, a summary of the SPH formulation is presented, together with the considered SPH schemes and a discussion on numerical consistency. Then, the *Results and Discussion* of numerical tests outline the performance of each method and their applicability to the problem of anisotropic dispersion for a variety of conditions likely to be found in aquifer models. Finally, the *Conclusions* of this work underline the advantages and challenges of the discussed SPH schemes for anisotropic transport through heterogeneous domains.

## 2. Methods

### 2.1. Advection Dispersion Equation (ADE)

Solute transport is represented by the Advection Dispersion Equation (ADE) with spatially varying aquifer properties described at the Darcy scale. The transport of an ideal non-reactive solute is assumed, meaning that concentration influences neither fluid density nor viscosity, without chemical reactions involved (Bear and Cheng, 2010). In Lagrangian form, the ADE under incompressible groundwater flow can be written as

$$ \frac{d\mathbf{r}}{dt} = \frac{\mathbf{q}}{\phi}, \qquad (1) $$

$$ \frac{d(\phi C)}{dt} = \nabla \cdot (\phi D \nabla C), \qquad (2) $$

where $d/dt(\cdot)$ is the total time derivative, $\mathbf{r}$ is the spatial coordinate vector for a particle's position, $\mathbf{q} = \phi \mathbf{v}$ is the Darcy velocity, $\phi$ is the medium porosity, $\mathbf{v}$ is the macroscopic flow velocity and $C$ is the solute concentration expressed as mass per unit water volume. The hydrodynamic dispersion tensor is given by (Bear and Cheng, 2010)

$$ D_{ij} = (\alpha_T |\mathbf{v}| + D_m)\delta_{ij} + (\alpha_L - \alpha_T) \frac{v_i v_j}{|\mathbf{v}|}, \qquad (3) $$
flow (e.g., Basser et al., 2017, 2019), the transport of bacterial concentrations by chemotaxis (e.g., Avesani et al., 2016), and applied to the simulation of solute transport through heterogeneous domains (e.g., Herrera et al., 2009; Boso et al., 2013). In particular, these last examples considered isotropic solute dispersion coefficients, and remarked the ability of SPH to provide accurate results for complex transport conditions, without the influence of numerical artifacts. However, macroscopic hydrodynamic dispersion is intrinsically anisotropic in porous media, with transverse components (both in the horizontal and vertical directions) much smaller than the longitudinal (Gelhar et al., 1992; Lichtner et al., 2002). The degree of dispersion anisotropy should be expected to cover a wide range of magnitudes in complex geological formations. Studies have discussed the scale- and time-dependence of dispersion parameters, and the influence on transport of hydraulic aspects of the system such as the spatial variability of the medium properties or the relative alignment of the flow direction with the preferential structures of the solid matrix (e.g., Pickens and Grisak, 1981; Silliman and Simpson, 1987; Lichtner et al., 2002; Fernández-Garcia et al., 2005). Furthermore, transverse dispersion plays a critical role in mixing, which ultimately determines the outcome and extent of reactive transport processes (e.g., Cirpka et al., 1999; Klenk and Grathwohl, 2002; Fernández-Garcia et al., 2008; Rolle et al., 2009; Hochstetler et al., 2013), hence, an accurate characterization of dispersion anisotropy is needed. This context underlines that if the SPH method is to be considered for practical field-scale studies, it should be verified that formulations of dispersion are able to properly handle scenarios with anisotropic coefficients for a wide range of anisotropy ratios. It is at this point where the method has encountered some challenges that may have hindered its adoption. Particularly, simple transport simulations in homogeneous porous media have reported until recently the occurrence of unphysical negative concentrations when considering anisotropic dispersion coefficients (e.g., Herrera and Beckie, 2013; Avesani et al., 2015; Alvarado-Rodríguez et al., 2019; Klapp et al., 2022). Therefore, the discussion has remained constrained to these results and no applications have been reported in the literature where the method is employed to simulate solute transport with more realistic porous media configurations, simultaneously considering an explicit representation of the medium heterogeneity and of dispersion anisotropy.

This context motivates a review of the SPH method in order to evaluate its potential for modeling anisotropic dispersion in heterogeneous domains. In this regard, the following major objectives have been established for this work: (i) to understand the formulation and numerical properties of SPH methods for dispersive transport, (ii) to address the origin of negative concentrations observed in scenarios of anisotropic dispersion, and (iii) to apply the SPH method to a problem of solute transport through heterogeneous media considering anisotropic transport coefficients. After an initial evaluation, three different SPH schemes have been selected for discussion: the *Español and Revenga* method (ER; *Español and Revenga*, 2003), Anisotropic SPH for Anisotropic Diffusion (ASPHAD; Tran-Duc et al., 2016) and Two First Derivatives (TFD; Biriukov and Price, 2018). In subsurface applications, the ER scheme has been the common choice for simulations with anisotropic dispersion and it is known to yield negative concentrations (e.g., Herrera and Beckie, 2013; Alvarado-Rodríguez et al., 2019; Klapp et al., 2022). Briefly, the scheme can reverse the physical direction of dispersive fluxes in some sections within the kernel, and the solute is then transported unphysically from low to high concentration particles (see Tran-Duc et al., 2016; Biriukov and Price, 2018). This effect occurs for ratios of anisotropy exceeding a relatively low threshold determined by the number of spatial dimensions, hence it is very likely that the scheme will lead to instabilities for most practical cases. Two novel alternative formulations, compatible with anisotropic transport coefficients, have been recently proposed in the literature after identifying the unphysical transport effect, reportedly leading to results without negative concentrations. These formulations have not

been discussed yet for the transport of solutes through heterogeneous porous media, and this work evaluates them for the first time for said purpose. By means of a change of the coordinates system, the ASPHAD scheme (Tran-Duc et al., 2016) transforms the anisotropic dispersion equation into an equivalent isotropic form, which is then integrated following a corrected version of the isotropic SPH method (Brookshaw, 1985). A different approach is taken in the TFD method (Biriukov and Price, 2018), where a double SPH integration is performed: first to compute the spatial concentration gradients and then to obtain the divergence of dispersive fluxes. Furthermore, these authors provide a comprehensive physical explanation for the origin of the instabilities displayed by the ER scheme.

Two numerical experiments are conducted in order to evaluate the performance of the selected SPH schemes, both of them exploring different scenarios of dispersion anisotropy. Firstly, the anisotropic dispersion of a Gaussian plume in a homogeneous medium is revisited (as in Herrera and Beckie, 2013; Avesani et al., 2015; Alvarado-Rodríguez et al., 2019; Klapp et al., 2022) and compared with the available analytical result. In the second numerical experiment, a solute is transported through a randomly heterogeneous aquifer. The focus is placed on the SPH integration of dispersive transport and particle velocities are simply interpolated from a pre-calculated steady, divergence-free groundwater flow (as in Herrera et al., 2009; Boso et al., 2013). SPH results from this test are compared with those from a high-resolution Random Walk Particle Tracking (RWPT) model of the same problem, which is assumed to represent the exact solution. The RWPT method has a well-established trajectory in hydrogeology as a particle formulation capable of providing results not influenced by numerical dispersion in heterogeneous media (e.g., LaBolle et al., 1996; Lichtner et al., 2002; Salamon et al., 2006; Fernández-Garcia et al., 2008; Le Borgne et al., 2010; Boso et al., 2013).

The structure of this article is as follows. In *Methods*, a summary of the SPH formulation is presented, together with the considered SPH schemes and a discussion on numerical consistency. Then, the *Results and Discussion* of numerical tests outline the performance of each method and their applicability to the problem of anisotropic dispersion for a variety of conditions likely to be found in aquifer models. Finally, the *Conclusions* of this work underline the advantages and challenges of the discussed SPH schemes for anisotropic transport through heterogeneous domains.

## 2. Methods

### 2.1. Advection Dispersion Equation (ADE)

Solute transport is represented by the Advection Dispersion Equation (ADE) with spatially varying aquifer properties described at the Darcy scale. The transport of an ideal non-reactive solute is assumed, meaning that concentration influences neither fluid density nor viscosity, without chemical reactions involved (Bear and Cheng, 2010). In Lagrangian form, the ADE under incompressible groundwater flow can be written as

$$ \frac{d\mathbf{r}}{dt} = \frac{\mathbf{q}}{\phi}, \qquad (1) $$

$$ \frac{d(\phi C)}{dt} = \nabla \cdot (\phi \mathbf{D}\nabla C), \qquad (2) $$

where $d/dt(\cdot)$ is the total time derivative, $\mathbf{r}$ is the spatial coordinate vector for a particle's position, $\mathbf{q} = \phi \mathbf{v}$ is the Darcy velocity, $\phi$ is the medium porosity, $\mathbf{v}$ is the macroscopic flow velocity and $C$ is the solute concentration expressed as mass per unit water volume. The hydrodynamic dispersion tensor is given by (Bear and Cheng, 2010)

$$ D_{ij} = (\alpha_T |\mathbf{v}| + D_m)\delta_{ij} + (\alpha_L - \alpha_T) \frac{v_i v_j}{|\mathbf{v}|}, \qquad (3) $$
where $\alpha_L$ and $\alpha_T$ are the longitudinal and transverse dispersivities, respectively, $D_m$ is the effective molecular diffusion (corrected for toruosity effects) and $\delta_{ij}$ is the Kronecker delta. For simplicity, transport is assumed in two dimensions for a single solute, without molecular diffusion and with unit porosity. In case of a spatially variable porosity, the effect should be accounted for while computing SPH quantities (e.g., Basser et al., 2019). Dispersion anisotropy is defined here as $\lambda_D = D_{xx}/D_{yy} \approx \alpha_L/\alpha_T$. In typical aquifer applications, $\lambda_D$ ranges between 3 and 20 for horizontal planes and could be even higher when considering the vertical transverse dispersivity (e.g., Gelhar et al., 1992).

## 2.2. Smoothed Particle Hydrodynamics (SPH)

The SPH method discretizes the fluid as particles, each one carrying physical properties such as mass, density and solute concentration. The formulation begins by expressing an arbitrary continuous field $A(\mathbf{r})$ as (Monaghan, 2005; Violeau, 2012)

$$A(\mathbf{r}) = \int_{\Omega} A(\mathbf{r}') \delta(\mathbf{r} - \mathbf{r}') d\mathbf{r}', \quad (4)$$

where $\Omega$ is the domain volume and $\delta$ is the Dirac delta function. The Dirac delta is then approximated by a positive continuous symmetric kernel function $W(|\mathbf{r} - \mathbf{r}'|, h)$, with $h$ being a smoothing length. Together with the number of particles, the smoothing length parameter controls the accuracy and spatial resolution of the approximation. In practical applications, the kernel function has a finite support volume characterized by a distance $H$, such that kernel density equals zero at any distance higher than $H$. Therefore, only a finite set of particles contribute to the interpolation at a given point. Besides the positivity and symmetry properties, the kernel function has a zeroth moment equal to one, a first order moment equal to zero (due to symmetry), and a non-zero second order moment. The latter is the source of the continuous interpolation error $O(h^2)$ (Violeau, 2012; Sigalotti et al., 2019). Continuous integrals are discretized on the set of particles, leading to (Monaghan, 2005)

$$A(\mathbf{r}_a) = \sum_{b}^{\mathcal{N}} \frac{m_b}{\rho_b} A_b W(|\mathbf{r}_a - \mathbf{r}_b|, h_a), \quad (5)$$

where subscripts a and b denote the center and neighbor particles respectively. Sums are performed over a finite set of $\mathcal{N}$ neighbors within the support volume. For the remainder of the article, SPH sums are assumed to always occur over the set of neighbors and $\mathcal{N}$ is not explicitly indicated. $m_b$ and $\rho_b$ represent the mass and density of particle b, so their quotient is a measure of the volume associated with that particle. Discretizing with particles introduces a second source of error that scales as $O(\mathcal{N}^{-\psi})$, where the exponent $\psi$ depends on the particle distribution and quantifies the rate of decrease in statistical error for an increasing number of neighbors (Zhu et al., 2015). For a uniform particle distribution, this error decreases faster ($\psi \approx 1$) than in the case of disordered particles ($\psi \approx 0.5$). Density is computed as (Monaghan, 2005)

$$\rho_a = \rho(\mathbf{r}_a) = \sum_b m_b W_{ab}(h_a), \quad (6)$$

where $W_{ab}(h_a)$ denotes the kernel function evaluated at distance $r_{ab} = |\mathbf{r}_a - \mathbf{r}_b|$ with a smoothing length $h_a$. In this paper, the particle mass $m_b$ is always assumed to be the same for all particles.

### 2.2.1. Kernel

The kernel is usually a piecewise continuous function (Cleary and Monaghan, 1999; Herrera and Beckie, 2013), being a common choice the family of B-Splines (Schoenberg, 1946). Smoothing functions have an associated residual error called the kernel bias, which is the fraction of error that does not decrease further after increasing the number of neighbor particles. The magnitude of the bias decreases with increasing order of the B-Spline function, thus higher order kernels allow

<figure><img src="image_5.png" /><figcaption>Fig. 1. Fifth order B-Spline kernel and characteristic distances. (a) Normalized shape function $w(q)$, and (b) normalized shape derivative $w'(q)$. $\sigma_W$ is the standard deviation of the kernel's shape.</figcaption></figure>

using a higher number of particles to improve accuracy (Biriukov and Price, 2018). This study employs the normalized kernel formulation by Dehnen and Aly (2012), defined as

$$W(\tilde{\mathbf{r}}, h) = \frac{c_v}{(\kappa_v h)^v} w(q), \quad (7)$$

where $\tilde{\mathbf{r}} = \mathbf{r} - \mathbf{r}'$ is the vector difference between the center and neighbor particle positions and $w(q)$ is a shape function normalized to a unit support radius, with $q = |\tilde{\mathbf{r}}|/(\kappa_v h)$ the normalized kernel coordinate. $c_v$ is a positive constant which makes the kernel integral unitary and whose value depends on the number of spatial dimensions $v$; and $\kappa_v = H/h$ is an aspect ratio specific to the kernel function. From Eq. (7) follows that kernel gradients are expressed as (Violeau, 2012)

$$\nabla W(\tilde{\mathbf{r}}, h) = F(q)\hat{\mathbf{r}}, \quad (8)$$

where $F(q)$ is defined as the kernel scalar derivative

$$F(q) = \frac{c_v}{(\kappa_v h)^{v+1}} w'(q), \quad (9)$$

which is negative for typical kernel functions (Monaghan, 2005), and $\hat{\mathbf{r}} = \tilde{\mathbf{r}}/|\tilde{\mathbf{r}}|$ is the unit vector between particles. In this article, a fifth order B-Spline kernel is used (Fig. 1)

$$w(q) = (1 - q)^5 + 6 \left(\frac{2}{3} - q\right)^5 + 15 \left(\frac{1}{3} - q\right)^5, \quad (10)$$

where the operator $x_+$ stands for $\max(0, x)$. In two dimensions, this kernel's constants have values $c_2 = 3^7 \cdot 7/478\pi$ and $\kappa_2 = 2.158131$ (Dehnen and Aly, 2012).

### 2.2.2. Errors and consistency

Both continuous and discrete interpolation errors are relevant for the accuracy of SPH estimates. Zhu et al. (2015) discussed this topic
where $\alpha_L$ and $\alpha_T$ are the longitudinal and transverse dispersivities, respectively, $D_m$ is the effective molecular diffusion (corrected for toruosity effects) and $\delta_{ij}$ is the Kronecker delta. For simplicity, transport is assumed in two dimensions for a single solute, without molecular diffusion and with unit porosity. In case of a spatially variable porosity, the effect should be accounted for while computing SPH quantities (e.g., Basser et al., 2019). Dispersion anisotropy is defined here as $\lambda_D = D_{xx}/D_{yy} \approx \alpha_L/\alpha_T$. In typical aquifer applications, $\lambda_D$ ranges between 3 and 20 for horizontal planes and could be even higher when considering the vertical transverse dispersivity (e.g., Gelhar et al., 1992).

## 2.2. Smoothed Particle Hydrodynamics (SPH)

The SPH method discretizes the fluid as particles, each one carrying physical properties such as mass, density and solute concentration. The formulation begins by expressing an arbitrary continuous field $A(\mathbf{r})$ as (Monaghan, 2005; Violeau, 2012)

$$A(\mathbf{r}) = \int_{\Omega} A(\mathbf{r}') \delta(\mathbf{r} - \mathbf{r}') d\mathbf{r}', \quad (4)$$

where $\Omega$ is the domain volume and $\delta$ is the Dirac delta function. The Dirac delta is then approximated by a positive continuous symmetric kernel function $W(|\mathbf{r} - \mathbf{r}'|, h)$, with $h$ being a smoothing length. Together with the number of particles, the smoothing length parameter controls the accuracy and spatial resolution of the approximation. In practical applications, the kernel function has a finite support volume characterized by a distance $H$, such that kernel density equals zero at any distance higher than $H$. Therefore, only a finite set of particles contribute to the interpolation at a given point. Besides the positivity and symmetry properties, the kernel function has a zeroth moment equal to one, a first order moment equal to zero (due to symmetry), and a non-zero second order moment. The latter is the source of the continuous interpolation error $O(h^2)$ (Violeau, 2012; Sigalotti et al., 2019). Continuous integrals are discretized on the set of particles, leading to (Monaghan, 2005)

$$A(\mathbf{r}_a) = \sum_{b}^{\mathcal{N}} \frac{m_b}{\rho_b} A_b W(|\mathbf{r}_a - \mathbf{r}_b|, h_a), \quad (5)$$

where subscripts a and b denote the center and neighbor particles respectively. Sums are performed over a finite set of $\mathcal{N}$ neighbors within the support volume. For the remainder of the article, SPH sums are assumed to always occur over the set of neighbors and $\mathcal{N}$ is not explicitly indicated. $m_b$ and $\rho_b$ represent the mass and density of particle b, so their quotient is a measure of the volume associated with that particle. Discretizing with particles introduces a second source of error that scales as $O(\mathcal{N}^{-\psi})$, where the exponent $\psi$ depends on the particle distribution and quantifies the rate of decrease in statistical error for an increasing number of neighbors (Zhu et al., 2015). For a uniform particle distribution, this error decreases faster ($\psi \approx 1$) than in the case of disordered particles ($\psi \approx 0.5$). Density is computed as (Monaghan, 2005)

$$\rho_a = \rho(\mathbf{r}_a) = \sum_b m_b W_{ab}(h_a), \quad (6)$$

where $W_{ab}(h_a)$ denotes the kernel function evaluated at distance $r_{ab} = |\mathbf{r}_a - \mathbf{r}_b|$ with a smoothing length $h_a$. In this paper, the particle mass $m_b$ is always assumed to be the same for all particles.

### 2.2.1. Kernel

The kernel is usually a piecewise continuous function (Cleary and Monaghan, 1999; Herrera and Beckie, 2013), being a common choice the family of B-Splines (Schoenberg, 1946). Smoothing functions have an associated residual error called the kernel bias, which is the fraction of error that does not decrease further after increasing the number of neighbor particles. The magnitude of the bias decreases with increasing order of the B-Spline function, thus higher order kernels allow

<figure><img src="image_4.png" /><figcaption>Fig. 1. Fifth order B-Spline kernel and characteristic distances. (a) Normalized shape function $w(q)$, and (b) normalized shape derivative $w'(q)$. $\sigma_W$ is the standard deviation of the kernel's shape.</figcaption></figure>

using a higher number of particles to improve accuracy (Biriukov and Price, 2018). This study employs the normalized kernel formulation by Dehnen and Aly (2012), defined as

$$W(\tilde{\mathbf{r}}, h) = \frac{c_v}{(\kappa_v h)^v} w(q), \quad (7)$$

where $\tilde{\mathbf{r}} = \mathbf{r} - \mathbf{r}'$ is the vector difference between the center and neighbor particle positions and $w(q)$ is a shape function normalized to a unit support radius, with $q = |\tilde{\mathbf{r}}|/(\kappa_v h)$ the normalized kernel coordinate. $c_v$ is a positive constant which makes the kernel integral unitary and whose value depends on the number of spatial dimensions $v$; and $\kappa_v = H/h$ is an aspect ratio specific to the kernel function. From Eq. (7) follows that kernel gradients are expressed as (Violeau, 2012)

$$\nabla W(\tilde{\mathbf{r}}, h) = F(q)\hat{\mathbf{r}}, \quad (8)$$

where $F(q)$ is defined as the kernel scalar derivative

$$F(q) = \frac{c_v}{(\kappa_v h)^{v+1}} w'(q), \quad (9)$$

which is negative for typical kernel functions (Monaghan, 2005), and $\hat{\mathbf{r}} = \tilde{\mathbf{r}}/|\tilde{\mathbf{r}}|$ is the unit vector between particles. In this article, a fifth order B-Spline kernel is used (Fig. 1)

$$w(q) = (1 - q)^5 + 6 \left( \frac{2}{3} - q \right)^5 + 15 \left( \frac{1}{3} - q \right)^5, \quad (10)$$

where the operator $x_+$ stands for $\max(0, x)$. In two dimensions, this kernel's constants have values $c_2 = 3^7 \cdot 7/478\pi$ and $\kappa_2 = 2.158131$ (Dehnen and Aly, 2012).

### 2.2.2. Errors and consistency

Both continuous and discrete interpolation errors are relevant for the accuracy of SPH estimates. Zhu et al. (2015) discussed this topic
considering density estimates. The authors argue that, as the total number of particles tends to infinity, $N \to \infty$, it should be satisfied simultaneously that $h \to 0$ and $N \to \infty$. Studies following said principles have shown that numerical error decreases while increasing the particle resolution (e.g., Sigalotti et al., 2016; Alvarado-Rodríguez et al., 2019). Still, it has been a common practice in SPH to set the smoothing length as proportional to the mean particle size $\delta_r = (m/\rho)^{1/\nu}$, that is, $h = \gamma\delta_r$. The parameter $\gamma$ controls the number of neighbors within the kernel, with a typically adopted value of $\gamma = 1.2$ (Cleary and Monaghan, 1999; Herrera and Beckie, 2013; Tran-Duc et al., 2016). It is worth noting that in order to achieve consistency (i.e., for the total interpolation error to approach zero as one adds more particles), $\gamma$ should actually increase as the total number of particles increases. Otherwise, the discrete error at the kernel scale remains constant and it is replicated more times as the total number of particles increases (Zhu et al., 2015). From a modeling perspective, definition of the smoothing length using a constant $\gamma$ is convenient due to its simplicity, but lacks theoretical rigor as $\gamma$ is not a function of the number of particles. In Appendix A it is shown that, by introducing consistency principles, the following expression for $\gamma$ as a function of the total number of SPH particles can be obtained

$$ \gamma(N) = \frac{1}{\kappa_v} \left( \frac{\mathcal{K}}{V_v} \right)^{1/\nu} N^{1/(\nu-1/\beta)}, \quad (11) $$

where $V_v$ is the volume of the unit sphere, $\mathcal{K}$ is a proportionality constant, and $\beta$ is a positive exponent relating the scaling of the smoothing length with the number of particles $h \propto N^{-1/\beta}$. For two dimensional problems, $\beta \ge 4$ has been established as a reference for a uniform particle distribution, ensuring the compliance of the triple limit for consistent scaling (for details refer to Appendix A).

## 2.3. SPH methods for anisotropic dispersion

### 2.3.1. Español & Revenga (ER)

This method has been the common choice in studies assessing anisotropic dispersion with SPH (Herrera and Beckie, 2013; Avesani et al., 2015; Alvarado-Rodríguez et al., 2019; Klapp et al., 2022). Negative concentrations have been reported while simulating the anisotropic dispersion of a Gaussian plume in a homogeneous domain. The scheme integrates dispersion as (Español and Revenga, 2003; Biriukov and Price, 2018)

$$ \frac{dC_a}{dt} = \sum_b \frac{m_b}{\bar{\rho}_{ab}} C_{ba} \bar{D}_{ij}^{ab} G_{ab}^{ij} \quad (12) $$

where $C_{ba} = C_b - C_a$, and $\bar{\rho}_{ab}$ is an averaged density. $\bar{D}_{ij}^{ab}$ is the $i, j$ component of the dispersion tensor averaged (typically harmonic) between particles $a$ and $b$. $G_{ab}^{ij}$ is a term describing the kernel gradient defined as

$$ G_{ab}^{ij} = - \left[ (\nu + 2) \frac{\hat{r}_{ab}^i \hat{r}_{ab}^j}{\delta_{ij}} - \frac{F_{ab}}{r_{ab}} \right] \quad (13) $$

where $\hat{r}_{ab}^i = r_{ab}^i / r_{ab}$ is the ith-component of the vector $\hat{r}_{ab}$ and $F_{ab} = F(q_{ab})$, with $q_{ab}$ the normalized kernel coordinate distance $q$ between particles $a$ and $b$. In Eq. (12), indexes $i, j$ follow the Einstein summation convention. The problem with this method lies in the term within squared brackets in Eq. (13), which can adopt negative values in some locations within the kernel, depending on the components of the unit vector $\hat{r}_{ab}$. This modifies the physical direction of dispersive flux which should be against the concentration gradient (Tran-Duc et al., 2016; Biriukov and Price, 2018). That is, assuming only one pair of particles $a, b$ for simplicity, for $C_a < C_b$ and under proper physical conditions, it is expected that $dC_a/dt > 0$, which is satisfied as long as the product between dispersion and the term $G_{ab}^{ij}$ is positive. In two dimensions, and assuming a diagonal dispersion tensor, such condition translates into the inequality

$$ \frac{D_{xx}}{D_{yy}} \left[ 1 - (\nu + 2)(\hat{r}_{ab}^x)^2 \right] \le (\nu + 2)(1 - (\hat{r}_{ab}^x)^2) - 1, \quad (14) $$

where it has been applied that $(\hat{r}_{ab}^x)^2 + (\hat{r}_{ab}^y)^2 = 1$. Limit cases can be obtained by considering the limit values of $(\hat{r}_{ab}^x)^2 \in \{0, 1\}$, which yields that the scheme preserves physical dispersive flux when dispersion anisotropy satisfies $\lambda_D \in [1/(\nu + 1), \nu + 1]$ (Tran-Duc et al., 2016; Biriukov and Price, 2018). In reported applications of the ER method for dispersion through porous media, the minimum (non-isotropic) considered degree of anisotropy is $\lambda_D = 10$, which falls outside the valid range for two dimensions. One question surrounding the results obtained with this method is related to the spatial distribution of negative concentrations. In applications to homogeneous domains, bands of negative values appear to be typically confined to regions far from the main solute plume. However, no reports are available about the impact that this artifact may have for smaller scales of concentration and gradients, particularly in heterogeneous domains.

### 2.3.2. Anisotropic SPH for anisotropic diffusion (ASPHAD)

This method by Tran-Duc et al. (2016) stems from a change of the coordinates system in the anisotropic dispersion equation, which allows rewriting the problem in isotropic form. For spatially heterogeneous dispersion, this transformation involves the assumption that differences in dispersion between neighboring particles are relatively small. This allows taking the dispersion tensor outside the divergence of fluxes in Eq. (2). A change of coordinates is performed, based on a transformation matrix defined from the dispersion coefficients. Briefly, the dispersion tensor is a positive definite symmetric matrix (Bear and Cheng, 2010) and can be decomposed as $D = BB^T$ (LaBolle et al., 1996), where the superscript $T$ stands for transpose. The matrix $B$ is well known in the context of Random Walk Particle Tracking (RWPT), where it is denoted as the displacement matrix. In two dimensions it is given by (Fernández-García et al., 2005; Salamon et al., 2006)

$$ B = \begin{bmatrix} \frac{v_x}{|v|} \sqrt{\alpha_L |v| + D_m} & - \frac{v_y}{|v|} \sqrt{\alpha_T |v| + D_m} \\ \frac{v_y}{|v|} \sqrt{\alpha_L |v| + D_m} & \frac{v_x}{|v|} \sqrt{\alpha_T |v| + D_m} \end{bmatrix}. \quad (15) $$

Furthermore, the inverse of this matrix can be computed analytically as

$$ B^{-1} = \begin{bmatrix} \frac{1}{\sqrt{\alpha_L |v| + D_m}} & \frac{v_x}{|v|} \\ -\frac{1}{\sqrt{\alpha_T |v| + D_m}} & \frac{v_y}{|v|} \end{bmatrix} \quad (16) $$

The coordinate transformation in ASPHAD is such that $X = B^{-1}X$, where $X$ is the position in the new frame of reference. This leads to the following SPH scheme for dispersion (Tran-Duc et al., 2016),

$$ \frac{dC_a}{dt} = 2 \sum_b \frac{m_b}{\rho_b} \frac{C_{ab}}{|\bar{B}_{ab}^{-1} \hat{r}_{ab}|^2} \frac{F_{ab}}{r_{ab}}, \quad (17) $$

where $\bar{B}_{ab}^{-1}$ is the inverse displacement matrix (Eq. (16)) averaged between particles $a$ and $b$, and $C_{ab} = C_a - C_b$. Tran-Duc et al. (2016) showed that this scheme may artificially increase dispersion in the direction of smaller dispersivity (transverse), and reduce it in the direction of larger dispersivity (longitudinal). However, a major advantage of this approach is that it is unconditionally stable (Brookshaw, 1985; Biriukov and Price, 2018).

### 2.3.3. Two first derivatives (TFD)

This method is presented by Biriukov and Price (2018) with application to heat transport, showing satisfactory results for an anisotropic conduction tensor. It differs from the previous methods in that the dispersive flux is first explicitly calculated and then used to compute the divergence of fluxes. It is thus a two-stage SPH integration. The jth-component of the gradient of solute concentrations is expressed as

$$ G_a^j = \nabla_a^j C_a = \frac{1}{\rho_a f_a} \sum_b m_b C_{ba} \nabla_a^j W_{ab}(h_a), \quad (18) $$

where

$$ f_a = 1 + \frac{h_a}{\nu \rho_a} \sum_b m_b \frac{\partial W_{ab}(h_a)}{\partial h_a} \quad (19) $$
considering density estimates. The authors argue that, as the total number of particles tends to infinity, $N \to \infty$, it should be satisfied simultaneously that $h \to 0$ and $N \to \infty$. Studies following said principles have shown that numerical error decreases while increasing the particle resolution (e.g., Sigalotti et al., 2016; Alvarado-Rodríguez et al., 2019). Still, it has been a common practice in SPH to set the smoothing length as proportional to the mean particle size $\delta_r = (m/\rho)^{1/\nu}$, that is, $h = \gamma\delta_r$. The parameter $\gamma$ controls the number of neighbors within the kernel, with a typically adopted value of $\gamma = 1.2$ (Cleary and Monaghan, 1999; Herrera and Beckie, 2013; Tran-Duc et al., 2016). It is worth noting that in order to achieve consistency (i.e., for the total interpolation error to approach zero as one adds more particles), $\gamma$ should actually increase as the total number of particles increases. Otherwise, the discrete error at the kernel scale remains constant and it is replicated more times as the total number of particles increases (Zhu et al., 2015). From a modeling perspective, definition of the smoothing length using a constant $\gamma$ is convenient due to its simplicity, but lacks theoretical rigor as $\gamma$ is not a function of the number of particles. In Appendix A it is shown that, by introducing consistency principles, the following expression for $\gamma$ as a function of the total number of SPH particles can be obtained

$$ \gamma(N) = \frac{1}{\kappa_v} \left( \frac{\mathcal{K}}{V_v} \right)^{1/\nu} N^{1/(\nu-1)/\beta}, \quad (11) $$

where $V_v$ is the volume of the unit sphere, $\mathcal{K}$ is a proportionality constant, and $\beta$ is a positive exponent relating the scaling of the smoothing length with the number of particles $h \propto N^{-1/\beta}$. For two dimensional problems, $\beta \ge 4$ has been established as a reference for a uniform particle distribution, ensuring the compliance of the triple limit for consistent scaling (for details refer to Appendix A).

## 2.3. SPH methods for anisotropic dispersion

### 2.3.1. Español & Revenga (ER)

This method has been the common choice in studies assessing anisotropic dispersion with SPH (Herrera and Beckie, 2013; Avesani et al., 2015; Alvarado-Rodríguez et al., 2019; Klapp et al., 2022). Negative concentrations have been reported while simulating the anisotropic dispersion of a Gaussian plume in a homogeneous domain. The scheme integrates dispersion as (Español and Revenga, 2003; Biriukov and Price, 2018)

$$ \frac{dC_a}{dt} = \sum_b \frac{m_b}{\bar{\rho}_{ab}} C_{ba} \bar{D}_{ij}^{ab} G_{ab}^{ij} \quad (12) $$

where $C_{ba} = C_b - C_a$, and $\bar{\rho}_{ab}$ is an averaged density. $\bar{D}_{ij}^{ab}$ is the $i, j$ component of the dispersion tensor averaged (typically harmonic) between particles $a$ and $b$. $G_{ab}^{ij}$ is a term describing the kernel gradient defined as

$$ G_{ab}^{ij} = - \left[ (\nu + 2) \frac{\hat{r}_{ab}^i \hat{r}_{ab}^j}{\delta_{ij}} - \frac{F_{ab}}{r_{ab}} \right], \quad (13) $$

where $\hat{r}_{ab}^i = r_{ab}^i / r_{ab}$ is the ith-component of the vector $\hat{r}_{ab}$ and $F_{ab} = F(q_{ab})$, with $q_{ab}$ the normalized kernel coordinate distance $q$ between particles $a$ and $b$. In Eq. (12), indexes $i, j$ follow the Einstein summation convention. The problem with this method lies in the term within squared brackets in Eq. (13), which can adopt negative values in some locations within the kernel, depending on the components of the unit vector $\hat{r}_{ab}$. This modifies the physical direction of dispersive flux which should be against the concentration gradient (Tran-Duc et al., 2016; Biriukov and Price, 2018). That is, assuming only one pair of particles $a, b$ for simplicity, for $C_a < C_b$ and under proper physical conditions, it is expected that $dC_a/dt > 0$, which is satisfied as long as the product between dispersion and the term $G_{ab}^{ij}$ is positive. In two dimensions, and assuming a diagonal dispersion tensor, such condition translates into the inequality

$$ \frac{D_{xx}}{D_{yy}} \left[ 1 - (\nu + 2)(\hat{r}_{ab}^x)^2 \right] \le (\nu + 2)(1 - (\hat{r}_{ab}^x)^2) - 1, \quad (14) $$

where it has been applied that $(\hat{r}_{ab}^x)^2 + (\hat{r}_{ab}^y)^2 = 1$. Limit cases can be obtained by considering the limit values of $(\hat{r}_{ab}^x)^2 \in \{0, 1\}$, which yields that the scheme preserves physical dispersive flux when dispersion anisotropy satisfies $\lambda_D \in [1/(\nu + 1), \nu + 1]$ (Tran-Duc et al., 2016; Biriukov and Price, 2018). In reported applications of the ER method for dispersion through porous media, the minimum (non-isotropic) considered degree of anisotropy is $\lambda_D = 10$, which falls outside the valid range for two dimensions. One question surrounding the results obtained with this method is related to the spatial distribution of negative concentrations. In applications to homogeneous domains, bands of negative values appear to be typically confined to regions far from the main solute plume. However, no reports are available about the impact that this artifact may have for smaller scales of concentration and gradients, particularly in heterogeneous domains.

### 2.3.2. Anisotropic SPH for anisotropic diffusion (ASPHAD)

This method by Tran-Duc et al. (2016) stems from a change of the coordinates system in the anisotropic dispersion equation, which allows rewriting the problem in isotropic form. For spatially heterogeneous dispersion, this transformation involves the assumption that differences in dispersion between neighboring particles are relatively small. This allows taking the dispersion tensor outside the divergence of fluxes in Eq. (2). A change of coordinates is performed, based on a transformation matrix defined from the dispersion coefficients. Briefly, the dispersion tensor is a positive definite symmetric matrix (Bear and Cheng, 2010) and can be decomposed as $D = BB^T$ (LaBolle et al., 1996), where the superscript $T$ stands for transpose. The matrix $B$ is well known in the context of Random Walk Particle Tracking (RWPT), where it is denoted as the displacement matrix. In two dimensions it is given by (Fernández-García et al., 2005; Salamon et al., 2006)

$$ B = \begin{bmatrix} \frac{v_x}{|v|} \sqrt{\alpha_L |v| + D_m} & - \frac{v_y}{|v|} \sqrt{\alpha_T |v| + D_m} \\ \frac{v_y}{|v|} \sqrt{\alpha_L |v| + D_m} & \frac{v_x}{|v|} \sqrt{\alpha_T |v| + D_m} \end{bmatrix}. \quad (15) $$

Furthermore, the inverse of this matrix can be computed analytically as

$$ B^{-1} = \begin{bmatrix} \frac{1}{\sqrt{\alpha_L |v| + D_m}} & \frac{v_x}{|v|} \\ -\frac{1}{\sqrt{\alpha_T |v| + D_m}} & \frac{v_y}{|v|} \end{bmatrix} \quad (16) $$

The coordinate transformation in ASPHAD is such that $X = B^{-1}X$, where $X$ is the position in the new frame of reference. This leads to the following SPH scheme for dispersion (Tran-Duc et al., 2016),

$$ \frac{dC_a}{dt} = 2 \sum_b \frac{m_b}{\rho_b} \frac{C_{ab}}{|\bar{B}_{ab}^{-1} \hat{r}_{ab}|^2} \frac{F_{ab}}{r_{ab}}, \quad (17) $$

where $\bar{B}_{ab}^{-1}$ is the inverse displacement matrix (Eq. (16)) averaged between particles $a$ and $b$, and $C_{ab} = C_a - C_b$. Tran-Duc et al. (2016) showed that this scheme may artificially increase dispersion in the direction of smaller dispersivity (transverse), and reduce it in the direction of larger dispersivity (longitudinal). However, a major advantage of this approach is that it is unconditionally stable (Brookshaw, 1985; Biriukov and Price, 2018).

### 2.3.3. Two first derivatives (TFD)

This method is presented by Biriukov and Price (2018) with application to heat transport, showing satisfactory results for an anisotropic conduction tensor. It differs from the previous methods in that the dispersive flux is first explicitly calculated and then used to compute the divergence of fluxes. It is thus a two-stage SPH integration. The jth-component of the gradient of solute concentrations is expressed as

$$ G_a^j = \nabla_a^j C_a = \frac{1}{\rho_a f_a} \sum_b m_b C_{ba} \nabla_a^j W_{ab}(h_a), \quad (18) $$

where

$$ f_a = 1 + \frac{h_a}{\nu \rho_a} \sum_b m_b \frac{\partial W_{ab}(h_a)}{\partial h_a}, \quad (19) $$
is a gradient correction which stems from considering a spatially variable smoothing length following the relation $h_a^v \propto 1/\rho_a$ (Price and Monaghan, 2007; Springel, 2010). The temporal variation of concentration at a particle is then calculated as the divergence of the product between the dispersion tensor and Eq. (18), which yields

$$ \frac{dC_a}{dt} = \rho_a \sum_b m_b \left[ \frac{D_{ij}^a G_a^i \nabla_a^j W_{ab}(h_a)}{f_a \rho_a^2} + \frac{D_{ij}^b G_b^i \nabla_a^j W_{ab}(h_b)}{f_b \rho_b^2} \right]. \quad (20) $$

As before, indexes i, j follow the Einstein convention. A formulation with a spatially variable smoothing length appears attractive in the context of heterogeneous systems. Nevertheless, in order to restrict the scope, in this work the kernel size is kept constant, and thus $f_a = 1$. The TFD scheme has been shown to be mass-conservative and stable. However, in the vicinity of sharp gradients, it may be prone to introduce extreme positive or negative values. Effects of this type are typically handled in SPH by artificial dissipation. That is, by artificial viscosity in hydrodynamics (Monaghan and Gingold, 1983) or by artificial conductivity in heat transport (Price, 2008; Biriukov and Price, 2018). The challenge lies in the fact that dissipation terms are only required in the vicinity of sharp gradients. Such corrections have not been discussed in the context of solute transport and potential analogies could be established with heat transport. Still, a definition of this kind is outside the scope of this article and it is important to stress that a formulation should be proposed with caution, in order to preserve the main feature of particle methods in heterogeneous systems, which is the absence of artificial dispersion.

## 2.4. Implementation

### 2.4.1. Non-uniform particle density

Some special considerations are needed while implementing the discussed SPH schemes for dispersion, in the context of particles emulating the advection through a heterogeneous domain with velocities linearly interpolated from a divergence-free, steady velocity field. Particularly, dispersive fluxes between neighbor particles should be anti-symmetric in order to satisfy mass conservation, this property being directly influenced by the estimation of local densities (see Herrera et al., 2009). Because the medium is heterogeneous, an initially uniform distribution of particles will most likely end up evolving into a seemingly quasi-random distribution (as in Herrera and Beckie, 2013). This means that the estimated particle density will be non-uniform and most likely evolving in time as particles are displaced, which is somewhat comparable to introducing artificial divergence in the flow (from the perspective of the SPH particles). Notice that the ER scheme (Eq. (12)) involves an averaged density between each pair of neighboring particles, hence the dispersive flux remains anti-symmetric. The same does not occur with default ASPHAD (Eq. (17)) and TFD (Eq. (20)) if the schemes are applied directly without any correction to the densities. In the implementation of these methods, densities are also replaced with the average between neighbors $\bar{\rho}_{ab}$, thus ensuring the anti-symmetry of dispersive fluxes.

### 2.4.2. Temporal integration

Temporal integration is performed with an explicit scheme. An adaptive time step based on dispersion, $\delta_t^D$, is computed as the minimum over all particles of the expression

$$ \delta_t^D \le C_T \frac{h^2}{\sum_l D_{il}}, \quad (21) $$

with $C_T = 0.1$ (Herrera and Beckie, 2013). For heterogeneous experiments, also the advection-limited time step $\delta_t^v$ is computed by imposing $CFL = 0.1$ and using the maximum spatial velocities, that is

$$ \delta_t^v \le \frac{CFL}{\max v_x / \Delta_x + \max v_y / \Delta_y}, \quad (22) $$

where $\Delta_x$, $\Delta_y$ are the flow grid sizes. Finally, time step is selected as the most restrictive, i.e.,

$$ \delta_t = \min\{\delta_t^D, \delta_t^v\}. \quad (23) $$

Solute concentrations are then integrated explicitly as

$$ C_a^{t+\delta_t} = C_a^t + \delta_t \left( \frac{dC_a}{dt} \right)^t, \quad (24) $$

where the time derivative is computed with any of the SPH schemes discussed previously.

### 2.4.3. Rounding error

The SPH implementation with floating point calculations introduces some rounding error. Sigalotti et al. (2019) showed that the relative error for SPH estimates $A$, following expressions of the form (5), is bounded by

$$ \frac{|A_a^{\text{fl}} - A_a|}{|A_a|} \le \frac{1}{2} \mathcal{N} \epsilon, \quad (25) $$

where $A_a^{\text{fl}}$ is the floating point representation of the SPH estimate and $\epsilon > 0$ is the machine epsilon. The latter is a property of the floating point system and is the maximum number that satisfies $1 = 1+\epsilon$. Models were implemented here using 64 bit precision with $\epsilon = 2.220446 \times 10^{-16}$. Expression (25) shows that the relative rounding error is proportional to the number of neighbor particles.

### 2.4.4. Grid interpolation

In the heterogeneous test case, results from SPH particles are interpolated into a reference grid for the purposes of analysis and comparison with the reference model. The interpolation towards the center of a cell $c$ is performed from the set of SPH particles as (Herrera et al., 2009)

$$ C_c = \frac{1}{\rho_c} \sum_b^{\mathcal{N}} m_b C_b W_{cb}, \quad (26) $$

where $\rho_c$ is the density computed at the cell center. Interpolations are performed with the same kernel size employed during simulations.

## 3. Results and discussion

### 3.1. Anisotropic dispersion in a homogeneous medium

This test considers the transport of a Gaussian injection through a homogeneous domain, already discussed in previous SPH literature (Herrera and Beckie, 2013; Avesani et al., 2015; Alvarado-Rodríguez et al., 2019; Klapp et al., 2022). The flow is uniform, without changes in the relative position of particles, implying that neighbors and density are computed only once. The analytical solution $C^*$ is given by

$$ \frac{C^*(\mathbf{r}, t)}{C_0} = \frac{d^2}{B_3} \exp \left[ \frac{-\bar{x}^2 A_1 - \bar{y}^2 A_2 + 4 \bar{x} \bar{y} A_3}{8t^2 B_1 + 4td^2 B_2 + 2d^4} \right], \quad (27) $$

where

$$ \begin{aligned} A_1 &= 2t D_{yy} + d^2 \\ A_2 &= 2t D_{xx} + d^2 \\ A_3 &= t D_{xy} \\ B_1 &= D_{xx} D_{yy} - D_{xy}^2 \\ B_2 &= D_{xx} + D_{yy} \\ B_3 &= (4t^2 B_1 + 2td^2 B_2 + d^4)^{1/2}. \end{aligned} \quad (28) $$

In the above expression, $C_0$ is the maximum initial concentration, which is employed to express concentrations in dimensionless form, $\bar{C} = C/C_0$, $(\bar{x}, \bar{y}) = (x-x_0, y-y_0)$ are the relative coordinates with respect to the plume center $(x_0, y_0) = (0, 0)$. The parameters follow values taken
is a gradient correction which stems from considering a spatially variable smoothing length following the relation $h_a^v \propto 1/\rho_a$ (Price and Monaghan, 2007; Springel, 2010). The temporal variation of concentration at a particle is then calculated as the divergence of the product between the dispersion tensor and Eq. (18), which yields

$$ \frac{dC_a}{dt} = \rho_a \sum_b m_b \left[ \frac{D_{ij}^a G_a^i \nabla_a^j W_{ab}(h_a)}{f_a \rho_a^2} + \frac{D_{ij}^b G_b^i \nabla_a^j W_{ab}(h_b)}{f_b \rho_b^2} \right]. \quad (20) $$

As before, indexes i, j follow the Einstein convention. A formulation with a spatially variable smoothing length appears attractive in the context of heterogeneous systems. Nevertheless, in order to restrict the scope, in this work the kernel size is kept constant, and thus $f_a = 1$. The TFD scheme has been shown to be mass-conservative and stable. However, in the vicinity of sharp gradients, it may be prone to introduce extreme positive or negative values. Effects of this type are typically handled in SPH by artificial dissipation. That is, by artificial viscosity in hydrodynamics (Monaghan and Gingold, 1983) or by artificial conductivity in heat transport (Price, 2008; Biriukov and Price, 2018). The challenge lies in the fact that dissipation terms are only required in the vicinity of sharp gradients. Such corrections have not been discussed in the context of solute transport and potential analogies could be established with heat transport. Still, a definition of this kind is outside the scope of this article and it is important to stress that a formulation should be proposed with caution, in order to preserve the main feature of particle methods in heterogeneous systems, which is the absence of artificial dispersion.

## 2.4. Implementation

### 2.4.1. Non-uniform particle density

Some special considerations are needed while implementing the discussed SPH schemes for dispersion, in the context of particles emulating the advection through a heterogeneous domain with velocities linearly interpolated from a divergence-free, steady velocity field. Particularly, dispersive fluxes between neighbor particles should be anti-symmetric in order to satisfy mass conservation, this property being directly influenced by the estimation of local densities (see Herrera et al., 2009). Because the medium is heterogeneous, an initially uniform distribution of particles will most likely end up evolving into a seemingly quasi-random distribution (as in Herrera and Beckie, 2013). This means that the estimated particle density will be non-uniform and most likely evolving in time as particles are displaced, which is somewhat comparable to introducing artificial divergence in the flow (from the perspective of the SPH particles). Notice that the ER scheme (Eq. (12)) involves an averaged density between each pair of neighboring particles, hence the dispersive flux remains anti-symmetric. The same does not occur with default ASPHAD (Eq. (17)) and TFD (Eq. (20)) if the schemes are applied directly without any correction to the densities. In the implementation of these methods, densities are also replaced with the average between neighbors $\bar{\rho}_{ab}$, thus ensuring the anti-symmetry of dispersive fluxes.

### 2.4.2. Temporal integration

Temporal integration is performed with an explicit scheme. An adaptive time step based on dispersion, $\delta_t^D$, is computed as the minimum over all particles of the expression

$$ \delta_t^D \le C_T \frac{h^2}{\sum_l D_{il}}, \quad (21) $$

with $C_T = 0.1$ (Herrera and Beckie, 2013). For heterogeneous experiments, also the advection-limited time step $\delta_t^v$ is computed by imposing $CFL = 0.1$ and using the maximum spatial velocities, that is

$$ \delta_t^v \le \frac{CFL}{\max v_x / \Delta_x + \max v_y / \Delta_y}, \quad (22) $$

where $\Delta_x$, $\Delta_y$ are the flow grid sizes. Finally, time step is selected as the most restrictive, i.e.,

$$ \delta_t = \min\{\delta_t^D, \delta_t^v\}. \quad (23) $$

Solute concentrations are then integrated explicitly as

$$ C_a^{t+\delta_t} = C_a^t + \delta_t \left( \frac{dC_a}{dt} \right)^t, \quad (24) $$

where the time derivative is computed with any of the SPH schemes discussed previously.

### 2.4.3. Rounding error

The SPH implementation with floating point calculations introduces some rounding error. Sigalotti et al. (2019) showed that the relative error for SPH estimates $A$, following expressions of the form (5), is bounded by

$$ \frac{|A_a^{\text{fl}} - A_a|}{|A_a|} \le \frac{1}{2} \mathcal{N} \epsilon, \quad (25) $$

where $A_a^{\text{fl}}$ is the floating point representation of the SPH estimate and $\epsilon > 0$ is the machine epsilon. The latter is a property of the floating point system and is the maximum number that satisfies $1 = 1+\epsilon$. Models were implemented here using 64 bit precision with $\epsilon = 2.220446 \times 10^{-16}$. Expression (25) shows that the relative rounding error is proportional to the number of neighbor particles.

### 2.4.4. Grid interpolation

In the heterogeneous test case, results from SPH particles are interpolated into a reference grid for the purposes of analysis and comparison with the reference model. The interpolation towards the center of a cell $c$ is performed from the set of SPH particles as (Herrera et al., 2009)

$$ C_c = \frac{1}{\rho_c} \sum_b^{\mathcal{N}} m_b C_b W_{cb}, \quad (26) $$

where $\rho_c$ is the density computed at the cell center. Interpolations are performed with the same kernel size employed during simulations.

## 3. Results and discussion

### 3.1. Anisotropic dispersion in a homogeneous medium

This test considers the transport of a Gaussian injection through a homogeneous domain, already discussed in previous SPH literature (Herrera and Beckie, 2013; Avesani et al., 2015; Alvarado-Rodríguez et al., 2019; Klapp et al., 2022). The flow is uniform, without changes in the relative position of particles, implying that neighbors and density are computed only once. The analytical solution $C^*$ is given by

$$ \frac{C^*(\mathbf{r}, t)}{C_0} = \frac{d^2}{B_3} \exp \left[ \frac{-\bar{x}^2 A_1 - \bar{y}^2 A_2 + 4\bar{x}\bar{y} A_3}{8t^2 B_1 + 4td^2 B_2 + 2d^4} \right], \quad (27) $$

where

$$ \begin{aligned} A_1 &= 2t D_{yy} + d^2 \\ A_2 &= 2t D_{xx} + d^2 \\ A_3 &= t D_{xy} \\ B_1 &= D_{xx} D_{yy} - D_{xy}^2 \\ B_2 &= D_{xx} + D_{yy} \\ B_3 &= (4t^2 B_1 + 2td^2 B_2 + d^4)^{1/2}. \end{aligned} \quad (28) $$

In the above expression, $C_0$ is the maximum initial concentration, which is employed to express concentrations in dimensionless form, $\bar{C} = C/C_0$. $(\bar{x}, \bar{y}) = (x-x_0, y-y_0)$ are the relative coordinates with respect to the plume center $(x_0, y_0) = (0, 0)$. The parameters follow values taken
<figure><img src="image_6.png" /><figcaption>Fig. 2. Numerical error $e_R$ with respect to the analytical solution as a function of the kernel size, for the anisotropic dispersion of a Gaussian injection through homogeneous medium at reference time $T = 300$[d]. (a) ER, (b) ASPHAD and (c) TFD. Simulations consider $N = 300 \times 300$ particles with uniform distribution, and flow aligned with the x-axis.</figcaption></figure>

from Herrera and Beckie (2013). The initial width of the Gaussian plume is $d = 44$[m] and the domain size is $L \times L = 2000 \times 2000$[m × m]. Flow velocity is constant, $|v| = 1$[m/d], following a vector that forms the angle $\theta$ with respect to the x-direction, although for most simulations is simply considered to be aligned with the x-axis. Dispersion neglects molecular diffusion, with longitudinal dispersivity $\alpha_L = 10$[m] and anisotropies $\lambda_D \in \{1, 10, 100\}$. The number of SPH particles is modified within the range $N \in [200 \times 200, 600 \times 600]$, and the $\gamma$ parameter for the smoothing length varies within the range $\gamma \in [1, 7]$. Transport simulations end at a reference time $T = 300$[d]. Results are mainly discussed with a focus on uniform and quasi-random particle distributions, and simulations with a random distribution are taken as a reference worst case scenario. The quasi-random distribution is generated by applying a random perturbation $\epsilon_r \sim U(-\delta_r/2, \delta_r/2)$ to the coordinates of the uniform distribution, with $\delta_r$ being the original spacing between particles. Differences between the numerical and analytical results are quantified by means of the Root Mean Squared Error (RMSE),

$$e_R = \sqrt{\frac{1}{N} \sum_a^N (C_a^* - C_a)^2}. \qquad (29)$$

### 3.1.1. Influence of dispersion anisotropy

Because each SPH scheme is formulated differently it is expected that their respective responses to anisotropic dispersion will be different. Fig. 2 depicts the error for all scenarios of anisotropy with varying kernel size, for a reference uniform distribution of particles, and flow aligned with the x-axis. In all cases, there is a smoothing that minimizes the error with some influence of anisotropy on the

<figure><img src="image_7.png" /><figcaption>Fig. 3. Concentrations for the problem of Gaussian dispersion with uniform distribution of particles ($N = 300 \times 300$), at the reference time $T = 300$[d], and flow aligned with the x-axis. First column group results for all SPH schemes considering $\lambda_D = 10$ (panels a,c,e,g) and second column for $\lambda_D = 100$ (panels b,d,f,h). Contour lines are drawn for $\tilde{C} \in \{0.4, 0.1, 10^{-3}, 10^{-5}\}$ (solid: analytical, dotted: numerical). White bands (if present) represent negative values.</figcaption></figure>

optimal size. The ER scheme's results are influenced by negative concentrations for the considered values of anisotropy, in concordance with the results reported in the literature (Fig. 3c, d). Both ASPHAD and TFD are free of negative concentrations for anisotropy $\lambda_D = 10$ (Fig. 3e, g respectively), with some differences in the small-scale spatial distribution of concentrations. Particularly, TFD appears to generate a more elliptical distribution for contours $\tilde{C} \sim 10^{-5}$, but in general both schemes are in good agreement with the analytical reference. Increasing the anisotropy to $\lambda_D = 100$ leads to some peculiar results. In particular, the concentration contours in the ER method (Fig. 3d) are now influenced by the bands of negative values generated by the instability, with a minimum concentration of $\tilde{C} = -4.6 \times 10^{-5}$. The TFD results (Fig. 3h) also present some negative values while moving away from the plume center, but these are of a different nature. The minimum concentration is $\tilde{C} = -9.8 \times 10^{-16}$, which is of the order of numerical precision. Furthermore, it has been noticed that reducing the kernel size (smaller number of neighbors in Eq. (25)) eliminates the effect. It has been detected that negative concentrations can be generated from very small differences (numerical precision) in the
<figure><img src="image_5.png" /><figcaption>Fig. 2. Numerical error $e_R$ with respect to the analytical solution as a function of the kernel size, for the anisotropic dispersion of a Gaussian injection through homogeneous medium at reference time $T = 300$[d]. (a) ER, (b) ASPHAD and (c) TFD. Simulations consider $N = 300 \times 300$ particles with uniform distribution, and flow aligned with the x-axis.</figcaption></figure>

from Herrera and Beckie (2013). The initial width of the Gaussian plume is $d = 44$[m] and the domain size is $L \times L = 2000 \times 2000$[m \times m]. Flow velocity is constant, $|v| = 1$[m/d], following a vector that forms the angle $\theta$ with respect to the x-direction, although for most simulations is simply considered to be aligned with the x-axis. Dispersion neglects molecular diffusion, with longitudinal dispersivity $\alpha_L = 10$[m] and anisotropies $\lambda_D \in \{1, 10, 100\}$. The number of SPH particles is modified within the range $N \in [200 \times 200, 600 \times 600]$, and the $\gamma$ parameter for the smoothing length varies within the range $\gamma \in [1, 7]$. Transport simulations end at a reference time $T = 300$[d]. Results are mainly discussed with a focus on uniform and quasi-random particle distributions, and simulations with a random distribution are taken as a reference worst case scenario. The quasi-random distribution is generated by applying a random perturbation $\epsilon_r \sim U(-\delta_r/2, \delta_r/2)$ to the coordinates of the uniform distribution, with $\delta_r$ being the original spacing between particles. Differences between the numerical and analytical results are quantified by means of the Root Mean Squared Error (RMSE),

$$e_R = \sqrt{\frac{1}{N} \sum_a^N (C_a^* - C_a)^2}. \qquad (29)$$

### 3.1.1. Influence of dispersion anisotropy

Because each SPH scheme is formulated differently it is expected that their respective responses to anisotropic dispersion will be different. Fig. 2 depicts the error for all scenarios of anisotropy with varying kernel size, for a reference uniform distribution of particles, and flow aligned with the x-axis. In all cases, there is a smoothing that minimizes the error with some influence of anisotropy on the

<figure><img src="image_6.png" /><figcaption>Fig. 3. Concentrations for the problem of Gaussian dispersion with uniform distribution of particles ($N = 300 \times 300$), at the reference time $T = 300$[d], and flow aligned with the x-axis. First column group results for all SPH schemes considering $\lambda_D = 10$ (panels a,c,e,g) and second column for $\lambda_D = 100$ (panels b,d,f,h). Contour lines are drawn for $\tilde{C} \in \{0.4, 0.1, 10^{-3}, 10^{-5}\}$ (solid: analytical, dotted: numerical). White bands (if present) represent negative values.</figcaption></figure>

optimal size. The ER scheme's results are influenced by negative concentrations for the considered values of anisotropy, in concordance with the results reported in the literature (Fig. 3c, d). Both ASPHAD and TFD are free of negative concentrations for anisotropy $\lambda_D = 10$ (Fig. 3e, g respectively), with some differences in the small-scale spatial distribution of concentrations. Particularly, TFD appears to generate a more elliptical distribution for contours $\tilde{C} \sim 10^{-5}$, but in general both schemes are in good agreement with the analytical reference. Increasing the anisotropy to $\lambda_D = 100$ leads to some peculiar results. In particular, the concentration contours in the ER method (Fig. 3d) are now influenced by the bands of negative values generated by the instability, with a minimum concentration of $\tilde{C} = -4.6 \times 10^{-5}$. The TFD results (Fig. 3h) also present some negative values while moving away from the plume center, but these are of a different nature. The minimum concentration is $\tilde{C} = -9.8 \times 10^{-16}$, which is of the order of numerical precision. Furthermore, it has been noticed that reducing the kernel size (smaller number of neighbors in Eq. (25)) eliminates the effect. It has been detected that negative concentrations can be generated from very small differences (numerical precision) in the
<caption>**Table 1**</caption>

Errors in homogeneous simulations with dispersion anisotropy, for different angles of alignment $\theta$ between the flow direction and the x-axis. Table shows the value of $\gamma$ that lead to the minimum error $e_R$, considering $N = 300 \times 300$ and uniform distribution of particles.

<table><thead><tr><th rowspan="2">λ<sub>D</sub></th><th rowspan="2">θ</th><th colspan="2">ER</th><th colspan="2">ASPHAD</th><th colspan="2">TFD</th></tr><tr><th>γ</th><th>e<sub>R</sub></th><th>γ</th><th>e<sub>R</sub></th><th>γ</th><th>e<sub>R</sub></th></tr></thead><tbody><tr><td rowspan="3">10</td><td>0</td><td>4.5</td><td>0.00013</td><td>1.4</td><td>0.00205</td><td>3.6</td><td>0.00062</td></tr><tr><td>π/6</td><td>4.5</td><td>0.00015</td><td>1.3</td><td>0.00420</td><td>3.6</td><td>0.00062</td></tr><tr><td>π/4</td><td>4.5</td><td>0.00015</td><td>1.4</td><td>0.00299</td><td>3.6</td><td>0.00062</td></tr><tr><td rowspan="3">100</td><td>0</td><td>4.2</td><td>0.00032</td><td>1.3</td><td>0.00063</td><td>3.6</td><td>0.00056</td></tr><tr><td>π/6</td><td>4.2</td><td>0.00035</td><td>2.4</td><td>0.01202</td><td>3.6</td><td>0.00056</td></tr><tr><td>π/4</td><td>4.2</td><td>0.00036</td><td>1.1</td><td>0.00434</td><td>3.6</td><td>0.00056</td></tr></tbody></table>

gradients following the direction of higher dispersivity near regions of zero concentration, only in scenarios of extreme anisotropy.

For this numerical problem, studies have discussed simulations with non-zero diagonal terms in the dispersion tensor (e.g., Herrera and Beckie, 2013; Avesani et al., 2015; Alvarado-Rodríguez et al., 2019; Klapp et al., 2022). These occur when the flow velocity vector is not aligned with the main axes of the reference coordinates (as seen in Eq. (3)), and is a condition to be expected in heterogeneous flow systems. At this point, it should be remarked that the relative alignment of the flow does not explain the occurrence of negative concentrations observed in results from the ER scheme while considering anisotropic dispersion. This was shown, for example, in the work of Alvarado-Rodríguez et al. (2019) where patterns of negative values were visible

while considering the flow aligned with the x-axis, which is consistent with the results previously shown in Figure (3c, d). To further expand on this topic, anisotropic simulations with all the considered SPH schemes were also performed with the velocity vector forming angles of $\theta \in \{\pi/6, \pi/4\}$ with respect to the x-axis (positive counterclockwise), and a uniform distribution of particles. The results from these simulations (refer to Table 1) show that the ASPHAD scheme is sensitive to the alignment between the flow and the particle lattice, more noticeably, in the cases with highest anisotropy. The worst case condition for this method is seen while considering $\theta = \pi/6$. In this particular case, the flow orientation is not aligned with any of the angles naturally formed by the uniform lattice of particles, as it occurs in the cases with $\theta = 0$ (particles aligned along the x-axis) or $\theta = \pi/4$ (diagonally aligned particles), leading to a noticeable increase in errors. This suggest that in cases of high anisotropy, the approximation of dispersion provided by ASPHAD will perform best in cases that the flow is aligned with the particles, although in any case a non-zero flow orientation increases the magnitude of errors obtained with this method (e.g., Fig. 4d). In contrast, anisotropic simulations performed with the ER (Fig. 4a, b) and TFD (Fig. 4e, f) schemes are robust to the relative alignment of the flow. With these methods, errors remain practically unaffected while modifying the flow orientation (Table 1) for both magnitudes of dispersion anisotropy, also preserving the previously observed agreement between the numerical and analytical concentration contours.

### 3.1.2. Errors and the particle distribution

The relative position of particles influences the rate at which the discrete error decreases while increasing the number of neighbors (Zhu et al., 2015). The case of a uniform distribution is an idealized scenario unlikely to occur in practical applications (e.g., Herrera and Beckie, 2013), thus understanding the performance of SPH schemes for different particle distributions is of relevance. The results indicate that all methods yield higher error in scenarios of particle disorder (Fig. 5). Errors obtained from ASPHAD show relatively large magnitudes for a

<figure><img src="image_8.png" /><figcaption>**Fig. 4.** Concentrations for the problem of Gaussian dispersion with uniform distribution of particles ($N = 300 \times 300$), at the reference time $T = 300[d]$, and diagonally aligned flow ($\theta = \pi/4$). First column group results for all SPH schemes considering $\lambda_D = 10$ (panels a, c, e) and second column for $\lambda_D = 100$ (panels b, d, f). Contour lines are drawn for $\tilde{C} \in \{0.4, 0.1, 10^{-3}, 10^{-5}\}$ (solid: analytical, dotted: numerical). White bands (if present) represent negative values.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 5.** Numerical error $e_R$ for different particle distributions, with respect to the kernel size, for the problem of Gaussian dispersion in a homogeneous medium. Simulations consider $N = 300 \times 300$, flow aligned with the x-axis, and errors computed at the reference time $T = 300[d]$. First column group results for isotropic dispersion (panels a, c, e) and second column for anisotropy $\lambda_D = 100$ (panels b, d, f). SPH method is indicated in each panel.</figcaption></figure>
<caption>**Table 1**</caption>

Errors in homogeneous simulations with dispersion anisotropy, for different angles of alignment $\theta$ between the flow direction and the x-axis. Table shows the value of $\gamma$ that lead to the minimum error $e_R$, considering $N = 300 \times 300$ and uniform distribution of particles.

<table><thead><tr><th rowspan="2">λ<sub>D</sub></th><th rowspan="2">θ</th><th colspan="2">ER</th><th colspan="2">ASPHAD</th><th colspan="2">TFD</th></tr><tr><th>γ</th><th>e<sub>R</sub></th><th>γ</th><th>e<sub>R</sub></th><th>γ</th><th>e<sub>R</sub></th></tr></thead><tbody><tr><td rowspan="3">10</td><td>0</td><td>4.5</td><td>0.00013</td><td>1.4</td><td>0.00205</td><td>3.6</td><td>0.00062</td></tr><tr><td>π/6</td><td>4.5</td><td>0.00015</td><td>1.3</td><td>0.00420</td><td>3.6</td><td>0.00062</td></tr><tr><td>π/4</td><td>4.5</td><td>0.00015</td><td>1.4</td><td>0.00299</td><td>3.6</td><td>0.00062</td></tr><tr><td rowspan="3">100</td><td>0</td><td>4.2</td><td>0.00032</td><td>1.3</td><td>0.00063</td><td>3.6</td><td>0.00056</td></tr><tr><td>π/6</td><td>4.2</td><td>0.00035</td><td>2.4</td><td>0.01202</td><td>3.6</td><td>0.00056</td></tr><tr><td>π/4</td><td>4.2</td><td>0.00036</td><td>1.1</td><td>0.00434</td><td>3.6</td><td>0.00056</td></tr></tbody></table>

gradients following the direction of higher dispersivity near regions of zero concentration, only in scenarios of extreme anisotropy.

For this numerical problem, studies have discussed simulations with non-zero diagonal terms in the dispersion tensor (e.g., Herrera and Beckie, 2013; Avesani et al., 2015; Alvarado-Rodríguez et al., 2019; Klapp et al., 2022). These occur when the flow velocity vector is not aligned with the main axes of the reference coordinates (as seen in Eq. (3)), and is a condition to be expected in heterogeneous flow systems. At this point, it should be remarked that the relative alignment of the flow does not explain the occurrence of negative concentrations observed in results from the ER scheme while considering anisotropic dispersion. This was shown, for example, in the work of Alvarado-Rodríguez et al. (2019) where patterns of negative values were visible

while considering the flow aligned with the x-axis, which is consistent with the results previously shown in Figure (3c, d). To further expand on this topic, anisotropic simulations with all the considered SPH schemes were also performed with the velocity vector forming angles of $\theta \in \{\pi/6, \pi/4\}$ with respect to the x-axis (positive counterclockwise), and a uniform distribution of particles. The results from these simulations (refer to Table 1) show that the ASPHAD scheme is sensitive to the alignment between the flow and the particle lattice, more noticeably, in the cases with highest anisotropy. The worst case condition for this method is seen while considering $\theta = \pi/6$. In this particular case, the flow orientation is not aligned with any of the angles naturally formed by the uniform lattice of particles, as it occurs in the cases with $\theta = 0$ (particles aligned along the x-axis) or $\theta = \pi/4$ (diagonally aligned particles), leading to a noticeable increase in errors. This suggest that in cases of high anisotropy, the approximation of dispersion provided by ASPHAD will perform best in cases that the flow is aligned with the particles, although in any case a non-zero flow orientation increases the magnitude of errors obtained with this method (e.g., Fig. 4d). In contrast, anisotropic simulations performed with the ER (Fig. 4a, b) and TFD (Fig. 4e, f) schemes are robust to the relative alignment of the flow. With these methods, errors remain practically unaffected while modifying the flow orientation (Table 1) for both magnitudes of dispersion anisotropy, also preserving the previously observed agreement between the numerical and analytical concentration contours.

### 3.1.2. Errors and the particle distribution

The relative position of particles influences the rate at which the discrete error decreases while increasing the number of neighbors (Zhu et al., 2015). The case of a uniform distribution is an idealized scenario unlikely to occur in practical applications (e.g., Herrera and Beckie, 2013), thus understanding the performance of SPH schemes for different particle distributions is of relevance. The results indicate that all methods yield higher error in scenarios of particle disorder (Fig. 5). Errors obtained from ASPHAD show relatively large magnitudes for a

<figure><img src="image_7.png" /><figcaption>**Fig. 4.** Concentrations for the problem of Gaussian dispersion with uniform distribution of particles ($N = 300 \times 300$), at the reference time $T = 300[d]$, and diagonally aligned flow ($\theta = \pi/4$). First column group results for all SPH schemes considering $\lambda_D = 10$ (panels a, c, e) and second column for $\lambda_D = 100$ (panels b, d, f). Contour lines are drawn for $\tilde{C} \in \{0.4, 0.1, 10^{-3}, 10^{-5}\}$ (solid: analytical, dotted: numerical). White bands (if present) represent negative values.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 5.** Numerical error $e_R$ for different particle distributions, with respect to the kernel size, for the problem of Gaussian dispersion in a homogeneous medium. Simulations consider $N = 300 \times 300$, flow aligned with the x-axis, and errors computed at the reference time $T = 300[d]$. First column group results for isotropic dispersion (panels a, c, e) and second column for anisotropy $\lambda_D = 100$ (panels b, d, f). SPH method is indicated in each panel.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Fig. 6.** Concentrations for anisotropic dispersion at $T = 300[d]$ with quasi-random distribution, $N = 300 \times 300$ particles, and flow aligned with the x-axis. First column group numerical results for each solver with $\lambda_D = 10$ and second column for $\lambda_D = 100$. White bands (if present) indicate negative values.</figcaption></figure>

quasi-random particle distribution (similar to the random case) for the highest dispersion anisotropy (Fig. 5d). The same does not occur with the ER (Fig. 5a, b) and TFD (Fig. 5e, f) methods which are able to preserve clear differences with respect to the worst case scenario regardless of anisotropy. Furthermore, the contour maps of anisotropic simulations with quasi-random particle distribution show that ER (Fig. 6a, b) and TFD (Fig. 6e, f) yield results which are consistent with the uniform case (Fig. 3). With this last scheme, the simulation of highest anisotropy now presents a more evident pattern of negative values, reaching a minimum concentration of $\tilde{C} = -4 \times 10^{-7}$, illustrating the influence of the particle disorder on error propagation. Regardless, it is clear that the method is robust with respect to the distribution of particles, leading to an accurate representation of concentration profiles (Fig. 7). Contours predicted by ASPHAD, in contrast, are highly impacted by the particle disorder (Fig. 6c, d). The effect was discussed in Tran-Duc et al. (2016), but only considering ratios of anisotropy of up to $\lambda_D = 10$. In this case, numerical peak concentration is ≈ 14% higher than the reference solution (Fig. 7a, b). Even more drastical discrepancies arise for anisotropy $\lambda_D = 100$, where the peak concentration is now ≈ 53% higher

<figure><img src="image_11.png" /><figcaption>**Fig. 7.** Numerical and analytical concentration profiles through the center of a Gaussian plume at time $T = 300[d]$, considering $N = 300 \times 300$, quasi-random distribution, and flow aligned with the x-axis. Panels (a, b) present the longitudinal (x) and transverse (y) profiles, respectively, for dispersion anisotropy $\lambda_D = 10$ and panels (c, d) the analogous for $\lambda_D = 100$.</figcaption></figure>

than the analytical model (Fig. 7c, d), also with major differences in the longitudinal extent of the plume (i.e., the actual anisotropy appears to be markedly lower than intended). These results indicate that ASPHAD will overpredict concentration gradients in scenarios of anisotropy and disordered particle distribution. Still, it remains to be seen whether this effect will be equally relevant in simulations where the main transport mechanism is the advection of particles, which is discussed in a later section.

### 3.1.3. Number of particles and consistency

The parameters of simulations that yielded the minimum error, considering the flow aligned with the x-axis, are discussed as a function of the total number of particles, in order to characterize the consistency of SPH schemes. As outlined previously, it should be satisfied that as the number of particles increases, kernel sizes decrease while simultaneously increasing the number of neighbors (Zhu et al., 2015; Sigalotti et al., 2016). Expression (11) relates the kernel size with these principles of consistent scaling. The trend of the minimum-error $\gamma$ for each SPH method, with respect to the number of particles is shown in Fig. 8, and compared with the reference scaling for two-dimensional domains ($\beta = 4$ in Eq. (11)). Overall, the ER and TFD schemes display an increasing value of $\gamma$ for increasing number of particles, almost independently from the dispersion anisotropy, which leads to a consistent decrease in error for both uniform and quasi-random distributions of particles (Fig. 8b and f, respectively). Errors

<caption>**Table 2**</caption>

Parameters of expression (11) obtained from fitting $\gamma$ and $N$ of the minimum error homogeneous simulations. Fitting is performed individually for each combination of dispersion anisotropy and particle distribution (Uni: uniform; Q.R.: quasi-random). $r^2$ is the coefficient of determination.

<table><thead><tr><th rowspan="2">Dist.</th><th rowspan="2">λ<sub>D</sub></th><th colspan="3">ER</th><th colspan="3">ASPHAD</th><th colspan="3">TFD</th></tr><tr><th>K</th><th>β</th><th>r<sup>2</sup></th><th>K</th><th>β</th><th>r<sup>2</sup></th><th>K</th><th>β</th><th>r<sup>2</sup></th></tr></thead><tbody><tr><td rowspan="3">Uni</td><td>1</td><td>2.72</td><td>3.52</td><td>0.91</td><td>2.70</td><td>3.61</td><td>0.99</td><td>0.46</td><td>4.39</td><td>0.99</td></tr><tr><td>10</td><td>2.03</td><td>3.64</td><td>0.98</td><td>32.9</td><td>2</td><td>-</td><td>0.79</td><td>3.85</td><td>0.99</td></tr><tr><td>100</td><td>0.56</td><td>4.26</td><td>0.97</td><td>24.7</td><td>2</td><td>-</td><td>0.93</td><td>3.77</td><td>0.99</td></tr><tr><td rowspan="3">Q.R.</td><td>1</td><td>284</td><td>2.08</td><td>0.01</td><td>0.04</td><td>8.11</td><td>0.82</td><td>0.31</td><td>4.87</td><td>0.95</td></tr><tr><td>10</td><td>43.6</td><td>2.48</td><td>0.41</td><td>2.36</td><td>2.68</td><td>0.33</td><td>1.08</td><td>3.80</td><td>0.93</td></tr><tr><td>100</td><td>266</td><td>2.05</td><td>0.01</td><td>3.62</td><td>3.04</td><td>0.80</td><td>0.33</td><td>4.79</td><td>0.98</td></tr></tbody></table>
<figure><img src="image_9.png" /><figcaption>**Fig. 6.** Concentrations for anisotropic dispersion at $T = 300[d]$ with quasi-random distribution, $N = 300 \times 300$ particles, and flow aligned with the x-axis. First column group numerical results for each solver with $\lambda_D = 10$ and second column for $\lambda_D = 100$. White bands (if present) indicate negative values.</figcaption></figure>

quasi-random particle distribution (similar to the random case) for the highest dispersion anisotropy (Fig. 5d). The same does not occur with the ER (Fig. 5a, b) and TFD (Fig. 5e, f) methods which are able to preserve clear differences with respect to the worst case scenario regardless of anisotropy. Furthermore, the contour maps of anisotropic simulations with quasi-random particle distribution show that ER (Fig. 6a, b) and TFD (Fig. 6e, f) yield results which are consistent with the uniform case (Fig. 3). With this last scheme, the simulation of highest anisotropy now presents a more evident pattern of negative values, reaching a minimum concentration of $\tilde{C} = -4 \times 10^{-7}$, illustrating the influence of the particle disorder on error propagation. Regardless, it is clear that the method is robust with respect to the distribution of particles, leading to an accurate representation of concentration profiles (Fig. 7). Contours predicted by ASPHAD, in contrast, are highly impacted by the particle disorder (Fig. 6c, d). The effect was discussed in Tran-Duc et al. (2016), but only considering ratios of anisotropy of up to $\lambda_D = 10$. In this case, numerical peak concentration is ≈ 14% higher than the reference solution (Fig. 7a, b). Even more drastical discrepancies arise for anisotropy $\lambda_D = 100$, where the peak concentration is now ≈ 53% higher

<figure><img src="image_10.png" /><figcaption>**Fig. 7.** Numerical and analytical concentration profiles through the center of a Gaussian plume at time $T = 300[d]$, considering $N = 300 \times 300$, quasi-random distribution, and flow aligned with the x-axis. Panels (a, b) present the longitudinal (x) and transverse (y) profiles, respectively, for dispersion anisotropy $\lambda_D = 10$ and panels (c, d) the analogous for $\lambda_D = 100$.</figcaption></figure>

than the analytical model (Fig. 7c, d), also with major differences in the longitudinal extent of the plume (i.e., the actual anisotropy appears to be markedly lower than intended). These results indicate that ASPHAD will overpredict concentration gradients in scenarios of anisotropy and disordered particle distribution. Still, it remains to be seen whether this effect will be equally relevant in simulations where the main transport mechanism is the advection of particles, which is discussed in a later section.

### 3.1.3. Number of particles and consistency

The parameters of simulations that yielded the minimum error, considering the flow aligned with the x-axis, are discussed as a function of the total number of particles, in order to characterize the consistency of SPH schemes. As outlined previously, it should be satisfied that as the number of particles increases, kernel sizes decrease while simultaneously increasing the number of neighbors (Zhu et al., 2015; Sigalotti et al., 2016). Expression (11) relates the kernel size with these principles of consistent scaling. The trend of the minimum-error $\gamma$ for each SPH method, with respect to the number of particles is shown in Fig. 8, and compared with the reference scaling for two-dimensional domains ($\beta = 4$ in Eq. (11)). Overall, the ER and TFD schemes display an increasing value of $\gamma$ for increasing number of particles, almost independently from the dispersion anisotropy, which leads to a consistent decrease in error for both uniform and quasi-random distributions of particles (Fig. 8b and f, respectively). Errors

<caption>**Table 2**</caption>

Parameters of expression (11) obtained from fitting $\gamma$ and $N$ of the minimum error homogeneous simulations. Fitting is performed individually for each combination of dispersion anisotropy and particle distribution (Uni: uniform; Q.R.: quasi-random). $r^2$ is the coefficient of determination.

<table><thead><tr><th rowspan="2">Dist.</th><th rowspan="2">λ<sub>D</sub></th><th colspan="3">ER</th><th colspan="3">ASPHAD</th><th colspan="3">TFD</th></tr><tr><th>K</th><th>β</th><th>r<sup>2</sup></th><th>K</th><th>β</th><th>r<sup>2</sup></th><th>K</th><th>β</th><th>r<sup>2</sup></th></tr></thead><tbody><tr><td rowspan="3">Uni</td><td>1</td><td>2.72</td><td>3.52</td><td>0.91</td><td>2.70</td><td>3.61</td><td>0.99</td><td>0.46</td><td>4.39</td><td>0.99</td></tr><tr><td>10</td><td>2.03</td><td>3.64</td><td>0.98</td><td>32.9</td><td>2</td><td>-</td><td>0.79</td><td>3.85</td><td>0.99</td></tr><tr><td>100</td><td>0.56</td><td>4.26</td><td>0.97</td><td>24.7</td><td>2</td><td>-</td><td>0.93</td><td>3.77</td><td>0.99</td></tr><tr><td rowspan="3">Q.R.</td><td>1</td><td>284</td><td>2.08</td><td>0.01</td><td>0.04</td><td>8.11</td><td>0.82</td><td>0.31</td><td>4.87</td><td>0.95</td></tr><tr><td>10</td><td>43.6</td><td>2.48</td><td>0.41</td><td>2.36</td><td>2.68</td><td>0.33</td><td>1.08</td><td>3.80</td><td>0.93</td></tr><tr><td>100</td><td>266</td><td>2.05</td><td>0.01</td><td>3.62</td><td>3.04</td><td>0.80</td><td>0.33</td><td>4.79</td><td>0.98</td></tr></tbody></table>
<figure><img src="image_12.png" /><figcaption>**Fig. 8.** Minimum error simulations with respect to the number of particles, for all scenarios of dispersion anisotropy, with uniform (solid) and quasi-random (dashed) distribution of particles. SPH method is indicated in each panel. Panels (a), (c), (e) present the $\gamma$ factor, and panels (b), (d), (f) the corresponding error. Reference curves consider $\beta = 4$.</figcaption></figure>

obtained from ASPHAD (Fig. 8d) decay while increasing the number of particles only in scenarios of isotropic dispersion. In anisotropic conditions, the minimum error is always obtained for the same value of $\gamma$ in cases with uniform particle distribution, regardless of the particle resolution (Fig. 8c). In the quasi-random case, this parameter displays more variability, although, in any case, the error remains insensitive to the resolution in anisotropic simulations (Fig. 8d). The $\gamma$ values from simulations with minimum error are fitted to expression (11) with coefficients shown in Table 2. The determination coefficients for simulations performed with TFD were always higher than $r^2 = 0.93$, confirming that the scheme follows the scaling from Eq. (11), with error decaying consistently with the proportionality $e_R \propto N^{-0.5}$. Although for the quasi-random distribution of particles the error magnitude increases, the results of this analysis underline the robustness of TFD with respect to the particle disorder, with fitted coefficients similar to the uniform case in all scenarios of anisotropy.

## 3.2. Transport through a heterogeneous medium

This test considers a two-dimensional aquifer with spatial variability of hydraulic properties (Fig. 9) composed of $N_c = 1500 \times 300$ unitary cells ($\Delta_x, \Delta_y$) = (1, 1)[m]. The natural logarithm of the hydraulic conductivity $Y(x)$ represents one individual realization of a sequential Gaussian simulation with zero mean and spherical variogram with correlation length $I_Y = 20\Delta_x$. Heterogeneities are controlled by log-conductivity variances $\sigma_Y^2 \in \{0.25, 1.00, 2.25\}$ through the relation $K(x) = \exp(\sigma_Y Y(x))$. For each hydraulic conductivity distribution, groundwater flow is obtained from solving on the numerical grid the steady-state Darcy's equation $\nabla \cdot \mathbf{q} = 0$, where $\mathbf{q} = -K(x)\nabla H$, and $H$ is the hydraulic head. Flow is induced by a unit-mean hydraulic head

<figure><img src="image_13.png" /><figcaption>**Fig. 9.** Initial condition for heterogeneous aquifer with high spatial variability of hydraulic conductivity. Groundwater flow is from left to right and reference streamlines are initially spaced by 66 flow cells.</figcaption></figure>

gradient along the x-axis, by prescribing the hydraulic head values at the upstream and downstream aquifer boundaries. North and south aquifer boundaries are defined as impermeable. Divergence-free flow is verified for all degrees of heterogeneity. As before, solute transport simulations are performed for different ratios of dispersion anisotropy, setting $\alpha_L = 0.05\Delta_x$ which leads to advection-dominated transport, with a field Péclet number of $Pe = I_Y/\alpha_L = 400$. A rectangular solute plume is initially released near the aquifer inlet with uniform concentration $C_0$. The plume size is $d_x \times d_y = 30\Delta_x \times 200\Delta_y$, with the lower left corner placed at $(x_0, y_0) = (15\Delta_x, 50\Delta_y)$. To minimize the singularity near extremely sharp gradients for the TFD scheme, a slightly smoothed version of the initial condition is used (as in Biriukov and Price, 2018, see Appendix B). All SPH simulations begin with a uniform distribution of particles. Groundwater flow velocities are known at the cell faces and linearly interpolated to the particle positions. Particles leaving the aquifer are reinjected upstream by resetting the x-coordinate to zero and defining the y-coordinate from the intersection of a random number $\xi_r \sim U(0, 1)$ with the normalized cumulative probability distribution of the inflow velocity, ensuring an entrance of particles that is statistically consistent with the non-uniform distribution of inflow velocities. Aquifer boundaries are surrounded by fixed SPH particles that compensate the loss of particle support while computing density near the borders. These fixed particles are defined as no-flux boundaries from the solute transport perspective (as in Wang et al., 2019). SPH results are compared with a high-resolution RWPT model, which in contrast, interprets particles as individual components of a solute cloud (e.g., Salamon et al., 2006). In RWPT, the number of particles determines the accuracy of concentrations. The individual particle mass for the RWPT simulations is $m_R = C_0V_M/N_R$, where $V_M$ is the volume occupied by the solute injection and $N_R = 10^8$ the number of particles. These parameters lead to a concentration resolution of $C_{min} = 6 \times 10^{-5}C_0$, estimated as the mass of one particle inside a flow-model cell. These parameters provide enough confidence to presume that the binning of solute particles is a reasonable proxy for the exact concentrations.

### 3.2.1. Particle resolution and peak concentrations

SPH simulations performed with different particle resolution reveal the influence of this parameter on the accuracy of peak concentrations. This is particularly relevant in cases of medium-to-low aquifer heterogeneity and high dispersion anisotropy, where strong concentration gradients can remain persistently throughout the simulation. Three degrees of particle resolution are considered, $N \in \{0.2M, 0.8M, 1.8M\}$, which translates into averaged densities of $\bar{\rho} \in \{0.4, 1.8, 4\}$ particles per flow-model cell, respectively. Herrera et al. (2009) reported results for a similar problem considering isotropic dispersion and a reference density of $\bar{\rho} = 8$ particles per cell. Notice that the particle resolution not only influences the accuracy of SPH schemes from a dispersion perspective (as in the homogeneous test case), but also determines the sampling of the groundwater flow velocities, the latter being independent from the selected dispersion scheme. To illustrate, the case of medium variability $\sigma_Y^2 = 1$, and the highest dispersion
<figure><img src="image_11.png" /><figcaption>**Fig. 8.** Minimum error simulations with respect to the number of particles, for all scenarios of dispersion anisotropy, with uniform (solid) and quasi-random (dashed) distribution of particles. SPH method is indicated in each panel. Panels (a), (c), (e) present the $\gamma$ factor, and panels (b), (d), (f) the corresponding error. Reference curves consider $\beta = 4$.</figcaption></figure>

obtained from ASPHAD (Fig. 8d) decay while increasing the number of particles only in scenarios of isotropic dispersion. In anisotropic conditions, the minimum error is always obtained for the same value of $\gamma$ in cases with uniform particle distribution, regardless of the particle resolution (Fig. 8c). In the quasi-random case, this parameter displays more variability, although, in any case, the error remains insensitive to the resolution in anisotropic simulations (Fig. 8d). The $\gamma$ values from simulations with minimum error are fitted to expression (11) with coefficients shown in Table 2. The determination coefficients for simulations performed with TFD were always higher than $r^2 = 0.93$, confirming that the scheme follows the scaling from Eq. (11), with error decaying consistently with the proportionality $e_R \propto N^{-0.5}$. Although for the quasi-random distribution of particles the error magnitude increases, the results of this analysis underline the robustness of TFD with respect to the particle disorder, with fitted coefficients similar to the uniform case in all scenarios of anisotropy.

## 3.2. Transport through a heterogeneous medium

This test considers a two-dimensional aquifer with spatial variability of hydraulic properties (Fig. 9) composed of $N_c = 1500 \times 300$ unitary cells ($\Delta_x, \Delta_y$) = (1, 1)[m]. The natural logarithm of the hydraulic conductivity $Y(x)$ represents one individual realization of a sequential Gaussian simulation with zero mean and spherical variogram with correlation length $I_Y = 20\Delta_x$. Heterogeneities are controlled by log-conductivity variances $\sigma_Y^2 \in \{0.25, 1.00, 2.25\}$ through the relation $K(x) = \exp(\sigma_Y Y(x))$. For each hydraulic conductivity distribution, groundwater flow is obtained from solving on the numerical grid the steady-state Darcy's equation $\nabla \cdot \mathbf{q} = 0$, where $\mathbf{q} = -K(x)\nabla H$, and $H$ is the hydraulic head. Flow is induced by a unit-mean hydraulic head

<figure><img src="image_12.png" /><figcaption>**Fig. 9.** Initial condition for heterogeneous aquifer with high spatial variability of hydraulic conductivity. Groundwater flow is from left to right and reference streamlines are initially spaced by 66 flow cells.</figcaption></figure>

gradient along the x-axis, by prescribing the hydraulic head values at the upstream and downstream aquifer boundaries. North and south aquifer boundaries are defined as impermeable. Divergence-free flow is verified for all degrees of heterogeneity. As before, solute transport simulations are performed for different ratios of dispersion anisotropy, setting $\alpha_L = 0.05\Delta_x$ which leads to advection-dominated transport, with a field Péclet number of $Pe = I_Y/\alpha_L = 400$. A rectangular solute plume is initially released near the aquifer inlet with uniform concentration $C_0$. The plume size is $d_x \times d_y = 30\Delta_x \times 200\Delta_y$, with the lower left corner placed at $(x_0, y_0) = (15\Delta_x, 50\Delta_y)$. To minimize the singularity near extremely sharp gradients for the TFD scheme, a slightly smoothed version of the initial condition is used (as in Biriukov and Price, 2018, see Appendix B). All SPH simulations begin with a uniform distribution of particles. Groundwater flow velocities are known at the cell faces and linearly interpolated to the particle positions. Particles leaving the aquifer are reinjected upstream by resetting the x-coordinate to zero and defining the y-coordinate from the intersection of a random number $\xi_r \sim U(0, 1)$ with the normalized cumulative probability distribution of the inflow velocity, ensuring an entrance of particles that is statistically consistent with the non-uniform distribution of inflow velocities. Aquifer boundaries are surrounded by fixed SPH particles that compensate the loss of particle support while computing density near the borders. These fixed particles are defined as no-flux boundaries from the solute transport perspective (as in Wang et al., 2019). SPH results are compared with a high-resolution RWPT model, which in contrast, interprets particles as individual components of a solute cloud (e.g., Salamon et al., 2006). In RWPT, the number of particles determines the accuracy of concentrations. The individual particle mass for the RWPT simulations is $m_R = C_0V_M/N_R$, where $V_M$ is the volume occupied by the solute injection and $N_R = 10^8$ the number of particles. These parameters lead to a concentration resolution of $C_{min} = 6 \times 10^{-5}C_0$, estimated as the mass of one particle inside a flow-model cell. These parameters provide enough confidence to presume that the binning of solute particles is a reasonable proxy for the exact concentrations.

### 3.2.1. Particle resolution and peak concentrations

SPH simulations performed with different particle resolution reveal the influence of this parameter on the accuracy of peak concentrations. This is particularly relevant in cases of medium-to-low aquifer heterogeneity and high dispersion anisotropy, where strong concentration gradients can remain persistently throughout the simulation. Three degrees of particle resolution are considered, $N \in \{0.2M, 0.8M, 1.8M\}$, which translates into averaged densities of $\bar{\rho} \in \{0.4, 1.8, 4\}$ particles per flow-model cell, respectively. Herrera et al. (2009) reported results for a similar problem considering isotropic dispersion and a reference density of $\bar{\rho} = 8$ particles per cell. Notice that the particle resolution not only influences the accuracy of SPH schemes from a dispersion perspective (as in the homogeneous test case), but also determines the sampling of the groundwater flow velocities, the latter being independent from the selected dispersion scheme. To illustrate, the case of medium variability $\sigma_Y^2 = 1$, and the highest dispersion
<figure><img src="image_14.png" /><figcaption>**Fig. 10.** Cell concentration profiles for different particle resolution in the aquifer with medium heterogeneity $\sigma_v^2 = 1$ and dispersion anisotropy $\lambda_D = 100$, at time $T = 300$[d]. Horizontal sectioned line at $\bar{C} = 0.1$ marks symmetric log scale towards negative values with linear threshold $10^{-6}$. Vertical grid in profiles correspond with grid in panel (a), spaced by $30\Delta_y$. Arrows indicate peak concentration values at $y = 240\Delta_y$ for each particle resolution.</figcaption></figure>

anisotropy ($\lambda_D = 100$) is considered for discussion (Fig. 10). The particle resolution determines the accuracy of the concentration profile. The lowest-resolution case leads to oversmoothed concentrations in comparison to the reference result, regardless of the selected SPH solver ($N_{0.2M}$ in Fig. 10b, c, d). Increasing the particle resolution above one particle per cell recovers most of the qualitative aspects of the concentration profile ($N_{0.8M}$ and $N_{1.8M}$ in Fig. 10b, c, d), and from here the magnitude of the peak concentrations appears to be determined by the specific properties of each SPH scheme. An aspect to note from the SPH results is that transitions towards low concentrations are always less abrupt than the reference RWPT model, due to the intrinsically continuous nature of SPH. Notice, for instance, in the shaded areas of Figure 10b, c, d), that all SPH solvers predicted non-zero concentrations regardless of the particle resolution, while simultaneously displaying a good agreement with the reference peak concentration in the cases with higher resolution. The exception to this remark is the ER scheme, which displays negative values within said region (Fig. 10b).

### 3.2.2. On the occurrence of negative concentrations

In general, SPH solvers displayed a reasonably good agreement with the reference RWPT models for all scenarios of heterogeneity and dispersion anisotropy (e.g., Fig. 11; Fig. 12). This is in itself a promising result for the method, taking into account that RWPT is known to be a robust approach for simulating transport through heterogeneous systems (e.g., LaBolle et al., 1996; Lichtner et al., 2002; Fernández-Garcia et al., 2005; Salamon et al., 2006). Still, some simulations, particularly those with anisotropic dispersion, displayed negative concentrations. The ER scheme's instabilities generate patterns of negative concentrations throughout the entire domain, mostly surrounding the solute plume and regardless of heterogeneity (e.g., Fig. 11e, f; Fig. 12b). In the simulations with the highest dispersion anisotropy the method predicted negative concentrations of significant magnitude (as high as nearly $\bar{C} = -10^{-1}$) in regions where every other approach yielded positive values (e.g., Fig. 10b, Fig. 13c). ASPHAD produced concentrations without negative values for all scenarios of anisotropy and heterogeneity (e.g., Fig. 11h, i; Fig. 12c). Although the scheme displayed a high sensitivity to the particle disorder, specifically in scenarios with anisotropic dispersion, the results from this test suggest in a broad sense that the impact of this flaw might be minor for mild dispersion anisotropies and advection-dominated transport. The method still overpredicts peak concentrations for the scenarios of highest anisotropy (Fig. 11i, Fig. 13c), but with discrepancies nowhere close to the effects observed in the homogeneous test. Nevertheless, it is worth reminding that ASPHAD might not actually produce the intended degree of local dispersion anisotropy under particle disorder, as seen in Section 3.1.2. Concentrations obtained with TFD are in agreement with the reference models, with some particles displaying negative values in scenarios of anisotropy (e.g., Fig. 11l; Fig. 12d). However, their occurrence appears as isolated events, not influencing the main solute plume. The origin of negative values in the TFD method can be two-fold. In the first place, homogeneous simulations exhibited low-magnitude negative concentrations for the case with highest dispersion anisotropy, rooted in the propagation of numerical rounding errors. The

<figure><img src="image_15.png" /><figcaption>**Fig. 11.** Cell concentrations for low heterogeneity aquifer ($\sigma_v^2 = 0.25$). Columns group results by dispersion anisotropy $\lambda_D$, for the reference model and all SPH solvers. Results correspond to $T = 800$[d] with $N = 0.8M$ particles and $\gamma = 2$. Black cells indicate negative values. Dashed white line marks the profile shown in Fig. 13. Concentrations are shown until a minimum of $\bar{C} = 10^{-5}$.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Fig. 10.** Cell concentration profiles for different particle resolution in the aquifer with medium heterogeneity $\sigma_v^2 = 1$ and dispersion anisotropy $\lambda_D = 100$, at time $T = 300$[d]. Horizontal sectioned line at $\bar{C} = 0.1$ marks symmetric log scale towards negative values with linear threshold $10^{-6}$. Vertical grid in profiles correspond with grid in panel (a), spaced by $30\Delta_y$. Arrows indicate peak concentration values at $y = 240\Delta_y$ for each particle resolution.</figcaption></figure>

anisotropy ($\lambda_D = 100$) is considered for discussion (Fig. 10). The particle resolution determines the accuracy of the concentration profile. The lowest-resolution case leads to oversmoothed concentrations in comparison to the reference result, regardless of the selected SPH solver ($N_{0.2M}$ in Fig. 10b, c, d). Increasing the particle resolution above one particle per cell recovers most of the qualitative aspects of the concentration profile ($N_{0.8M}$ and $N_{1.8M}$ in Fig. 10b, c, d), and from here the magnitude of the peak concentrations appears to be determined by the specific properties of each SPH scheme. An aspect to note from the SPH results is that transitions towards low concentrations are always less abrupt than the reference RWPT model, due to the intrinsically continuous nature of SPH. Notice, for instance, in the shaded areas of Figure 10b, c, d), that all SPH solvers predicted non-zero concentrations regardless of the particle resolution, while simultaneously displaying a good agreement with the reference peak concentration in the cases with higher resolution. The exception to this remark is the ER scheme, which displays negative values within said region (Fig. 10b).

### 3.2.2. On the occurrence of negative concentrations

In general, SPH solvers displayed a reasonably good agreement with the reference RWPT models for all scenarios of heterogeneity and dispersion anisotropy (e.g., Fig. 11; Fig. 12). This is in itself a promising result for the method, taking into account that RWPT is known to be a robust approach for simulating transport through heterogeneous systems (e.g., LaBolle et al., 1996; Lichtner et al., 2002; Fernández-Garcia et al., 2005; Salamon et al., 2006). Still, some simulations, particularly those with anisotropic dispersion, displayed negative concentrations. The ER scheme's instabilities generate patterns of negative concentrations throughout the entire domain, mostly surrounding the solute plume and regardless of heterogeneity (e.g., Fig. 11e, f; Fig. 12b). In the simulations with the highest dispersion anisotropy the method predicted negative concentrations of significant magnitude (as high as nearly $\bar{C} = -10^{-1}$) in regions where every other approach yielded positive values (e.g., Fig. 10b, Fig. 13c). ASPHAD produced concentrations without negative values for all scenarios of anisotropy and heterogeneity (e.g., Fig. 11h, i; Fig. 12c). Although the scheme displayed a high sensitivity to the particle disorder, specifically in scenarios with anisotropic dispersion, the results from this test suggest in a broad sense that the impact of this flaw might be minor for mild dispersion anisotropies and advection-dominated transport. The method still overpredicts peak concentrations for the scenarios of highest anisotropy (Fig. 11i, Fig. 13c), but with discrepancies nowhere close to the effects observed in the homogeneous test. Nevertheless, it is worth reminding that ASPHAD might not actually produce the intended degree of local dispersion anisotropy under particle disorder, as seen in Section 3.1.2. Concentrations obtained with TFD are in agreement with the reference models, with some particles displaying negative values in scenarios of anisotropy (e.g., Fig. 11l; Fig. 12d). However, their occurrence appears as isolated events, not influencing the main solute plume. The origin of negative values in the TFD method can be two-fold. In the first place, homogeneous simulations exhibited low-magnitude negative concentrations for the case with highest dispersion anisotropy, rooted in the propagation of numerical rounding errors. The

<figure><img src="image_14.png" /><figcaption>**Fig. 11.** Cell concentrations for low heterogeneity aquifer ($\sigma_v^2 = 0.25$). Columns group results by dispersion anisotropy $\lambda_D$, for the reference model and all SPH solvers. Results correspond to $T = 800$[d] with $N = 0.8M$ particles and $\gamma = 2$. Black cells indicate negative values. Dashed white line marks the profile shown in Fig. 13. Concentrations are shown until a minimum of $\bar{C} = 10^{-5}$.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Fig. 12. Cell concentrations for high heterogeneity aquifer at time $T = 800[\text{d}]$. Dispersion anisotropy is $\lambda_D = 10$ and SPH models consider $N = 1.8\text{M}$ particles and $\gamma = 2$. Black cells indicate negative values. Concentrations are shown until a minimum of $\tilde{C} = 10^{-5}$.</figcaption></figure>

second source of negative concentration values is the known response to sharp concentration gradients. Depending on the medium heterogeneity, fluid stretching can rapidly enhance concentration gradients. A particle with high concentration and traveling on a fast streamline, could be suddenly surrounded by low concentration particles, creating the conditions for the singular behavior. This can occur independently of the dispersion anisotropy and some particles with negative values could be expected even in scenarios of isotropic dispersion. Anisotropy plays a role because it ultimately controls the magnitude of gradients in the transverse direction. In any case, this problem differs both by cause and by significance from the case of ER, where negative concentration values originate from a known unphysical behavior of the scheme. In simulations with this method considering anisotropic dispersion, approximately 40% of the particles presented negative values persistently in time. The TFD scheme is fundamentally stable for anisotropic dispersion, with a reported singularity near sharp gradients which needs further research for proper handling. Similar issues exist in other SPH applications ultimately requiring some form of dissipation to regulate the abrupt response (e.g., Monaghan and Gingold, 1983; Price, 2008). The medium heterogeneity modulates the sudden generation of strong concentration gradients eventually leading to the occurrence of isolated particles with negative values.

### 3.2.3. Potential for mixing-controlled reactive transport

A direct potential application for the SPH schemes assessed in this paper is reactive transport modeling. Under certain assumptions, a reactive transport problem can be conveniently reformulated as an equivalent conservative transport simulation (e.g., Saaltink et al., 1998; Cirpka and Valocchi, 2007; De Simoni et al., 2007). Reaction rates can then be locally calculated from the mixing rate of the conservative component and a term characterizing the chemical system. The domain-integrated mixing rate is also known as the scalar dissipation rate (Le Borgne et al., 2010; Engdahl et al., 2013),

$$ \chi(t) = -\frac{1}{2} \frac{d}{dt} \int_{\Omega} C^2 d\Omega = \int_{\Omega} \nabla^T C D \nabla C d\Omega, \quad (30) $$

<figure><img src="image_17.png" /><figcaption>Fig. 13. Grid concentration profiles at $x = 800\lambda_x$ for aquifer of low spatial variability ($\sigma_{\tilde{\nu}} = 0.25$), $N = 0.8\text{M}$ and $\gamma = 2$. Horizontal dashed line at $\tilde{C} = 0.1$ marks the beginning of a symmetric log-scale plot towards negative values with linear threshold $10^{-6}$. Scenario of dispersion anisotropy is indicated in each respective panel.</figcaption></figure>

which can be computed from the temporal derivative of the domain-integrated squared concentrations (left-hand side in Eq. (30)), or as the domain-integrated mixing rate (right-hand side in Eq. (30)). This index can be used as a proxy representation of the integrated reaction rates under mixing-limited transport conditions. The relation shown in Eq. (30) remains valid for cases with zero solute flux through the domain boundaries (refer to Le Borgne et al., 2010, for a discussion), as is the case in the previously shown heterogeneous simulations. The index $\chi(t)$ is employed here while comparing SPH results with the reference RWPT models, in order to assess from an integral perspective the potential of SPH for simulating mixing-controlled reactive transport in heterogeneous systems. A reference RWPT scalar dissipation rate is obtained for each combination of dispersion anisotropy and heterogeneity by numerically evaluating the temporal derivative of the domain-integrated squared concentrations. From the SPH results, the index is obtained by evaluating the domain-integrated mixing rate at different instants of the simulation. The scalar dissipation rate is discussed in non-dimensional form, $\tilde{\chi} = \chi\tau/\Omega$, where $\tau = I_y/\tilde{v}$ is the characteristic advection timescale of the porous medium, with $\tilde{v}$ the domain-averaged flow velocity. Concentrations in Eq. (30) are considered dimensionless. The results (Fig. 14) indicate that for isotropic dispersion, regardless of heterogeneity, all SPH models would correctly predict integrated reaction rates. Visible differences between schemes arise for anisotropic dispersion coefficients. In particular, the ER method provides accurate results up to the anisotropy ratio $\lambda_D = 10$, supporting the observation that for these cases concentration gradients are not strongly influenced by the instability. In contrast, the scheme deviates considerably from the reference when $\lambda_D = 100$, with overprediction for the case of low heterogeneity and underprediction for high heterogeneity. The AS-PHAD method consistently overpredicted the scalar dissipation rate for
<figure><img src="image_15.png" /><figcaption>Fig. 12. Cell concentrations for high heterogeneity aquifer at time $T = 800[d]$. Dispersion anisotropy is $\lambda_D = 10$ and SPH models consider $N = 1.8M$ particles and $\gamma = 2$. Black cells indicate negative values. Concentrations are shown until a minimum of $\tilde{C} = 10^{-5}$.</figcaption></figure>

second source of negative concentration values is the known response to sharp concentration gradients. Depending on the medium heterogeneity, fluid stretching can rapidly enhance concentration gradients. A particle with high concentration and traveling on a fast streamline, could be suddenly surrounded by low concentration particles, creating the conditions for the singular behavior. This can occur independently of the dispersion anisotropy and some particles with negative values could be expected even in scenarios of isotropic dispersion. Anisotropy plays a role because it ultimately controls the magnitude of gradients in the transverse direction. In any case, this problem differs both by cause and by significance from the case of ER, where negative concentration values originate from a known unphysical behavior of the scheme. In simulations with this method considering anisotropic dispersion, approximately 40% of the particles presented negative values persistently in time. The TFD scheme is fundamentally stable for anisotropic dispersion, with a reported singularity near sharp gradients which needs further research for proper handling. Similar issues exist in other SPH applications ultimately requiring some form of dissipation to regulate the abrupt response (e.g., Monaghan and Gingold, 1983; Price, 2008). The medium heterogeneity modulates the sudden generation of strong concentration gradients eventually leading to the occurrence of isolated particles with negative values.

### 3.2.3. Potential for mixing-controlled reactive transport

A direct potential application for the SPH schemes assessed in this paper is reactive transport modeling. Under certain assumptions, a reactive transport problem can be conveniently reformulated as an equivalent conservative transport simulation (e.g., Saaltink et al., 1998; Cirpka and Valocchi, 2007; De Simoni et al., 2007). Reaction rates can then be locally calculated from the mixing rate of the conservative component and a term characterizing the chemical system. The domain-integrated mixing rate is also known as the scalar dissipation rate (Le Borgne et al., 2010; Engdahl et al., 2013),

$$ \chi(t) = -\frac{1}{2} \frac{d}{dt} \int_{\Omega} C^2 d\Omega = \int_{\Omega} \nabla^T C D \nabla C d\Omega, \quad (30) $$

<figure><img src="image_16.png" /><figcaption>Fig. 13. Grid concentration profiles at $x = 800A_x$ for aquifer of low spatial variability ($\sigma_{\tilde{v}} = 0.25$), $N = 0.8M$ and $\gamma = 2$. Horizontal dashed line at $\tilde{C} = 0.1$ marks the beginning of a symmetric log-scale plot towards negative values with linear threshold $10^{-6}$. Scenario of dispersion anisotropy is indicated in each respective panel.</figcaption></figure>

which can be computed from the temporal derivative of the domain-integrated squared concentrations (left-hand side in Eq. (30)), or as the domain-integrated mixing rate (right-hand side in Eq. (30)). This index can be used as a proxy representation of the integrated reaction rates under mixing-limited transport conditions. The relation shown in Eq. (30) remains valid for cases with zero solute flux through the domain boundaries (refer to Le Borgne et al., 2010, for a discussion), as is the case in the previously shown heterogeneous simulations. The index $\chi(t)$ is employed here while comparing SPH results with the reference RWPT models, in order to assess from an integral perspective the potential of SPH for simulating mixing-controlled reactive transport in heterogeneous systems. A reference RWPT scalar dissipation rate is obtained for each combination of dispersion anisotropy and heterogeneity by numerically evaluating the temporal derivative of the domain-integrated squared concentrations. From the SPH results, the index is obtained by evaluating the domain-integrated mixing rate at different instants of the simulation. The scalar dissipation rate is discussed in non-dimensional form, $\tilde{\chi} = \chi\tau/\Omega$, where $\tau = I_y/\tilde{v}$ is the characteristic advection timescale of the porous medium, with $\tilde{v}$ the domain-averaged flow velocity. Concentrations in Eq. (30) are considered dimensionless. The results (Fig. 14) indicate that for isotropic dispersion, regardless of heterogeneity, all SPH models would correctly predict integrated reaction rates. Visible differences between schemes arise for anisotropic dispersion coefficients. In particular, the ER method provides accurate results up to the anisotropy ratio $\lambda_D = 10$, supporting the observation that for these cases concentration gradients are not strongly influenced by the instability. In contrast, the scheme deviates considerably from the reference when $\lambda_D = 100$, with overprediction for the case of low heterogeneity and underprediction for high heterogeneity. The AS-PHAD method consistently overpredicted the scalar dissipation rate for
<figure><img src="image_18.png" /><figcaption>**Fig. 14.** Comparison of scalar dissipation rates obtained from the reference RWPT model (solid lines) and the SPH results (scatter). Simulations were performed with $N = 0.8M$ particles and $\gamma = 2$. Each panel group the results for a given degree of heterogeneity and color codes group the degree of dispersion anisotropy.</figcaption></figure>

the highest anisotropy ($\lambda_D = 100$). This stems from the overestimation of concentration gradients under scenarios of particle disorder. The TFD method, on the other hand, consistently exhibits the best agreement with the reference curves for all of the simulation scenarios, underlining the robustness of this method when faced with different conditions of dispersion anisotropy and aquifer heterogeneity. Overall, and in agreement with results from previous sections as well, the TFD scheme appears as a suitable alternative for modeling mixing-limited reactive transport, still with some challenges and improvements to be addressed in future developments.

## 4. Conclusions

This article reviewed three different SPH dispersion schemes compatible with anisotropic dispersion coefficients, focusing on their applicability to simulate solute transport through heterogeneous porous media. The analysis was motivated by a recurrent discussion in the literature where negative concentrations are reported in conservative simulations under anisotropic dispersion. This result poses a potential problem for the application of SPH to solute transport in porous media, where hydrodynamic dispersion is intrinsically anisotropic. Furthermore, the accurate representation of anisotropy is necessary for reactive transport modeling purposes because transverse mixing has been shown to control the amount and extent of reactions. In one of the integration methods (ER) the occurrence of negative concentrations is explained by an unphysical solute transfer from particles of low concentration to particles of high concentration. This occurs for values of anisotropy above a relatively low threshold so it is likely for instabilities to

develop in most practical scenarios. The artifact can influence local concentration gradients near plume limits, particularly in scenarios of high anisotropy and medium-to-low aquifer heterogeneity, where large transverse concentration gradients are likely to be preserved throughout simulations. Two novel SPH formulations compatible with anisotropic dispersion coefficients were applied for the first time in heterogeneous porous media. A dispersion scheme based on a modified form of the isotropic SPH interpolator (ASPHAD) was shown to be strongly influenced by the particles' disorder, specifically in scenarios with anisotropy. This causes an overprediction of peak concentrations and underprediction of the longitudinal extent of solute plumes. This effect was highly visible in the homogeneous medium test. While simulating solute transport through a heterogeneous medium under advection-dominated conditions, in particular for scenarios of mild anisotropy, the shortcomings of the scheme were less obvious than in the homogeneous case when simply comparing concentration profiles, but the analysis of the domain-integrated mixing rate revealed that the overprediction of longitudinal concentration gradients produced a consistent overprediction of mixing, regardless of the degree of aquifer heterogeneity. Still, the scheme did not generate negative concentrations for any of the discussed conditions and it can potentially be considered for some porous media applications with mid-low dispersion anisotropy. A second alternative SPH interpolator based on a two-stage integration of dispersion (TFD) exhibits desirable numerical properties with a robust behavior when faced to the different conditions of particle distribution and dispersion anisotropy. The scheme always provided results in agreement with the reference solutions, although some challenges remain to be addressed. In particular, the method is known to be affected by a singularity near sharp concentration gradients. This has the potential to produce negative concentration values regardless of dispersion anisotropy. Effects of this kind exist in other SPH applications, and are usually handled by the introduction of artificial dissipation, which was not addressed in this study. Regardless, particles with negative concentration were rather isolated and of low magnitude, not influencing the main characteristics of the solute plume and not affecting the mixing rate calculations, which displayed accurate and robust results for the various conditions of anisotropy and heterogeneity. In summary, the results of this work contribute to the advance of SPH as an alternative for modeling solute transport with anisotropic dispersion through heterogeneous porous media.

### CRediT authorship contribution statement

**Rodrigo Pérez-Illanes:** Formal analysis, Software, Writing – original draft. **Guillem Sole-Mari:** Methodology, Writing – original draft. **Daniel Fernández-García:** Supervision, Writing – original draft.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Data availability

Data will be made available on request.

### Acknowledgments

To Escola Tècnica Superior d'Enginyers de Camins, Canals i Ports de Barcelona for granting computational resources in the cluster Titani, which contributed to the results presented in this article. The research leading to these results has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement no. 814066 (Managed Aquifer Recharge Solutions Training Network - MARSoluT) and project GRA-DIENT (reference no. PID2021-127911OB-I00) from the Ministry of Science and Innovation of Spain.
<figure><img src="image_17.png" /><figcaption>**Fig. 14.** Comparison of scalar dissipation rates obtained from the reference RWPT model (solid lines) and the SPH results (scatter). Simulations were performed with $N = 0.8M$ particles and $\gamma = 2$. Each panel group the results for a given degree of heterogeneity and color codes group the degree of dispersion anisotropy.</figcaption></figure>

the highest anisotropy ($\lambda_D$ = 100). This stems from the overestimation of concentration gradients under scenarios of particle disorder. The TFD method, on the other hand, consistently exhibits the best agreement with the reference curves for all of the simulation scenarios, underlining the robustness of this method when faced with different conditions of dispersion anisotropy and aquifer heterogeneity. Overall, and in agreement with results from previous sections as well, the TFD scheme appears as a suitable alternative for modeling mixing-limited reactive transport, still with some challenges and improvements to be addressed in future developments.

## 4. Conclusions

This article reviewed three different SPH dispersion schemes compatible with anisotropic dispersion coefficients, focusing on their applicability to simulate solute transport through heterogeneous porous media. The analysis was motivated by a recurrent discussion in the literature where negative concentrations are reported in conservative simulations under anisotropic dispersion. This result poses a potential problem for the application of SPH to solute transport in porous media, where hydrodynamic dispersion is intrinsically anisotropic. Furthermore, the accurate representation of anisotropy is necessary for reactive transport modeling purposes because transverse mixing has been shown to control the amount and extent of reactions. In one of the integration methods (ER) the occurrence of negative concentrations is explained by an unphysical solute transfer from particles of low concentration to particles of high concentration. This occurs for values of anisotropy above a relatively low threshold so it is likely for instabilities to

develop in most practical scenarios. The artifact can influence local concentration gradients near plume limits, particularly in scenarios of high anisotropy and medium-to-low aquifer heterogeneity, where large transverse concentration gradients are likely to be preserved throughout simulations. Two novel SPH formulations compatible with anisotropic dispersion coefficients were applied for the first time in heterogeneous porous media. A dispersion scheme based on a modified form of the isotropic SPH interpolator (ASPHAD) was shown to be strongly influenced by the particles' disorder, specifically in scenarios with anisotropy. This causes an overprediction of peak concentrations and underprediction of the longitudinal extent of solute plumes. This effect was highly visible in the homogeneous medium test. While simulating solute transport through a heterogeneous medium under advection-dominated conditions, in particular for scenarios of mild anisotropy, the shortcomings of the scheme were less obvious than in the homogeneous case when simply comparing concentration profiles, but the analysis of the domain-integrated mixing rate revealed that the overprediction of longitudinal concentration gradients produced a consistent overprediction of mixing, regardless of the degree of aquifer heterogeneity. Still, the scheme did not generate negative concentrations for any of the discussed conditions and it can potentially be considered for some porous media applications with mid-low dispersion anisotropy. A second alternative SPH interpolator based on a two-stage integration of dispersion (TFD) exhibits desirable numerical properties with a robust behavior when faced to the different conditions of particle distribution and dispersion anisotropy. The scheme always provided results in agreement with the reference solutions, although some challenges remain to be addressed. In particular, the method is known to be affected by a singularity near sharp concentration gradients. This has the potential to produce negative concentration values regardless of dispersion anisotropy. Effects of this kind exist in other SPH applications, and are usually handled by the introduction of artificial dissipation, which was not addressed in this study. Regardless, particles with negative concentration were rather isolated and of low magnitude, not influencing the main characteristics of the solute plume and not affecting the mixing rate calculations, which displayed accurate and robust results for the various conditions of anisotropy and heterogeneity. In summary, the results of this work contribute to the advance of SPH as an alternative for modeling solute transport with anisotropic dispersion through heterogeneous porous media.

### CRediT authorship contribution statement

**Rodrigo Pérez-Illanes:** Formal analysis, Software, Writing – original draft. **Guillem Sole-Mari:** Methodology, Writing – original draft. **Daniel Fernández-García:** Supervision, Writing – original draft.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Data availability

Data will be made available on request.

### Acknowledgments

To Escola Tècnica Superior d'Enginyers de Camins, Canals i Ports de Barcelona for granting computational resources in the cluster Titani, which contributed to the results presented in this article. The research leading to these results has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement no. 814066 (Managed Aquifer Recharge Solutions Training Network - MARSoluT) and project GRA-DIENT (reference no. PID2021-127911OB-I00) from the Ministry of Science and Innovation of Spain.
<figure><img src="image_19.png" /><figcaption>Fig. B.1. Concentrations from TFD one time step after injection for sharp (a) and smoothed (b) initial conditions. Dispersion anisotropy is $\lambda_D = 10$, $N = 0.8M$ and $\gamma = 2$. In panel (a) negative concentrations are shown in black. Minimum value is $\bar{C} = -6.3 \times 10^{-5}$.</figcaption></figure>

## Appendix A. Scaling relation for kernel size

The number of neighbor particles within the smoothing length distance is given by (Dehnen and Aly, 2012)

$$ \mathcal{N}_h = V_v h^\nu \left(\frac{\rho}{m}\right) \qquad (A.1) $$

where $V_v$ is the volume of the unit sphere in $v$ dimensions and the mean particle size satisfies $\delta_r = (m/\rho)^{1/\nu}$, being $m$ the particle mass and $\rho$ the SPH density estimate. Zhu et al. (2015) introduced the proportionality $h \propto N^{-1/\beta}$, with $\beta > 0$ a parameter that determines the rate at which the smoothing length decreases as the total number of particles increases. Further, by considering that the ratio between the neighbor particles and $h^\nu$ is proportional to the ratio of total particles per domain volume, and introducing the aforementioned scaling for the smoothing length, it is obtained

$$ \mathcal{N} \propto N^{1-\nu/\beta} \qquad (A.2) $$

where $\beta$ should satisfy $\beta > \nu$ in order to increase the number of neighbors as the total number of particles increases. Considering that both the continuous and discrete SPH interpolation errors are approximately of the same order and introducing (A.2) yields $\beta \sim \nu + 2/\psi$, where $\psi$ characterizes the influence of the particle distribution on the discrete interpolation error (Zhu et al., 2015). For the best case scenario of uniform particle distribution $\psi \approx 1$, which leads to the limit condition $\beta \ge \nu + 2$. Introducing a parameter $\mathcal{K}$ to the scaling in Eq. (A.2) yields

$$ \mathcal{N} = \mathcal{K} N^{1-\nu/\beta}, \qquad (A.3) $$

where neither the value of $\mathcal{K}$ nor its exact dependence is known beforehand. Using Eq. (A.1), the parameter $\gamma = h/\delta_r$ is obtained as

$$ \gamma = h \left(\frac{\rho}{m}\right)^{1/\nu} = \left(\frac{\mathcal{N}_h}{V_v}\right)^{1/\nu}, \qquad (A.4) $$

which can be expressed in terms of $\mathcal{N}$ by considering that $\mathcal{N}_h = (h/H)^\nu \mathcal{N}$ (Dehnen and Aly, 2012). This yields

$$ \gamma = \frac{h}{H} \left( \frac{\mathcal{N}}{V_v} \right)^{1/\nu}. \qquad (A.5) $$

Introducing Eq. (A.3), it is finally obtained

$$ \gamma(N) = \frac{1}{\kappa_v} \left(\frac{\mathcal{K}}{V_v}\right)^{1/\nu} N^{1-\nu-1/\beta}, \qquad (A.6) $$

where the quotient between the kernel characteristic distances has been replaced by the inherent kernel aspect ratio $\kappa_v$.

## Appendix B. Smoothed initial condition

The rectangular initial condition for the heterogeneous test case is smoothed in both directions as the product of error functions

$$ \bar{C}(\mathbf{x}, t_0) = \frac{1}{4} \left[ \operatorname{erf}\left(\frac{\mathbf{x}}{l_x}\right) - \operatorname{erf}\left(\frac{\mathbf{x} - d_x}{l_x}\right) \right] \cdot \\ \left[ \operatorname{erf}\left(\frac{\mathbf{y}}{l_y}\right) - \operatorname{erf}\left(\frac{\mathbf{y} - d_y}{l_y}\right) \right], \qquad (B.1) $$

where $l_x, l_y$ are smoothing distances in the $x$ and $y$ directions, respectively. An example of the singularity near sharp gradients of the TFD scheme is shown in Fig. B.1.

## References

* Alvarado-Rodríguez, C.E., Sigalotti, L.D.G., Klapp, J., 2019. Anisotropic dispersion with a consistent smoothed particle hydrodynamics scheme. Adv. Water Resour. 131, 103374. http://dx.doi.org/10.1016/j.advwatres.2019.07.004.

* Avesani, D., Dumbser, M., Chiogna, G., Bellin, A., 2016. An alternative smooth particle hydrodynamics formulation to simulate chemotaxis in porous media. J. Math. Biol. 74 (5), 1037–1058. http://dx.doi.org/10.1007/s00285-016-1049-6.

* Avesani, D., Herrera, P., Chiogna, G., Bellin, A., Dumbser, M., 2015. Smooth particle hydrodynamics with nonlinear moving-least-squares WENO reconstruction to model anisotropic dispersion in porous media. Adv. Water Resour. 80, 43–59. http://dx.doi.org/10.1016/j.advwatres.2015.03.007.

* Bandara, U., Tartakovsky, A., Oostrom, M., Palmer, B., Grate, J., Zhang, C., 2013. Smoothed particle hydrodynamics pore-scale simulations of unstable immiscible flow in porous media. Adv. Water Resour. 62, 356–369. http://dx.doi.org/10.1016/j.advwatres.2013.09.014.

* Basser, H., Rudman, M., Daly, E., 2017. SPH modelling of multi-fluid lock-exchange over and within porous media. Adv. Water Resour. 108, 15–28. http://dx.doi.org/10.1016/j.advwatres.2017.07.011.

* Basser, H., Rudman, M., Daly, E., 2019. Smoothed particle hydrodynamics modelling of fresh and salt water dynamics in porous media. J. Hydrol. 576, 370–380. http://dx.doi.org/10.1016/j.jhydrol.2019.06.048.

* Bear, J., Cheng, A.H.D., 2010. Modeling Groundwater Flow and Contaminant Transport. Springer Netherlands, http://dx.doi.org/10.1007/978-1-4020-6682-5.

* Benson, D.A., Aquino, T., Bolster, D., Engdahl, N., Henri, C.V., Fernández-Garcia, D., 2017. A comparison of Eulerian and Lagrangian transport and non-linear reaction algorithms. Adv. Water Resour. 99, 15–37. http://dx.doi.org/10.1016/j.advwatres.2016.11.003.

* Biriukov, S., Price, D.J., 2018. Stable anisotropic heat conduction in smoothed particle hydrodynamics. Mon. Not. R. Astron. Soc. 483 (4), 4901–4909. http://dx.doi.org/10.1093/mnras/sty3413.

* Boso, F., Bellin, A., Dumbser, M., 2013. Numerical simulations of solute transport in highly heterogeneous formations: A comparison of alternative numerical schemes. Adv. Water Resour. 52, 178–189. http://dx.doi.org/10.1016/j.advwatres.2012.08.006.

* Brookshaw, L., 1985. A method of calculating radiative heat diffusion in particle simulations. Proc. Astron. Soc. Aust. 6 (2), 207–210. http://dx.doi.org/10.1017/S1323358000018117.

* Cirkpa, O.A., Frind, E.O., Helmig, R., 1999. Numerical simulation of biodegradation controlled by transverse mixing. J. Contam. Hydrol. 40 (2), 159–182. http://dx.doi.org/10.1016/s0169-7722(99)00044-3.

* Cirkpa, O.A., Valocchi, A.J., 2007. Two-dimensional concentration distribution for mixing-controlled bioreactive transport in steady state. Adv. Water Resour. 30 (6–7), 1668–1679. http://dx.doi.org/10.1016/j.advwatres.2006.05.022.

* Cleary, P.W., Monaghan, J.J., 1999. Conduction modelling using smoothed particle hydrodynamics. J. Comput. Phys. 148 (1), 227–264. http://dx.doi.org/10.1006/jcph.1998.6118.

* De Simoni, M., Sanchez-Vila, X., Carrera, J., Saaltink, M.W., 2007. A mixing ratios-based formulation for multicomponent reactive transport. Water Resour. Res. 43 (7), http://dx.doi.org/10.1029/2006wr005256.
<figure><img src="image_18.png" /><figcaption>Fig. B.1. Concentrations from TFD one time step after injection for sharp (a) and smoothed (b) initial conditions. Dispersion anisotropy is $\lambda_D = 10$, $N = 0.8M$ and $\gamma = 2$. In panel (a) negative concentrations are shown in black. Minimum value is $\bar{C} = -6.3 \times 10^{-5}$.</figcaption></figure>

## Appendix A. Scaling relation for kernel size

The number of neighbor particles within the smoothing length distance is given by (Dehnen and Aly, 2012)

$$ \mathcal{N}_h = V_v h^\nu \left(\frac{\rho}{m}\right) \qquad (A.1) $$

where $V_v$ is the volume of the unit sphere in $v$ dimensions and the mean particle size satisfies $\delta_r = (m/\rho)^{1/\nu}$, being $m$ the particle mass and $\rho$ the SPH density estimate. Zhu et al. (2015) introduced the proportionality $h \propto N^{-1/\beta}$, with $\beta > 0$ a parameter that determines the rate at which the smoothing length decreases as the total number of particles increases. Further, by considering that the ratio between the neighbor particles and $h^\nu$ is proportional to the ratio of total particles per domain volume, and introducing the aforementioned scaling for the smoothing length, it is obtained

$$ \mathcal{N} \propto N^{1-\nu/\beta} \qquad (A.2) $$

where $\beta$ should satisfy $\beta > \nu$ in order to increase the number of neighbors as the total number of particles increases. Considering that both the continuous and discrete SPH interpolation errors are approximately of the same order and introducing (A.2) yields $\beta \sim \nu + 2/\psi$, where $\psi$ characterizes the influence of the particle distribution on the discrete interpolation error (Zhu et al., 2015). For the best case scenario of uniform particle distribution $\psi \approx 1$, which leads to the limit condition $\beta \ge \nu + 2$. Introducing a parameter $\mathcal{K}$ to the scaling in Eq. (A.2) yields

$$ \mathcal{N} = \mathcal{K} N^{1-\nu/\beta}, \qquad (A.3) $$

where neither the value of $\mathcal{K}$ nor its exact dependence is known beforehand. Using Eq. (A.1), the parameter $\gamma = h/\delta_r$ is obtained as

$$ \gamma = h \left(\frac{\rho}{m}\right)^{1/\nu} = \left(\frac{\mathcal{N}_h}{V_v}\right)^{1/\nu}, \qquad (A.4) $$

which can be expressed in terms of $\mathcal{N}$ by considering that $\mathcal{N}_h = (h/H)^\nu \mathcal{N}$ (Dehnen and Aly, 2012). This yields

$$ \gamma = \frac{h}{H} \left( \frac{\mathcal{N}}{V_v} \right)^{1/\nu}. \qquad (A.5) $$

Introducing Eq. (A.3), it is finally obtained

$$ \gamma(N) = \frac{1}{\kappa_v} \left(\frac{\mathcal{K}}{V_v}\right)^{1/\nu} N^{1-\nu-1/\beta}, \qquad (A.6) $$

where the quotient between the kernel characteristic distances has been replaced by the inherent kernel aspect ratio $\kappa_v$.

## Appendix B. Smoothed initial condition

The rectangular initial condition for the heterogeneous test case is smoothed in both directions as the product of error functions

$$ \bar{C}(\mathbf{x}, t_0) = \frac{1}{4} \left[ \operatorname{erf}\left(\frac{\mathbf{x}}{l_x}\right) - \operatorname{erf}\left(\frac{\mathbf{x} - d_x}{l_x}\right) \right] \cdot \\ \left[ \operatorname{erf}\left(\frac{\mathbf{y}}{l_y}\right) - \operatorname{erf}\left(\frac{\mathbf{y} - d_y}{l_y}\right) \right], \qquad (B.1) $$

where $l_x, l_y$ are smoothing distances in the $x$ and $y$ directions, respectively. An example of the singularity near sharp gradients of the TFD scheme is shown in Fig. B.1.

## References

* Alvarado-Rodríguez, C.E., Sigalotti, L.D.G., Klapp, J., 2019. Anisotropic dispersion with a consistent smoothed particle hydrodynamics scheme. Adv. Water Resour. 131, 103374. http://dx.doi.org/10.1016/j.advwatres.2019.07.004.

* Avesani, D., Dumbser, M., Chiogna, G., Bellin, A., 2016. An alternative smooth particle hydrodynamics formulation to simulate chemotaxis in porous media. J. Math. Biol. 74 (5), 1037–1058. http://dx.doi.org/10.1007/s00285-016-1049-6.

* Avesani, D., Herrera, P., Chiogna, G., Bellin, A., Dumbser, M., 2015. Smooth particle hydrodynamics with nonlinear moving-least-squares WENO reconstruction to model anisotropic dispersion in porous media. Adv. Water Resour. 80, 43–59. http://dx.doi.org/10.1016/j.advwatres.2015.03.007.

* Bandara, U., Tartakovsky, A., Oostrom, M., Palmer, B., Grate, J., Zhang, C., 2013. Smoothed particle hydrodynamics pore-scale simulations of unstable immiscible flow in porous media. Adv. Water Resour. 62, 356–369. http://dx.doi.org/10.1016/j.advwatres.2013.09.014.

* Basser, H., Rudman, M., Daly, E., 2017. SPH modelling of multi-fluid lock-exchange over and within porous media. Adv. Water Resour. 108, 15–28. http://dx.doi.org/10.1016/j.advwatres.2017.07.011.

* Basser, H., Rudman, M., Daly, E., 2019. Smoothed particle hydrodynamics modelling of fresh and salt water dynamics in porous media. J. Hydrol. 576, 370–380. http://dx.doi.org/10.1016/j.jhydrol.2019.06.048.

* Bear, J., Cheng, A.H.D., 2010. Modeling Groundwater Flow and Contaminant Transport. Springer Netherlands, http://dx.doi.org/10.1007/978-1-4020-6682-5.

* Benson, D.A., Aquino, T., Bolster, D., Engdahl, N., Henri, C.V., Fernández-Garcia, D., 2017. A comparison of Eulerian and Lagrangian transport and non-linear reaction algorithms. Adv. Water Resour. 99, 15–37. http://dx.doi.org/10.1016/j.advwatres.2016.11.003.

* Biriukov, S., Price, D.J., 2018. Stable anisotropic heat conduction in smoothed particle hydrodynamics. Mon. Not. R. Astron. Soc. 483 (4), 4901–4909. http://dx.doi.org/10.1093/mnras/sty3413.

* Boso, F., Bellin, A., Dumbser, M., 2013. Numerical simulations of solute transport in highly heterogeneous formations: A comparison of alternative numerical schemes. Adv. Water Resour. 52, 178–189. http://dx.doi.org/10.1016/j.advwatres.2012.08.006.

* Brookshaw, L., 1985. A method of calculating radiative heat diffusion in particle simulations. Proc. Astron. Soc. Aust. 6 (2), 207–210. http://dx.doi.org/10.1017/S1323358000018117.

* Cirkpa, O.A., Frind, E.O., Helmig, R., 1999. Numerical simulation of biodegradation controlled by transverse mixing. J. Contam. Hydrol. 40 (2), 159–182. http://dx.doi.org/10.1016/s0169-7722(99)00044-3.

* Cirkpa, O.A., Valocchi, A.J., 2007. Two-dimensional concentration distribution for mixing-controlled bioreactive transport in steady state. Adv. Water Resour. 30 (6–7), 1668–1679. http://dx.doi.org/10.1016/j.advwatres.2006.05.022.

* Cleary, P.W., Monaghan, J.J., 1999. Conduction modelling using smoothed particle hydrodynamics. J. Comput. Phys. 148 (1), 227–264. http://dx.doi.org/10.1006/jcph.1998.6118.

* De Simoni, M., Sanchez-Vila, X., Carrera, J., Saaltink, M.W., 2007. A mixing ratios-based formulation for multicomponent reactive transport. Water Resour. Res. 43 (7), http://dx.doi.org/10.1029/2006wr005256.
Dehnen, W., Aly, H., 2012. Improving convergence in smoothed particle hydrodynamics simulations without pairing instability. Mon. Not. R. Astron. Soc. 425 (2), 1068-1082. http://dx.doi.org/10.1111/j.1365-2966.2012.21439.x

Engdahl, N.B., Ginn, T.R., Fogg, G.E., 2013. Scalar dissipation rates in non-conservative transport systems. J. Contam. Hydrol. 149, 46-60. http://dx.doi.org/10.1016/j.jconhyd.2013.03.003.

Español, P., Revenga, M., 2003. Smoothed dissipative particle dynamics. Phys. Rev. E 67 (2), http://dx.doi.org/10.1103/physreve.67.026705.

Fernández-Garcia, D., Illangasekare, T.H., Rajaram, H., 2005. Differences in the scale-dependence of dispersivity estimated from temporal and spatial moments in chemically and physically heterogeneous porous media. Adv. Water Resour. 28 (7), 745-759. http://dx.doi.org/10.1016/j.adwatres.2004.12.011.

Fernández-Garcia, D., Sanchez-Vila, X., 2011. Optimal reconstruction of concentrations, gradients and reaction rates from particle distributions. J. Contam. Hydrol. 120-121, 99-114. http://dx.doi.org/10.1016/j.jconhyd.2010.05.001, URL: https://www.sciencedirect.com/science/article/pii/S0169772210000483.

Fernández-Garcia, D., Sanchez-Vila, X., Guadagnini, A., 2008. Reaction rates and effective parameters in stratified aquifers. Adv. Water Resour. 31 (10), 1364-1376. http://dx.doi.org/10.1016/j.adwatres.2008.07.001.

Gelhar, L.W., Welty, C., Rehfeldt, K.R., 1992. A critical review of data on field-scale dispersion in aquifers. Water Resour. Res. 28 (7), 1955-1974. http://dx.doi.org/10.1029/92WR00607, URL: https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/92WR00607.

Gingold, R.A., Monaghan, J.J., 1977. Smoothed particle hydrodynamics: theory and application to non-spherical stars. Mon. Not. R. Astron. Soc. 181 (3), 375-389. http://dx.doi.org/10.1093/mnras/181.3.375.

Herrera, P.A., Beckie, R.D., 2013. An assessment of particle methods for approximating anisotropic dispersion. Internat. J. Numer. Methods Fluids 71 (5), 634-651. http://dx.doi.org/10.1002/fld.3676.

Herrera, P.A., Massabó, M., Beckie, R.D., 2009. A meshless method to simulate solute transport in heterogeneous porous media. Adv. Water Resour. 32 (3), 413-429. http://dx.doi.org/10.1016/j.adwatres.2008.12.005.

Hochstetler, D.L., Rolle, M., Chiogna, G., Haberer, C.M., Grathwohl, P., Kitanidis, P.K., 2013. Effects of compound-specific transverse mixing on steady-state reactive plumes: Insights from pore-scale simulations and Darcy-scale experiments. Adv. Water Resour. 54, 1-10. http://dx.doi.org/10.1016/j.adwatres.2012.12.007.

Kitanidis, P.K., 1994. The concept of the dilution index. Water Resour. Res. 30 (7), 2011-2026. http://dx.doi.org/10.1029/94wr00762.

Klapp, J., Sigalotti, L.D.G., Alvarado-Rodríguez, C.E., Rendón, O., Díaz-Damacillo, L., 2022. Approximately consistent SPH simulations of the anisotropic dispersion of a contaminant plume. Comput. Part. Mech. 9 (5), 987-1002. http://dx.doi.org/10.1007/s40571-022-00461-1.

Klenk, I., Grathwohl, P., 2002. Transverse vertical dispersion in groundwater and the capillary fringe. J. Contam. Hydrol. 58 (1-2), 111-128. http://dx.doi.org/10.1016/j.adwatres.2010.08.006.

Kunz, P., Zarikos, I.M., Karadimitriou, N.K., Huber, M., Nieken, U., Hassanizadeh, S.M., 2015. Study of multi-phase flow in porous media: Comparison of SPH simulations with micro-model experiments. Transp. Porous Media 114 (2), 581-600. http://dx.doi.org/10.1007/s11242-015-0599-1.

LaBolle, E.M., Fogg, G.E., Tompson, A.F.B., 1996. Random-walk simulation of transport in heterogeneous porous media: Local mass-conservation problem and implementation methods. Water Resour. Res. 32 (3), 583-593. http://dx.doi.org/10.1029/95wr03528.

Le Borgne, T., Dentz, M., Bolster, D., Carrera, J., de Dreuzy, J.-R., Davy, P., 2010. Non-Fickian mixing: Temporal evolution of the scalar dissipation rate in heterogeneous porous media. Adv. Water Resour. 33 (12), 1468-1475. http://dx.doi.org/10.1016/j.adwatres.2010.08.006.

Lichtner, P.C., Kelkar, S., Robinson, B., 2002. New form of dispersion tensor for axisymmetric porous media with implementation in particle tracking. Water Resour. Res. 38 (8), 21-21-21-16. http://dx.doi.org/10.1029/2000wr000100.

Lucy, L., 1977. A numerical approach to the testing of the fission hypothesis. Astron. J. 82, 1013. http://dx.doi.org/10.1086/112164.

Monaghan, J., 2005. Smoothed particle hydrodynamics. Rep. Progr. Phys. 68 (8), 1703-1759. http://dx.doi.org/10.1088/0034-4885/68/8/r01.

Monaghan, J., 2012. Smoothed particle hydrodynamics and its diverse applications. Annu. Rev. Fluid Mech. 44 (1), 323-346. http://dx.doi.org/10.1146/annurev-fluid-120710-101220.

Monaghan, J.J., Gingold, R.A., 1983. Shock simulation by the particle method SPH. J. Comput. Phys. 52 (2), 374-389. http://dx.doi.org/10.1016/0021-9991(83)90036-0.

Pickens, J.F., Grisak, G.E., 1981. Scale-dependent dispersion in a stratified granular aquifer. Water Resour. Res. 17 (4), 1191-1211. http://dx.doi.org/10.1029/wr017i004p01191.

Price, D., 2008. Modelling discontinuities and Kelvin-Helmholtz instabilities in SPH. J. Comput. Phys. 227 (24), 10040-10057. http://dx.doi.org/10.1016/j.jcp.2008.08.011.

Price, D.J., Monaghan, J.J., 2007. An energy-conserving formalism for adaptive gravitational force softening in smoothed particle hydrodynamics and N-body codes. Mon. Not. R. Astron. Soc. 374 (4), 1347-1358. http://dx.doi.org/10.1111/j.1365-2966.2006.11241.x.

Rolle, M., Eberhardt, C., Chiogna, G., Cirpka, O.A., Grathwohl, P., 2009. Enhancement of dilution and transverse reactive mixing in porous media: Experiments and model-based interpretation. J. Contam. Hydrol. 110 (3-4), 130-142. http://dx.doi.org/10.1016/j.jconhyd.2009.10.003.

Saaltink, M.W., Ayora, C., Carrera, J., 1998. A mathematical formulation for reactive transport that eliminates mineral concentrations. Water Resour. Res. 34 (7), 1649-1656. http://dx.doi.org/10.1029/98wr00552.

Salamon, P., Fernández-Garcia, D., Gómez-Hernández, J.J., 2006. A review and numerical assessment of the random walk particle tracking method. J. Contam. Hydrol. 87 (3-4), 277-305. http://dx.doi.org/10.1016/j.jconhyd.2006.05.005.

Schoenberg, I., 1946. Contributions to the problem of approximation of equidistant data by analytic functions. Part a: On the problem of smoothing or graduation. A first class of analytic approximation formulae. Quart. Appl. Math. 4 (2), 45. http://dx.doi.org/10.1007/978-1-4899-0433-1_1.

Sigalotti, L.D.G., Klapp, J., Rendón, O., Vargas, C.A., Peña-Polo, F., 2016. On the kernel and particle consistency in smoothed particle hydrodynamics. Appl. Numer. Math. 108, 242-255. http://dx.doi.org/10.1016/j.apnum.2016.05.007.

Sigalotti, L.D.G., Rendón, O., Klapp, J., Vargas, C.A., Cruz, F., 2019. A new insight into the consistency of the SPH interpolation formula. Appl. Math. Comput. 356, 50-73. http://dx.doi.org/10.1016/j.amc.2019.03.018.

Silliman, S.E., Simpson, E.S., 1987. Laboratory evidence of the scale effect in dispersion of solutes in porous media. Water Resour. Res. 23 (8), 1667-1673. http://dx.doi.org/10.1029/wr023i008p01667.

Springel, V., 2010. Smoothed particle hydrodynamics in astrophysics. Annu. Rev. Astron. Astrophys. 48 (1), 391-430. http://dx.doi.org/10.1146/annurev-astro-081309-130914.

Tartakovsky, A.M., Meakin, P., Scheibe, T.D., Wood, B.D., 2007. A smoothed particle hydrodynamics model for reactive transport and mineral precipitation in porous and fractured porous media. Water Resour. Res. 43 (5), http://dx.doi.org/10.1029/2005wr004770.

Tartakovsky, A.M., Redden, G., Lichtner, P.C., Scheibe, T.D., Meakin, P., 2008. Mixed-induced precipitation: Experimental study and multiscale numerical analysis. Water Resour. Res. 44 (6), http://dx.doi.org/10.1029/2006wr005725.

Tartakovsky, A.M., Scheibe, T.D., Meakin, P., 2009. Pore-scale model for reactive transport and biomass growth. J. Porous Media 12 (5), 417-434. http://dx.doi.org/10.1615/jpormedia.v12.i5.30.

Tartakovsky, A.M., Trask, N., Pan, K., Jones, B., Pan, W., Williams, J.R., 2015. Smoothed particle hydrodynamics and its applications for multiphase flow and reactive transport in porous media. Comput. Geosci. 20 (4), 807-834. http://dx.doi.org/10.1007/s10596-015-9468-9.

Tran-Duc, T., Bertevas, E., Phan-Thien, N., Khoo, B.C., 2016. Simulation of anisotropic diffusion processes in fluids with smoothed particle hydrodynamics. Internat. J. Numer. Methods Fluids 82 (11), 730-747. http://dx.doi.org/10.1002/fld.4238.

Vacondió, R., Altomare, C., Leffe, M.D., Hu, X., Touzé, D.L., Lind, S., Marongiu, J.-C., Marrone, S., Rogers, B.D., Souto-Iglesias, A., 2021. Grand challenges for smoothed particle hydrodynamics numerical schemes. Comput. Part. Mech. 8 (3), 575-588. http://dx.doi.org/10.1007/s40571-020-00354-1.

Vacondió, R., Rogers, B.D., Stansby, P.K., Mignosa, P., 2012. SPH modeling of shallow flow with open boundaries for practical flood simulation. J. Hydraul. Eng. 138 (6), 530-541. http://dx.doi.org/10.1061/(asce)hy.1943-7900.0000543.

Violeau, D., 2012. Fluid Mechanics and the SPH Method: Theory and Applications. Oxford University Press, Oxford.

Wang, J., Hu, W., Zhang, X., Pan, W., 2019. Modeling heat transfer subject to inhomogeneous Neumann boundary conditions by smoothed particle hydrodynamics and peridynamics. Int. J. Heat Mass Transfer 139, 948-962. http://dx.doi.org/10.1016/j.ijheatmasstransfer.2019.05.054.

Ye, T., Pan, D., Huang, C., Liu, M., 2019. Smoothed particle hydrodynamics (SPH) for complex fluid flows: Recent developments in methodology and applications. Phys. Fluids 31 (1), http://dx.doi.org/10.1063/1.5068697.

Zhu, Y., Fox, P.J., 2001. Smoothed particle hydrodynamics model for diffusion through porous media. Transp. Porous Media 43 (3), 441-471. http://dx.doi.org/10.1023/a:1010769915901.

Zhu, Y., Fox, P.J., 2002. Simulation of pore-scale dispersion in periodic porous media using smoothed particle hydrodynamics. J. Comput. Phys. 182 (2), 622-645. http://dx.doi.org/10.1006/japh.2002.7189.

Zhu, Q., Hernquist, L., Li, Y., 2015. Numerical convergence in smoothed particle hydrodynamics. Astrophys. J. 800 (1), 6. http://dx.doi.org/10.1088/0004-637x/800/1/6.
Dehnen, W., Aly, H., 2012. Improving convergence in smoothed particle hydrodynamics simulations without pairing instability. Mon. Not. R. Astron. Soc. 425 (2), 1068-1082. http://dx.doi.org/10.1111/j.1365-2966.2012.21439.x

Engdahl, N.B., Ginn, T.R., Fogg, G.E., 2013. Scalar dissipation rates in non-conservative transport systems. J. Contam. Hydrol. 149, 46-60. http://dx.doi.org/10.1016/j.jconhyd.2013.03.003.

Español, P., Revenga, M., 2003. Smoothed dissipative particle dynamics. Phys. Rev. E 67 (2), http://dx.doi.org/10.1103/physreve.67.026705.

Fernández-García, D., Illangasekare, T.H., Rajaram, H., 2005. Differences in the scale-dependence of dispersivity estimated from temporal and spatial moments in chemically and physically heterogeneous porous media. Adv. Water Resour. 28 (7), 745-759. http://dx.doi.org/10.1016/j.adwatres.2004.12.011.

Fernández-García, D., Sanchez-Vila, X., 2011. Optimal reconstruction of concentrations, gradients and reaction rates from particle distributions. J. Contam. Hydrol. 120-121, 99-114. http://dx.doi.org/10.1016/j.jconhyd.2010.05.001, URL: https://www.sciencedirect.com/science/article/pii/S0169772210000483.

Fernández-García, D., Sanchez-Vila, X., Guadagnini, A., 2008. Reaction rates and effective parameters in stratified aquifers. Adv. Water Resour. 31 (10), 1364-1376. http://dx.doi.org/10.1016/j.adwatres.2008.07.001.

Gelhar, L.W., Welty, C., Rehfeldt, K.R., 1992. A critical review of data on field-scale dispersion in aquifers. Water Resour. Res. 28 (7), 1955-1974. http://dx.doi.org/10.1029/92WR00607, URL: https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/92WR00607.

Gingold, R.A., Monaghan, J.J., 1977. Smoothed particle hydrodynamics: theory and application to non-spherical stars. Mon. Not. R. Astron. Soc. 181 (3), 375-389. http://dx.doi.org/10.1093/mnras/181.3.375.

Herrera, P.A., Beckie, R.D., 2013. An assessment of particle methods for approximating anisotropic dispersion. Internat. J. Numer. Methods Fluids 71 (5), 634-651. http://dx.doi.org/10.1002/fld.3676.

Herrera, P.A., Massabó, M., Beckie, R.D., 2009. A meshless method to simulate solute transport in heterogeneous porous media. Adv. Water Resour. 32 (3), 413-429. http://dx.doi.org/10.1016/j.adwatres.2008.12.005.

Hochstetler, D.L., Rolle, M., Chiogna, G., Haberer, C.M., Grathwohl, P., Kitanidis, P.K., 2013. Effects of compound-specific transverse mixing on steady-state reactive plumes: Insights from pore-scale simulations and Darcy-scale experiments. Adv. Water Resour. 54, 1-10. http://dx.doi.org/10.1016/j.adwatres.2012.12.007.

Kitanidis, P.K., 1994. The concept of the dilution index. Water Resour. Res. 30 (7), 2011-2026. http://dx.doi.org/10.1029/94wr00762.

Klapp, J., Sigalotti, L.D.G., Alvarado-Rodríguez, C.E., Rendón, O., Díaz-Damacillo, L., 2022. Approximately consistent SPH simulations of the anisotropic dispersion of a contaminant plume. Comput. Part. Mech. 9 (5), 987-1002. http://dx.doi.org/10.1007/s40571-022-00461-1.

Klenk, I., Grathwohl, P., 2002. Transverse vertical dispersion in groundwater and the capillary fringe. J. Contam. Hydrol. 58 (1-2), 111-128. http://dx.doi.org/10.1016/j.adwatres.2010.08.006.

Kunz, P., Zarikos, I.M., Karadimitriou, N.K., Huber, M., Nieken, U., Hassanizadeh, S.M., 2015. Study of multi-phase flow in porous media: Comparison of SPH simulations with micro-model experiments. Transp. Porous Media 114 (2), 581-600. http://dx.doi.org/10.1007/s11242-015-0599-1.

LaBolle, E.M., Fogg, G.E., Tompson, A.F.B., 1996. Random-walk simulation of transport in heterogeneous porous media: Local mass-conservation problem and implementation methods. Water Resour. Res. 32 (3), 583-593. http://dx.doi.org/10.1029/95wr03528.

Le Borgne, T., Dentz, M., Bolster, D., Carrera, J., de Dreuzy, J.-R., Davy, P., 2010. Non-Fickian mixing: Temporal evolution of the scalar dissipation rate in heterogeneous porous media. Adv. Water Resour. 33 (12), 1468-1475. http://dx.doi.org/10.1016/j.adwatres.2010.08.006.

Lichtner, P.C., Kelkar, S., Robinson, B., 2002. New form of dispersion tensor for axisymmetric porous media with implementation in particle tracking. Water Resour. Res. 38 (8), 21-21-21-16. http://dx.doi.org/10.1029/2000wr000100.

Lucy, L., 1977. A numerical approach to the testing of the fission hypothesis. Astron. J. 82, 1013. http://dx.doi.org/10.1086/112164.

Monaghan, J., 2005. Smoothed particle hydrodynamics. Rep. Progr. Phys. 68 (8), 1703-1759. http://dx.doi.org/10.1088/0034-4885/68/8/r01.

Monaghan, J., 2012. Smoothed particle hydrodynamics and its diverse applications. Annu. Rev. Fluid Mech. 44 (1), 323-346. http://dx.doi.org/10.1146/annurev-fluid-120710-101220.

Monaghan, J.J., Gingold, R.A., 1983. Shock simulation by the particle method SPH. J. Comput. Phys. 52 (2), 374-389. http://dx.doi.org/10.1016/0021-9991(83)90036-0.

Pickens, J.F., Grisak, G.E., 1981. Scale-dependent dispersion in a stratified granular aquifer. Water Resour. Res. 17 (4), 1191-1211. http://dx.doi.org/10.1029/wr017i004p01191.

Price, D., 2008. Modelling discontinuities and Kelvin-Helmholtz instabilities in SPH. J. Comput. Phys. 227 (24), 10040-10057. http://dx.doi.org/10.1016/j.jcp.2008.08.011.

Price, D.J., Monaghan, J.J., 2007. An energy-conserving formalism for adaptive gravitational force softening in smoothed particle hydrodynamics and N-body codes. Mon. Not. R. Astron. Soc. 374 (4), 1347-1358. http://dx.doi.org/10.1111/j.1365-2966.2006.11241.x.

Rolle, M., Eberhardt, C., Chiogna, G., Cirpka, O.A., Grathwohl, P., 2009. Enhancement of dilution and transverse reactive mixing in porous media: Experiments and model-based interpretation. J. Contam. Hydrol. 110 (3-4), 130-142. http://dx.doi.org/10.1016/j.jconhyd.2009.10.003.

Saaltink, M.W., Ayora, C., Carrera, J., 1998. A mathematical formulation for reactive transport that eliminates mineral concentrations. Water Resour. Res. 34 (7), 1649-1656. http://dx.doi.org/10.1029/98wr00552.

Salamon, P., Fernández-García, D., Gómez-Hernández, J.J., 2006. A review and numerical assessment of the random walk particle tracking method. J. Contam. Hydrol. 87 (3-4), 277-305. http://dx.doi.org/10.1016/j.jconhyd.2006.05.005.

Schoenberg, I., 1946. Contributions to the problem of approximation of equidistant data by analytic functions. Part a: On the problem of smoothing or graduation. A first class of analytic approximation formulae. Quart. Appl. Math. 4 (2), 45. http://dx.doi.org/10.1007/978-1-4899-0433-1_1.

Sigalotti, L.D.G., Klapp, J., Rendón, O., Vargas, C.A., Peña-Polo, F., 2016. On the kernel and particle consistency in smoothed particle hydrodynamics. Appl. Numer. Math. 108, 242-255. http://dx.doi.org/10.1016/j.apnum.2016.05.007.

Sigalotti, L.D.G., Rendón, O., Klapp, J., Vargas, C.A., Cruz, F., 2019. A new insight into the consistency of the SPH interpolation formula. Appl. Math. Comput. 356, 50-73. http://dx.doi.org/10.1016/j.amc.2019.03.018.

Silliman, S.E., Simpson, E.S., 1987. Laboratory evidence of the scale effect in dispersion of solutes in porous media. Water Resour. Res. 23 (8), 1667-1673. http://dx.doi.org/10.1029/wr023i008p01667.

Springel, V., 2010. Smoothed particle hydrodynamics in astrophysics. Annu. Rev. Astron. Astrophys. 48 (1), 391-430. http://dx.doi.org/10.1146/annurev-astro-081309-130914.

Tartakovsky, A.M., Meakin, P., Scheibe, T.D., Wood, B.D., 2007. A smoothed particle hydrodynamics model for reactive transport and mineral precipitation in porous and fractured porous media. Water Resour. Res. 43 (5), http://dx.doi.org/10.1029/2005wr004770.

Tartakovsky, A.M., Redden, G., Lichtner, P.C., Scheibe, T.D., Meakin, P., 2008. Mixed-induced precipitation: Experimental study and multiscale numerical analysis. Water Resour. Res. 44 (6), http://dx.doi.org/10.1029/2006wr005725.

Tartakovsky, A.M., Scheibe, T.D., Meakin, P., 2009. Pore-scale model for reactive transport and biomass growth. J. Porous Media 12 (5), 417-434. http://dx.doi.org/10.1615/jpormedia.v12.i5.30.

Tartakovsky, A.M., Trask, N., Pan, K., Jones, B., Pan, W., Williams, J.R., 2015. Smoothed particle hydrodynamics and its applications for multiphase flow and reactive transport in porous media. Comput. Geosci. 20 (4), 807-834. http://dx.doi.org/10.1007/s10596-015-9468-9.

Tran-Duc, T., Bertevas, E., Phan-Thien, N., Khoo, B.C., 2016. Simulation of anisotropic diffusion processes in fluids with smoothed particle hydrodynamics. Internat. J. Numer. Methods Fluids 82 (11), 730-747. http://dx.doi.org/10.1002/fld.4238.

Vacondió, R., Altomare, C., Leffe, M.D., Hu, X., Touzé, D.L., Lind, S., Marongiu, J.-C., Marrone, S., Rogers, B.D., Souto-Iglesias, A., 2021. Grand challenges for smoothed particle hydrodynamics numerical schemes. Comput. Part. Mech. 8 (3), 575-588. http://dx.doi.org/10.1007/s40571-020-00354-1.

Vacondió, R., Rogers, B.D., Stansby, P.K., Mignosa, P., 2012. SPH modeling of shallow flow with open boundaries for practical flood simulation. J. Hydraul. Eng. 138 (6), 530-541. http://dx.doi.org/10.1061/(asce)hy.1943-7900.0000543.

Violeau, D., 2012. Fluid Mechanics and the SPH Method: Theory and Applications. Oxford University Press, Oxford.

Wang, J., Hu, W., Zhang, X., Pan, W., 2019. Modeling heat transfer subject to inhomogeneous Neumann boundary conditions by smoothed particle hydrodynamics and peridynamics. Int. J. Heat Mass Transfer 139, 948-962. http://dx.doi.org/10.1016/j.ijheatmasstransfer.2019.05.054.

Ye, T., Pan, D., Huang, C., Liu, M., 2019. Smoothed particle hydrodynamics (SPH) for complex fluid flows: Recent developments in methodology and applications. Phys. Fluids 31 (1), http://dx.doi.org/10.1063/1.5068697.

Zhu, Y., Fox, P.J., 2001. Smoothed particle hydrodynamics model for diffusion through porous media. Transp. Porous Media 43 (3), 441-471. http://dx.doi.org/10.1023/a:1010769915901.

Zhu, Y., Fox, P.J., 2002. Simulation of pore-scale dispersion in periodic porous media using smoothed particle hydrodynamics. J. Comput. Phys. 182 (2), 622-645. http://dx.doi.org/10.1006/japh.2002.7189.

Zhu, Q., Hernquist, L., Li, Y., 2015. Numerical convergence in smoothed particle hydrodynamics. Astrophys. J. 800 (1), 6. http://dx.doi.org/10.1088/0004-637x/800/1/6.