A Lagrangian Method for Reactive Transport with Solid/Aqueous Chemical Phase Interaction

A significant drawback of Lagrangian (particle-tracking) reactive transport models has been their inability to properly simulate interactions between solid and liquid chemical phases, such as dissolution and precipitation reactions. This work addresses that problem by implementing a mass-transfer algorithm between mobile and immobile sets of particles that allows aqueous species of reactant that are undergoing transport to interact with stationary solid species. This mass-transfer algorithm is demonstrated to solve the diffusion equation and thus does not introduce any spurious mixing. The algorithm is capable of simulating an arbitrarily small level of diffusion, and can be combined with diffusive random walks to simulate the desired level of diffusion in a reactive transport system.


Introduction
Reactive transport models fall into two main categories: Eulerian (gridbased) models (e.g., finite difference, finite volume, and finite element) or Lagrangian models (e.g., particle-tracking (PT)). For a comparison between Eulerian and Lagrangian reactive transport methods, in terms of accuracy and computational gains, see [1]. Eulerian models, as used in research and industry, are far more common, but, for certain problems, they are known to suffer from important drawbacks due to their use of a static grid. Namely, concentration fluctuations occurring at a length scale smaller than that of the grid are not captured. This imposes a cutoff point for the model, such that heterogeneities and fluctuations may be resolved between grid elements, but concentrations are perfectly-mixed within each grid element. Similarly, if an Eulerian model user wishes to capture smaller-scale fluctuations, then a more finely discretized grid must be employed, leading to longer computation times. This discretizationdependent homogeneity assumption has important implications for the modeling of chemical reactions undergoing transport via diffusion and/or advection. The traditional way to handle this has been by empirical adjustments to reaction rate constants [2,3,4]. Another known shortcoming of Eulerian methods is their introduction of numerical diffusion in the simulation of advection [5,6,7], a phenomenon which does not occur in PT methods.
In contrast, Lagrangian reactive transport models do not use a static grid that reactants move through. Instead, concentrations (or masses) of reactant are carried with the elements (particles) as the particles, themselves, move through the computational domain. For this reason, particle methods impose no cutoff resolution and can capture heterogeneities and fluctuations on a fully continuous scale. This may be achieved both by allowing for mass differences between particles [e.g., 8,9], and letting the spatial distribution of particles to change with time. The heterogeneity captured by inter-particle spacings is related to the level of discretization (i.e., particle number), and this method of representing spatial fluctuations of the concentrations is not possible with Eulerian methods. Another valuable property of these PT algorithms is that systems with a higher level of concentration heterogeneity (i.e., high spatial variance in concentration) may be modeled using fewer particles than would be required for a more homogeneous (or "smoother") system, as the number of particles is directly tied to a system's degree of spatial concentration variability [10].
However, this efficiency in capturing heterogeneity can also be a disadvantage when one needs to model a system including a very smooth, or homogeneous, concentration field. In order to properly resolve a very smooth field, classic particle methods would require a very large number of particles, leading to high computational burden (number of particles is analogous to number of grid points in an Eulerian model). Two methods for addressing this problem have recently been proposed. First, a spatial kernel with nonzero support may be employed (a Gaussian kernel is the standard choice, due to its natural compatibility with the diffusion process), in contrast to the classic Dirac delta kernel. This allows for significant reduction in particle number while introducing a controllably low level of error [11]. Alternatively, for complex reactions involving many chemical species, particle number may be reduced by allowing each particle to carry an arbitrary number of different species (rather than a single species per particle) that are transferred between/among particles according to a diffusive masstransfer algorithm [9,12,13].
Another shortcoming of Lagrangian PT methods, which we address in this paper, is their inability to model interactions between solid and aqueous chem-ical phases (e.g., dissolution or precipitation reactions, wherein flowing liquid carrying aqueous chemicals interacts with the surrounding rock). The difficulty of this problem lies at the core of the Lagrangian framework because the particles in a PT simulation move, but stationary solids must not. A naïve approach to this problem may employ a "reaction grid" that particles carrying aqueous species move through, and any particles within a given grid element are considered to be a well-mixed reaction volume. However, this type of approach imposes the previously-described over-mixing problem. Thus, we propose a framework in which mobile particles (carrying aqueous chemical phases) interact with immobile particles (carrying solid chemical phases) through a diffusive mass-transfer (MT) process that imposes no spurious mixing.
First, in Section 2, we define the reactive transport problem considered herein and describe the mobile-immobile reactive particle-tracking (miRPT) algorithm used to model the problem. In Section 2.2, the mass-transfer algorithm that is central to the miRPT model is derived in detail, and computational considerations necessary for simulating such a model are presented. We derive a "diffusion operator" in Section 2.3 that is a discretized Green's function solution and use the solutions given by this diffusion operator for comparison with the miRPT MT algorithm throughout Section 3. In Section 3.1, we first examine the behavior of the miRPT MT algorithm in isolation (without reaction) and demonstrate that it solves the diffusion equation with a controllable level of accuracy. We also examine the effect of the various discretization parameters (number of mobile or immobile particles and time step) on the accuracy of the miRPT MT algorithm and describe a stability condition that should be taken into consideration when implementing the algorithm. Finally, in Section 3.2, we discuss the results of using the miRPT algorithm to model a chemicallycomplex reactive transport model that involves calcite dissolution and dolomite precipitation/dissolution. We describe the process of choosing the correct discretization parameters for this system in Appendix B.1, and we give some important computational details of the model in Appendix B.2.

Governing equations
We consider a reactive system of arbitrary chemical complexity that is undergoing transport by advection and diffusion. For an infinite, d-dimensional domain (d = 1, 2, 3), the governing (system of) equation(s) for this reactive transport system is the advection-diffusion-reaction equation (ADRE) where the concentration of species i is given by , and reaction rate constant , respectively, for species i and reaction channel j = 1, . . . , m. The reaction term, r i (·, . . . , ·) is some, possibly nonlinear, function assumed to be governed by the law of mass action. We note that for a solid species (not undergoing transport) v = D ≡ 0, as concentrations will only be altered by (dissolution or precipitation) reactions. As well, for the purposes of this paper, we assume that, for aqueous species, v and D are constant in both space and time.

Particle tracking model
Benson and Meerschaert [14] developed their reactive particle-tracking (RPT) method in order to simulate chemical reaction in diffusive media. The RPT algorithm models reactions by calculating the probability of reaction between particle pairs that is conditioned upon the event of particles co-locating by molecular diffusion. Under this paradigm, each particle carries a given mass of only one type of reactant. When a pair of particles react they are "killed," converting the mass that they carry to product, and the mass of the corresponding species are thus altered. This approach has been extended beyond the particlekilling approach to assign continuously-varying masses to the particles that are reduced in proportion to their co-location probability [15]. It has also been extended to allow individual particles to carry multiple species of reactant (obviating the need for a distinct, large set of particles for each individual species) and exchange mass among one another, thus allowing reactions to occur within a given particle that is treated as a well-mixed reaction volume [9].
It is this final, arbitrarily-complex reactive particle-tracking (aRPT) algorithm that we base our work on in this paper. It was recently shown that the mass-transfer (MT) process employed by aRPT solves the diffusion equation with O(∆t) accuracy [13]. Thus, the total diffusion in a system may be partitioned between the random walk and MT numerical processes with a controllable level of error. Benson and Bolster [9] attempted to simulate mobile-immobile particle interaction with a naïve implementation of the aRPT algorithm that resulted in large errors, attributable to excess numerical diffusion. However, we modify this MT algorithm to handle interactions between aqueous and solid (mobile/immobile) phases of reactant and thus formulate our approach to a mobile-immobile RPT (miRPT) algorithm in the following section. Additionally, we show that the altered MT algorithm also solves the diffusion equation with minimal error.

Mobile-immobile reactive particle tracking
The essential motivation behind the miRPT model is that the algorithm employs two sets of particles (mobile and immobile) because aqueous and solid species (henceforth referred to as mobile/immobile mass, respectively, within the context of the algorithm) must interact with one another, but solid species must not be transported. As a result, all chemical reaction calculations are performed by treating the immobile particles as a well-mixed volume and allowing the mobile mass to be transferred between the mobile and immobile particles within Step in time a time step according to an MT algorithm that is shown to be diffusive. The size of this well-mixed volume may be defined by the user with the chemical reactions in mind, as opposed to, in the Eulerian paradigm, choosing the wellmixed volume size based on constraints imposed by the transport algorithm (i.e., bounds on the Péclet or Courant number). A step-by-step outline of the miRPT algorithm is given by the flowchart in Figure 1.
We specify here that, for the miRPT algorithm described in this paper, all of the mobile mass is transfered to the immobile particles for reaction and back to the mobile particles for transport (steps (c) and (g) in Figure 1). Presumably, this could be modified by performing reaction calculations for interactions between mobile and immobile particles; however, this would require precise control over both the reaction calculations and over which reactions are calculated between certain particles. For this paper, we wish to consider the most general case, wherein all reactions are calculated simultaneously among all species of reactant. This treatment allows for the use of geochemical solvers such as PHREEQC (PhreeqcRM), CrunchFlow, or Reaktoro [16,17,18].

Algorithm details
Particle-tracking algorithms represent concentrations by employing a spatial kernel. The most common choice of kernel is a Dirac delta, though other choices have been shown to be effective, particularly if a smaller number of particles is desired for a given initial condition [11,19]. We note here the distinction that, due to the infinitesimal support of a Dirac delta, particles typically carry masses, as opposed to the concentrations tracked at grid locations of analogous Eulerian methods. However, concentrations may be easily generated by grouping particles into "bins" of fixed volume (in this case, the immobile particles represent a fixed volume).
Concentrations, in phase p = M, I (mobile or immobile), are represented in space in the following way where Ω is the d-dimensional domain (d = 1, 2, 3), m  j (t) are the mass and position of the j th particle in phase p (the time-dependence of each of these is suppressed above), N p is the number of particles of phase p, φ is the spatial kernel (assumed to be symmetric and integrate to unity), and δ is a Dirac delta. For the specific choice of Dirac delta spatial kernel (as in this work), (2) becomes ( The positions of mobile particles evolve according to the stochastic Langevin equation, as derived from Fick's Law where ξ j is a d-dimensional (d = 1, 2, 3) random vector, and each component is an independent standard normally-distributed (N (0, 1)) random variable. We formulate the miRPT MT algorithm as a generalization of the work of Benson and Bolster [9] such that transfer of mobile mass from one phase to another may be formulated as a weighted sum of mass exchanges is the weighting function for the mass transfer from the j th particle in phase p to the k th particle in phase p and need not be the same for p = I or M . Since all mobile mass is transferred from one phase to the other by the miRPT algorithm, this implies that the weights sum to 1 over the p index, i.e., or, We see that the form of the sum in (7) implies that this may be written in matrix-vector form as where the entries of the matrix and vectors are [ The specific weighting function we use for miRPT is based on the probability density for the co-location of a mobile-immobile particle pair by diffusion within a timestep of length ∆t [14] v κ (s) = 1 where d = 1, 2, 3, is the number of spatial dimensions, D is the diffusion coefficient, |s| is the (Euclidean) distance between the pair of particles, and κ (p) is a scaling constant on the portion of diffusion to be simulated by a single mass transfer that may be different for p = M, I. Now, in order to enforce (6), we select w (p) to be The reason for this choice is that we would like the miRPT MT process to be diffusive, as was shown for a similar MT algorithm in [13], and division by the sum term in brackets ensures that the weights sum to unity and all mass is transferred from one phase to the other (computational details for this process are discussed in Appendix B.2.1). Choosing κ = 0.5 implies that half of this diffusive mixing should occur with each mass transfer "direction" (this choice is examined in Appendix A). Thus, the composition of MT mixing and macrodispersive random walks may additively simulate the entirety of the diffusion in the system. As a side note, the reader may notice that v 2 (s) would be the co-location probability density for a mobile-mobile collision. Additionally, if, in (8), we treat phase p as mobile particles at t = t 0 and p as mobile particles at t = t 0 +∆t (this requires that N p ≡ N p ), we may use this density to define the weight matrix where , diag(x) denotes a diagonal matrix with the entries of the vector x on the main diagonal, and 1 is a 1 × N vector of ones. This is the matrix associated with the "explicit matrix" method discussed in [13]. That is to say, the miRPT method recovers the classical mobile-mobile mass transfer algorithm.

Computational details
In this section, we discuss details that are not essential to the reactive transport portion of the algorithm but are nonetheless necessary from a computational/coding standpoint. The first concern to address is the (typically) computationally expensive nature of chemistry calculations. The miRPT algorithm has a natural ability to handle this, as all of the chemistry calculations are performed "on" the immobile particles, and the number of immobile particles (N I ) is a parameter that may be adjusted according to the desired balance between error and computational expense. As we show in Section 3.1.1, the error in the miRPT algorithm as a diffusive operator is at its lowest when the number of immobile particles is greater than or equal to the number of mobile particles (N I ≥ N M ); however, N I may be adjusted downward (e.g., N I = N M /2) in favor of conducting fewer chemistry calculations, and the aforementioned error may still be made relatively small for appropriate choices of other simulation parameters (∆t, N M ).
From a theoretical standpoint, according to the miRPT MT algorithm outlined in Section 2.2.2, every particle exchanges mass with every other particle, regardless of separation distance. However, because the co-location probability density given in (9) is Gaussian, we know that particles with a separation distance larger than three standard deviations σ = √ κ (p) 4D∆t have 0.3% probability of co-locating and exchanging mass. For this reason, a cutoff distance for MT interactions is typically imposed, so that particles with a separation distance greater than 3σ do not interact. If this cutoff distance approach is combined with a linear algebra package that takes advantage of sparse matrix structure, this can provide substantial computational speedup and memory savings. As well, imposing this cutoff distance lends the problem to highly efficient fixed-radius search algorithms O (N log N ) , as compared to O N 2 for a naïve search , where the standard choice is a kD Tree search algorithm [20] such as MATLAB's rangesearch() [21] or an open-source version available in C ++ or Fortran 95 [22]. Finally, all of the above computational strategies become essential for models with higher than one spatial dimension, as dense matrices of the required size often will not fit in memory, and a naïve distance search quickly becomes computationally intractable.

Derivation of diffusion operator
In order to analyze the accuracy of the miRPT MT algorithm of Section 3.1, we wish to compare it to another numerical method that solves the continuum diffusion equation given a particle-like spatial discretization. We outline the derivation of such a method here.
In an infinite d-dimensional domain, the solution to (12) at time T can be exactly expressed as a convolution of the initial condition (IC) with the associated Green's function. Given C 0 (x) := C(t = 0, x), and the Green's function for the diffusion equation where denotes convolution. If we discretize time, such that T = k∆t, we may equivalently find our solution in the following way (employing pseudo-code) If we wish to use this same method to compute an exact solution in the case where our initial condition is composed of N Dirac deltas, each with position x i and mass m i (an IC that corresponds to a particle-tracking simulation), we may insert the Dirac particle representation for C(t = 0, x) from (3) into (13) Note, however, that this is a continuous exact solution, given a particle-like IC. If we wish our final solution to be represented by N discrete particles, it must have the form where However, we see that this space-discretized diffusion operator, D T , does not conserve mass because each column is simply N function evaluations from a Gaussian distribution function. Thus, we must normalize the columns of D T to sum to 1, which may be equivalently viewed as constructing a discrete probability mass function from the continuous distribution defined by G. So, we define where, as previously used, 1 is a 1 × N vector of ones, and diag −1 (x) is a square matrix with the element-wise reciprocal of the vector x on its main diagonal. So, substituting D T , above, into (16), we have a representation for MT by a particle diffusion operator where the error in the approximation above is due to discretization and becomes equivalent to (15) as N → ∞. As discussed in [13], the simulated solutions generated by this discretized matrix diffusion operator (henceforth referred to as the diffusion operator ) are nearly error-free when there is no influence from boundary conditions. However, because the derivation assumes an infinite domain, error is introduced in a finite computational domain when particles near the boundaries have non-negligible mass. This is due to the fact that the normalization in (17) alters the "shape" of the Green's function near the boundary. Nonetheless, the diffusion operator serves as a good algorithmic benchmark for the miRPT MT algorithm if boundary effects are minimized.

Results
In this section, we show the results of testing the accuracy of the miRPT MT algorithm in solving the diffusion equation (Section 3.1). As well, we consider the results of using the full miRPT reactive transport algorithm to model a calcite/dolomite geochemical system in which CO 2 -saturated brine is injected into a porous domain with solids composed of calcite and quartz (Section 3.2). This problem is of interest because it involves both the dissolution of calcite and precipitation of dolomite [23], making it an ideal test problem for the miRPT algorithm.
Most numerical simulations in this section were conducted in MATLAB, using a MacBook Pro with a 2.9 GHz Intel Core i5 processor and 8 GB of RAM. The more computationally-intensive MT-analysis simulations (2D and 100-member ensemble) were also conducted in MATLAB but were conducted on a 16-core node with two 2.7 GHz Intel Xeon E5 processors, with 8 cores each, and 64 GB of total RAM. The full reactive transport simulation of the calcite/dolomite system was also run on the 16-core node, but code was written in Fortran 90 and employed the PhreeqcRM geochemical solver [16] and kD Tree code from [22]. Additionally, we make all code used in this section available at doi.org/10.5281/zenodo.1244245 [24].

Analysis of mobile-immobile mass-transfer algorithm
We wish to gain an understanding of how the miRPT MT algorithm behaves in isolation and whether it solves the diffusion equation with a sufficient level of accuracy. Therefore, we follow the approach of Schmidt et al. [13] and consider only the MT algorithm, without other transport (random walks and advection) or reaction (i.e., simulating only steps (c) and (e) in Figure 1). Thus, we simulate a transfer of mass from the mobile phase to the immobile phase and then back to the mobile phase. Numerically, using (8), this is represented by where To asses the performance of the operator over the time period T = k∆t, we may apply the operator, W I W M , k times to the initial mass vector as follows We first consider the miRPT MT algorithm with equally-spaced particle positions and employ initial conditions (ICs) that have known analytic solutions, so that we may compare the results to these analytic solutions. In addition, we compare these results to simulated solutions generated by the diffusion operator discussed in Section 2.3 that has been shown to yield near-exact solutions in the absence of influence from spatial boundaries [13]. The two ICs chosen for this purpose are a domain-centered Gaussian (to avoid finite-domain boundary effects) and a Heaviside function (to investigate both the infinite gradient at a discontinuous interface and performance when boundary effects come into play). We examine the sensitivity of the error between the analytic reference solution and the simulated solution generated by both miRPT MT and the diffusion operator. This is done by successively refining the resolution provided by the parameter under consideration (∆t, N M , N I ) and tracking the resulting error, with the desired effect being that finer discretization will provide more accurate results. Further, having found evidence of a stability condition for the miRPT MT numerical scheme, we investigate the influence of that condition. Also, we investigate whether the stability condition holds for a 2D simulation.
Lastly, we consider the miRPT MT algorithm's performance for randomly assigned particle positions. In the case of the Gaussian IC, we first conduct a 100-member ensemble run and average the results. We also conduct a single run and quantify the error by tracking the spatial variance in mass (given in Equation (24)) in order to verify that it grows at the proper rate. We also examine the variance growth for a non-analytic "noisy box" IC in the same manner.
We specify here that we include the diffusion operator approach as a comparison tool for the miRPT MT algorithm. As is seen in the following sections, the diffusion operator approach achieves consistently lower error for many of the Gaussian IC cases, which is unsurprising since the diffusion operator is essentially interpolating a Gaussian solution using N point-evaluations for each of N Gaussian basis functions, a task for which it is naturally suited. However, the diffusion operator is solving a different problem than the miRPT MT algorithm. The diffusion operator approach is simulating diffusion with a discretized Green's function solution, using only the mobile particles, i.e., (using (18)) whereas, the miRPT MT algorithm is simulating diffusion by performing the alternating mobile-immobile and immobile-mobile transfers (given in (20)) that are required for the more complicated reactive transport problem we are ultimately interested in simulating. Thus, converging to the low error of the diffusion operator approach for a Gaussian IC is exactly the desired outcome, and outperforming the diffusion operator approach for the Heaviside IC, the 2D case, and all simulations with randomly-spaced particles is highly encouraging. Unless specified otherwise, all simulations are conducted on a 1-dimensional domain, Ω = [0, 1], for 1 second of simulation time and a diffusion coefficient of D = 1.0 × 10 −3 m 2 /s. The parameter ∆t is chosen in order to properly explore the stability condition that is discussed below but typically varies between 1.0 × 10 −3 and 1.0 × 10 −1 . Further, the cutoff-distance approach described in Section 2.2.3 is not employed here (other than in the 2D case when memory becomes an issue), as computation time was not a concern for these simulations.
The measure of error used in these analyses will be root-mean-square error (RMSE), defined as follows for the RMSE between a given simulated solution C and reference solution C * (both vectors of length N )

Mass-transfer analysis for equally-spaced particles
We first look to demonstrate that the miRPT MT algorithm solves the diffusion equation in the most idealized case, wherein particles are assigned fixed, equally-spaced positions in order to remove the randomness introduced by the diffusive random walks. We stress here that using such a particle simulation to model the diffusion equation in practice would be non-optimal, as we have essentially recreated an Eulerian grid but are using an N -point stencil. In this case, a finite difference (FD) or finite element method with a much narrower stencil would likely be a more computationally-efficient choice. However, to re-iterate, we do wish to show that the miRPT MT algorithm is diffusive so that we may partition a system's total diffusion between random walks and this MT process with a controllable level of numerical error.
In this section, convergence of the miRPT MT-simulated solutions to the analytic solution is investigated as it relates to all discretization parameters (∆t, N M , N I ). To do this, the parameter of interest is refined by successive halves to increase the resolution of the simulation, and we seek to verify that error is decreasing with the increased resolution.
Relationship between error and ∆t. Figures 2 and 3 show the results of examining the sensitivity of error to the time discretization parameter, ∆t. Both figures demonstrate that the miRPT MT algorithm is relatively insensitive to changes in ∆t (as long as the stability condition discussed in a subsequent section is not violated). For these simulations, N I = 500, N M = 1000. We choose this ratio because, in a reactive transport model, it provides a desirable trade-off between minimizing the number of chemistry calculations to be performed (on the immobile particles) and minimizing the error due to the MT algorithm. The effects of altering this ratio are explored in a following section.
In Figure 2 (Gaussian IC), we consider the results when finite-boundary effects are minimized by employing the Gaussian IC and we see that error does, in fact, decrease with ∆t. As well, we see that the error in the miRPT MT algorithm appears to approach a minimal level of error that is achieved by the diffusion operator (whose error grows marginally as ∆t decreases). In Figure 3 (Heaviside IC), we see a different story. As discussed in [13], the diffusion operator performs poorly when particles near the boundary have non-negligible mass, and that is clearly displayed in Figure 3(a) where we see a vaguely Gaussianshaped "bump" near the left-hand boundary, when those particles should all have a mass of exactly 1. However, the miRPT MT algorithm attains an error many orders of magnitude smaller in this case, and hardly changes with refinements in ∆t.
Relationship between error and particle number. Figures 4 and 5 show the results of examining the sensitivity of error to the spatial discretization parameters, N M and N I , the number of mobile and immobile particles, respectively. For these simulations, ∆t = 1.0 × 10 −1 , N I = N M /2 is held at a constant ratio (the reason for this choice of ratio is discussed in the previous section), and N M is successively doubled. In Figure 4 (Gaussian IC) we see similar, but more dramatic, results as in the previous section. As the choice of N M drives the miRPT MT algorithm into the stable region of the parameter space (this stability condition is discussed in a subsequent section), we see the error in the miRPT solution decrease by approximately six orders of magnitude and then level off to the minimal value attained by the diffusion operator. In Figure 5 (Heaviside IC) we see that miRPT MT initially has larger error than the diffusion operator but, the error decreases steadily as particle number increases and becomes many orders of magnitude lower than that achieved by the diffusion operator.     Relationship between error and particle number ratio. Figure 6 shows the results of exploring the sensitivity of error to the ratio between the spatial discretization parameters, N M and N I , the number of mobile and immobile particles, respectively. For these simulations, only a Gaussian IC is considered, ∆t = 1.0 × 10 −3 , N M = 1000 is held constant, and N I is successively doubled. Considering Figure 6, we see the error in the miRPT MT algorithm decreasing steadily as N I increases and the ratio N I /N M approaches and exceeds 1. We also see that the error approaches the minimal level attained by the diffusion operator, which is constant because there are no immobile particles in the diffusion operator approach. One might have expected the minimal level of error to be attained when N I = N M = 1000, but for the chosen parameters this is not the case.
To explore this, we added an extra data point (N I = 1200) to our doubling values of N I , and we see in Figure 6 that the minimal level of error is achieved relatively soon after we move into the stable region of the parameter space and the ratio N I /N M exceeds 1.
Stability condition. We see in Figures 2(b), 4(b) and 6 convincing evidence of something like a stability condition at work in the miRPT MT algorithm. We formulate this stability condition (in 1D) in the form of an empirical bound on the dimensionless number where L is the length of the domain, Ω. The form of η is reminiscent of a corresponding von Neumann stability condition, as applied to FD schemes (in fact, it is the reciprocal), where the term L/ min(N I , N M ) takes the place of ∆x and can be thought of as the minimum spatial resolution provided by the miRPT MT algorithm (or, alternatively, as the maximum average inter-particle spacing). We acknowledge that the stability condition proposed here does not operate in exactly the same manner as a von Neumann stability condition for FD schemes. However, there is numerical intuition behind why this MT stability condition operates, and it lies in the form of the weighting function (9) used for the miRPT mass transfers (i.e., the co-location probability density for a mobileimmobile particle pair). The argument of the exponential in (9) contains the quantity −s 2 /(D∆t), where s is the distance between a given particle pair. So, the ratio of squared average inter-particle spacing L/ min(N I , N M ) 2 must stay proportional to D∆t. Otherwise, as this ratio (η) becomes much larger than 1, the magnitude of mass transfers (or probability of particle interactions) rapidly becomes, on average, indistinguishable from 0, to machine precision. We see, in all of the aforementioned figures, a distinct drop in the error (sometimes several orders of magnitude) when η crosses below the threshold value of 1, where, afterward, the error appears to converge to a minimal level. As well, we note here that for the results depicted in Figure 3(b), there is in fact a distinct, relative jump in error for the final data point, where η > 1; however, the scaling of the figure does not make this apparent. The results shown in Figure 6 inform us that there may be something slightly more complex occurring in regards to the stability condition. In Figure 6, we do see a marked drop in error (roughly 4 orders of magnitude) between data points 4 and 5 (N I = 500, 1000, respectively), when the stable region of parameter space is entered. However, the minimal level of error is not attained until the value of N I is increased to 1200. The form of η also provides insight as to why the error in the miRPT MT algorithm is relatively insensitive to changes in ∆t yet highly sensitive to varying particle number (whether mobile or immobile) because the stability condition varies linearly with ∆t and quadratically with the minimum particle discretization.

Stability in 2D
We see the results of the stability analysis for a 2D simulation depicted in Figure 7, where the form of the stability condition becomes where L 1 and L 2 are the lengths of the sides of the rectangular 2D domain, Ω. The changes in the numerator of (23), as compared to (22), serve to, again, allow this quantity to represent the minimum spatial resolution provided by the simulation parameters (corresponding to ∆x in an Eulerian framework). For the results shown in Figure 7, L 1 = L 2 = 1, , ∆t = 1.0 × 10 −1 , N I = N M /2, and √ N M was successively doubled. We see in Figure 7, that for the properly formulated stability condition, there is the anticipated drop in error (≈ 5 orders of magnitude) when crossing from the unstable to stable region of parameter space. However, since the form of η is changed for the 2D case, we no longer N M have a quadratic decrease in η with particle number but only a linear relation. As a result, driving the algorithm into the stable region requires much larger increases in particle number than it does in the 1D case.

Mass-transfer analysis for randomly-spaced particles
Having verified that the miRPT MT algorithm solves the diffusion equation satisfactorily in the idealized, equally-spaced particle case, we must also consider the case that corresponds to mobile particles moving via random walk. For this analysis, immobile particles are still spaced evenly across the domain, but the mobile particles are assigned positions according to independent draws from a Uniform, U(0, 1), distribution to simulate their random motion. We perform a convergence analysis for the randomly-spaced particles in two cases. First, we conduct an ensemble run of 100 simulations for each set of parameters and average the results to see whether the expectation of these stochastically perturbed simulations converges to the expected value of the analytic solution ( Figure 8). Next, we examine the results of a single randomly-spaced particle simulation from a statistical standpoint to confirm that the spatial variance in mass grows at the expected rate of 2Dt. This second analysis is conducted both for a Gaussian IC (Figure 9) with an analytic solution that can be visually verified and also for a non-analytic "noisy box" IC ( Figures 10 and 11).
Ensemble run for Gaussian initial condition. The results of the 100-member ensemble run for a Gaussian IC are shown in Figure 8. For this simulation, ∆t = 1.0 × 10 −3 , N I = N M , and N M was successively doubled. The results of this ensemble run were then averaged by first placing particles into 250 spatial "bins" and averaging the masses in each bin across the ensemble to generate an ensemble-averaged concentration for each bin. This binning was necessary because, given randomly assigned particle positions from a continuous probability distribution, a certain particle will never be in the same position for any two simulations, so a particle-wise average does not make sense.
While strict convergence is difficult to prove for a stochastic algorithm, we do see favorable results in Figure 8. First, in Figure 8(a), for the final refinement in particle number (N I = N M = 12800), we see good visual agreement between the analytic solution and both the miRPT MT algorithm and the diffusion operator approach. However, as would be expected, the error shown in Figure  8(b) does not drop nearly as quickly with increases in particle number as it did for the equally-spaced case. Additionally, we see miPRT MT outperforming the diffusion operator approach for all tested parameters. Further, there is evidence to suggest that the stability condition may become more restrictive with randomly-spaced particles. In Figure 8(b), the largest drop in error is seen between points 4 and 5 (N M = 1600 and 3200, respectively), both of which are associated with values of η < 1, suggesting that η < 0.25 may be better practice when particle positions are not equally-spaced. This follows logically, since simulations with randomly-spaced particles will inevitably include some larger gaps between particles (physically, representing a poorly mixed area of the domain). As a result, the quantity L/ min(N I , N M ) (representing maximum average particle spacing) must be made smaller than in the equally-spaced case.
Variance analysis. The results of variance growth analysis for a single run with randomly-spaced mobile particles are shown in Figures 9-11. The spatial variance in mass is calculated as where m (I) and m is the first spatial moment, or center of mass. Here we note that the center of mass should be 0.5 for all results in this section; however, we do not plot results for mean position, as spatial mean was nearly exact for the miRPT MT algorithm and gained a minimal amount of error with time for the diffusion operator approach. These results may be visually inferred by the shape of the final mass distributions in Figures 9(a)-11(a). For ease in visual comparison, we depict spatial variance growth in Figures 9-11 that is computed by subtracting the spatial variance of the IC from σ 2 S . In addition, the single-realization results, as shown in this section, will vary based on initial mobile particle positions, but the results depicted in the figures are typical for the given parameters and ICs.
For these simulations, N I = N M = 1600, ∆t = 1.0 × 10 −1 for Figures 9 and 10, and ∆t = 1.0 × 10 −2 for Figure 11. In Figure 9(b), we see that for the Gaussian IC, the miRPT MT algorithm is virtually error-free in terms of variance growth and that a small amount of error is introduced by the diffusion operator approach. In Figures 10 and 11 we consider a non-analytic "noisy box" IC, wherein the middle 20% of particles are assigned their initial mass according  to independent draws from a U(0, 1) distribution. This is done so as to rule out finite boundary effects while eliminating the smooth gradients provided by the Gaussian IC. Considering Figure 10, we see that both miRPT MT and the diffusion operator initially introduce extra spatial variance; although, since this error stays relatively constant over the 1 second simulation, the error actually decreases, in a relative sense, with time. A somewhat different result can be seen in Figure 11, where the diffusion operator at first increases the spatial variance by nearly and order of magnitude, and miRPT MT to a lesser degree. Both algorithms, however, closely approach the analytic variance after a certain point in time, with miRPT MT becoming nearly error-free by the end of the simulation. This notion of "accuracy as a function of repeated operation" is discussed in [13].

Reactive transport model
Having established that the miRPT MT algorithm is capable of solving the diffusion equation with a sufficient level of accuracy, we now wish to apply the full miRPT algorithm (all of the steps depicted in Figure 1) to a chemicallycomplex reactive transport problem. The test problem we choose is based on one considered by Leal et al. [23]. In this system, which we refer to as the calcite-dolomite reactive transport (CDRT) system, we consider a high temperature and pressure (60°C, 100 bar) 1-dimensional domain (0.5 m in length) composed of calcite (CaCO 3 ) and quartz (SiO 2 ) (2% and 98%, respectively, at uniform 50% porosity) into which a constant concentration of CO 2 -saturated brine composed of various salts (≈ 0.88 mol/L NaCl, 0.05 mol/L MgCl, 0.01 mol/L CaCl 2 ) is injected at x = 0. Injection of an acidic brine with Mg concentrations higher than in the initial domain causes the fluid to become undersaturated with respect to calcite and oversaturated with respect to dolomite, resulting in the co-located dissolution of calcite and precipitation of dolomite as a reaction front that advances through the domain in the direction of advection in time. Calcite dissolution buffers the pH at the leading edge of this reaction front, while dolomite dissolution buffers the pH within the reaction front. Once both calcite and dolomite fully dissolve at any given location, the pH is is no longer buffered and is equal to the pH of the inlet brine (≈ 3.0), defining the trailing edge of the reaction front. For a comparable schematic diagram of this IC and boundary condition (BC), see Figure 1 in [23]. This system is governed by the ADRE given in (1), where v = 2.4 × 10 −5 m/s and D = 1.2 × 10 −7 m 2 /s. Reactions are considered to be instantaneous, and these calculations are handled by an operator-splitting approach using the PhreeqcRM geochemical solver [16] and the phreeqc.dat thermodynamic database [25].
This CDRT system is an ideal test bed for miRPT for a few reasons. First, calcite and dolomite are solid phases that are not transported but must interact with the flowing CO 2 and salts in the injection brine that are being transported by advection and diffusion. Thus, calcite and dolomite concentrations will always be stored on immobile particles. The masses of the aqueous ions composing the injection brine (10 total species) will be stored on mobile    particles for transport (Steps (a)-(b) in Figure 1) and then transferred to the immobile particles using the miRPT MT algorithm (Step (c) in Figure 1). Once there, reaction mass-balance calculations are conducted before the resulting aqueous ions are transferred back, again according to miRPT MT (Steps (d)-(h) in Figure 1). Another reason this problem was chosen is that the low level of diffusion O 10 −7 m 2 /s will present a challenge for miRPT MT, due to the stability condition given in (22). As discussed in Section 3.1.2, when particles are randomly-spaced (as will be induced by the diffusive random walks in this model), the stability condition becomes more restrictive than imposing η < 1, possibly requiring η < 0.25. For that reason, we devote Appendix B.1 to choosing the proper model parameters to ensure a stable simulation. We see the results for the miRPT model of the CDRT system in Figure 12. The miRPT results are compared to an analogous FD simulation that uses operator splitting to solve the ADRE (1) using explicit upwinding for advection, explicit centered differencing for diffusion, and PhreeqcRM for reaction calculations. As in [23], the spatial distributions of pH and concentrations of calcite and dolomite are recorded and analyzed to ensure both the FD and miRPT solutions are displaying proper behavior. We see near exact fit between the two model solutions for the times depicted in Figure 12 (t = 50, 250, 1000 min), with some slight discrepancy in the position of the dissolution and precipitation fronts for calcite and dolomite at t = 1000. This, however, is not particularly concerning because the difference is small, and the "true" solution is often unclear for systems of this type. This is because various modeling decisions may generate solution variations on the same order of magnitude, for example, the choice of thermodynamic database, the level of discretization, or the order in which the transport and reaction operations are performed [23,26].

Summary and conclusions
In this paper, we present a Lagrangian (particle-tracking) method for modeling reactive transport with interactions between solid and aqueous chemical phases using mobile and immobile particles. These particles transfer mass between phases (mobile to immobile, or vice versa) using an algorithm that is demonstrated to solve the diffusion equation with a controllable level of error, meaning that a high degree of accuracy may be attained for appropriate parameter choices. This mobile-immobile reactive particle-tracking mass-transfer (miRPT MT) algorithm is physically-motivated and based on collision probabilities between particle pairs. Mathematically, it has the form of a finitedimensional convolution with a Green's function, which is the basis for the diffusion operator approach, derived in Section 2.3. However, the key difference is the double application of this convolution operator that is applied to two separate "particle grids" in the miRPT MT algorithm. The regularly-spaced grid represented by the immobile particles allows miRPT MT to overcome the sometimes irregular spacing of the mobile particles (which can represent heterogeneity or reactant segregation), a problem which degraded the accuracy of the diffusion operator approach. Additionally, we show that miRPT produces the proper concentration spatial variance growth in simulations with random mobile particle spacings, when strict convergence to an analytic solution would not be expected or an analytic solution may not be available. This is important because, in real-world systems, global convergence is not a realistic metric, but as long as the proper dynamics are being simulated at the local (particle) level, the algorithm remains valid. We know that the miRPT MT algorithm is simulating local diffusion dynamics because each mass transfer is a discretized convolution of the particle's spatial kernel (Dirac delta, for this work) with the diffusion Green's function, and, in the infinite particle limit, this solution is exact.
To summarize, we present the following conclusions: • The mobile-immobile reactive particle-tracking (miRPT) algorithm extends the work of Benson and Bolster [9] that enables particle-tracking methods to simulate arbitrarily complex chemical reactions by allowing particles to carry many species of reactant and transfer mass between particles.
• The miRPT mass-transfer (MT) algorithm allows for mass transfers between aqueous and solid (mobile and immobile) phases, enabling it to model chemical reactions such as dissolution and precipitation.
• The miRPT MT algorithm is demonstrated to solve the diffusion equation with a controllable level of error, depending on the level of discretization.
-For this reason, the total diffusion of a system may be partitioned between diffusive random walks and miRPT MT, allowing for the necessary mass transfers to occur without introducing spurious transport effects. -Because the miRPT MT algorithm takes the form of a convolution with the diffusion Green's function, it can simulate any length of time step (∆t), as long as attention is given to the stability condition, discussed throughout Section 3. -The miRPT MT algorithm displays rapid convergence as the number of mobile and/or immobile particles (N M , N I ) in a simulation is increased, and the solution approaches the exact solution in the limiting case.
• In the context of a full reactive transport simulation (calcite/dolomite dissolution/precipitation is considered in this work), the miRPT algorithm is demonstrated to perform effectively and generate similar solutions to an analogous Eulerian finite difference model.

Acknowledgments
Appendix A. Relationship between error and κ (p) In this section, we consider the effect of the scaling constant, κ (p) , p = M, I in (9), on the behavior of the miRPT MT algorithm. For these simulations, ∆t = 1.0 × 10 −1 , N M = 1000, the total diffusion in the simulation is partitioned between κ (M ) and κ (I) such that κ (M ) = 1 − κ (I) , and a Gaussian IC is employed. We consider the parameter range κ (p) ∈ {0.1, 0.2, . . . , 0.9} and test these parameters both for equally-spaced (as in Section 3.1.1) and randomlyspaced (as in Section 3.1.2) mobile particles. We conduct these tests for the range of values N I ∈ {100, 200, 500, 1000}.
For these results, we consider a generalization of the stability condition given by (22)  , we see that the error is mostly insensitive to the values of κ (p) , until the arithmetic stability condition is violated (η > 1), and we see a jump in error of approximately two orders of magnitude. Figure A.13(d) shows similar results with one key difference. We see the characteristic large jump in error when η > 1; however, we also see a sizable jump in error on the low side, before we attain the minimal level of error in the range κ (M ) ∈ {0.3, 0.4}. This suggests that there is something more nuanced occurring and that perhaps the stability condition formulated above is not quite restrictive enough when the ratio N I /N M becomes small. We propose one possibility that involves introducing the quantity which is the geometric mean of η κ (M ) and η κ (I) (denoted geometric stability condition (SC) in Figures A.13-A.14). The geometric mean is a suitable choice for averaging these two quantities because it does not scale linearly, like the arithmetic mean, and so is able to capture the average of quantities that may have a significant difference in scale. Using this, the stability region appears to be of the form {η 1} ∩ { η 0.25}, as there does not appear to be a strict cutoff when considering η and η together. In Figure A.14, we show the results for the randomly-spaced mobile particle test. For these simulations, mobile particle positions were assigned according to independent draws from a U(0, 1) distribution, as in Section 3.1.2. In this analysis, the magnitude of the errors is relatively large, as no ensemble-averaging (as in Section 3.1.2) is conducted; however, the trend is clear in that there are obvious minima at κ (M ) = κ (I) = 0.5 for both the N I = N M = 1000 and N I = N M /10 = 100 cases (Figures A.14(a) and (b), respectively).
Additionally, as in the ensemble-averaged results given in Section 3.1.2 the stability condition has a different behavior for randomly-spaced particles than it does for equally-spaced. The important result, though, is that the minimum error reliably occurs at the point where η is also minimized, and this result holds reliably for all tested ratios of N I /N M between 0.1 and 1.0. For this reason, the parameters κ (M ) = κ (I) = 0.5 will be held constant for all tests conducted in Section 3.

Appendix B. Supporting information for reactive transport system
Appendix B.1. Selection of model parameters As a starting point, we first ensure that the FD model replicates the results of the CDRT system presented in [23]. This led to choosing a time step of ∆t = 1.0 seconds and a spatial discretization of ∆x = 1.0 × 10 −3 m, which results in a Courant number of C F D := (v∆t)/∆x = 0.024 and a grid Péclet number of P e F D := (v∆x)/D = 0.2. We note here that the choice of ∆t requires consideration also be given to the timescale of the chemical reactions being modeled, as the equilibrium reactions are considered to be instantaneous. However, that is not the primary focus of this text, so we merely verify that the FD model reproduces the original results of Leal et al. and proceeded to use that value of ∆t in our miRPT model as well.
The next step is to choose the proper level of spatial discretization (i.e., the number of mobile and immobile particles, N M and N I ), given the the stability constraint provided by the miRPT MT algorithm and the defined parameters of the CDRT system, D = 1.2 × 10 −7 , Ω = 0.5, along with the choice of ∆t = 1.0. For our miRPT model of the CDRT system, we choose to simulate half of the total diffusion by random walks and half using the miRPT MT algorithm (for this reason we impose D M T = 0.6×10 −7 for the following miRPT MT analysis). This 50/50 partitioning of diffusion is somewhat arbitrary in this case, though the analysis would be similar regardless of this modeling choice because the value of the stability parameter, η, stays on the same order of magnitude for any partition that sums to unity. A discussion of employing this partition between random walk and MT diffusion in order to differentiate between macro-scale dispersion and micro-scale mixing may be found in [13].
The first step we take is to verify that the miRPT MT algorithm is capable of simulating the low level of diffusion in the CDRT for stationary mobile particles   in both the equally-and randomly-spaced cases. The results of this analysis are not shown here, as they are positive and mirror those seen in Section 3.1. Next, we consider spatial variance growth in the case where the mobile particles are simulating diffusion via random walks. These conditions capture the total action (in terms of variance growth) of the conservative transport portion of the CDRT system, as advection simulation should not induce variance growth. This is, of course, ideally speaking, as FD schemes for advection simulation are known to induce numerical diffusion, a feature which is not present in particle-tracking methods. The results shown in Figure  . The maximum and final error in both cases is nearly two orders of magnitude smaller than the value of the variance, itself. This reinforces the previously mentioned assertion that the random walk diffusion for mobile particles works to smooth the effects of any irregular spatial position distributions. An important result here is that a larger number of particles (both mobile and immobile) is required to adequately model the desired level of diffusion than in the cases without random walks (satisfactory results were achieved using half as many particles for some of the stationary particle cases), indicating that a good rule of thumb may be to keep η < 0.5. Additionally, for both N I = 4000, N M = 5000 and N I = 3000, N M = 6000 there is a significant increase in error if N I is reduced much below these values or if N M is reduced much below 5000. From an accuracy perspective, either of the parameter pair choices from the random walk/miRPT MT simulations, shown in Figure B.15, would be appropriate. As well, from a computational perspective, if we are considering a conservative system (i.e., without chemical reactions), the choice of N I = 4000, N M = 5000 may be superior. This is because, while both choices result in 9000 total particles (leading to approximately the same number of total masstransfer calculations), choosing to employ fewer mobile particles will lead to fewer transport calculations and a faster simulation. However, the primary computational cost, by far, in the CDRT system is due to the chemical reaction calculations that are performed by PhreeqcRM. So, we would like to minimize the number of these calculations that are performed by choosing the parameter pair with the minimum appropriate number of immobile particles. For that reason, the miRPT CDRT model will use N I = 3000, N M = 6000.

Appendix B.2. Computational details for the reactive transport system
In the preceding sections, we have laid the theoretical groundwork necessary for developing a reactive transport model using the miRPT algorithm. However, from an application standpoint, actually writing the code for such a model involves some computational choices that may not be readily apparent. For this reason, we intend this section to act as something of a high-level "user's guide" for the miRPT model. Additionally, we make our CDRT FD and miRPT code available at doi.org/10.5281/zenodo.1244245 [24].  We begin the CDRT simulation by initializing particle positions. Immobile particles are equally-spaced on the interval Ω = [0.0, 0.5] (this may be thought of as creating an Eulerian grid with ∆x = 1/N I ), and mobile particles are assigned positions according to draws from a U(0.0, 0.5) distribution. Once the geochemical solver, PhreeqcRM, is initialized, we pass our IC (calcite and quartz within the domain) and BC concentrations (CO 2 and salts to be injected) to it, so that we may assign initial concentrations of our solid phases and aqueous ions to the appropriate particles. All of the immobile particles are assigned the concentrations generated from the IC. The BC is slightly more complicated in the particle context. In the FD context, the "constant concentration injection" BC may be interpreted as a constant-flux, or Neumann BC. To match this type of BC, concentrations are assigned to the lower-boundary immobile particle (x = 0.0). Once time-stepping and transport begin, negative fluxes at the lower boundary are prevented by imposing a "reflecting" BC on the diffusive random walks of the mobile particles (advection only moves in the positive direction for this problem). So, a particle that randomly walks across the lower boundary and attains a negative x-position is re-assigned the absolute value of that x-position.
With the simulation initialized, we convert the concentrations (mol/L) of the aqueous species on the immobile particles (the species undergoing transport) to masses (mol/particle) by multiplying by the representative volume of an immobile particle (V 0 = (0.5 m × 1000 L/m) /N I ) (step (f) in Figure 1). These masses are then transferred to the mobile particles according to the miRPT MT algorithm (step (g) in Figure 1, details described in the following section). Here, time-stepping begins, and random-walk diffusion and advection are simulated sequentially by altering mobile particle positions according to (4) (step (b) in Figure 1). A minor computational note regarding the random-walk algorithm is that, as this code was written in Fortran, there is not a built-in normallydistributed random number generator, so we use the Box-Muller transformation [27] to generate N (0, 1) random numbers using U(0, 1) random numbers (for example code, see [24,28]). Next, BCs are enforced according to the previously described reflecting lower boundary, as well as an "absorbing" upper boundary (x = 0.5) that can be interpreted as a homogeneous Dirichlet BC in the FD context. Computationally, this absorbing boundary is enforced by reassigning any mobile particle that crosses the upper boundary to have position x = 0.0 and carry no mass. This "wraparound" reassignment to the lower boundary allows the simulation to maintain a constant particle number. Next, the masses of the aqueous species on the mobile particles are transferred to the immobile particles via the miRPT MT algorithm and converted to concentrations (steps (c) and (d) in Figure 1). These concentrations are passed to PhreeqcRM for the chemical reaction calculations (step (e) in Figure 1). After reaction, we apply the injection BC again, and then proceed to steps (f) and (g) in Figure  1 and repeat the steps in Figure 1 until we reach the final simulation time. We note here that the BC injection scheme we employ allows for a constant flux of aqueous species mass through the boundary even if no mobile particles leave the upper boundary (and then enter the lower boundary) because the BC is injected into the lower-boundary immobile particle in every time step. That is to say, assigning the BC mass to an immobile particle ensures the proper mass flux in each time step, as injecting the mass on mobile particles (that may or may not enter the lower boundary in a given time step) would introduce an undesirable random variability to the BC.
Appendix B.2.1. Computational details for the mass-transfer process In order to conduct the miRPT MT calculation, two pieces of information are required: the positions and masses of aqueous species on all mobile and immobile particles (x (p) and m (p) , respectively, for p = M, I). Using the position vectors, x (M ) and x (I) , we construct a pairwise distance matrix, ∆, where . These entries of ∆ are calculated using a kD tree algorithm [20,21,22] to perform a fixed-radius search and find the nearest neighbors within a specified cutoff distance (this is discussed in further detail in Section 2.2.3). It is worth noting that ∆ may be constructed in sparse fashion, which can often lead to significant computational speedup and memory savings if sparse linear algebra methods are employed. We then use ∆ to construct the relevant MT matrix (W p , p = M, I, depending on the direction of transfer) in the following way, using the weighting function given in (9)  where exponentiation is considered to be an element-wise operation, and D M T is the diffusion constant corresponding to the portion of diffusion to be simulated by miRPT MT. Also, 1 is a 1 × N vector of ones, and diag −1 (x) is a square matrix with the element-wise reciprocal of the vector x on its main diagonal, as employed previously. We also note that (B.1) and (B.3) do not include multiplication by the constant term in (9) because it would divide out due to the column-normalization step in lines (B.2) and (B.4). Finally, the mass-transfers are computed according to (8), where for a immobile to mobile transfer, the formula would be W I m (I) = m (M ) .