{'fraction_non_alphanumeric': 0.039048926346506425, 'fraction_numerical': 0.012920318304424201, 'mean_word_length': 4.905112799002381, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 4, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Reactive Transport: A Review of Basic Concepts with Emphasis on Biochemical Processes

Reactive transport (RT) couples bio-geo-chemical reactions and transport. RT is important to understand numerous scientific questions and solve some engineering problems. RT is highly multidisciplinary, which hinders the development of a body of knowledge shared by RT modelers and developers. The goal of this paper is to review the basic conceptual issues shared by all RT problems, so as to facilitate advance along the current frontier: biochemical reactions. To this end, we review the basic equations to point that chemical systems are controlled by the set of equilibrium reactions, which are easy to model, but whose rate is controlled by mixing. Since mixing is not properly represented by the standard advection-dispersion equation (ADE), we conclude that this equation is poor for RT. This leads us to review alternative transport formulations, and the methods to solve RT problems using both the ADE and alternative equations. Since equilibrium is easy, difficulties arise for kinetic reactions, which is especially true for biochemistry, where numerous frontiers are open (how to represent microbial communities, impact of genomics, effect of biofilms on flow and transport, etc.). We conclude with the basic 10 issues that we consider fundamental for any conceptually sound RT effort.


Introduction
The term reactive transport (RT) describes the simultaneous occurrence of bio-geo-chemical reactions and transport of biochemical species (Error! Reference source not found.). RT is important to understand numerous large number of Earth Sciences problems, ranging from engineering applications (geologic carbon storage, enhanced geothermal systems, or groundwater remediation, etc.), as well as geological studies (like sediment diagenesis, seismic activities, hydrothermal ore deposition, etc.) to watershed or global issues [1].
The importance we attribute to RT contrasts with its impact in the literature. It is frequent to start a review with a bibliometric analysis to highlight the importance of the field. We searched for "Reactive transport" in the Science Citation Index Expanded (SCIE) bibliographic database. Results are quite sobering (Error! Reference source not found.): only 4880 papers, mostly coming from Environmental and Earth Sciences (very few from Health Sciences). Worse, the most cited papers are not "strictly" RT papers. The first and third most cited papers are quite general papers in Soil Science [3] and Seawater Intrusion [4]. The third one refers to non-Fickian transport [5]. The fourth one describes a specific surface reaction [6]. The four are relevant but not really RT papers by themselves. It is clear that this kind of search does not portray an accurate picture of the field. In fact, the list does not include the early works of Rubin and James [7], Grove and Wood [8] or Lichtner [9], who set the basis for the current paradigm of RT. But the list is more complete for recent years, which may reflect that the reactive transport term only became common in the 1990s. Since then, the growth has been steady, albeit only linear, which is modest compared to the standard exponential growth of sister fields (e.g. Cao et al., [10]). A second issue that emerges from the search is the large number of review papers (174, or 3.6% of the RT papers identified in the SCIE). We attribute the large number of reviews not so much to the difficulty of RT as to its diversity. RT reviews are driven by the nature of the chemical system, by the type of medium, by the types of reactions, or by the solution approach, etc. Several reviews emphasize geochemical issues [11]. Some are centered on problems related to CO storage [12,13], which generates acidity, so it is somewhat related to acid mine drainage [14], more formally revised by Amos et al. [15]. Nuclear waste disposal has motivated reviews on Uranium mobility [16], validity of existing models [17], cement-clay interaction [18] or oxygen transport at depth, which motivated the thorough review of MacQuarrie and Mayer [19] on RT in fractured rock. Global issues motivate the reviews of Dignac et al. [20], Li et al. [21], Arndt at al. [22] and Regnier et al. [23] concentrate on degradation of organic carbon in marine sediments. Wetlands are reviewed by Langergraber [24]. Energy related subsurface problems are starting to consider RT [25][26][27][28].
Numerous recent works address the multiscale nature of the RT, especially in evolving porous media [29][30][31][32][33], pore scale and pore network models [34,35]. Finally, many papers review and/or compare methods and codes [36][37][38]. We feel that the modest publication record and parallel growth on reviews may reflect a limited interest. That is, scientists are interested, but funding agencies not so much (you may reach similar conclusions at a far lower effort). Beyond difficulty (either in understanding and/or publishing), multidisciplinarity is another difficulty with RT, which involves specialists in transport and in reactions. The former typically involves quantitatively oriented hydrologists, engineers and physicists. The latter might sound like referring to chemists, but it refers to a broad and heterogeneous group which includes geochemists (originally geologists), chemists, molecular biologists and many others. They are all needed but their languages are so different that they need to work together. Multidisciplinarity is often hallowed, but interdisciplinary work is difficult, as it requires a determined effort in developing shared language and understanding. Worse, in practice it is penalized by evaluation panels.
The goal of this paper is to frame and review basic conceptual issues, which we highlight, with special emphasis on the fast-growing biochemical field. Most of the conceptual issues are widely accepted but might go unnoticed to newcomers into the RT field. Some of them are sufficiently old to be treated in textbooks (e.g., geochemical zonation), but some are not (e.g., importance and impacts of mixing, or geochemical localization). We will not enter into the fascinating world of parameter estimation [39-41] nor on the numerous application of RT, which are discussed in the field specific reviews mentioned above, except for illustrating conceptual issues.
We start with the basic equations for transport (Advection Dispersion Equation, or ADE) and reactions (Section 0). These allow us to frame simple reactive transport paradigms (Section 0) that illustrate what we consider the main conceptual issues of RT. In turn, the importance of properly representing mixing motivates Section 0, where we analyze transport equations other than the RT. Biochemical concepts are reviewed in Section 0 and solution methods in Section 0.

Basic equations
The first conceptual issue is that RT is multidisciplinary and requires interdisciplinarity. Thus, we start by reviewing the basic equations for transport and reactions.

Transport.
Fluid flow needs to be solved prior to transport, in order to know the time and space variability of fluid fluxes. Solving for flow requires establishing conservation of fluid mass and momentum. Fluid mass conservation is expressed as: [ / ], which is equal to the resident fluid density for sinks (i.e., when < 0). Note that, given the possibility of confusion, we find it useful to express units using not only the first name (e.g., L, T or M), but also the "last name" (i.e., subindex for water, actually fluid in general, or for medium). Mass conservation needs to be coupled to momentum conservation, given by the Navier-Stokes equation. The latter can be simplified to Stokes equation when flow is laminar, which is the usual assumption in pore-scale analyses (e.g., [34]), or further simplified to Darcy's Law when the porous medium is treated as a continuum [42].
Transport is usually solved with the Advection-Dispersion Equation (ADE), which expresses the mass balance of each aqueous species when transport is driven by advection and dispersion. Explicit ] includes the contributions of chemical reactions to the mass balance of species , expressed per unit volume of porous medium. This equation is written in various forms depending on the context. For example, it is frequent to write in terms of velocities instead of fluxes for historical reasons, but it is inconvenient in porous media. Also, we have written concentration as mass fraction (e.g., ppm) to ensure that dispersion does not cause advection and that adding Equation (2)(2) for all aqueous species yields the fluid mass balance Equation (1) (e.g., [43]). However, reactions usually express the solute mass [ ] in moles, so that can express moles of species per unit mass of solution (still a mass fraction, for practical purposes). When fluid density is approximately constant then it is valid to write transport in terms of (molarity). But, as we shall see, chemical data bases, usually require working with molality (moles solute per kg of , or = / , where is the mass fraction of ). The conservative form of the ADE (Equations (2) and (3)) is robust, but lacks an explicit representation of velocity, which is needed for Lagrangian methods. We can deduce a form, in which velocity is explicit, by subtracting Equation (1), multiplied by , from Equations (2)  ], which is why the most frequent option is to assume water incompressible and work with molarities (i.e., moles per unit volume of fluid). A somewhat easier to understand form results from dividing Equation Error! Reference source not found. by , while neglecting its spatial variations for the dispersive term, which yields , = ∇ ⋅ · ∇ , − ∇ , + , − , + , = 1, ⋯ , 5) where , = / is the contribution of reactions, but expressed per unit mass of fluid, is similar to , but expressed in terms of velocity instead of fluid fluxes, and = ⁄ is the mean water velocity. This form is preferred by many because it is identical to the transport equation used in the much broader community of fluid continuum transport, so that it facilitates recycling of codes developed in sister sciences. It is also the form preferred by the pore scale simulation community. Furthermore, it leads naturally to a Lagrangian formulations of transport (e.g., random walk methods) , = ∇ ⋅ · ∇ , + , − , + , = 1, ⋯ , 6) where ⁄ = ⁄ + ∇ is the material derivative. It is important to point that all these advective forms (Equations Error! Reference source not found., (5) and (6)) are equivalent to Equations (2) and (3). We have spelled them out to illustrate that they can be found in numerous forms. In practice, however, all of them require being extremely careful in the application of numerical methods. For example, very few authors use the reverse time weighting [45] for transient flow cases, which may lead to severe mass balance errors. Still, for simplicity, we will write Note that, written like this, the concentration of immobile species is expressed per unit volume of liquid, which is unusual. Normally, immobile species concentrations are expressed per unit mass of solid (or volume of porous medium). We will keep it as written for the sake of simplicity.
Some insight into the nature of this equation can be obtained by writing it in dimensionless form and from Error! Reference source not found.. The figure shows that the basic behavior of the ADE is defined by the dispersion length = √ , and the advection length = ⁄ = . Error! Reference source not found. displays the time evolution of these two lengths for different values of , , and . It is clear that, for short times and distances, diffusion (or dispersion) dominates, but advection always dominates at long time (and space scales). The characteristic time of the ADE is chosen as the one that equates these two lengths (i.e. = ⁄ ), which leads to a characteristic length = ⁄ and the dimensionless ADE is a characteristic concentration, chosen on the basis of the chemical system, often a boundary or initial concentration, here = ), = / 5 of 46 (i.e., the fraction of water flux, , that enters over ), and = / . This dimensionless form is unusual. The most frequent choice to select as a length of interest and then define the Peclet number as = ⁄ , which appears in front of the advection term. This way, > 1 is the same as saying > . We prefer using Equation Error! Reference source not found. because it makes apparent that the Peclet number is not an attribute of the ADE, but of the length of interest, which is relevant because any problem may have numerous spatial lengths of interest. For example, under normal conditions the ADE will be diffusion dominated at the pore scale, but advection dominated at the field and laboratory scales. This discussion leads to Conceptual Issue No. 2: Diffusion controls microscales; Advection controls macroscales. Other dimensional forms are possible and, especially, many other transport equations can be chosen. We will analyze them in Section 0.

Reactions.
The term represent the contributions of all reactions to the mass balance of species , expressed per unit volume of porous medium. This section is devoted to the evaluation of . Stated this way, it may sound irreverent, given the complexity of the concepts behind (reactions include all biogeochemical processes from mineral dissolution to hormones production). But, from a RT perspective, evaluating is all one needs. The first step is to define the chemical system, that is, the set of species we are going to deal with and their reactions. In essence the chemical system is defined by: is the chemical formula of the i-th species, is the number of species, is the number of reactions, and is the stoichiometric coefficient of species in reaction . Then, is given by where is the reaction rate (i.e., moles of reactants that evolve into products, per unit time and per unit volume of medium).

·
= 0 ( 11) RT usually involves many reactions and species, which makes it convenient to write these equations in a matrix-vector notation for concise writing. Thus, we rewrite Equations (9) and (10) as = ( 12) where is the × stoichiometric matrix, whose rows contain the stoichiometric coefficients in each reaction, and a vector of size , containing rates of all reactions, and a vector of size , containing the chemical formula of all species. To illustrate these concepts, let us consider a simple 6 of 46 calcite dissolution problem, which one might define as consisting of three species = ( ( ) , , ) and one reaction ( ( ) − − = 0). Then, the stoichiometric coefficients would be = (1, −1, −1), so that = , = − , and = − . Note, however, that in most groundwater problems, this definition would be insufficient because most of the carbonate evolves to bicarbonate ( ) under neutral pH conditions, so that both and would have to be included in the set of species and a new reaction ( − − = 0) in the set of reactions. Let this simple example serve to illustrate conceptual issue No. 3: "it is important to select the appropriate set of species and reactions". While this may sound trivial, it is not. There are few absolute objective rules to define an "appropriate" chemical system. We will come back to this issue. But for now, let us simply state that, in practice, it requires the active participation of a specialist (a geochemist for inorganic reactions, a biochemist for biochemical reactions, often both), a reminder of conceptual issue No. 1.
Mathematically oriented readers may use these definitions to view the chemical system as a vector space of dimension , and reactions as vectors of the changes in concentrations due to chemical reactions. This view is extensively discussed by Rubin [2] and Friedly and Rubin [46], who attribute it to Prigogine and Defay [47].
For the purpose of computing reaction rates, , it is convenient to distinguish between kinetic ("slow" in Error! Reference source not found.) and equilibrium ("fast" in Error! Reference source not found.) reactions. Obviously, the concept of slow or fast is relative to transport time or to the characteristic time of other reactions. For now, let us simply state that kinetic reactions are those whose rate is solely controlled by the concentrations at the point where the reaction is occurring, that is = ( ), where we have introduced the subindex to indicate that, in general, only a subset of all reactions are kinetic. ( ) are empirical functions. Lasaga [48] covers them extensively for inorganic reactions. Biochemical reactions will be reviewed in Section 0. It should be noted the apparent simplicity of this expression hides the complexity of most Earth processes, including life, of course.
Equilibrium reactions are governed by the Mass Action Law (MAL), which states that consists of the rows of corresponding to the equilibrium reactions, is the vector of their equilibrium constants and is the vector of their activities. Activity is a dimensionless quantity, related to the chemical potential of the species in solution, that quantifies its "effective" concentration (i.e., its availability to react). The activity of aqueous species should be equal to the molar fraction of the species for ideal, low salinity, solutions. Unfortunately, standard chemical data bases have been introduced for molalities, so that, ideally, the activity is equal to the molality. As salinity increases, the "effective" concentration decreases, so that the activity becomes = , where is its "activity coefficient". For moderate salinities, the molality is very close to the molarity (i.e., ≈ ). In general, however, one would need to write = / . MAL also applies for sorption to solid surfaces, where activity is expressed in mol per unit mass of solid. As we express concentrations of solid species in mol per unit mass of water (Section 0), the activity of adsorbed species becomes = /( ), where and are the density and volume fraction of the solid. For cation exchange the activity is an equivalent fraction, that is, where is the charge of the solute and is the Cation Exchange Capacity expressed in equivalents per unit of mass of solid. Moreover, in some models the electrical charge of the solid is assumed to be different from that of the liquid and an electrical double layer is assumed consisting of an electrical potential that decreases rapidly away from the charged solid surface. The equilibrium constant in the MAL (Equation (13)) is corrected for the free energy associated with the electrical potential [38]. As we will see in Section 0, model representation of sorption can range from very simple to very complex (Appelo and Postma, [49], chapter 6). For minerals one has to understand that the solid comprises of various mineral phases. Generally, it is assumed that they are pure phases and their activity is equal to 1.0. It is important to realize that the set of present minerals (or mineral zone), and hence the MAL, may change in space and time. The activity of gases is called fugacity, but is usually expressed in terms of pressure (in bars), so that under ideal conditions, the fugacity of a gas is equal to its partial pressure. For high pressures and temperatures, the behavior of gases is not ideal and their partial pressure needs to be multiplied by a fugacity coefficient, analogous in concept to the activity coefficient. Note that numer-7 of 46 ous species (pure minerals, constant pressure gases, ideal ) display constant activity. Therefore, from the chemical point of view, they are not unknowns. The number of unknowns in Equation (13) is − , where is the number of constant activity species, a feature we will profit from later. The equilibrium constant, K, is a characteristic of every chemical reaction that depends (a little) on pressure and (possibly a lot) on temperature. In practice, it is obtained by minimizing the Gibbs free energy of the solution or from chemical data bases. The selection of the database is a topic in itself. The most widely used databases for geochemical modeling are EQ3/EQ6 [50] and PHREEQE [51]. Specific data bases are often available for specific purposes. For example, [52] review data for CO2 storage. There are databases in the framework of cement/clay interactions [53] and the NEA (Nuclear Energy Agency) has produced thermochemical databases intended for the radioactive waste management community [54]. Since databases contain stoichiometric coefficients, as well as data to compute activity coefficients, and they are well tested, we recommend using them, at the very least to reduce errors during input data preparation (this can also be done by means of a good preprocessor).
Databases are also available for kinetic rate laws, especially for dissolution of minerals [55,56]. In general, they consist of rather empirical data from laboratory experiments. Therefore, their use at the field scale is not straightforward, as rates at the field scale are known to be different from the laboratory scale [57]. Databases are also available for kinetic biochemical reactions. Their use is less frequent because reaction rates can be written in numerous forms. A good one is the NIST Chemical Kinetics Database (https://kinetics.nist.gov/kinetics). Things get more complicated for biochemical reactions. PubChem (https://pubchem.ncbi.nlm.nih.gov) is a "popular chemical information resource" [58]. The problem is that using it requires some expertise. For example, searching for "trichloroethylene" yields 18 compounds (the "standard" 1,1,2-trichloroethene, but also 1,1,2-trichloro-2-deuterioethene, and many others that we would never call TCE) and 185 substances. But if you search for "trichloroethene", another synonym of TCE for us, pubchem yields 167 compounds, but only 80 substances.
The MAL describes equilibrium conditions that species must satisfy. Therefore, Equation (13) can be viewed as a set of constraints on valid values of concentrations. Or, returning to the vector space view of chemical systems, equilibrium reactions reduce the dimensionality of the concentrations space leading to Conceptual Issue No. 4. Equilibrium reactions reduce the size RT problems. The issue is illustrated in Error! Reference source not found. for a simple gypsum equilibrium example. The chemical system consists of three species (Gypsum, , and ). The stoichiometric matrix contains only one row describing the single reaction of the system (Error! Reference source not found., left). Concentrations (center) and activities (right) must lie on the equilibrium line. This line is straight (i.e., linear) in log-activities (Equation (13)), but non-linear in concentrations because of the logs and the activity coefficients. Note that reactions (red line) occur only when the system is driven away from the equilibrium line. They are given by Equation (12). Therefore, they are orthogonal to the equilibrium line in log-activities, recall that the coefficients of a (hyper) plane equation are the components of its orthogonal vector. The concentration space is convenient for mass balance, which is given by the transport equations. The activity space is convenient for chemical calculations. In this example, knowing one activity ( = , in Error! Reference source not found., right) allows easy calculation of the other ( = ).  The concepts of Error! Reference source not found. are general. Each equilibrium reactions reduces the dimension of the concentrations space by 1. Therefore, the number of degrees of freedom is − and one can choose a subset of primary species as basis of the subspace of "valid" concentrations. The remaining species are called secondary. Their activities can be calculated from those of primary species through the MAL. We can write the mass action laws in terms of primary and secondary species by splitting matrix and vector into two parts ( = ( , ); = ( , )), where contains the activities of the ( − ) primary species and the remaining ( ), secondary, species and where contains the stoichiometric coefficients for the primary and those for the secondary species. Substituting these definitions in Equation (13) and multiplying by yields: where * = −( ) and * = ( ) . In this way we obtain a new description of the same chemical system, which is now defined by equilibrium constants * instead of and of stoichiometric matrix for the primary species * instead of . RT computations are simplified by introducing components, defined here as chemical entities that are not affected by equilibrium reactions. To this end, it is convenient to define the components matrix as the ( − ) × kernel matrix (or null space) of . That is, = . In vector space nomenclature, this matrix defines the orthogonal complement to the equilibrium subspace of valid concentrations. Matrix is not unique. In fact, the set of possible component matrices is in itself a vector space of dimension − . Hence, there are numerous methods to obtain it. Molins et al. (2004) describe some of them [59]. In particular, they define approaches to ensure that facilitates reactive transport computations. Still, the most widely used method is based on Gauss-Jordan elimination [60], which consists of defining the component matrix as The identity matrix in the front portion of the component matrix (i.e., the one coinciding with the primary species) ensures that the components are somehow associated to the primary species. This association is sometimes weak, but examining the components often provides some insight into the chemical system. Since some species are immobile, = ( ) , components may also contain an immobile part. This can be formalized by splitting = ( ), so that = = + , where = , and = . As mentioned before constant activity species are not really unknowns during chemical calculations and it is profitable to eliminate them. This is precisely what Saaltink et al. (1998) did [61]. They used Equation (15), but then multiplied the components matrix by an elimination matrix, thus reducing the number of primary species. The resulting method was the first one consistent with Gibbs rule. Yet, another possibility is to eliminate constant activity species from the outset, which leads to a reduced stoichiometric matrix (red square in Error! Reference source not found.), prior to using Equation (15) [62]. In practice, this implies reducing the unknowns from to − . These concepts are discussed further in Error! Reference source not found., where two pure mineral (calcite and gypsum) equilibrium reactions are considered. The system consists now of 5 species and 2 reactions, but its dimension is reduced to 3 by neglecting constant activity species. The component space is 1D, so that any specie may serve as primary. In the full concentration space, the single component is represented as a plane orthogonal to the equilibrium (and, thus, to both the calcite and gypsum equilibrium subspaces) reactions (red arrows are contained in this plane). Speciation calculations are needed to obtain the vector of concentrations of all variable activity species from or . In the first case, Equation (14) yields directly . Concentrations are then obtained from the activity coefficients as = ⁄ . Iterations may still be needed because = ( ). But this dependence is weak. A more problematic case arises when speciation needs to be performed from components. In this case, the concentrations are obtained from the MAL of the equilibrium reactions (Equation (12)) and the − definitions of components (Equation Error! Reference source not found.). The solution of these equations can be difficult and may require special solutions methods [63][64][65]. A number of codes are available to perform speciation calculations, such as PHREEQC [66], CHEPROO [67] and GEMSK [68]. Still, since these calculations need to be repeated numerous times (number of nodes times number of iterations times number of time steps), it makes them ideal for artificial neural networks [69], which is possibly what the future holds for RT.

Coupled Transport and Reactions
Reactive transport is the phenomenon resulting from the interaction and coupling of solute transport, Section 0, and chemical reactions, explained in the previous section. Thus, reactive transport equations simply represent the mass balance of chemical species subject to both transport processes and chemical reactions. This statement suggests posing the problem as solving the system of equations: This formulation of reactive transport contains + unknowns. These are matched by the N mass balances of all species (Equation (17)), plus the MAL (Equation Error! Reference source not found.), and the N definitions of kinetic reaction rates. The latter can be eliminated by substitution of Equation (18) into (17). Still, the large number of coupled non-linear equations makes the problem intractable. Therefore, various strategies have been proposed to simplify the problem [70,71]. But most of them are based on the equilibrium reaction rates by multiplying Equation (17) by a component matrix, which leads to: Care must be taken in acknowledging the phase to which species belong. Only the aqueous fraction of components, = , is transported. In fact, the transport operator does not apply to the solid (immobile) fraction, = . An important issue arises here for problems where the solid phase consists of several ideal (constant activity) minerals, depending on how components have been derived. If components were derived using all reactions in the system, then the transport equations will be the same throughout the domain, which sounds advantageous. However, if components were defined after eliminating constant activity species (reduced stoichiometric and components matrices) the transport problem would be different for every mineral zone. While this sounds disadvantageous, the reality is that the problems are indeed different, as discussed by Rubin [2]. As a result, sharp concentration fronts develop at the interface between different mineral zones because water with the conditions of the upstream chemical system must adapt instantaneously to the chemical system downstream (Error! Reference source not found.). Figure 6. Transport of pulse injected at x=-8 m and t=-18 d in a medium consisting of two minerals (similar to the one in Error! Reference source not found.). A front occurs for A3 at x>50 m, where it equilibrates with mineral 2, causing its sharp precipitation and an apparent discontinuity in concentrations. The analytical solution (A_r1) [62] is displayed for the precipitation rate of Min. 1.
The apparent discontinuity in Error! Reference source not found. simply reflects that the set of equations in each mineral zone will be different, thus yielding different solutions. This leads to Conceptual issue 5: Variations in mineral zone, leads to sharp fronts in aqueous chemistry. The transport equation does not allow discontinuities. The solution will display a sharp front with a width comparable to dispersivity (0.3 m in the example of Error! Reference source not found.). In practice, this will cause convergence problems to numerical methods that do not account explicitly for chemical heterogeneity and force modelers to use large dispersion coefficients, which may cause other problems, as we discuss in Section Error! Reference source not found.. Ways to overcome the problem numerically are described in Section 0.

Basic reactive transport paradigms
In the previous section, we have analyzed the ADE for transport, the formulation of general chemical reactions, and how to couple them to write a RT model. We now discuss, how the coupling leads to vastly different regimes, as a first step to understand how the mathematically "innocent-looking" reactions term, , can produce the rich biochemical diversity of Nature.

Sorption. Conceptual issue 6: Each species travels at a different velocity
Arguably, the simplest reactive transport model is the one of adsorption of a solute onto the solid surface (A + X ⇌ X − A, where A is a cation in solution, X is a negatively charged "sorption site", and X − A is the adsorbed cation occupying a sorption site). The mass action law for this reaction is = . However, equilibrium is usually obtained experimentally from batch tests where (moles of A sorbed per unit mass of solid) is measured for different values of . Results ( vs plots) often display a straight line through the origin (its slope is called distribution coefficient), that is Note that the units of are weird (Vw/Ms) and unrelated to its meaning. Instead, sometimes the reported value is the dimensionless partition coefficient = ⁄ , which relates the sorbed to the aqueous solute mass in a volume of aquifer [72]. Including sorption into the transport equation simply implies realizing that every mol of that sorbs onto the solid disappears from solution. Therefore, sorption can be viewed as a sink term equal to = − ⁄ , where is the bulk density, mass of solid per unit volume of aquifer. Including this term into the solute transport equation (Equation (17)), while using Equation (21) for leads to where is the "retardation factor", given by = 1 + ⁄ = 1 + . The impact of can be deduced from the velocity of (coefficient of the gradient divided by the coefficient of the time derivative), which equals to ⁄ , compared to the velocity of water, ⁄ (i.e., is indeed retarded by a factor with respect to water). Another way to look at it relies on noticing that the left hand side of Equation (22) can be rewritten as ( / ) ⁄ , which cause Equation (22) to be identical to that of a conservative solute, if is substituted by / (i.e., the time scale is slowed down by a factor , hence the name retardation). This is illustrated in Error! Reference source not found. for a case with constant = 2. While linear adsorption is a nice simple concept, several remarks can be made. First, recall that, despite its empirical nature, Equation (21) is a rewrite of the mass action law. That is, includes the effect of expressing the concentration of − as per unit mass of solid (i.e., = / ). Therefore, a constant requires that the activity coefficients are constant and that the concentration of sorption sites is constant (i.e., very large compared to A). In practice, varies in response to changes in salinity (which changes activity coefficients and competes for sorption sites) or temperature (changes ) or pH (protons also compete). A simple example illustrates the effect of varying . Assume that salinity, which reduces , enters into an aquifer. This type of problem has been studied in depth for Radium, which is widely used as a tracer of submarine groundwater discharge [73]. Ra tends to sorb strongly and its retardation depends on salinity (Michael et al., [72], report values of ranging between 10 and 10 5 for fresh water and between 1.5 and 10 4 for salt water). Under these conditions, Equation (22) is not valid because is not constant. Thus, the left hand side must be changed to ( Since depends on salinity, the last term becomes a Ra source term when salinity increases, which explains the sharp increase in at the salinity front (dashed lines in Error! Reference source not found. phenomenon is well-known in the Ra literature, Equation (22) and retardation are so deeply ingrained in the hydrological literature that Equation (22) remains the standard. Sorbed ions may be displaced by other solutes [74]. Charge balance requires that when one ion sorbs, another one must desorb, which leads to competition for sorption sites and to ion exchange reactions. Cation Exchange (CE) is nicely discussed by Valocchi et al. [75], who explains the chromatographic effect (separation of cations depending on their selectivity). Details on the relationship between adsorption and ion exchange are discussed by Charbeneau [76] and Appelo [77,78]. In essence, CE should be considered always, but it is not easy to handle because selectivity coefficients depend on numerous parameters and because access to exchange sites may be delayed in time (think of diffusion within clay layers). As a result, specific models are needed to simulate CE [79][80][81]. This complexity explain why numerous empirical models have been developed to estimate CE parameters (e.g., [82]), including Artificial Neural Networks [83]. But, even if the CE model is known, spatial variability of minerals hinders their rigorous application.
The summary of this short discussion on cation exchange is two-fold. First, the retardation representation of sorption displays numerous limitations. But given the complexity of more accurate exchange models, it must be considered a valid option when simplicity is needed (simplicity is usually associated to robustness). This is why safety assessment of nuclear waste repositories typically relies on retardation and seeks empirical approximations to deal with its limitations (e.g., [84][85][86][87]) The second remark is that different species will travel at different velocities, which places a heavy demand of numerical models, as transport may be diffusion dominated for highly retarded species and advection dominated for (nearly) conservative species.

Kinetic reactions. Conceptual issue 7: kinetics are important and complex.
Kinetic reactions are those whose rate is solely controlled by concentrations at the point where the reaction is occurring. The simplest kinetic reaction is the first order build-up or decay reaction ( = ), which may be appropriate for biodegradation when the concentration of the reacting species is the rate limiting factor (e.g., when the concentration of oxygen is constant or when the concentrations of catalyzers remain constant). For the purpose of discussion, let us consider -th order reactions; whose rate is given by In the absence of transport, the integration of Equation (23) yields the evolution of concentration The characteristic time for these reactions is = 1 ⁄ . To find out the relevance of a kinetic reaction, one should compare the reaction time with the transport time. This comparison is performed with the Damkhöler number, defined as: where is the characteristic transport time. The Damkhöler number can also be defined as the reaction rate divided by the transport rate, which leads to equivalent expressions. Note that this expression raises two issues. First, the reaction time depends on for this simple kinetic model. In general, it may depend on numerous other factors. Specifically, biochemical reactions can be very fast or slow depending on the type and development of bacterial communities (this issue will be addressed in Section 0). Second, choosing an appropriate transport time is a relevant conceptual issue. As characteristic transport time it is frequent to adopt the advection time. However, when interested in pore scale reactions, it may be more appropriate to adopt the diffusion time (recall C.I. No. 2: Diffusion controls microscales). The resulting number would be the product of (defined with the advective time) and the Peclet number.
The Damkhöler number is used for assessing whether a given kinetic reaction needs to be taken into account. When is small (much smaller than 1, say < 0.1), then one may presume that the reaction has barely modified the concentration. That is, transport is so fast that the reaction can be neglected. Reversely, if is large (much larger than 1, say > 10), then one may guess that the solute has been exhausted and that a new equilibrium has been reached, so that the reaction may be assumed to be an equilibrium reaction. This analysis is easy to confirm for first order reactions (Error! Reference source not found.). The analysis is more complex, and the usefulness of the Damkhöler number is more limited for non-linear kinetics. The time evolution of solutes undergoing kinetic reactions of order 1, 2 and 3 are also displayed in Error! Reference source not found.. It is clear that close to 1 implies a sizeable reaction in all cases. However, it is not true that a large implies that the solution is close to equilibrium. is not necessarily true for non-linear kinetics (right). Da=1 implies a sizeable reaction, but note that significant amounts of solute remain for Da=100, if n=4. Note also the abrupt approach to zero for n<1, which may cause convergence problems.
One might be tempted to use kinetic laws of the type in Equation (23) in an attempt to simplify complex models (recall, again, that simplicity helps robustness). However, Reigner et al. [23] point out a six orders of magnitude range in reported rate constants for anaerobic oxidation of methane. Such a broad range suggests that an important factor may be missing. They suggest that microbial concentration may be the missing factor. Similarly, White and Peterson [88] reviewed mineral dissolution rates and found laboratory rates to be 2-4 orders of magnitude larger than observed field rates. They attributed this difference to differences in access to the minerals (see also Maher et al. [89]). Both examples point out that too simple models or direct extrapolation of laboratory rates to field conditions may not be appropriate. Local conditions are what define actual reaction rates. Valhondo et al. [90] found three orders of magnitude ranges on the removal rates of emerging organic compounds (EOCs).
And, yet, the rates they observed in the field turned out to be even faster, an additional order of magnitude for some EOCs, which they attributed to the diversity of microbial populations in their enhanced Soil Aquifer Treatment design [91].
The summary of this discussion is that kinetic reactions, while numerically a simple sink/source term, turn out to be the main source of both numerical and conceptual problems in practice. From the numerical point of view, they can cause convergence problems when non-linear and fast. From the RT conceptual point of view, they often imply coupling and, when fast, they may lead to chemical zonation similar to equilibrium reactions [92]. From the chemical point of view, they may be affected by numerous factors (local variability, access to reaction sites, and presence of catalyzers). This is especially true in biochemical processes, which we review in Section 0. In this context, the Damkhöler number should be viewed as a rough preliminary indicator of the relative importance of kinetic reactions, but cannot be used blindly.
3.3. Equilibrium. Conceptual Issue No. 8 The rate of fast reactions is controlled by mixing and access to reaction sites .
An easy paradigm for equilibrium reactions is the one of two ions, and , that precipitate as an immobile solid mineral in equilibrium ( + ⇄ ). The solution to this problem illustrates the methodology for solving reacting transport [62]. In essence, assuming ideal conditions (pure mineral and constant activity coefficients, so that they can be incorporated in the equilibrium constant), the problem is formulated as Equations (17) and (18): is the concentration of , and is the concentration of the mineral. These four equations contain four unknowns: , , , and . They are coupled, and there is no explicit expression for the reaction rate, , which results from the solution. The problem, formulated like this, is highly non-linear and coupled. Solving it requires eliminating some variables, which is done by the use of components, and decoupling transport and chemistry. The procedure can be formalized as consisting of five steps.
Step 1: Definition of components. The stoichiometric and component matrices for this case are given in Error! Reference source not found.. Multiplying the vector of the first two equations times the component matrix is equivalent to subtracting Equation (26) from (25), which yields − · · + = 0 ( 29) where = − , which is conservative indeed a component, and easy to solve. We have reduced the problem from 4 coupled non-linear PDEs to 1 linear PDE (recall conceptual issue No. 4, Equilibrium reactions reduce the size of RT problems, although this example is extreme).
Step 2: Solving for components. In our example, this simply requires solving Equation (29) for . In general, where kinetic reactions may also occur, and where minerals may not be present throughout the domain, this step may require iterations. For now, let us assume that Equation (29) can be solved for .
Step 3: Speciation calculations. Once components have been obtained, their definition is used together with the mass action law for obtaining the concentrations of aqueous species. These calculations can be complex (recall discussion on speciation codes in Section 0). But in this case, speciation is easy, because it consists of solving the system of equations given by Equation (26)  Step 4: Evaluation of the reaction rate. Substitution of the concentration of the secondary species, , into its transport equation leads to [62] = · · ( 30) where we have used bold-face to emphasize that the equation is valid for the case of multiple components and equilibrium reactions. In fact, this equation may also include temperature by including thermal conduction in . Beyond its basic simplicity, this equation illustrates that under our ideal conditions the reaction rate is always positive (i.e., precipitation occurs because ⁄ is positive), as discussed by Rubin [2]. This observation could be deduced from Error! Reference source not found., which also shows that the slope of the reaction amount (red arrow) is equal to the ratio of the stoichiometric coefficients (1 in our example). Note, however, that this is a result of having assumed constant activity coefficients. In reality, equilibrium is controlled by numerous other factors, such as temperature or ion strength [93], which affect ⁄ . As a result, dissolution may occur. For example, changes in ion strength produce dissolution of calcite when calcite saturated fresh and salt water mix, which explains the development of coastal karst [94,95].
Beyond the limitations of the example, Equation (30)  also on transport processes, controlling the gradient of . This may not be surprising but implies that the two effects are separable, so that after all, chemistry and transport are not as tightly linked as one might have thought. We will come back to this is issue when discussing the "water mixing approach" to transport in Section 0. Another paradoxical feature of Equation (30) is that reaction does not necessarily take place where concentrations attain their maximum values. In fact, the reaction rate equals zero when is maximum or minimum. It is also interesting to notice that, in the absence of the first contribution on the right hand side of Equation (30), all terms are proportional to the dispersion tensor, , thus strengthening the relevance of mixing processes to the development of fast reactions. In particular, the term is a measurement of the mixing rate, which is consistent with the concept of dilution index, as defined by Kitanidis [96] on the basis of entropy arguments. This result also suggests that evaluating mixing ratios may help to properly identify not only the sources of water [97], but also the geochemical processes occurring in the system. In fact, De Simoni et al. [62] produced a similar expression in terms of the mixing ratios of end member waters.
Step 5: Computation of constant activity species. Once reaction rates have been obtained, the mass balances of constant activity species are used to obtain their concentrations. In our case, the concentration of mineral is obtained using in Equation (28), which had not been used until now. The procedure described here can be considered unrealistically simple, because it does not account for kinetic or exchange reactions, which are universal. Full solution methods, accounting for all the complexities of real systems will be described in Section 0. These equations were used in the early part of Error! Reference source not found. to compute mineral precipitation in response to a pulse injection, for which the transport Equation (29) can be solved analytically.
The summary of this section is that mixing, is fundamental for RT modeling. As we will see in the next section, a problem with the ADE is that it does not simulate mixing correctly, which has led to numerous efforts to seek alternative transport formulations.

Heterogeneity
The advection-dispersion equation (ADE) (Equation (2)) is the standard formulation for solute transport in all kinds of media, not only porous media. However, it does not provide a good representation of transport in natural materials. Tracer experiments in the field rarely display the symmetric spatial distribution predicted by the ADE and breakthrough curves (BTCs) always display tailing, inconsistent with the ADE [98,99]. Calibrated dispersivity grow with the experiment scale [100,101], and calibrated advective porosities often depend on the residence time [102] or on the flow direction [103,104]. This type of behavior is called anomalous, which we find ironic because it is the one always observed in nature, not only in porous media. Since an adequate transport equation is critical for RT, we devote this section to discuss alternative transport formulations. We must state from the outset that (1) no universal convincingly adequate equation has been found, (2) the problem largely arises from broad range of scales involved in RT (Error! Reference source not found.), and (3) heterogeneity and preferential flow paths occur at all scales. Therefore, here we outline the evolution of transport concepts, as related to RT, in the hope of conveying the problem and preventing error as the ones we have made in the recent past.  Efforts to quantify anomalous transport started during the last part of the 20 th century. They concentrated on explaining the scale dependence of dispersivity and were largely successful in demonstrating that dispersivity increases with scale due to the scale dependence of heterogeneity [105][106][107]. By the turn of the century, it became apparent that scale dependence of dispersivity was not sufficient, as it did not help in addressing the other limitations of the transport equation (tailing, time dependence of apparent porosity, etc.).

The search for an effective transport equation
After establishing that the ADE is a poor representation of transport, which conform our conceptual issue No. 9, researchers focused on finding "effective equations" in the hope to find a better representation than the ADE. Numerous formulations have been proposed at different scales to address anomalous transport [108]. We make here a short detour into methods to solve the ADE and, in general solute transport, so as to properly frame the numerous equations that were proposed to overcome the limitations of the ADE. They are schematically described in Error! Reference source not found., which is somewhat crowded, but whose goal is to show that, despite their mathematical difficulty, they are not that different (in fact, they are virtually identical!). Solution methods for the ADE can be divided into Eulerian and Lagrangian (top row of Error! Reference source not found.). Eulerian methods discretize the equation in space and time over a fixed grid. Lagrangian methods solve the equation over a grid that moves with the fluid velocity and it is used to solve for dispersion. Lagrangian grids can become highly deformed if the velocity field is spatially variable, which is the rule. This motivated particle and Random Walk (RW) methods, based on Einstein's [109] concept that Fick's Law results from the mean behavior of particles subject to Brownian motion [110,111]. In RW methods, the goal is to approximate the solute by a cloud of particles that move with an advection component and a random (Brownian) component. In fact, if the random component displays a finite variance, the resulting behavior of the cloud is identical to the ADE. Particle distributions can be transformed into concentrations by Binning (Error! Reference source not found.d) or, far better, by associating a kernel to each particle (Bells method, Error! Reference source not found.g) [112,113]. Still, the usual goal of RW methods is not so much reconstructing concentrations as to obtain transport information from the distribution of particles (e.g., evolution of their mean location or their spatial variance).
RW methods require very large numbers of particles especially for pore scale simulations (recall that Lagrangian methods are not well suited for microscales). At these scales, RW require fractions of Avogadro number particles, which are best approximated by molecular dynamics models. Two sets of methods that approximate molecular dynamics and have been used successfully for RT are the smoothed particle hydrodynamics (SPH) and the lattice-Boltzman methods. SPH is similar to RW, except that particles are "smoothed" with a finite support volume (much like in the "bells" method) with the added advantage that they also allow mechanical particle interactions, which implies solving for momentum conservation. SPH allowed modelers to perform highly resolved simulations that yield insights into pore scale transport and reactions [114,115]. Lattize-Botzman method is similar in treating particle interactions, but on a fixed grid [116]. A description of methods to approximate molecular behavior is presented by Meakin and Tartakovsky [34].
The search for effective transport equations, alternative to the ADE, started with the Continuous Time Random Walk (CTRW) which, in its basic form, can be viewed as a RW in which the time-steps for velocity transitions are random [117]. Variable time steps allow the modeler to represent non-fickian behavior (both waiting in immobile zones and fast travelling through preferential flow paths). By selecting appropriately the distribution of times between velocity transitions, a broad range of transport behaviors can be obtained. CTRW became very popular because it could explain a large number of non-trivial observations, such as tailing in field and laboratory experiments [5], proper scaling of variable time push-pull tests [118] and especially pore network model outputs [119]. A second option is to substitute the time and/or space derivatives in the ADE by fractional order derivatives [120,121]. While using fractional order derivatives sounds awkward, the concept dates back to Leibnitz himself and has been the subject of much mathematical work [122]. Its use in hydrology emerges from a RW, much like the CTRW, but where the random component arises from a fractional Levy distribution or, a bit more generally a fractional Brownian motion, which allows simulating long jumps. The approach became popular because it can be viewed as emerging from the popular fractal view of Nature and because fractional behavior has been frequently observed in well hydraulics [123]. Lagrangian methods on (b) moving grid and (c) particle tracking (or Random Walk); recovering concentrations from particles (d) binning, or (g) the «bells» method; Effective transport equations can be viewed as (e) a multicontinuum formulation, or (f) a random walk with «anomalous» noise.
A third approach was introduced at about the same time by considering the medium as consisting of a mobile region, representing fast flow paths, superimposed to an immobile region. Both exchange mass by means of a memory function, , that represents the delay of slowly flowing portions of the medium and the diffusive exchange between them [124]. As it turns out, by appropriately choosing the memory function, the three methods are identical [125,126]. From a RW perspective, the memory function can be viewed as the residence time distribution in immobile water. Therefore, these methods are directly to the Multi-Rate Mass Transfer (MRMT) method [127]. This method considers a discrete number of well mixed immobile zones that exchange mass with transfer coefficients, inverse of mean residence times, ( = 1, ⋯ , , number of immobile zones). The probability of each zone, , can be viewed as the fraction of immobile porosity it occupies. Therefore, it can be seen as a discrete revision of the other effective equations. Numerous variants of these methods are possible [128][129][130]. A variation of MRMT are the hybrid models in which the set of immobile zones is substituted by a pore scale model, which allows reproducing reactions at the pore scale where they really take place [131][132][133][134][135][136][137]. The different methods for this coupling are described by Molins and Knaber [30].

Lessons learned, formal upscaling, and on-going efforts
A problem with these effective equations is that they are largely empirical. The CTRW and memory function approaches were originally designed to simulate transport in fractured media. The MRMT was devised to handle variable retention times for non-instantaneous adsorption. As such, they suffer two fundamental problems: (1) their parameters do not result from descriptions of heterogeneity and geological understanding, but need to be calibrated (virtually impossible in most cases, where transport may take years), and (2) they were developed before transport researchers realized the importance of mixing. Therefore, while they do reproduce many of the anomalous features of transport though natural media, they fail to adequately distinguish mixing from spreading (Error! Reference source not found.). Worse, with the exception of MRMT, they do not provide a value of concentration in its geological setting, which is required for RT. For this reason, we have preferred MRMT. We have met with some successes, such as reproducing reaction rates in a complex heterogeneous system with a simple MRMT model [138] or geochemical localization, which stands for the occurrence of reactions that would not be possible in a locally equilibrated medium [139,140]. Still, much of the tailing occurs due to very slow advection and does not necessarily imply mixing [141,142]. The need to distinguish between mixing and spreading has prompted Soler-Sagarra et al. [143] to rewrite the dispersion term as ∇( ∇ ), where is the dispersion flux (Error! Reference source not found.), instead of ∇( ∇ ), which implies mixing. Therefore. It is not surprising that MRMT failed to meet the criteria that an effective equation for anomalous transport should meet: proper reproduction of advection, dispersion and mixing, hopefully linked to geological understanding [144,145] as shown in Error! Reference source not found.. To us, this is a relevant conceptual issue (No. 9): Solutes are displaced by advection, spread by dispersion, diluted by mixing. Figure 11. The ADE and effective equations fail to distinguish spreading and mixing, thus overestimating mixing and fast reactions, even if spatial (vertically integrated during sampling) and BTCs look fine.
By 2010, many started to believe that effective equations were failing, which prompted Neuman and Tartakovsky [146] to recall formal upscaling. Upscaling is the process of finding large scale behavior from small-scale principles. The problem is that processes governing small scale behavior (e.g., molecular diffusion) may become irrelevant at large scales, whereas new processes (e.g., dispersion) emerge in response to the change in scale. Recent reviews [147,148] summarize upscaling procedures. But the basic problem is that such procedures are well-defined to evaluate how a parameter (e.g., dispersivity) changes with scale. There is no well-defined approach to identify how an equation evolves, which explains why effective equations were largely empirical.
Therefore, recent efforts on characterizing transport through heterogeneous media have concentrated two avenues. The first one is to follow the traditional stochastic upscaling methodology that consists of (1) generate realizations of heterogeneous fields of hydraulic conductivity (and/or chemical properties of immobile phases, pore geometries, etc.), (2) simulate explicitly RT through each realization, and (3) extract relevant statistics and/or conceptual lessons from the simulation results. The three steps display uncertainties. Proper generation of geological variability remains a challenge [147,149,150]. RT transport difficulties are the subject of this paper. And much of the stochastic transport research on dispersion upscaling turned out to be of little relevance to RT because of computing ensemble statistics instead of analyzing the behavior of each realization [151,152]. Even if the procedure is properly followed, formal upscaling remains difficult, because physical and chemical heterogeneity cannot be upscaled separately [153]: upscaling reactions depends on physical heterogeneity and upscaling transport is affected by chemical heterogeneity. Therefore, heterogeneity of chemical micro-environments must be acknowledged when reactions are upscaled from pore to continuum scale [11,140,154].
Still, significant advances have been made on the third step. For example, advances have been made on how transverse mixing controls mixing rate and reaction extent by controlling the interfacial area between reacting water bodies or through flow focusing [155,156]. This brings renewed interest on the role of transverse dispersion, which has not received sufficient attention in the past [157,158]. But, from an effective transport perspective, the most important finding is that for RW methods to reproduce transport through heterogeneous media it is best to produce velocity transitions not at fixed intervals in time, but at fixed intervals in space [159]. Recall that RW methods were originally devised to reproduce molecular diffusion. Anomalous transport caused researchers to attribute random steps not only to molecular diffusion, but also to velocity variability ( in Error! Reference source not found.). While diffusion steps are well reproduced by sequential time steps, the impact of velocity fluctuations is not (if a particle is moving at 1 mm/century today, chances are that it will stay in the same low permeability region tomorrow). This prompted Le Borgne et al. [160] to propose an alternative CTRW approach and has led to a revival in the search for effective transport equations.
Efforts still continue to find an effective equation. Thus, two new steps are required: (i) the velocity distribution of the domain has to be discretized [161] (ii) and a matrix transition has to be calculated and used for the transitions [162][163][164]. This lead to a new phase space formulation [165], which means that velocity is considered as another variable of solute, like space and time. An immediate extension from the above discussion is to distinguish between velocity transitions due to heterogeneity, which are markovian in space and do not produce mixing, and transitions due to molecular diffusion, which are markovian in time and cause mixing [143]. Sherman et al. [148] review spatial Markov approaches. At this stage, it is not clear if an effective transport equation will ever be found. But we are convinced that in order to properly reproduce advection, dispersion, and the mixing enhancement caused by dispersion, the effective equation will be markovian in space, which will facilitate accommodating geological understanding (i.e., the fact that improved geological knowledge reduces uncertainty but not variability).

Biochemical reactions. Conceptual issue No. 10: Biochemical reactions involve microorganisms.
Biochemical processes are those involving living organisms, which thrive on the energy liberated by the reactions. Biochemical reactions typically involve exchange of electrons between a species that releases them (i.e., an electron donor, which is oxidized) and one that receives them (an electron acceptor, which is reduced). As stated by Thullner and Regnier [166], "… environmental conditions can often be traced back directly or indirectly to the oxidation of reduced carbon molecules… which… control, among others, the recycling of inorganic carbon and nutrients, the precipitation/dissolution of carbonate and sulfide minerals, the pH of (pore)-waters…". Therefore, we start by reviewing the basic OC (Organic Carbon) decay and redox concepts, prior to discussing the concepts related to kinetics and biofilms.

Organic carbon decay and redox sequence
Decay of organic carbon (OC, any C with valence less than +4) involves transforming it into another organic species with higher valence or full mineralization to inorganic carbon (valence +4).  168]. Error! Reference source not found. illustrates the redox sequence and its effect on groundwater chemistry and mineralogy for a case without transport and all reactions in equilibrium. As water flows, this temporal sequence tends to display a spatial zonation [169]. These reactions are catalyzed by microorganisms and, hence, imply biomass growth (cell synthesis) and decay that one may want to take into account. A distinction is made between the growth of biomass, and endogenous decay. The full reactions of endogenous decay are simply reactions R1+R3 to R1+R8, depending on the TEA. Cell synthesis reactions result from combining the oxidation of organic matter by some TEA which yields the energy, and the formation of the cell material itself, which consumes the energy. The yield factor ( ) is the amount of cell material synthesized per amount of oxidized OC. It can be expressed in mass fraction or electron equivalent fraction [167]. For example, aerobic degradation of CH2O can be written as: y(R1−R2) + (1−y)(R1+R3): 0.25CH2O + 0.05yNH4 + + 0.25yO2 → 0.05yC5H7O2N(bio) + 0.25(1−y)HCO3 − + 0.15yH2O + (0.25−0.2y)H + where the chemical formula of biomass is assumed to be C5H7O2N. For other TEAPs, one can construct cell synthesis reactions in a similar way. In theory, the yield factor depends on the free Gibbs energy of the reactions involved and energy losses due to inefficiencies in the cell and, thus, depends on chemical composition of groundwater. In practice, yield factors are usually assumed constant parameters. Since yield depends on the released energy, the yield factor drops along the redox sequence. Values for yield factors range from less than 0.1 for sulphate reduction and methanogenesis to 0.6 or 0.7 for aerobic reactions ( [167], table 3.1). Alternatives to the yield concept are the Carbon Use and the Microbial Growth Efficiencies, defined as the ratio of growth to overall C use (i.e., growth plus maintenance) which are widely used in soil science [170].
These reactions display characteristic times ranging from less than 1 day (aerobic oxidation) to more than 1 year (sulfate reduction or methanogenesis), depending on the local conditions. Thus they may be considered slow for many field problems and, definitely slow for laboratory scale problems, thus requiring kinetic rate laws (see discussion of Cirpka and Valocchi, [171] for semianalytical solutions under "fast" conditions). Microorganisms in the subsurface can catalyze these reactions through numerous pathways, steps and biocatalyzers or enzymes. The process is as complex as life, literally. Therefore, simple kinetic models have been proposed to handle redox reactions.

Traditional kinetic models (Monod and Michaelis-Menten)
A very simple enzyme kinetic rate law was developed by Michaelis and Menten [172]. They assumed that an enzyme, , catalyzes the transformation of a substrate, , to a product, , according to two reactions ( + ↔ → + , where is an enzyme substrate complex of and ). Assuming that the first reaction occurs in equilibrium (i.e., the MAL is valid), that the total is fixed and much larger than , and that the second is kinetic with its rate proportional to , it is easy to show that the rate of production of is is the kinetic rate constant of the second reaction and , is the equilibrium constant of the first reaction (also called the Michaelis-Menten coefficient or half-saturation constant). This rate law was also derived empirically by Monod [173]. This equation expresses the dependence of on limiting factors in a simple way: grows linearly on , with slope ⁄ , when is small (approximately up to = 0.5 , when = , which is why it is called half-saturation constant), but tends to the maximum rate, , when is much larger than . The elegance and simplicity of this dependence on the limiting factor explains why the Michaelis-Menten (Monod) factor has become widely used far beyond its original context.
An interesting extension of the Michaelis-Menten equation is its use to quantify the rate of cell synthesis [167,174] is the biomass concentration, which can be considered proportional to the term + (total enzyme plus enzyme complex). The rate of cell synthesis may be limited not only by the substrate, but also by TEAs, such as or , which is qualitatively reproduced by a multiplicative Monod equation [167,175]: = 33) where is the TEA. For the rate of endogenous decay, one often assumes a linear relation with the biomass concentration.
= 34) This type of approach is widely used to model organic carbon oxidation rates in the redox sequence. However, instead of including biomass and endogenous decay (presumably, one would have to include a different type of biomass for each TEAP), the standard option is to include inhibiting factors ( = , / , + , where , is the inhibition constant and subscript refers to inhibiting electron acceptor in the redox sequence, = 1 for , = 2 for , etc.). These factors display a behavior opposite to the Monod factors (constant and equal to 1 when is too small to inhibit the process and tends to 0, thus inhibiting it, when becomes larger than , ). Several multiplicative factors can be added. The OC decay law becomes (e.g., [176]):

35)
where subscript refers to the TEAP in the redox sequence. This equation replicates the preferential use of TEAs in the redox sequence. Numerous variants of equation (35) have been used. Instead of a Monod factor for , a linear relation can be used [177]. Also, instead of the inhibition factors, one can use switch functions that are either one or zero, depending on the concentration of an electron acceptor [178]. As equation (35) does not consider Gibbs free energy, it could give rates that are thermodynamically unfavorable [179]. Therefore, factors may have to be included to correct this (e.g., [180]). This equation can also be modified to allow for negative concentrations that may arise due to instabilities during transport solution [181].
An alternative approach for organic carbon decay in the redox sequence is the Partial Equilibrium Approach [182]. Rather than formulating rate laws for each TEAP, in equation (35), it assumes the oxidation of organic carbon (R1) to be the rate-determining step and controlled by a kinetic rate law, which could be first order, Monod or multiplicative Monod. The other redox reactions (R3 to R8) are assumed in equilibrium. The advantage is its simplicity and the need of less parameters. In addition, it guarantees that electron acceptors will be reduced only when it is thermodynamically favorable. However, one should be careful with the assumption on equilibrium. For instance, in some cases it may produce spurious reactions, such as oxidation of N2 by O2 (Saaltink and Rodríguez-Escales, 2021, submitted).
These approaches have been used widely. General reviews have been written by [183] and [166]. One important application is the modelling of redox sequences in aquifers with contamination plumes. To this both the multiplicative Monod equation (35) is applied (e.g., [178,184,185]) as well as the Partial Equilibrium Approach (e.g., [186]). Another important application, though at a much smaller scale, is the redox sequence in marine or lake sediments (e.g., [187][188][189]). Reviews exists for the modelling of this particular phenomenon [22,190]. Modelling of carbon oxidation and redox sequence has also been done for domains closer to the surface, such as natural wetlands (e.g., [191,192]) and managed aquifer recharge (e.g., [193][194][195]).

Microbial based kinetic models
Models that use equation (35) and, though to a less extent, the partial equilibrium approach, require many parameters (for each TEAP one needs a value for , , , , and several , ). Usually, they are obtained from the literature or by fitting the model to measured concentrations. If we compare equation (33) and (35), we can see that the maximum rate of a particular TEAP, , , is related to the biomass, , capable of using the particular electron acceptor . Therefore, it might seem reasonable to use the biomass concentration to control degradation rates. In fact, there are several methods to identify and quantify biomass and its activity [169,174], including counting colony-forming units in Petri dishes, measuring the end product for the electron acceptor of concern in tubes with a defined media, and measuring intracellular components, such as RNA, DNA, and proteins. These data can be used to better assess , and other rate law parameters. In short, one may model the dynamic behavior of microbes.
Microbial (community) growth is controlled by both extrinsic and intrinsic characteristics. The former refers to the access of microbes to substrates, including processes such as solutes diffusion and advection, or microbial sorption. The latter refer to microbial traits such as substrate preference or metabolic capabilities [196].
The basic mass balance for a specific microbe, assumed immobile, simply results from Equation (17), which reads:

36
) where superindex refer to the different substrates that may sustain the growth of the specific microbe . The simplicity of this expression hides some of the complexities of microbial metabolism. For example, separation between the energy required to generate new biomass and that required for maintenance (e.g., respiration).
As an increasing number of factors are taken into account, the simplicity of the Michaelis-Menten terms has come under question. For example, the basic derivation of the Michaelis-Menten equation (31) was based on the assumption that substrate concentration was not a limiting factor, but this assumption must be revised. Aquifers are oligotrophic and we are starting to address situations where OC is low. A reverse Michaelis-Menten factor has been proposed to address these situations. Both Michaelis-Menten and reverse Michaelis-Menten factors are special cases of the "equilibrium chemical approximation" [197].
Another surprising feature of Monod-based equations is the linear growth behavior of biomass during early stages. One would expect a Malthusian exponential growth, possibly limited by some internal factor, which would lead to a logistic or S-shaped curve behavior. This behavior could be reproduced by increasing the pool of "factors" with a logistic factor that limits the maximum concentration (density) of biomass .
= 1 − 37) This term is needed to prevent a boundless increase in when neither substrate nor TEAs are limited. In fact, terms with this form were proposed by Kindred and Celia [198] and Schäfer et al. [181]. They limited biomass growth by the available pore space. However, as we shall see in Section 0, we view microorganisms as species belonging to a specific phase (the biofilm), where cannot be limited by the porosity. We conjecture that will be described by some internal, microbial intrinsic, factor. This, together with the fact that different microbial groups may live together in a biofilm, suggests an expression like is the maximum total microbial concentration. Equation (38) also contains factors to limit growth by the availability of substrate and electron acceptors, similar to Equation (35), which could be added or not. Inhibition factors are not necessary, because the inhibition by some TEA should come implicitly from the larger growth of biomass of the corresponding TEAP. Nevertheless, rate laws with inhibition factors in combination with biomass can be found in the literature [199,200]. Xu [201] and Tang [197] studied Monod ( = ), logistic ( = ) and hybrid rate laws ( = ) and argue in favor of the hybrid rate laws.
Many other factors must be acknowledged when simulating microbial growth. Stolpovsky et al. [202] account for the possibility of biomass entering into a "dormant" state during starvation periods, and then reactivate when substrate concentrations recover. This is modeled for a batch system which does not consider transport. Since energy required for reactivation is lower than for synthesis, this strategy may be advantageous. Dormant species still consume OC for maintenance, which they get from death cells, but at a very slow rate. Dormancy is especially relevant in soils, where it may be provoked not only by OC, but also by water scarcity [203]. Actually, while microbial communities may feed on different substrates, a lag period may be needed for the microbes to "get used to" (i.e., to develop the adequate enzymes for) new substrates when they change. Other issues include biomass mobility and pathogen predator functions.
Reactive transport models that consider biomass, as in equation (36) and (38), are mostly applied to constructed wetlands, which are an alternative for treating waste water. Reviews exist on this type of modeling [204,205]. Arguably, they contain the most complex and complete biochemical systems, used for porous media up to now. An example is CMW1 (Constructed Wetland Model No1) [199], a biochemical system designed for constructed wetlands. It consists of species (inorganic electron acceptors, slowly and readily degradable organic carbon, biomass of bacteria for different electron accepting reactions) and stoichiometry and kinetic rate laws and parameters for hydrolysis of slowly to readily degradable organic carbon, and cell synthesis and decay of the biomass of the different bacteria. Most rate laws contain Monod and inhibition factors, such as in equation (35), but, they do not contain logistic factors, such as in equation (38). CMW1 has been coupled to flow and transport codes, such as Retraso [206,207] and Hydrus [208,209]. There are a few modeling studies in marine and lake sediments that also consider biomass (e.g., [200]), although most use the traditional rate laws of Section 0. They more or less follow similar rate laws as for the constructed wetlands with Monod and inhibition factors but without logistic factors.
We feel that the lack of consensus on empirical models of microbial metabolism, and the large number of controlling factors (and parameters) reflect the developing stage of these models. It also explains the wide use of the simplified models of Section 0. But we also take it as indicative of the need for further research.

Genomics based metabolic models
The fact that full genome of microbial populations is becoming increasingly available, makes it irresistible to try and model the growth of microbial communities. Ideally, it should be possible to model the metabolism from the genetic information available from sequencing, if one knew the functions of the genes and how they evolve. We do not know them. But, even if we did, the problem is so complex (the large number of genes, of reaction paths and possible interactions) as to become non-computable, not to mention the characterization of the involved parameters (rate constants, inhibition factors, etc.). Therefore, simplifications are needed. The most immediate is to use sequencing data as indicators of the microbial community behavior much in the way CFUs (Colony-Forming Units) have been used to calibrate models of microbial growth Reed et al. [210] and Louca et al. [211] went further and applied a complex 'gene-centric' approach using marker genes as experimental ref- erences, the production of which was linked to biogeochemical reactions in the ocean. Pagel et al [212] linked genetic information on abundances of bacteria, fungi and pesticide degrading specific microorganisms to the biogeochemical dynamics of DOC, and a pesticide in soil by multiobjective calibration. Instead, Scheibe and coworkers propose an alternative approach: constraint-based modelling [213][214][215][216]. Genome sequence data are used to define a metabolic reaction network within the cells. This may consist of several hundreds of metabolites and reactions and exchange rates between the cell and its environment [217] and is represented, similar to Equation (17), by a stoichiometric matrix and metabolic and exchange rates, but neglecting the transport terms and assuming steady state within the cells, yielding = 0. Given concentrations of chemical species of the cell's environment, the rates, , are optimized for maximum biomass production, ATP (Adenosine Triphosphate) synthesis, or minimize nutrient uptake, subject to thermodynamic constraints, upper and lower bounds of the rates and mass balances of metabolites (i.e., = 0). In essence, exchange rates between the cell and its environment are calculated from concentrations outside the cell, using relatively complex models of metabolic reaction network within the cells. Coupling this kind of models to RT has been restricted to single microbial species metal redox problems [217], which again suggests that much remains to be done.

Biofilm Modeling
Microorganisms are mostly concentrated in the biofilm, a biological substance, that cover grain surfaces and/or span through the pores in web-like structures. Biofilms mainly consist of water. Of the dry mass cells of microorganisms account for less than 10% and more than 90% is composed of extracellular polymeric substances or EPS [218]. EPS is a fibrous gel-type matrix composed of polysaccharides, lipids, proteins, and DNA materials excreted by the microbes [219,220]. The importance of the biofilm lays in the fact that almost all biochemical reactions, such as organic carbon decay, described in Section 0, take place in the biofilm. As the hydraulic conductivity of biofilm is very low [221], exchange of solutes between pore water and biofilm is controlled by diffusion. As a result, rates of biochemical reactions may become affected by the size and thickness of the biofilm [222].
Thus, biofilms can be seen as a particular form of heterogeneity. Hence, for its modeling, one can adopt a similar approach as for models for heterogeneity, such as the MRMT (multirate mass transfer), explained in Section 0. The biofilm is represented by the immobile zones, which can be of different size and thickness. Alternatively, one can think of the biofilm as another immobile phase, but different from minerals or surface phases. Exchange of mass between the pore water and the biofilm becomes a "chemical" reaction, for which a kinetic rate law can be defined. Either way, the need for the exchange between biofilm and porewater has been acknowledged for a long time [181]. In fact, its inclusion in continuum formulations has been formalized mathematically [223]. The biofilm-porewater exchange is assumed proportional to the difference between the concentration in the biofilm and porewater ( ( − ), where is an exchange coefficient and and are the concentrations in the porewater and biofilm), which is identical to the exchange rate between mobile and immobile zone for MRMT in Figure 10e.
As the biofilm grows, it can alter the properties of the porous medium. This is particularly important for cases with much organic carbon decay, such as constructed wetlands and managed aquifer recharge. The effect on porosity and permeability (or clogging) has received a lot of attention (e.g., [224][225][226][227]). Also, its effect on the retention curve has been studied and modeled [228][229][230]. Maybe the most interesting is the effect of biofilm growth on transport properties, such as, dispersivity and parameters for the porewater-biofilm exchange, which in their turn effect biochemical reactions and biofilm growth [231]. Experimental results have demonstrated that, indeed, biofilm growth induces heterogeneities that affect the transport in porous media [232,233]. Moreover, it forms preferential flow paths and stagnation zones [230]. Therefore, to understand the reactive transport related to biofilm growth, we need challenging reactive transport models that are capable of simulating dynamic biofilm growth and the increase of the heterogeneity of porous media simultaneously. Beyond the conceptual difficulties introduced in Sections 0 and 0, RT modeling is challenging because of the large number of non-linear coupled PDEs (Equations (17), (18) and (19)) and AEs (Equation (12) or (13)) that have to be solved. The challenge is worsened by the need of high spatial resolution needed to address the ever present chemical and physical heterogeneity, discussed in Section 0. Since analytical solutions are restricted to very simple conditions (we discussed a few of them in Section 0), numerical methods are required. Traditionally, solution methods were divided in two groups: The Sequential Iteration Approach (SIA), based on Picard method, and the Direct Substitution Approach (DSA), based on Newton-Raphson. Here, we discuss three other possible groups: methods based on generic PDE solvers, the Water Mixing Approach (WMA), and RW methods. While the SIA and DSA groups require full coupling and iterations on both transport and chemistry, the WMA and RW only require iterating on the chemical step (Error! Reference source not found.). All these methods accept numerous variants. Here we only describe the main features. Specific features are discussed in detail in methodological reviews [36,234].
The SIA consists of two steps at each iteration: (1) solve the transport equations of the aqueous component ( ) and (2) compute the sink source term R, which takes into account the concentrations of solid species and kinetic reactions of a previous iteration (see Error! Reference source not found.a).
For the first step we rewrite the transport Equations (17), (18) and (19) as: where superscript n refers to iteration number. By using the source-sink term ( ) of the previous iteration, we effectively decouple the numerical solution of the various components. Transport of each component can be treated as a "conservative" solute with an additional term for chemical reactions. Therefore, transport is solved numerically by any appropriate method (typically Eulerian finite element or differences, with various techniques to control numerical dispersion and instability). One of the problems with SIA is that splitting transport from chemistry requires using the same components throughout the domain, which is not convenient in problems with several mineral zones (recall Conceptual Issue No. 5).
The second step of SIA consists of calculating the concentrations of all chemical species ( ) from the aqueous components (ua) to evaluate the rates of kinetic reactions. For this we impose: which has to be solved together with the mass action law for equilibrium reactions, equations (12) or (13). This leads to a standard speciation calculation from components, briefly discussed in Section 0. A second problem with SIA is that these equations are non-linear, which requires an additional loop of iterations within the loop of the transport equations (39) (Error! Reference source not found.a). These calculations are performed for each node separately, which facilitates parallelization. After calculation of all concentrations, the sink-source term ( ) can be updated easily using equation(40).
The advantage of SIA is the relative ease of building a reactive transport code by coupling a code for speciation to a code for conservative transport to which a chemical source-sink term (f) has to be added. This may explain the large number of codes that use SIA. PHREEQC can solve 1D reactive transport problems [235] or be used as a speciation library to solve the chemical calculations [236]. Some examples of transport codes, that have been coupled to PHREEQC, are MODFLOW/MT3D [237], COMSOL [238], or PHWAT [239]. HYDRUS-1D [240] also solves 1D reactive transport (RT) problems along with heat transport and density dependent flow cases. HYDROGEOCHEM [241] solves coupled Thermo-Hydro-Mechanical-Chemical (THMC) 3D models. TOUGHREACT is another THMC and multiphase flow code [242], largely based on CORE-2D [243]. Finally, Van Der Lee [244] coupled the transport code HYTEC with speciation code CHESS. Some codes avoid iteration between the two steps, a methodology termed Sequential Non-Iteration Approach (SNIA). Examples of codes using SNIA are the Pore Network code coupled to PHREEQC, PoreFlow [245] and OpenGeoSys [246], which is also a THMC code.
The DSA consists of doing one single step (see Error! Reference source not found.b) per iteration by substituting chemical equations into the transport equations and applying method of the New- ton-Raphson. Concentrations of secondary species can be written as a function of primary concentrations by writing the mass action law (Equation (13)) as log = * log + log * + log ( 42) which can be substituted into the transport Equations (17), (18) and (19). The size of the system to be solved at each iteration is equal to the number of components times the number of nodes or cells, much larger than that of the SIA. One advantage of the DSA is that it facilitates the elimination of constant activity species, even with different mineral zones each with its own reduced component matrix [61]. In that case, in discretized form through finite differences or finite elements, each node or cell can belong to a different mineral zone. This implies that the number of degrees of freedom will be different for every node. As a result, more programming is involved. In fact, since the coupled problem is (highly) non-linear, DSA requires using Newton-Raphson, which is considered good for convergence, but implies an added programming difficulty: need to code not only speciation calculations, but also their derivatives. However, there is a high degree of freedom in the choice of unknowns. In particular, primary species can be the unknowns, which facilitates speciation calculations (they are restricted to updating activity coefficients, since is an explicit function of in Equation (42)). Therefore, comparison of methods is not straight-forwards. SIA is simpler to code and leads to smaller systems of equations than DSA, but DSA should converge faster. Initial comparisons favored SIA [36]. But more recent comparisons favor DSA in terms of performance [59,234,247]. Still, the main advantage of DSA over SIA is that it is more robust than the SIA, in the sense that its convergence is less sensitive to time step size, and that it simulates better sharp fronts [37].
Some examples of DSA codes are MIN3P [248], PFLOTRAN [249] and RETRASO [45]. CRUNCHFLOW [250] is a code that can use both SIA and DSA. Some codes have been adapted to high performance computing, such as MIN3P-HPC [251] or ParCRUNCHFLOW [252]. They speed the performances by implementing parallelization programming and compute the models on multiple processors at the same time.
It could be argued that RT modelers do not need to be experts in the solution of PDEs and, instead, adopt general PDE solvers. In particular, coding the Jacobian is difficult and error-prone. Therefore, some researchers have opted for generic, including Jacobian-free, solvers [64,253,254]. In particular, the Method Of Lines [255], a general method for solving coupled PDEs, has been used with notable success to solve RT problems [256,257].
The fact that transport equations are virtually identical for all mobile species prompted Soler-Sagarra [258]to propose an alternative concept, the Water Mixing Approach (WMA). The basic idea is to understand transport as water advection and mixing (being solute just an attribute of this water) instead of transport of each individual solute. The concept is very convenient because any transport master equation (and, thus, any form of the transport equation with any solution method) can be viewed as a mixing equation [258]. The concept can be extended to reactive transport by simply adding the reaction term to the mixing equation, which leads to rewrite Equation (17) as: where is the mixing ratio, that is the volume (or mass, depending on the units of c) of water (actually fluid) that transits from cell (or node, or particle, etc.) to during the time step, expressed as a fraction of the water volume (or mass) in cell ; is the number of cells connected to , and is the vector of reactive contributions to all species in the cell. This equation is complemented by a similar one (only without the mixing term) for immobile species. Note that this equation is similar to the one used by Pelizardi et al. [259] to compute end-member mixing ratios for the interpretation of concentration data. Here, the goal is the opposite; given the mixing ratios, compute reaction amounts ∆ and the updated concentrations. This calculation is relatively easy to implement in existing speciation codes. Soler-Sagarra et al. [258] details the procedure and implemented it on CHEPRO [67]. Obviously, the main advantage of the WMA lies in decoupling transport (restricted to the calculation of mixing ratios) and chemical calculations, so that no transport iterations are needed (see Error! Reference source not found.c). The WMA method looks promising but requires further testing. The vast majority of SIA and DSA codes presented above use Eulerian methods, such as Finite Differences or the Finite Element Method, to model transport by using the classic ADE. However, the proper modeling of transport is crucial for RT problems, as discussed in Section 0. Lagrangian formulations (Equation (43)) are more adequate for advection dominated problems. This can be done with mixed Eulerian-Lagrangian methods using variants of the method of characteristics [260][261][262][263], Localized Adjoint Methods [264,265] or stream-based grids [266]. Most of these methods risk suffering fluctuations or numerical dispersion, which motivates purely Lagrangian methods (Error! Reference source not found.b and c). Error! Reference source not found. and Error! Reference source not found. were obtained with a Lagrangian method (isochronal grids) coupled to the WMA. Error! Reference source not found. illustrates one difficulty to be expected with Lagrangian methods. They usually assume solutes to move at the velocity of water, which is not true when sorption reactions take place (recall Section 0 and Conceptual Issue No. 6). If the actual solute velocity is different from the one assumed by the Lagrangian method, some numerical dispersion may occur (Error! Reference source not found.).
The most widely used Lagrangian methods (and the exclusive ones in recent transport research) belong to the family of RW methods. As mentioned in Section 0, RW methods do not simulate mixing or, worse, concentrations explicitly, which makes it hard to assess kinetic reactions and direct application of RW to RT. Some researchers seek to overcome this problem by getting around the need for concentrations and instead trying to find the probability of particles encountering each other to interact as a function of their distribution in space and time [267][268][269][270][271]. The probability of particle collisions and interactions, and transformations may be made to depend on the scale [272]. This implies treating particles as individual species, and facilitates assigning each particle its specific transport properties, thus getting around the above problem with sorption [273]. The approach has been successful in reproducing bi-molecular reactions in anomalous transport settings [274][275][276], including biomass growth [272]. While it may appear unrealistic for chemical systems with numerous species and reactions, advances in computational capabilities make it an active field of research.
In most of the above methods, particles represent molecules, so that they get annihilated upon reaction, which is not practical, as it requires a huge number of particles. Therefore, many RW researchers have opted by assigning them an initial mass, so that only the reaction fraction is used upon reaction [277] assign each particle a fraction of the total mass of a given species, Still, the Kernels introduced in Section 0 as part of the method to recover concentration fields from particles (Error! Reference source not found.g) tend to be used to define a support volume to facilitate the evaluation of particles interaction [113,278,279]. A final twist on the method is to extend the concept to total concentrations [280], in the hope that it allows describing transport accurately (as discussed in Sections 0 and 0, this hope is vain because each species in the component will move at its own velocity and because mobility is non local, i.e., it depends on other particles). Still, this last direction converges on WMA concepts, as it becomes conceptually similar to the Multi-Advective Water Mixing Approach [143], and suggests that the RW community may eventually blend with the rest of the RT community.

Discussion and conclusion.
This review consists of two parts. In the first part (Sections 0 and 0) we have revised RT concepts of traditional RT, largely "geochemical", models. In the second, we have reviewed emerging developments on transport, biochemistry, and RT solution methods.
The first part suggests that RT has reached a significant level of maturity, which means that basic concepts are broadly accepted by the community. We have summarized our views on these basic concepts as a set of conceptual issues (Table 1). Still, broad acceptance does not mean unanimous acceptance. In fact, the first and probably least controversial issue (multidisciplinarity) suggests diverse views. Therefore, our goal in setting these issues is two-fold. On the one hand, given the diversity of views, we hope to promote a debate on what it is that we agree on. On the other hand, we hope to facilitate the integration of newcomers into RT. Coming from hydrology, it took us some time to understand the importance of these "basic" concepts. We hope that confronting them upfront will save some of this time to others. Some of the issues may sound trivial. It is obvious that numerous disciplines are involved in RT (C.I. 1) or that appropriate reactions and species must be included (C.I. 2). We have listed them to point that they deserve further thought. It is not sufficient to bring together different specialists, you need them to work together (interdisciplinarity) and to understand each other. For example, when, why, and which aqueous complexes need to be included in a RT model may not be evident to a solute transport specialist. As a result, some published papers appear to ignore it. Reversely, the fact that RW methods are used at the pore scale suggests that diffusion dominance at the pore scale is not widely acknowledged yet. In fact, the scale dependence of RT processes is broadly mentioned but not so broadly accounted for.
The fact that the field appears to have reached a significant level of maturity does not imply either that RT is well understood. We have reviewed the literature and pointed to issues where we feel much remains to be done. The field is extremely active in several directions.
There are several formulations for solute transport. The ADE remains, by far, the most widely used, but the fact that it equates spreading and mixing, makes it a poor representation. Numerous alternatives have been proposed to improve on this. They improve on mixing and overcome many of the limitations of the transport equation, such as scale dependence and tailing. Some of these methods are hard to apply directly to RT. MRMT, which is the simplest one, seems to be the most suitable, because reactions terms can be added in a straightforward way. We found very few cases where RT is modeled in combination with effective non-fickian transport formulations. However, we feel that these formulations are particularly relevant for RT, because mixing (not spreading) enhances reaction rates, especially of equilibrium reactions.
Regarding the formulation of reactions, the mass action law is the standard for equilibrium reactions. It is used for calculations in conventional codes for geochemical calculations. However, for RT the use of artificial neural networks may become appealing, because of the high number of such calculations that are required.
RT models are performed by combining transport and reaction formulations and require the numerical solution of large number of nonlinear equations. The two traditional methods are DSA (or Newton-Raphson) and SIA (or operator splitting). DSA is considered the most effective, but its complexity discourages broad application and SIA remains the most widely used. A new method is the WMA (Water Mixing Approach), which may become the method of choice, because of its accuracy and simplicity. Biochemical reactions are widely modeled with Monod-like kinetic laws. However, it is becoming increasingly apparent that they do not suffice and that one also should consider microbial communities, that can grow and live in biofilms. Genome based formulations sound like the right approach. However, their complexity and the difficulties inherent to assigning functionality to geochemical reactions are far from resolved.
As a final conclusion, it can be argued that the field is mature. But one should not confuse maturity with uniformity or, much less, senescence.

No. Conceptual issue 1
RT is multidisciplinary and requires interdisciplinarity RT requires expertise on groundwater flow, hydrogeochemistry, biology, transport in porous heterogeneous media, numerical solution. 2 Diffusion controls microscales; Advection controls macroscales This follows from the characteristic times of diffusion and advection and typical values for diffusion coefficients and water velocity (section 3.2). 3 It is important to select the appropriate set of species and reactions Selection of abiotic geochemical reactions may be done relatively easy by using speciation codes (section 2.2). Selection of biological reaction is more difficult, though at least they should be thermodynamically favorable (section 5.1). 4 Equilibrium reactions reduce the size of RT problems The number of components, or the number of transport equations to be solved, equals the number of species minus the number of equilibrium reactions (section 2.2 and 2.3). 5 Variations in mineral zone, leads to sharp fronts in aqueous chemistry The set of equilibrium minerals present in some domain (or mineral zone) determines the set of mass action laws that are imposed and, thus, effects the chemistry of the water (section 3.3). 6 Each species travels at a different velocity Solid species do not move at all, whereas aqueous species do. Species that may adsorb are retarded with respect to the water (section 3.1). Moreover, it makes RT more difficult to solve (section 6), because it produces a usually non-linear relation between total concentration of a component ( ) and its aqueous part ( ) in the transport equation (20). 7 Kinetics are important and complex Equilibrium constants are well-known, the gist of RT often lies on the adopted kinetic laws and their parameters (section 2.2 and 5). Since kinetic rates are often variable, beware of the Damkhöler number (section 5.2). 8 The rate of fast reactions is controlled by mixing and access to reaction sites For cases with only aqueous and mineral reaction in equilibrium reaction rates can be expressed by equation (30). Solutes are displaced by advection, spread by dispersion, diluted by mixing The ADE is a poor representation of transport as it does not distinguish dispersion from mixing. This has been acknowledged for conservative transport, but is even more relevant for RT 10 Biochemical reactions involve microorganisms Microorganisms mean complex kinetic rate laws. Moreover, they can create biofilms that alter heterogeneity and flow and transport characteristics of the porous media.