Smoothed Particle Hydrodynamics for anisotropic dispersion in heterogeneous porous media

Advances in Water Resources 183 (2024) 104601 Available online 7 December 2023 0309-1708/© 2023 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). Contents lists available at ScienceDirect Advances in Water Resources journal homepage: www.elsevier.com/locate/advwatres Smoothed Particle Hydrodynamics for anisotropic dispersion in heterogeneous porous media Rodrigo Pérez-Illanes a,b,∗, Guillem Sole-Mari c, Daniel Fernàndez-Garcia a,b aDepartment of Civil and Environmental Engineering (DECA), Universitat Politècnica de Catalunya, Barcelona, Spain bHydrogeology Group (UPC-CSIC), Barcelona, Spain cGéosciences Rennes, Université de Rennes, Rennes, France 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 ∗Corresponding author at: Department of Civil and Environmental Engineering (DECA), Universitat Politècnica de Catalunya, Barcelona, Spain. E-mail address: [email protected] (R. Pérez-Illanes). 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 (Bandara 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 https://doi.org/10.1016/j.advwatres.2023.104601 Received 28 August 2023; Received in revised form 26 November 2023; Accepted 6 December 2023 Advances in Water Resources 183 (2024) 104601 2 R. Pérez-Illanes et al. 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 𝑑𝐫 𝑑𝑡 =𝐪 𝜙,(1) 𝑑(𝜙𝐶) 𝑑𝑡 = ∇ ⋅(𝜙𝐃∇𝐶),(2) where 𝑑∕𝑑𝑡(⋅)is the total time derivative, 𝐫is the spatial coordinate vector for a particle’s position, 𝐪=𝜙𝐯is the Darcy velocity, 𝜙is the medium porosity, 𝐯is the macroscopic flow velocity and 𝐶is the solute concentration expressed as mass per unit water volume. The hydrodynamic dispersion tensor is given by (Bear and Cheng,2010) 𝐷𝑖𝑗 = (𝛼𝑇|𝐯|+𝑚)𝛿𝑖𝑗 +(𝛼𝐿−𝛼𝑇)𝑣𝑖𝑣𝑗 |𝐯|,(3) Advances in Water Resources 183 (2024) 104601 3 R. Pérez-Illanes et al. where 𝛼𝐿and 𝛼𝑇are the longitudinal and transverse dispersivities, respectively, 𝑚is the effective molecular diffusion (corrected for tortuosity effects) and 𝛿𝑖𝑗 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 𝜆𝐃=𝐷𝑥𝑥∕𝐷𝑦𝑦 ≈𝛼𝐿∕𝛼𝑇. In typical aquifer applications, 𝜆𝐃 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 𝐴(𝐫) as (Monaghan,2005;Violeau,2012) 𝐴(𝐫) = ∫𝛺 𝐴(𝐫′)𝛿(𝐫−𝐫′)𝑑𝐫′,(4) where 𝛺is the domain volume and 𝛿is the Dirac delta function. The Dirac delta is then approximated by a positive continuous symmetric kernel function 𝑊(|𝐫−𝐫′|, ℎ), with ℎ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 𝐻, such that kernel density equals zero at any distance higher than 𝐻. 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 (ℎ2)(Violeau,2012;Sigalotti et al., 2019). Continuous integrals are discretized on the set of particles, leading to (Monaghan,2005) 𝐴(𝐫𝑎) =  ∑ 𝑏 𝑚𝑏 𝜌𝑏 𝐴𝑏𝑊(|𝐫𝑎−𝐫𝑏|, ℎ𝑎),(5) where subscripts 𝑎and 𝑏denote the center and neighbor particles respectively. Sums are performed over a finite set of neighbors within the support volume. For the remainder of the article, SPH sums are assumed to always occur over the set of neighbors and is not explicitly indicated. 𝑚𝑏and 𝜌𝑏represent the mass and density of particle 𝑏, 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 (−𝜓), where the exponent 𝜓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 (𝜓≈ 1) than in the case of disordered particles (𝜓≈ 0.5). Density is computed as (Monaghan, 2005) 𝜌𝑎=𝜌(𝐫𝑎) = ∑ 𝑏 𝑚𝑏𝑊𝑎𝑏(ℎ𝑎),(6) where 𝑊𝑎𝑏(ℎ𝑎)denotes the kernel function evaluated at distance 𝑟𝑎𝑏 = |𝐫𝑎−𝐫𝑏|with a smoothing length ℎ𝑎. In this paper, the particle mass 𝑚𝑏 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 Fig. 1. Fifth order B-Spline kernel and characteristic distances. (𝑎)Normalized shape function 𝑤(𝑞), and (𝑏)normalized shape derivative 𝑤′(𝑞).𝜎𝑊is the standard deviation of the kernel’s shape. 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 𝑊( 𝐫, ℎ) = 𝑐𝜈 (𝜅𝜈ℎ)𝜈𝑤(𝑞),(7) where  𝐫=𝐫−𝐫′is the vector difference between the center and neighbor particle positions and 𝑤(𝑞)is a shape function normalized to a unit support radius, with 𝑞=| 𝐫|∕(𝜅𝜈ℎ)the normalized kernel coordinate. 𝑐𝜈is a positive constant which makes the kernel integral unitary and whose value depends on the number of spatial dimensions 𝜈; and 𝜅𝜈=𝐻∕ℎis an aspect ratio specific to the kernel function. From Eq. (7) follows that kernel gradients are expressed as (Violeau, 2012) ∇𝑊( 𝐫, ℎ) = 𝐹(𝑞) 𝐫,(8) where 𝐹(𝑞)is defined as the kernel scalar derivative 𝐹(𝑞) = 𝑐𝜈 (𝜅𝜈ℎ)𝜈+1 𝑤′(𝑞),(9) which is negative for typical kernel functions (Monaghan,2005), and  𝐫= 𝐫∕| 𝐫|is the unit vector between particles. In this article, a fifth order B-Spline kernel is used (Fig. 1) 𝑤(𝑞) = (1 − 𝑞)5 +− 6 (2 3−𝑞)5 ++ 15 (1 3−𝑞)5 + ,(10) where the operator 𝑥+stands for max(0, 𝑥). In two dimensions, this kernel’s constants have values 𝑐2= 37⋅7∕478𝜋and 𝜅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 Advances in Water Resources 183 (2024) 104601 4 R. Pérez-Illanes et al. considering density estimates. The authors argue that, as the total number of particles tends to infinity, 𝑁→∞, it should be satisfied simultaneously that ℎ→0and →∞. 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 𝛿𝑟= (𝑚∕𝜌)1∕𝜈, that is, ℎ=𝛾𝛿𝑟. The parameter 𝛾controls the number of neighbors within the kernel, with a typically adopted value of 𝛾= 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), 𝛾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 𝛾is convenient due to its simplicity, but lacks theoretical rigor as 𝛾is not a function of the number of particles. In Appendix A it is shown that, by introducing consistency principles, the following expression for 𝛾as a function of the total number of SPH particles can be obtained 𝛾(𝑁) = 1 𝜅𝜈( 𝑉𝜈)1∕𝜈 𝑁1∕𝜈−1∕𝛽,(11) where 𝑉𝜈is the volume of the unit sphere, is a proportionality constant, and 𝛽is a positive exponent relating the scaling of the smoothing length with the number of particles ℎ∝𝑁−1∕𝛽. For two dimensional problems, 𝛽≥4has 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) 𝑑𝐶𝑎 𝑑𝑡 =∑ 𝑏 𝑚𝑏 𝜌𝑎𝑏 𝐶𝑏𝑎  𝐷𝑎𝑏 𝑖𝑗 𝐺𝑖𝑗 𝑎𝑏 (12) where 𝐶𝑏𝑎 =𝐶𝑏−𝐶𝑎, and 𝜌𝑎𝑏 is an averaged density.  𝐷𝑎𝑏 𝑖𝑗 is the 𝑖, 𝑗 component of the dispersion tensor averaged (typically harmonic) between particles 𝑎and 𝑏.𝐺𝑖𝑗 𝑎𝑏 is a term describing the kernel gradient defined as 𝐺𝑖𝑗 𝑎𝑏 = − [(𝜈+ 2)𝑟𝑖 𝑎𝑏 𝑟𝑗 𝑎𝑏 −𝛿𝑖𝑗 ]𝐹𝑎𝑏 𝑟𝑎𝑏 ,(13) where 𝑟𝑖 𝑎𝑏 =𝑟𝑖 𝑎𝑏∕𝑟𝑎𝑏 is the 𝑖th-component of the vector  𝐫𝑎𝑏 and 𝐹𝑎𝑏 = 𝐹(𝑞𝑎𝑏), with 𝑞𝑎𝑏 the normalized kernel coordinate distance 𝑞between particles 𝑎and 𝑏. In Eq. (12), indexes 𝑖, 𝑗 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  𝐫𝑎𝑏. 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 𝑎, 𝑏 for simplicity, for 𝐶𝑎< 𝐶𝑏and under proper physical conditions, it is expected that 𝑑𝐶𝑎∕𝑑𝑡 > 0, which is satisfied as long as the product between dispersion and the term 𝐺𝑖𝑗 𝑎𝑏 is positive. In two dimensions, and assuming a diagonal dispersion tensor, such condition translates into the inequality 𝐷𝑥𝑥 𝐷𝑦𝑦 [1−(𝜈+ 2)(𝑟𝑥 𝑎𝑏)2]≤(𝜈+ 2)(1 − (𝑟𝑥 𝑎𝑏)2)−1,(14) where it has been applied that (𝑟𝑥 𝑎𝑏)2+ (𝑟𝑦 𝑎𝑏)2= 1. Limit cases can be obtained by considering the limit values of (𝑟𝑥 𝑎𝑏)2∈ {0,1}, which yields that the scheme preserves physical dispersive flux when dispersion anisotropy satisfies 𝜆𝐃∈ [1∕(𝜈+ 1), 𝜈 + 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 𝜆𝐃= 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 𝐃=𝐁𝐁𝑇(LaBolle et al., 1996), where the superscript 𝑇stands for transpose. The matrix 𝐁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-Garcia et al.,2005;Salamon et al.,2006) 𝐁=[𝑣𝑥 |𝐯|√𝛼𝐿|𝐯|+𝑚−𝑣𝑦 |𝐯|√𝛼𝑇|𝐯|+𝑚 𝑣𝑦 |𝐯|√𝛼𝐿|𝐯|+𝑚 𝑣𝑥 |𝐯|√𝛼𝑇|𝐯|+𝑚].(15) Furthermore, the inverse of this matrix can be computed analytically as 𝐁−1 =⎡⎢⎢⎣ 1 √𝛼𝐿|𝐯|+𝑚 𝑣𝑥 |𝐯| 1 √𝛼𝐿|𝐯|+𝑚 𝑣𝑦 |𝐯| −1 √𝛼𝑇|𝐯|+𝑚 𝑣𝑦 |𝐯| 1 √𝛼𝑇|𝐯|+𝑚 𝑣𝑥 |𝐯|⎤⎥⎥⎦ .(16) The coordinate transformation in ASPHAD is such that 𝐗=𝐁−1𝐱, where 𝐗is the position in the new frame of reference. This leads to the following SPH scheme for dispersion (Tran-Duc et al.,2016), 𝑑𝐶𝑎 𝑑𝑡 = 2 ∑ 𝑏 𝑚𝑏 𝜌𝑏 𝐶𝑎𝑏 |||𝐁−1 𝑎𝑏  𝐫𝑎𝑏||| 2 𝐹𝑎𝑏 𝑟𝑎𝑏 ,(17) where 𝐁−1 𝑎𝑏 is the inverse displacement matrix (Eq. (16)) averaged between particles 𝑎and 𝑏, and 𝐶𝑎𝑏 =𝐶𝑎−𝐶𝑏.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 𝑗th-component of the gradient of solute concentrations is expressed as 𝐺𝑗 𝑎= ∇𝑗 𝑎𝐶𝑎=1 𝜌𝑎𝑓𝑎∑ 𝑏 𝑚𝑏𝐶𝑏𝑎∇𝑗 𝑎𝑊𝑎𝑏(ℎ𝑎),(18) where 𝑓𝑎= 1 + ℎ𝑎 𝜈𝜌𝑎∑ 𝑏 𝑚𝑏 𝜕𝑊𝑎𝑏(ℎ𝑎) 𝜕ℎ𝑎 (19) Advances in Water Resources 183 (2024) 104601 5 R. Pérez-Illanes et al. is a gradient correction which stems from considering a spatially variable smoothing length following the relation ℎ𝜈 𝑎∝ 1∕𝜌𝑎(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 𝑑𝐶𝑎 𝑑𝑡 =𝜌𝑎∑ 𝑏 𝑚𝑏[𝐷𝑎 𝑖𝑗 𝐺𝑖 𝑎∇𝑗 𝑎𝑊𝑎𝑏(ℎ𝑎) 𝑓𝑎𝜌2 𝑎 + 𝐷𝑏 𝑖𝑗 𝐺𝑖 𝑏∇𝑗 𝑎𝑊𝑎𝑏(ℎ𝑏) 𝑓𝑏𝜌2 𝑏].(20) As before, indexes 𝑖, 𝑗 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 𝑓𝑎= 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 quasirandom 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 𝜌𝑎𝑏, 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, 𝛿𝐃 𝑡, is computed as the minimum over all particles of the expression 𝛿𝐃 𝑡≤𝐶𝑇 ℎ2 ∑𝑖𝐷𝑖𝑖 ,(21) with 𝐶𝑇= 0.1(Herrera and Beckie,2013). For heterogeneous experiments, also the advection-limited time step 𝛿𝐯 𝑡is computed by imposing 𝐶𝐹 𝐿 = 0.1and using the maximum spatial velocities, that is 𝛿𝐯 𝑡≤𝐶𝐹 𝐿 max 𝑣𝑥∕𝛥𝑥+ max 𝑣𝑦∕𝛥𝑦 ,(22) where 𝛥𝑥, 𝛥𝑦are the flow grid sizes. Finally, time step is selected as the most restrictive, i.e., 𝛿𝑡= min{𝛿𝐃 𝑡, 𝛿𝐯 𝑡}.(23) Solute concentrations are then integrated explicitly as 𝐶𝑡+𝛿𝑡 𝑎=𝐶𝑡 𝑎+𝛿𝑡(𝑑𝐶𝑎 𝑑𝑡 )𝑡 ,(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 𝐴, following expressions of the form (5), is bounded by |||𝐴f l 𝑎−𝐴𝑎||| |𝐴𝑎|≤1 2𝜖, (25) where 𝐴f l 𝑎is the floating point representation of the SPH estimate and 𝜖 > 0is the machine epsilon. The latter is a property of the floating point system and is the maximum number that satisfies 1 = 1+𝜖. Models were implemented here using 64 bit precision with 𝜖= 2.220446 × 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 𝑐is performed from the set of SPH particles as (Herrera et al., 2009) 𝐶𝑐=1 𝜌𝑐  ∑ 𝑏 𝑚𝑏𝐶𝑏𝑊𝑐𝑏,(26) where 𝜌𝑐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 𝐶∗is given by 𝐶∗(𝐫, 𝑡) 𝐶0 =𝑑2 3 exp [−𝑥21−𝑦22+ 4𝑥 𝑦3 8𝑡21+ 4𝑡𝑑22+ 2𝑑4],(27) where 1= 2𝑡𝐷𝑦𝑦 +𝑑2 2= 2𝑡𝐷𝑥𝑥 +𝑑2 3=𝑡𝐷𝑥𝑦 1=𝐷𝑥𝑥𝐷𝑦𝑦 −𝐷2 𝑥𝑦 2=𝐷𝑥𝑥 +𝐷𝑦𝑦 3=(4𝑡21+ 2𝑡𝑑22+𝑑4)1∕2 . (28) In the above expression, 𝐶0is the maximum initial concentration, which is employed to express concentrations in dimensionless form,  𝐶=𝐶∕𝐶0.(𝑥, 𝑦)=(𝑥−𝑥0, 𝑦−𝑦0)are the relative coordinates with respect to the plume center (𝑥0, 𝑦0) = (0,0). The parameters follow values taken Advances in Water Resources 183 (2024) 104601 6 R. Pérez-Illanes et al. Fig. 2. Numerical error 𝑒𝑅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 𝑇= 300[d].(𝑎)ER, (𝑏)ASPHAD and (𝑐)TFD. Simulations consider 𝑁= 300 × 300 particles with uniform distribution, and flow aligned with the 𝑥-axis. from Herrera and Beckie (2013). The initial width of the Gaussian plume is 𝑑= 44[m] and the domain size is 𝐿×𝐿= 2000 × 2000[m × m]. Flow velocity is constant, |𝐯|= 1[m∕d], following a vector that forms the angle 𝜃with respect to the 𝑥-direction, although for most simulations is simply considered to be aligned with the 𝑥-axis. Dispersion neglects molecular diffusion, with longitudinal dispersivity 𝛼𝐿= 10[m] and anisotropies 𝜆𝐃∈ {1,10,100}. The number of SPH particles is modified within the range 𝑁∈ [200 × 200,600 × 600], and the 𝛾parameter for the smoothing length varies within the range 𝛾∈ [1,7]. Transport simulations end at a reference time 𝑇= 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 𝜉𝑟∼𝑈(−𝛿𝑟∕2, 𝛿𝑟∕2) to the coordinates of the uniform distribution, with 𝛿𝑟being the original spacing between particles. Differences between the numerical and analytical results are quantified by means of the Root Mean Squared Error (RMSE), 𝑒𝑅=√ √ √ √1 𝑁 𝑁 ∑ 𝑎(𝐶∗ 𝑎−𝐶𝑎)2.(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 𝑥-axis. In all cases, there is a smoothing that minimizes the error with some influence of anisotropy on the Fig. 3. Concentrations for the problem of Gaussian dispersion with uniform distribution of particles (𝑁= 300 × 300), at the reference time 𝑇= 300[d], and flow aligned with the 𝑥-axis. First column group results for all SPH schemes considering 𝜆𝐃= 10 (panels 𝑎, 𝑐, 𝑒, 𝑔) and second column for 𝜆𝐃= 100 (panels 𝑏, 𝑑, 𝑓 , ℎ). Contour lines are drawn for  𝐶∈ {0.4,0.1,10−3,10−5}(solid: analytical, dotted: numerical). White bands (if present) represent negative values. 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. 3𝑐, 𝑑). Both ASPHAD and TFD are free of negative concentrations for anisotropy 𝜆𝐃= 10 (Fig. 3𝑒, 𝑔 respectively), with some differences in the small-scale spatial distribution of concentrations. Particularly, TFD appears to generate a more elliptical distribution for contours  𝐶∼ 10−5, but in general both schemes are in good agreement with the analytical reference. Increasing the anisotropy to 𝜆𝐃= 100 leads to some peculiar results. In particular, the concentration contours in the ER method (Fig. 3𝑑) are now influenced by the bands of negative values generated by the instability, with a minimum concentration of  𝐶= −4.6 × 10−5. The TFD results (Fig. 3ℎ) also present some negative values while moving away from the plume center, but these are of a different nature. The minimum concentration is  𝐶= −9.8 × 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 Advances in Water Resources 183 (2024) 104601 7 R. Pérez-Illanes et al. Table 1 Errors in homogeneous simulations with dispersion anisotropy, for different angles of alignment 𝜃between the flow direction and the 𝑥-axis. Table shows the value of 𝛾that lead to the minimum error 𝑒𝑅, considering 𝑁= 300 × 300 and uniform distribution of particles. 𝜆𝐃𝜃ER ASPHAD TFD 𝛾 𝑒𝑅𝛾 𝑒𝑅𝛾 𝑒𝑅 10 0 4.5 0.00013 1.4 0.00205 3.6 0.00062 𝜋∕6 4.5 0.00015 1.3 0.00420 3.6 0.00062 𝜋∕4 4.5 0.00015 1.4 0.00299 3.6 0.00062 100 0 4.2 0.00032 1.3 0.00063 3.6 0.00056 𝜋∕6 4.2 0.00035 2.4 0.01202 3.6 0.00056 𝜋∕4 4.2 0.00036 1.1 0.00434 3.6 0.00056 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 Fig. 4. Concentrations for the problem of Gaussian dispersion with uniform distribution of particles (𝑁= 300 × 300), at the reference time 𝑇= 300[d], and diagonally aligned flow (𝜃=𝜋∕4). First column group results for all SPH schemes considering 𝜆𝐃= 10 (panels 𝑎, 𝑐, 𝑒) and second column for 𝜆𝐃= 100 (panels 𝑏, 𝑑, 𝑓). Contour lines are drawn for  𝐶∈ {0.4,0.1,10−3,10−5}(solid: analytical, dotted: numerical). White bands (if present) represent negative values. while considering the flow aligned with the 𝑥-axis, which is consistent with the results previously shown in Figure (3𝑐, 𝑑). To further expand on this topic, anisotropic simulations with all the considered SPH schemes were also performed with the velocity vector forming angles of 𝜃∈ {𝜋∕6, 𝜋∕4} with respect to the 𝑥-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 𝜃=𝜋∕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 𝜃= 0 (particles aligned along the 𝑥-axis) or 𝜃=𝜋∕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. 4𝑑). In contrast, anisotropic simulations performed with the ER (Fig. 4𝑎, 𝑏) and TFD (Fig. 4𝑒, 𝑓) 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 Fig. 5. Numerical error 𝑒𝑅for different particle distributions, with respect to the kernel size, for the problem of Gaussian dispersion in a homogeneous medium. Simulations consider 𝑁= 300 × 300, flow aligned with the 𝑥-axis, and errors computed at the reference time 𝑇= 300[d]. First column group results for isotropic dispersion (panels 𝑎, 𝑐, 𝑒) and second column for anisotropy 𝜆𝐃= 100 (panels 𝑏, 𝑑, 𝑓). SPH method is indicated in each panel. Advances in Water Resources 183 (2024) 104601 8 R. Pérez-Illanes et al. Fig. 6. Concentrations for anisotropic dispersion at 𝑇= 300[d] with quasi-random distribution, 𝑁= 300 × 300 particles, and flow aligned with the 𝑥-axis. First column group numerical results for each solver with 𝜆𝐃= 10 and second column for 𝜆𝐃= 100. White bands (if present) indicate negative values. quasi-random particle distribution (similar to the random case) for the highest dispersion anisotropy (Fig. 5𝑑). The same does not occur with the ER (Fig. 5𝑎, 𝑏) and TFD (Fig. 5𝑒, 𝑓) 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. 6𝑎, 𝑏) and TFD (Fig. 6𝑒, 𝑓) 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  𝐶= −4 × 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. 6𝑐, 𝑑). The effect was discussed in Tran-Duc et al. (2016), but only considering ratios of anisotropy of up to 𝜆𝐃= 10. In this case, numerical peak concentration is ≈ 14% higher than the reference solution (Fig. 7𝑎, 𝑏). Even more drastical discrepancies arise for anisotropy 𝜆𝐃= 100, where the peak concentration is now ≈ 53% higher Fig. 7. Numerical and analytical concentration profiles through the center of a Gaussian plume at time 𝑇= 300[d], considering 𝑁= 300 × 300, quasi-random distribution, and flow aligned with the 𝑥-axis. Panels (𝑎, 𝑏)present the longitudinal (𝑥) and transverse (𝑦) profiles, respectively, for dispersion anisotropy 𝜆𝐃= 10 and panels (𝑐, 𝑑)the analogous for 𝜆𝐃= 100. than the analytical model (Fig. 7𝑐, 𝑑), 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 𝑥-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 𝛾for each SPH method, with respect to the number of particles is shown in Fig. 8, and compared with the reference scaling for twodimensional domains (𝛽= 4 in Eq. (11)). Overall, the ER and TFD schemes display an increasing value of 𝛾for increasing number of particles, almost independently from the dispersion anisotropy, which leads to a consistent decrease in error for both uniform and quasirandom distributions of particles (Fig. 8𝑏and 𝑓, respectively). Errors Table 2 Parameters of expression (11) obtained from fitting 𝛾and 𝑁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). 𝑟2is the coefficient of determination. Dist. 𝜆𝐃ER ASPHAD TFD 𝛽 𝑟2𝛽 𝑟2𝛽 𝑟2 Uni 1 2.72 3.52 0.91 2.70 3.61 0.99 0.46 4.39 0.99 10 2.03 3.64 0.98 32.9 2 – 0.79 3.85 0.99 100 0.56 4.26 0.97 24.7 2 – 0.93 3.77 0.99 Q.R. 1 284 2.08 0.01 0.04 8.11 0.82 0.31 4.87 0.95 10 43.6 2.48 0.41 2.36 2.68 0.33 1.08 3.80 0.93 100 266 2.05 0.01 3.62 3.04 0.80 0.33 4.79 0.98 Advances in Water Resources 183 (2024) 104601 9 R. Pérez-Illanes et al. 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 (𝑎),(𝑐),(𝑒) present the 𝛾factor, and panels (𝑏),(𝑑),(𝑓)the corresponding error. Reference curves consider 𝛽= 4. obtained from ASPHAD (Fig. 8𝑑) 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 𝛾in cases with uniform particle distribution, regardless of the particle resolution (Fig. 8𝑐). 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. 8𝑑). The 𝛾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 𝑟2= 0.93, confirming that the scheme follows the scaling from Eq. (11), with error decaying consistently with the proportionality 𝑒𝑅∝𝑁−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 𝑁𝑐= 1500 × 300 unitary cells (𝛥𝑥, 𝛥𝑦) = (1,1)[m]. The natural logarithm of the hydraulic conductivity 𝑌(𝐱)represents one individual realization of a sequential Gaussian simulation with zero mean and spherical variogram with correlation length 𝐼𝑌= 20𝛥𝑥. Heterogeneities are controlled by logconductivity variances 𝜎2 𝑌∈ {0.25,1.00,2.25} through the relation 𝐾(𝐱) = exp (𝜎𝑌𝑌(𝐱)). For each hydraulic conductivity distribution, groundwater flow is obtained from solving on the numerical grid the steady-state Darcy’s equation ∇⋅𝐪= 0, where 𝐪= −𝐾(𝐱)∇, and  is the hydraulic head. Flow is induced by a unit-mean hydraulic head 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. gradient along the 𝑥-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 𝛼𝐿= 0.05𝛥𝑥which leads to advection-dominated transport, with a field Péclet number of Pe = 𝐼𝑌∕𝛼𝐿= 400. A rectangular solute plume is initially released near the aquifer inlet with uniform concentration 𝐶0. The plume size is 𝑑𝑥×𝑑𝑦= 30𝛥𝑥× 200𝛥𝑦, with the lower left corner placed at (𝑥0, 𝑦0) = (15𝛥𝑥,50𝛥𝑦). 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 𝑥-coordinate to zero and defining the 𝑦-coordinate from the intersection of a random number 𝜉𝑟∼𝑈(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 𝑚𝑅=𝐶0𝑉𝑀∕𝑁𝑅, where 𝑉𝑀is the volume occupied by the solute injection and 𝑁𝑅= 108the number of particles. These parameters lead to a concentration resolution of 𝐶min = 6 × 10−5𝐶0, estimated as the mass of one particle inside a flowmodel 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, 𝑁∈ {0.2M,0.8M,1.8M}, which translates into averaged densities of 𝜌 ∈ {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 𝜌 = 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 𝜎2 𝑌= 1, and the highest dispersion