MODPATH‐RW: A Random Walk Particle Tracking Code for Solute Transport in Heterogeneous Aquifers

Random walk particle tracking (RWPT) is a discrete particle method that offers several advantages for simulating solute transport in heterogeneous geological systems. The formulation is a discrete solution to the advection‐dispersion equation, yielding results that are not influenced by grid‐related numerical dispersion. Numerical dispersion impacts the magnitude of concentrations and gradients obtained from classical grid‐based solvers in advection‐dominated problems with relatively large grid Péclet numbers. Accurate predictions of concentrations are crucial for reactive transport studies, and RWPT has been recognized for its potential benefits for this kind of application. This highlights the need for a solute transport program based on RWPT that can be seamlessly integrated with industry‐standard groundwater flow models. This article presents a solute transport code that implements the RWPT method by extension of the particle tracking model MODPATH, which provides the base infrastructure for interacting with several variants of MODFLOW groundwater flow models. The implementation is achieved by developing a method for determining the exact cell‐exit position of a particle undergoing simultaneous advection and dispersion, allowing for the sequential transfer of particles between flow model cells. The program is compatible with rectangular unstructured grids and integrates a module for the smoothed reconstruction of concentrations. In addition, the program incorporates parallel processing of particles using the OpenMP library, enabling faster simulations of solute transport in heterogeneous systems. Numerical test cases involving different applications in hydrogeology benchmark the RWPT model with well‐known transport codes.


Introduction
Solute transport models based on particles offer advantages over classical grid-based formulations in advection-dominated transport problems with relatively large grid Péclet numbers.These conditions are known to be challenging for grid-based solvers in models with a detailed representation of heterogeneity, with results influenced by artificial effects such as numerical dispersion and oscillations (e.g., Cirpka et al. 1999;Herrera et al. 2009;Boso et al. 2013;Sanchez-Vila and Fernàndez-Garcia 2016;Benson et al. 2017).In contrast, simulations with particles are free of these effects, leading to an improved representation of spatial solute gradients, concentrations and dilution metrics.These quantities are notably relevant for reactive transport in heterogeneous systems, where solute plumes adopt complex spatial distributions controlled by the magnitude and distribution of aquifer heterogeneities.Random walk particle tracking (RWPT) is a formulation where particles represent individual components of the solute plume undergoing simultaneous advection and dispersion.Its applicability and challenges for hydrogeological systems has been discussed in the literature since a relatively long time (Ahlstrom et al. 1977;Prickett et al. 1981;Kinzelbach 1988;Tompson et al. 1988;Uffink 1988;Tompson and Dougherty 1992;Kitanidis 1994;LaBolle et al. 1996).The formulation is based on an equivalent Fokker-Planck form of the Advection-Dispersion Equation (ADE), with concentration related to the probability density function of finding a particle at a given location and time.
Despite recent and novel developments, solute transport codes based on RWPT have not permeated to practical field applications at the same rate than Eulerian grid-based alternatives.This is in part explained by some intrinsic challenges of the method, like the application of consistent boundary conditions and the need to use a large number of particles to attenuate the statistical fluctuations in concentrations.Furthermore, reactive transport simulations based on RWPT can easily involve the displacement of multiple particle groups (for different chemical species), potentially leading to high computational demand in scenarios where these are interacting via chemical reactions.Still, some codes implementing RWPT have been addressed in the literature: RWHet (LaBolle 2006) with reported applications based on MODFLOW, EcoSLIM (Maxwell et al. 2018;Yang et al. 2021) a parallel code compatible with ParFlow, RW3D (Fernàndez-Garcia et al. 2005;Henri and Fernàndez-Garcia 2015) a Fortran implementation compatible with MODFLOW-2000, PARTRACE (Neuendorf 1997;Bechtold et al. 2011) employed in discussions of the reflection method for displacing particles, and more recently PAR 2 (Rizzo et al. 2019) a GPU implementation able to interpret the flow-transport link (FTL) file that connects MODFLOW with MT3DMS (Zheng and Wang 1999).
A groundwater flow model widely employed in hydrogeology is MODFLOW, developed by the U.S. Geological Survey (McDonald and Harbaugh 2003;Harbaugh 2005;Langevin et al. 2017), with such a significance in research and engineering that developments in solute transport have pursued a straightforward interaction with the different versions of this program (e.g., Konikow et al. 1996;Zheng and Wang 1999;Prommer et al. 2003;Langevin et al. 2022).Along this line, it is then reasonable to consider the implementation of a RWPT code capable to be easily integrated with MODFLOW in order to facilitate its adoption for practical field studies.A particle tracking engine from the same family of programs is MODPATH (Pollock 2016), which provides the base infrastructure for interacting with groundwater flow models based on MODFLOW-2005 (Harbaugh 2005), MODFLOW-USG (Panday et al. 2013;Pollock 2015) and MODFLOW 6 (Langevin et al. 2017).Among others, the program is classically employed in aquifers to delineate well capture zones and streamlines (e.g., Bair et al. 1991;Buxton et al. 1991;Riva et al. 2006;Horn and Harter 2009).Hence, is a strong alternative to be considered as initial platform from where to start the development of a robust RWPT code.However, MODPATH displaces the equivalent to a fluid particle (only by advection) implementing the semi-analytical algorithm of Pollock (1988), which moves a particle from its current position within a flowcell to the exit face with a single step.MODPATH source code is built around this idea and particles are transfered sequentially between flow model cells until one of several possible stopping conditions is met.In RWPT, knowing the exact cell-exit position of a particle is not straightforward because on their way toward an exit face, particles undergo random displacements due to hydrodynamic dispersion.
The main objective of this work is to extend the current version of MODPATH source code (Pollock and Provost 2017), in order to integrate a RWPT solver for the ADE, while preserving the base functionalities provided by the program.An specific objective is to implement the new tracking methodology minimizing interventions to the base code structure, although inevitably some major changes need to be introduced.In this regard, this article presents a method for computing the exact cell-exit position of a particle displaced by simultaneous advection and dispersion, which allows to preserve the sequential transfer of particles between flow model cells.Furthermore, a feature already integrated into MOD-PATH is the compatibility with rectangular unstructured grids and preserving this characteristic for the RWPT implementation is also an specific objective of this work.Besides the RWPT solver, an additional Fortran module for the smoothed reconstruction of concentrations with Grid Projected Kernel Density Estimation (GPKDE) is integrated into the program, which allows to obtain continuous concentration maps in simulations with noninfinite number of particles (as in Sole-Mari et al. 2019).The release of the code is complemented with Python utilities extended from the interface FloPy (Bakker et al. 2016), assisting the writing of model input files.These and other program features lead to a set of tools integrating an important part of the workflow required for applications of RWPT to solute transport studies in aquifers.
The structure of this article is as follows.In methods, specifics about the integration of the RWPT solver into MODPATH are presented, followed by results and discussion where numerical test cases are employed to benchmark the code with grid-based solute transport models, considering aquifer applications with homogeneous and heterogeneous hydraulic properties.In both these sections, some limitations of the current implementation are addressed, motivating further development pathways for the program.The article finalizes with the conclusions of this development and the presentation of MODPATH-RW.

Methods
The integration of RWPT into MODPATH was performed in two main stages.The first stage considered the implementation of parallel processing for particles' displacement (Pérez-Illanes and Fernàndez-Garcia 2023) with the OpenMP library (OpenMP 2020).A RWPT solver inherently requires more calculations than pure advection, hence, parallel processing is a necessary feature for these programs (e.g., Maxwell et al. 2018;Rizzo et al. 2019).The second stage is the integration of the new solver and related solute transport functionalities.This section begins with an overview of RWPT and continues with specific aspects of the implementation.

Random Walk Particle Tracking
It is considered that solute transport through porous media can be described by the Advection-Dispersion Equation (ADE), which is written similar to the MOD-FLOW 6 framework as (Bear and Cheng 2010;Langevin et al. 2022) where R s (x, t) is the retardation factor (dimensionless), S w (x, t) is the water saturation (dimensionless), φ(x, t) is the medium porosity (dimensionless), q(x, t) the specific discharge (L/T), D(x, t) is the hydrodynamic dispersion tensor (L 2 /T), f is the term comprising external solute sources (M/L 3 T), and c is the dissolved solute concentration (M/L 3 ).The dispersion tensor is considered to satisfy where α T and α L are the transverse and longitudinal dispersivities, D m is the apparent molecular diffusion coefficient (corrected by tortuosity), and I d the identity matrix.Besides this form, the program also permits to alternatively choose the axisymmetric dispersion tensor from Lichtner et al. (2002), with the vertical axis taken as the axis of symmetry.RWPT simulates the ADE by partitioning the solute mass into a large number of representative particles, which are transported according to simple physical relationships derived from stochastic integral theory (Tompson et al. 1988;Risken 1989) or by the method of spatial moments (Kitanidis 1994).The latter is perhaps more intuitive and directly provides a physical interpretation to a complex problem.In this case, each particle is seen as a very small plume of the total concentration ρ = R s S w φc, the motion of which is determined by its spatial moments.The mass of a particle represents the total mass of the chemical compound, considering both the sorbed and the aqueous phase (as in Tompson 1993).Working with the total concentration ρ allows the method to properly represent the changes in space and time of R s , S w , and φ.Following Kitanidis (1994), integration by parts of (1) yields the following expression for the drift velocity A (L/T) of the center of mass of the small plume and its rate of spreading, The tensor B quantifies the strength of spreading (L/T 1/2 ), known in the literature as the displacement matrix (LaBolle et al. 1996;Lichtner et al. 2002;Salamon et al. 2006a).RWPT moves the particles in a given time interval t following the drift velocity (Equation 3) plus a Brownian motion with covariance matrix BB T t, which is written as where X p (t) is a random variable describing the position of the pth particle at time t, t is a time step of particle displacements, and ξ(t) is a vector of independent random numbers following a standard normal distribution ξ(t) ∼ N (0, 1).For large number of particles N → ∞ and small time steps t → 0, the particle density distribution p(x, t), defined as the probability per unit volume 1/L 3 of finding a particle at a given coordinate x and time t, fulfills the Fokker-Planck equation (Tompson et al. 1988;Risken 1989;Salamon et al. 2006a) Considering that p is proportional to the total mass density ρ, the Fokker-Planck equation obtained from (3), (4) and (5), and the ADE are exactly the same.Estimates of particle densities can be determined by counting the mass of the particles located within a given volume.However, this estimation only converges to the true value when the number of particles tends to infinity.Consequently, statistical fluctuations develop when a non-infinite number of particles is used.To avoid this effect, a KDE of the probability function can be used, in which case the total mass concentration is mathematically expressed as (Tompson 1993) where m p is the mass of a particle (M) located at X p and W (1/L 3 ) is a symmetric weighting function known as the kernel.This function is of unitary volume (resembling a Dirac delta), expressing that points close to the particle position are well known, and a typical choice is the Gaussian function.The kernel contribution of particles NGWA.orgR. Pérez-Illanes and D. Fernàndez-Garcia Groundwater 62, no.4: 617-634 far away from the position of interest is negligible and is common practice to consider that only particles within a finite limit distance explain most of the density.The parameter H p is the kernel bandwidth, which defines the spatial influence of a particle and controls the degree of smoothing.Earlier KDE approaches to estimate concentrations from particles considered a fixed bandwidth based on the flow grid discretization (e.g., Tompson and Dougherty 1992;Tompson 1993).However, the choice of this parameter is crucial for an accurate estimation of concentrations (Fernàndez-Garcia and Sanchez-Vila 2011).Intuitively, if the kernel size is large, small scale features of the particle distribution might be lost due to oversmoothing, and if the kernel is too small, the estimated concentrations can display a strong statistical variability.An optimal bandwidth balancing these two limits exists, obtained from minimizing the asymptotic mean integrated squared error (AMISE) of the estimation.This optimal bandwidth depends on both the number of particles used and the shape of the distribution.Recently, a locally adaptive optimal bandwidth selection was discussed in the literature (e.g., Sole-Mari et al. 2019), allowing kernels to evolve in conjuction with the solute.
In RWPT, the resolution of a simulation is controlled by the magnitude of the particle mass m p , which directly determines the number of particles.In general, simulations of higher particle resolution (lower m p ) will lead to kernels of smaller size, and a lower uncertainty on the estimated concentrations.It needs to be remarked that the displacement of particles is not influenced by this uncertainty.This means that even in simulations of low particle resolution (high m p ), the cumulative displacement of a particle followed the appropriate advection-dispersion equation, without numerical dispersion.In said case, however, the estimation of concentrations will display a higher statistical uncertainty.Once the total mass density ρ is estimated through kernels, the dissolved concentrations can then be calculated as Particles at Cell Interfaces MODPATH displaces particles to the cell-exit face in one step applying the semi-analytical algorithm developed by Pollock (1988).The source code design is closely linked to this method, transferring particles sequentially between cells.Tracking-relevant decisions are performed once the cell-exit face is known, like the analysis of the stopping conditions or the internal particle transfer for models with unstructured grids.In RWPT, knowing the exact cell-exit position of a particle is not as straightforward.Due to the random nature of displacements, even when particles are moved with very small time steps, there is a low probability that particles land on the interfaces.This article proposes a method for displacing particles exactly to the interface by computing the required time step to land on this location.This allows the reuse of most of the already implemented cell-transfer procedures.To illustrate, consider the case in which particles are moved with very small time steps such that the displacement function (Equation 5) for the x -coordinate can be written as where x p (t), A x , and (B • ξ(t)) x represent the x -components of the particle position, drift velocity and random displacement, respectively.Once the particle is close to the interface, it is possible to reinterpret the displacement equation and impose that the distance to be traveled is exactly the distance to the interface x * , and calculate the time step t * x required to satisfy this condition.This is possible by solving Equation 9 as a quadratic polynomial for t * x , with solutions if A x = 0 (10) The limit case A x = 0 occurs for example in scenarios that groundwater discharge is zero, combined with spatially uniform effective molecular diffusion, water saturation and porosity.RWPT displacements are performed incrementally inside a cell.The interfacelanding procedure is only applied once the particle passes through the interface.Near a cell corner it is possible that a particle travels through more than one interface (Figure 1), in which case the interface-landing procedure is applied for each spatial coordinate.The final time step is selected as the minimum between the crossing events.An important aspect to remark is that the direction of the displacement vector is a non-linear function of the time step, meaning that the final location where the particle leaves the cell is not aligned with the initial displacement vector.Random numbers from the initial displacement are preserved while computing the final step.
Although this method requires detecting the moment when a particle leaves the cell, it enables the reuse of cell-transfer procedures, and recomputing the final displacement is a straightforward process, thanks to Equation 10.It has been noticed that consistent solutions to the quadratic polynomial should satisfy t * x > 0. From Equation 10, one may be tempted to select the final time step as the minimum square of the two solutions (±) to the quadratic polynomial, however, in the case that t * x < 0 the particle does not land on the interface.Notice that the final time step will be always smaller than the initial one, which takes the particle outside the cell.In this context, the interface-landing procedure can be seen as a refinement of the transition near interfaces.

Hybrid Flow Interpolation Method
The RWPT form used in this development is called the hybrid scheme, which combines linear interpolation of face velocities to calculate the velocity at the particle's position, along with trilinear interpolation for calculating the divergence of local dispersion (e.g., Rizzo et al. 2019).This has been shown to provide continuity of mass fluxes at cell faces (LaBolle et al. 1996;Salamon et al. 2006a).To compute the drift displacement (Equation 3) it is first necessary to obtain the cell corner velocities, leading to a notable difference with respect to the advective transport of particles; besides the current cell velocities, the dispersion model necessitates the face-flows from neighboring cells.Traditionally, MODFLOW cells store information about surrounding cells connected through faces (Figure 2a), which are denoted as direct cell connections.To interpolate the corner values, it is necessary to extract information from cells surrounding the corner that are not directly connected to the center cell.These cells are referred to as indirect cell connections because their data can only be accessed by first initializing a direct connection.In three-dimensional problems, the center cell has access to the face-flows of 18 neighbors (refer to Figure 2b).To obtain the specific discharge in a given direction at the corner, the program adds the contribution of four surrounding cells and divides by an equivalent area (see Figure 2c).The same procedure is applied for models with rectangular unstructured grids, in which case the flows surrounding the corner are obtained differently depending on the local grid refinement conditions (refer to Supporting Information: Unstructured grid ).Notice that flow model layers are allowed to have different height while using this method.Interpolated discharges are then employed to evaluate the dispersion tensor at the corner.The porosity in Equation 2 is estimated as the volume weighted-average of the cell-porosities surrounding the corner.From here, the divergence of local dispersion (in Equation 3) is evaluated at the particle position by trilinear interpolation.

Boundary Conditions
MODPATH-RW introduces various types of particle boundary conditions that can be applied to simulate the interactions between solute transport and the surrounding environment: impermeable boundaries, time-dependent particle injection, and particles removal.As demonstrated in the study by Szymczak and Ladd (2003), the best method to simulate an impermeable boundary is the Specular Reflection Method (or elastic rebound).By default, MODPATH-RW applies this condition to flow model cell faces marked as inactive and domain limits.However, users have the option to modify this behavior and choose to fallback to the removal of particles considered by MODPATH.For the impermeable boundary, the interface-landing method is applied to detect the point of rebound, and from here there are two possible after-rebound scenarios: (i) the particle remains inside the origin cell, or (ii) the particle crosses toward a neighboring cell (see Figure 3a).In the first case, the displacement after rebound continues as usual, and in the second case, the interface-landing procedure is applied.A time-dependent particle injection can be assigned to the source cells of the flow model (e.g., an injection well, losing stream, prescribed heads with positive flow rates) to specify the mass per unit of time J in entering into the system at that location.The question here is how to configure a release of particles that ensures a consistent boundary model.One may consider that a given number of particles is released with a certain frequency, and in such a case the relation to determine the particles' mass is m p = J in δ R t /N R , where δ R t is the release interval and N R the number of particles released each time.However, the difficulty of this approach lies in the determination of a release interval δ R t that ensures continuity of the injection, and is also consistent with the temporal variability of the boundary condition.Szymczak and Ladd (2004) suggested that a constant mass flux can be simulated by an stochastic selection of the individual release time of particles, taken from a uniform distribution between zero and the simulation time step.However, in practical applications, it is most likely that the mass flux is variable in time.To overcome this, the release time of particles in MODPATH-RW is determined from the cumulative distribution function (CDF) of the injected mass (Figure 3b).The user controls the simulation resolution with an input scale of particles' mass m p , and a release template which determines the relative positions of particles inside the release cell.Briefly, a cell injecting a total mass M will require a number of particles N = M /m p to be represented.Considering a release template with N R particles, means that N i = N /N R injections are required, which can be seen as differential mass increments sequentially intersected with the CDF.The release time is then calculated by using the CDF inverse transform sampling technique.The release template specifies where particles should be initially placed during the injection.Particles can be distributed internally inside the cell or at the cell faces.Both N and N i are likely to be non-integer quantities, so an integer approximation is involved, after which the boundary model is made consistent by updating the particles mass as where t i , t f is the time interval of the injection.In MODPATH-RW, users specify this boundary condition by providing a concentration c in (t) to a source cell.The mass flux entering to the system is then estimated as the product between the flow rate and the given concentration at a particular time J in (t) = Q(t)c in (t).The advantage of this approach is that temporal changes in flow rates or concentrations are naturally incorporated by means of the release time selection, sharing some similarities with the method proposed in Szymczak and Ladd (2004).It is important to remark, that the value of m * p (Equation 11) is specific for a boundary, hence, models considering multiple injections with different characteristics (flow rates, concentrations, temporal variability), can lead to different values of m * p given the same input mass scale m p .The last solute transport boundary allows to define the removal of particles in cells with sink flow rates.Similar to classical MODPATH, users are able to specify the fate of particles upon landing on cells where groundwater is being extracted from the system.At this stage, this kind of boundary applies from a global perspective, meaning, is a model configuration affecting all of the cells with extraction flow rates.In case the user decides to remove the particles, then these cells constitute a strong sink (e.g., an extraction well, drain, prescribed heads with negative flow rates).Simply stated, this boundary imposes that the solute mass flux entering the cell is the same as the flux leaving the groundwater system.
At this stage, constant concentration and reactive boundaries are not implemented into MODPATH-RW and should be considered for future work (refer to Supporting Information: Boundary conditions).

Initial Conditions
Users have the option to specify a distribution of concentrations at the initial simulation time, which needs to be transformed to an equivalent distribution of mass particles.In MODPATH-RW, this can be performed with different approaches.Here, the method based on a uniform distribution of particles inside cells with non-zero concentration is outlined.Given an initial particle mass m p , the initial concentration c 0 is transformed into an estimated number of particles per cell as N c = c 0 φR s V c S w /m p , where V c is the cell volume.The cell aspect ratio is employed to define a release template that distributes these particles inside the cell.In direction i , the cell aspect ratio is given by s i = i /(S w V c ) 1/d , where i is the cell size and d is the number of dimensions.For the vertical direction, the aspect ratio is obtained with the cell size also accounting for the water saturation.Then, the number of template sub-divisions in direction i is obtained as c , which undergoes an integer approximation.Similar to the case of prescribed-flux boundaries, the final value of particles' mass is updated to ensure consistency between the total mass represented by particles and the conceptual initial condition as Attention should be paid to the selected particles' mass, because it will determine both the total number of particles and the agreement between the discrete (particles) and initial continuous concentration distribution.Other implemented approaches for defining the initial condition allow to release particles placed randomly or quasi-randomly inside each cell.Also, mass consistency can be achieved by recalculating the mass of particles on a per-cell basis (similar to Equation 12), eventually leading to a set of particles with non-uniform mass.

Concentration Reconstruction
A major challenge in RWPT models is the interpretation of the particle distribution as a continuous concentration.Methods for mixing-controlled reactive transport are strongly dependent on the estimation of solute gradients (e.g., Cirpka and Valocchi 2007;De Simoni et al. 2007), and the statistical variability of simple concentration estimates (histogram) for simulations with insufficiently large number of particles can have a negative impact on reaction related calculations (Fernàndez-Garcia and Sanchez-Vila 2011).The effect can be attenuated with smoothed concentration reconstruction, which is integrated into MODPATH-RW by the incorporation of an external Fortran module that performs GPKDE.The latter employs cell-averaged kernel functions with automatic bandwidth selection, combining the speed of the histogram reconstruction with the ability to smooth concentrations using locally adaptive KDE (as in Sole-Mari et al. 2019).GPKDE is employed for both spatial and timeseries reconstruction of concentrations, taking into account kernels corrected by reflection when they intersect the domain limits.This functionality allows the program to report output results as smoothed concentrations and not necessarily as individual particle positions (like in traditional MODPATH).
Furthermore, the smoothed reconstruction reduces the statistical variability of simple concentration estimates (histogram) reported since early applications of RWPT (e.g., Ahlstrom et al. 1977;Tompson et al. 1988).

Observation Cells
In MODPATH-RW, the user can specify cells or group of cells within the numerical groundwater flow model where the simulated particle information and corresponding concentrations are monitored during the simulation.Two types of observations are implemented in the program.The first type of observation monitors resident concentrations (see Parker and van Genuchten 1984;Kreft and Zuber 1986).These concentrations, denoted as c r , are calculated as R s c r = M obs /(S w φV) obs , where M obs is the total particle mass contained inside the group of cells assigned to the observation, and (S w φV) obs is the corresponding total water volume.The second type of observation estimates flux concentrations and is specifically applicable to sink flow cells that remove particles (strong sinks) from the system.In this case, the program records the first arrival time of particles and their particle mass at observation cells during the simulation.From here, the program performs a timeseries reconstruction of the mass per unit of time leaving the system J out at that location using the GPKDE module.In combination with the total extraction flow rate Q out , the mass flux concentration is then estimated as c f = J out /Q out .In this sense, the immediate removal of particles can sometimes overestimate the mass outflux in model cells with discrepancies on their size with respect to the actual sink size (e.g., a large cell simulating a small well).For these cases, studies have suggested to include the uncertainty in the decision of removing a particle (e.g., Visser et al. 2009;Abrams et al. 2012), which is aligned with the consideration that solute particles may re-enter the domain due to their random motion (e.g., Tompson et al. 1988;Lattanzi et al. 2019).

Backward Tracking
Backward tracking is also implemented in MODPATH-RW.This is achieved by reversing the direction of face-flow velocities at each cell, while preserving the same dispersion terms as in forward tracking (Bagtzoglou et al. 1992;Neupauer and Wilson 1999).In the context of solute transport, random walk backward tracking allows to determine the probability density functions of the origin of a contaminant source, both from a spatial or temporal perspective (Neupauer andWilson 1999, 2004).The primary distinction introduced by backward tracking in the solute transport components of the program pertains to the particle injection in source cells and the removal of particles in sink cells.In the first case, cells that withdraw water from the aquifer can now be designated as particle injection cells, whereas cells injecting water into the aquifer are defined as sink cells.Thus, in backward tracking, a consistent sink cell should take into account that particles are removed at cells with source flows.It is crucial to emphasize that the NGWA.orgR. Pérez-Illanes and D. Fernàndez-Garcia Groundwater 62, no.4: 617-634 key output during backward tracking is the Probability Density Function (PDF) of the concentration, rather than the estimated magnitudes.Consequently, any concentration distribution obtained from backward tracking must be normalized to integrate to one to facilitate interpretation.

Time Step for Displacements
The default setting in MODPATH-RW entails an independent time step for each particle, which is calculated on a per-cell basis.One of the criteria for selecting the time step relies on advection, aiming to maintain a predefined Courant-Friedrichs-Lewy (CFL) number, , where v i is the maximum absolute groundwater velocity obtained from the cell faces oriented perpendicular to the i th direction, and i is the cell size.The second time step selection method is based on a dispersion criteria, determined as where C T is a prescribed constant with typical values of C T ≈ 0.1, and D ii is the diagonal component of dispersion in the i th direction.The program considers the previous time step selection methods independently or selecting the minimum between them, according to the user input.Another option available to users is the ability to manually specify a constant time step for all particles.

Results and Discussion
This section presents the results of numerical test cases that cover a variety of solute transport applications in one, two, and three dimensions.MODPATH-RW is benchmarked with equivalent simulations carried out with the solute transport module of MODFLOW@6.4.1, simply referred to as GWT (Langevin et al. 2022).Additionally, some tests are compared with simulations using MT3D-USGS@1.1.0(Bedekar et al. 2016), simply referred to as MT3D.For all simulations involving GWT and MT3D, the total variation diminishing (TVD) advection scheme is consistently employed.

Pulse Injection into Porous Medium
The first test considers a semi-infinite onedimensional homogeneous porous medium column with constant prescribed mass flux and zero initial concentration (Lindstrom et al. 1967;Parker and van Genuchten 1984).Injection has a duration of t 0 , with inflow rate Q and concentration c 0 .During the injection, c 0 and the resident concentration c r at the transition zone are related by the expression.with D being the coefficient of hydrodynamic dispersion, v the flow velocity, and the ratio D/v the apparent dispersivity.The analytical solution for this problem is given in Appendix (Equation A.1). Numerical simulations consider a constant velocity and different values of the grid-cell Péclet number Pe = x v/D, well below advection-dominated transport (Pe < 0.5).For these conditions, it is expected that Eulerian grid-based models outperform the random walk particle method.
The injection time and the final simulation time T satisfy respectively that vt 0 = 0.5 and vT = 1.Numerical results from both GWT and MODPATH-RW present a good agreement with the analytical solution (Figure 4).For the random walk particle model, the time step selection plays a relevant role and values C T , CFL ≈ 0.05 show a satisfactory agreement for concentration magnitudes and arrival times.Increasing the particle resolution leads to a consistent decrease in error with respect to the analytical solution (Figure 4b), supporting the suitability of the method for modeling the prescribed mass flux boundary condition.However, in all cases, errors obtained with the random walk model where higher than those with GWT, approaching similar values with an increasing number of particles.Notice that the GPKDE consistently leads to smaller statistical fluctuations than the base histogram (counting particles inside cells, RW-HIST in Figure 4b).However, as the number of particles increases, these fluctuations become closer in magnitude, eventually reaching a point where the resolution provided by particles is so high that reconstruction might not be needed (i.e., infinite particles).Still, it has been observed that although close to the analytical solution, histogram reconstruction with a large number of particles presents some local fluctuations smoothed by kernels.An aspect to remark from these results, is the characterization of concentrations near the domain origin (Figure 4a).Here, boundary conditions from both the particle model and kernel reconstruction are playing a role on the values estimated near the interface.
A second similar test is analyzed using the same onedimensional column setup.Now the inflow concentration is variable in time following the expression where c a and c b are concentration scales (taken as c a = c b = 0.5), and λ is a constant that modulates the increase/decrease of concentrations in time.The analytical solution for this problem is presented in van Genuchten and Alves (1982), also given in Appendix (Equation A.3).This problem illustrates a more practical application with a variable in time prescribed mass flux injection.Simulations with MODPATH-RW for different magnitudes of λ (in Equation 14), illustrating different time scales of the inflow concentration decaying function, are in good agreement with the analytical solution (Figure 5).

Solute Transport in Heterogeneous Aquifer
This test considers an instantaneous release of a rectangular unit concentration into a two-dimensional heterogeneous aquifer subjected to a steady-state mean flow oriented along the x-coordinate.The parameters adopted for this test are summarized in Table 1.The spatial distribution of the hydraulic conductivity K( x   which quantifies the magnitude of the flux over which the solute is distributed.This metric is computed through the specific discharge of the main flow direction q x as with units of inverse flow rate.The index is normalized by the total flow rate ẼQ = E Q /Q and its temporal trend is discussed with respect to the non-dimensional time t = tv/I Y , where v is the mean flow velocity.
Simulations are executed until a total time of T ≈ 44t.For all models, a Courant-Friedrichs-Lewy number of CFL = 0.1 is imposed.In GWT, this is done by defining the number of steps, calculating the time step from the maximum spatial velocity.
Transport conditions are particularly challenging for GWT and MT3D due to the heterogeneous nature of the random conductivity field, and some influence of numerical dispersion is expected, particularly in the case of high heterogeneity (large spatial variability of velocities).Finite-differences codes should perform best in dispersion-dominated problems with low grid-cell Péclet numbers.So, in this setup, it is expected that the random walk particle model outperforms the Eulerian grid-based programs, provided that simulations are carried out with enough particle resolution.In this regard, particle simulations are performed with two different magnitudes of resolution (N 1 = 0.28M and N 2 = 7M).
Results indicate that MODPATH-RW provides the smaller magnitudes of dilution for all test conditions (Figure 6).For low degrees of heterogeneity, the random walk model is generally in agreement with MT3D, both in the magnitude and temporal trend of the dilution index (Figure 6a and 6c).Flow variability appears to be the key factor controlling the results obtained from GWT, particularly for high degrees of heterogeneity (Figure 6b  and 6d), where the temporal trend of dilution presents low sensitivity with respect to the considered grid-Péclet numbers.In this last case, simulations with MODPATH-RW present significantly lower magnitudes of dilution, and the particle resolution mostly plays a role in determining the short-term temporal variability of the index.MT3D is able to provide metric values similar to MODPATH-RW for high degrees of heterogeneity and Pe = 5.However, important differences appear after increasing the grid-Péclet number with dilution predicted by GWT being ≈3.9 times larger than the obtained from MODPATH-RW at the moment of highest difference (vertical line in Figure 6d), and at the same instant, the discrepancy of MT3D is ≈2.1 times.A visual representation of these effects is presented in Figure 7. Notice that results from GWT appear to be influenced by numerical dispersion with lower peak concentrations and small-scale values surrounding the main solute branches.To a lesser extent, simulation with MT3D is also influenced by similar effects, which are not observed in the concentration distributions from the random walk particle model.The last activity performed with this test is a comparison of runtimes between the different transport codes, for one of the simulated scenarios.The runtime of simulations on the aquifer of high heterogeneity (σ 2 Y = 2.5), and with transport strongly dominated by advection (Pe = 50), is considered.Simulations are performed on a desktop computer equipped with an Intel ® Core™ i7-9700 CPU @ 3.00GHz processor.All model runs are considered to be equivalent in terms of the transport problem, also considering the same number of concentration maps written to the output files (150 equispaced snapshots).MODPATH-RW is compiled with gfortran@11.4.0 and simulations are executed in parallel modifying the number of threads N thr ∈ {1,2,4,8}, for both magnitudes of resolution (N 1 = 0.28M and N 2 = 7M).Concentration reconstruction is performed by pre-allocating a database of kernels, which is of benefit for models requiring a frequent reconstruction (Sole-Mari et al. 2019).Simulations with GWT and MT3D are executed with a single processor on the same desktop computer, and in both cases the runtime considers only the execution of the transport model.For GWT, groundwater flow is assigned using the Flow Model Interface (FMI) package.In all cases, the concentration maps shown in Figure 7 are compared against a high resolution distribution generated with random walk (N HR = 50M), via the root mean squared error e R .Results are shown in Figure 8.For both magnitudes of particle resolution (N 1 , N 2 ), the random walk model displays competitive runtimes with respect to the reference codes.Interestingly, this occurs even for the simulations performed with a single thread (N thr = 1), and runtimes can be further reduced by taking advantage of the parallel implementation.For further comparison, simulations with the finite-differences codes were also executed with CFL = 1, without significant changes in the magnitude of differences with respect to the reference solution (e R ).Overall, the random walk model always provided the smaller values of e R , which illustrates the potential of the program to generate accurate results in aquifer models with spatial variability of the flow velocities and advection-dominated transport.Furthermore, these can be generated with runtimes that are competitive with respect to well-established Eulerian alternatives.

Transport near Low Permeability Zone
This test is based on a well-known benchmark solution from MT3DMS (Zheng and Wang 1999).The transport problem deals with a constant injection into a two-dimensional aquifer, characterized by a main hydraulic conductivity denoted as K m , along with an embedded low permeability zone of hydraulic conductivity K l .The low permeability zone partially intersects the solute plume as it moves downstream toward an extraction well (refer to Figure 9).The problem can be seen as a simplified heterogeneous system, with a single hydraulic conductivity contrast.The flow model grid for the base test is structured with a uniform cell size, and it is further modified with refinement features to demonstrate the suitability of MODPATH-RW for unstructured grids (Figure 9b).To ensure conceptual consistency between both aquifer models, multiple cells are assigned for injection/extraction wells in the unstructured case.Groundwater flow is induced by prescribed head cells in the northern and southern domain boundaries.All flow model boundaries remain constant during the simulation.Groundwater flow is assumed to be at steady-state.Flow and transport parameters are summarized in Table 2.
A relevant output of this problem is the concentration breakthrough curve reported by models at the outflow well W out .The MT3DMS concentration breakthrough curve is obtained from Zheng and Wang (1999), with models solved on the original structured grid.In this sense, the random walk particle model considers an observation at the extraction well reporting flux concentrations.In contrast, both GWT and MT3DMS consider by default that the concentration at the well represents the value predicted for the corresponding aquifer cell (Zheng and Wang 1999;Langevin et al. 2022), which holds more resemblance to a resident concentration.In MODPATH-RW, it is difficult to conceive a resident concentration at the extraction well since no particles reside in this area at any given moment.It is important to highlight NGWA.orgR. Pérez-Illanes and D. Fernàndez-Garcia Groundwater 62, no.4: 617-634  this topic as it seems to be the primary factor explaining the significant difference in concentration magnitudes obtained from the random walk particle model versus the grid-based programs (see Figure 10).Additionally, this difference may also be influenced by the immediate removal of particles at the well in the random walk model.Despite these considerations, it is noteworthy that the estimated concentrations using the particle model remained independent of the flow model discretization, showing consistent arrival and peak concentrations.
A final exercise is performed with this test case to illustrate the applicability of backward tracking in random walk models.As discussed previously, this feature can be used to estimate the probability that a contaminant source originates at a certain location and time (e.g., Neupauer and Wilson 2004).To demonstrate this, the simulated concentration breakthrough curve at the outflow well (W out , original structured grid) is now given as input concentration at the same location and particles are tracked considering groundwater flow in reversed (backward) direction.The injection well (W in ) now becomes a sink observation and the corresponding computed concentrations are interpreted as the probability density function of the injection p c (Figure 11).Results indicate that the probability density function is higher for times between t ∈ [0, 1] years, which is consistent with the real problem.Furthermore, within this range of times, the curve exhibits an extended plateau, indicating a sustained injection over time.Of course, this is a simple problem involving only a single injection well.However, in more intricate setups involving multiple injections and an unknown source, a comprehensive analysis must be conducted with observations from all existing injection wells to ascertain the most probable source of concentrations.

Field Scale Aquifer Remediation
This last test case is also based on a well-known example of MT3DMS (Zheng and Wang 1999).The test case exemplifies the application of a remediation scheme in a three-dimensional field scale aquifer.The setup represents a more realistic application and exposes the program to conditions potentially found in real-world problems.Particularly, the use of non-regular grids on a quasi two-dimensional system, and the specification of a non-uniform initial condition for a solute subjected to sorption effects.The model aquifer is composed of four layers of the same thickness.The two upper layers represent a sandy material of medium-to-fine grain sizes with hydraulic conductivity K f , and the two lower layers consist of coarser grains and hydraulic conductivity K c .The flow model grid is structured, with cells of the same horizontal extent inside a perimeter where the solute plume has been detected.The cell size increases as one moves toward the region boundaries.The hydraulic head is prescribed at the domain boundaries, and the system experiences a homogeneous vertical recharge.All flow model boundaries remain constant during the simulation for T = 1000 [d].The remediation scheme is implemented by means of 8 extraction wells with screen opened at layer 3 (from top to bottom), where a non-uniform distribution of 1,2-dicloroethane (DCA) was detected (Figure 12).Higher concentrations are found in layer 3. Layer 2 has the same   3.
While studying this test, several aspects of the program were evaluated.First and foremost, the smoothed reconstruction of concentrations for three-dimensional domains was assessed.Particularly, the domain of this problem is quasi two-dimensional, and the reconstruction module limits the size of some isotropic kernels employed during the automatic bandwidth selection to avoid issues related to boundary corrections.This means that when performing a purely three-dimensional reconstruction, some of these kernels will most likely adopt smaller values than might be required by the horizontal characteristics of the particle distribution, limited by the vertical extent of the domain.The latter can lead to an inconsistent reconstruction.To overcome this issue, the module was extended to handle the reconstruction of quasi twodimensional problems on a per-layer basis.
Another important aspect is the influence of the particles' mass on transport problems with non-uniform initial conditions.This parameter controls the closeness between the nominal concentration map and the particle distribution.Large values of mass (low resolution) can lead to inconsistent concentration maps (contours).To demonstrate this, particle simulations were performed with different magnitudes of particles' mass.Concentration maps for layer 3 obtained at time t = 500[d] were compared with simulation results from GWT and MT3D, the latter employing the full form of dispersion (Figure 13).Results show that maximum resident concentrations obtained from high resolution MODPATH-RW runs (Figure 13c and 13d) are smaller in magnitude than those obtained from GWT and MT3D.Remarkably, the contours of low concentrations (≈ 5[ppb]) exhibit more similarities with results from MT3D than with GWT.The effect of particles' mass is clear: larger particle mass values lead to lower number of particles and eventually the contours of concentration appear oversmoothed (Figure 13e and 13f), consequence of the higher uncertainty on the reconstruction of concentrations.This indicates that the choice of the particle mass should also be consistent with the expected resolution of concentrations, meaning that this selection is not completely arbitrary.To illustrate this point, Table 4 provides an estimated resolution for the resident concentration, obtained from considering a single mass particle inside one of the cells within the perimeter with regular sizes (min c r = m p /R s φV c ), for various magnitudes of particles' mass.For this particular problem, mass values below m p = 1000[mg] would guarantee a resolution in resident concentrations below 1 [ppb].This elucidates the discrepancies in contours observed in simulations with low particle resolution.
The fact that the maximum concentrations in the random walk particle simulations of high resolution are smaller than those obtained with the finite-differences codes is, at first sight, counter-intuitive when considering results from previous test cases.To further investigate this, the solute breakthrough curve at one reference well (W 4 ) is analyzed (Figure 14a).It is observed that, during early simulation times (t ≈ 5[d]), concentrations with MODPATH-RW are large (>100[ppb]) in comparison    to the other models.This suggests that once the simulation begins, almost instantly, a very large number of particles is removed from those cells were wells are located.Similar to the previous test case, flow model cells are considerably larger in extent than a real-world well.However, at late simulation times, breakthrough concentrations from the random walk particle model overcome the magnitude of those obtained from GWT and MT3D.For comparison purposes, simulations are also executed without the influence of sorption (R s = 1, Figure 14b) and in both cases the qualitative behavior of breakthrough curves is practically the same.

Conclusions
This article discusses the integration of RWPT into MODPATH, achieved through the extension of the program's source code.The main motivation for adopting this approach (instead of developing from the ground-up) was the inherent interoperability of MODPATH with a variety of MODFLOW models.Some interesting features already implemented in the base program were extended and adapted for random walk, most notably, the compatibility with rectangular unstructured grids (for models whose layers possess the same horizontal structure of cells) and the backward tracking of particles.By implementing a methodology to compute the exact cell-exit position of a particle displaced by advection and dispersion, the random walk functions were successfully integrated, making use of most of the cell-transfer procedures upon which the base program is built.The program implements solute related boundary conditions including the impermeable, the time-dependent solute mass flux and the removal of particles at sink cells.Prescribed concentration and reactive boundaries are not yet implemented into the code.Additional functionalities like the observation cells and the smoothed reconstruction of concentrations, contribute significantly to the necessary workflow for applying the RWPT method to reactive transport studies.Numerical test cases have illustrated the applicability of the code, benchmarked with analytical solutions and well-known Eulerian grid-based programs.In problems with an explicit representation of heterogeneities, the new model provided solute distributions without the influence of numerical dispersion and similar dilution metrics for different magnitudes of particle resolution.Observation cells illustrated the estimation of mass flux concentrations at extraction wells and the estimation of the injection probability density function when employing backward tracking.One aspect to consider when removing particles immediately at sink cells, is that the magnitude of concentrations can be overestimated in models where cells representing wells do not accurately reflect the real well size.Future research should consider additional approaches to determine the fate of particles arriving at extraction cells.In reactive transport, one of the main (if not the main) components controlling the accuracy of results, is the solute transport model.Real aquifer systems are inherently heterogeneous, posing a wellknown challenge for grid-based programs.To contribute in addressing these issues, this article has introduced the MODPATH-RW program as a particle-based alternative for modeling solute transport in heterogeneous systems.

Figure 1 .
Figure 1.Diagram of the transition of particles from the current cell to a neighbor.Gray arrows indicate the initial displacement obtained with the uncorrected time step ( t 0 ), and black arrows the actual displacement after forcing particles to the interface.Particles p 1 and p 3 cross a single interface whereas p 2 crosses two interfaces and thus the final time step must be selected as the minimum between the two crossing events.

Figure 2 .
Figure 2. Cell diagrams relevant to the RWPT implementation.(a) MODPATH cell indicating face and corner id's, (b) neighbor cells required for interpolating specific discharge to the corners of the MODPATH cell.Black cells are direct connections and the rest indirect.(c) Schematic of the interpolation procedure for corner c 101 .Each surrounding cell contributes a quarter of the total flow rate through the face and specific discharge q c 101 is computed as the sum of contributions divided by the area A c 101 .

Figure 3 .
Figure 3. (a) After rebound scenarios for impermeable walls.Particle p 1 remains within the current cell and p 2 is transferred with modified time step t * * .Difference in exit point coordinates is exaggerated.(b) Diagram of cumulative mass function for determination of release times in prescribed-flux boundaries.

Figure 4 .
Figure 4. Results for the problem of constant injection into porous media.(a) Normalized resident concentration c = c r /c 0 for different apparent dispersivities.(b) Root mean squared error of numerical results with respect to the analytical solution.RW-HIST and RW-GPKDE represent the error curves obtained with concentrations estimated from a simple histogram and kernels, respectively.Horizontal lines are estimated errors from GWT.
) is obtained from one realization of a sequential Gaussian simulation Y (x) following an isotropic exponential variogram with correlation lengthI Y = 10[m].The problem is discussed for two degrees of aquifer heterogeneities (low and high) controlled by log-conductivity variances σ 2 Y ∈ {0.25,2.5}.The hydraulic conductivity distribution is obtained as K(x) = exp(σ Y Y (x)).In both cases, groundwater flow is solved with MODFLOW by imposing a constant mean-unit head gradient.The grid-cell Péclet number is set to Pe ∈ {5, 50} by considering two scenarios of longitudinal dispersivity α L ∈ {0.2,0.02}[m].The anisotropy of dispersion, defined in terms of the relative magnitude of dispersivities, is always taken to be α L /α T = 10.Concentration maps obtained from the random walk model are compared with equivalent simulations performed with GWT and MT3D.The performance of simulations is discussed in terms of the flux-related dilution index (Rolle and Kitanidis 2014),

Figure 5 .
Figure 5. Results for injection with variable inflow concentration until time vT = 1 and Pe = 0.5.(a) Inflow concentration for different values of λT.(b) Comparison between analytical solution and profiles obtained with RW-GPKDE.Simulations consider N = 0.61M particles.

)
NGWA.org R. Pérez-Illanes and D. Fernàndez-Garcia Groundwater 62, no.4: 617-634 for a reference coordinate x R = 200[m] (∼ 17I Y downstream the solute release area), and integrals performed over the cross-sectional area .The index has units of flow rate, and p Q (x, t) is a probability density function of the flux, given by

Figure 6 .
Figure 6.Comparison of flux-related dilution index at reference section x R = 200[m], obtained with GWT, MT3D, and MODPATH-RW employing N 1 = 0.28M (light-gray) and N 2 = 7M particles (dark-gray).Index obtained from histogram reconstruction is shown in dotted lines and the corresponding GPKDE is shown in solid format (-). Figure columns group results for a given degree of aquifer heterogeneity and the Péclet number is indicated in each respective panel.Vertical line in panel (d) marks the time of snapshots shown in Figure 7.

Figure 7 .
Figure 7. Solute distribution obtained with (a) GWT, (b) MT3D, and MODPATH-RW considering the smoothed concentration reconstruction from (c) N 1 = 0.28M and (d) N 2 = 7M particles.Dashed vertical line at x R = 200[m] indicates the cross-section where the flux-related dilution index is computed.Rectangular perimeter illustrates the initial injection.Concentrations are shown until a minimum value of c = 1 × 10 −4 .

Figure 8 .
Figure 8.Comparison of computational performance in two-dimensional aquifer of high heterogeneity (σ 2 Y = 2.5) and Pe = 50.Runtimes (T ) are shown in the left-axis, and the root mean squared error (e R ) of reference concentrations maps (Figure 7) with respect to a high resolution MODPATH-RW run (N HR = 50M) is shown in the rightaxis.The runtime of the finite-differences codes are shown for CFL ∈ {0.1, 1}, and for the random walk particle model for different number of threads N thr .e R corresponds to simulations with CFL = 0.1.

Figure 9 .
Figure 9. Two-dimensional aquifer for the problem of transport near a low permeability zone.(a) Original structured grid, (b) modified unstructured.Colormap indicates the steady state head distribution.W in and W out are the well cells for solute injection and extraction respectively.Main medium conductivity is K m and low conductivity is K l .White lines illustrate advective streamlines released from the injection well, and scatter points are illustrative solute particles at the reference time T in = 1[yr].

Figure 10 .
Figure 10.Concentrations at the extraction well W out .Data from MT3DMS is partially transparent, obtained from Zheng and Wang (1999), representing solutions obtained with the hybrid method of characteristics (HMOC), third-order TVD (ULTIMATE), upstream finitedifferences (UPSTREAM-FD) and central finite-differences (CENTRAL-FD).

Figure 11 .
Figure 11.Probability density function of the injection at well W in (p c , left-axis).Curve is obtained from backward tracking of the flux concentration at W out from the original forward model (structured grid, right-axis).Shaded area between zero and T in marks the period of injection in the forward model.

Figure 12 .
Figure 12.Initial condition of concentrations at layer 3 for the remediation in three-dimensional aquifer test.Extraction wells are indicated with black circles W i .Contour lines are equispaced every 30[ppb].Dashed perimeter indicates the area with regular size cells, where DCA was detected.
m p is the input scale, m * p is the corrected mass that ensures consistency with the nominal initial condition, N the number of particles and min c r the resident concentration of one particle inside a cell within the perimeter with regular sizes (1[ppb] = 1 mg/m 3 ).

Figure 14 .
Figure 14.Breakthrough curves at extraction well W 4 obtained from the different models.Panel (a) show results while considering sorption and (b) without sorption.Results from MODPATH-RW are obtained considering m p = 100[mg].Vertical dashed line indicates the reference time t = 500[d].

Table 3 Model Parameters for Test Case Exemplifying Remediation in Field Scale Aquifer
s = 1 + ρ b κ d /φ ≈ 2,where ρ b is the aquifer bulk density, κ d the distribution coefficient, and φ is the aquifer porosity.Model parameters for this problem are summarized in Table