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 scientiﬁc 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 advancement along the current frontier: biochemical reactions. To this end, we review the basic equations to indicate 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, difﬁculties arise for kinetic reactions, which is especially true for biochemistry, where numerous challenges are open (how to represent microbial communities, impact of genomics, effect of bioﬁlms on ﬂow and transport, etc.). We conclude with the basic eleven conceptual issues that we consider fundamental for any conceptually sound RT effort.


Introduction
The term reactive transport (RT) describes the simultaneous occurrence of (bio-and/or geo-) chemical reactions and transport of biochemical species (Figure 1).RT is important to understand numerous Earth Sciences problems, ranging from engineering applications (geologic carbon storage, enhanced geothermal systems, or groundwater remediation, etc.), as well as geological studies (such as sediment diagenesis, seismic activities, hydrothermal ore deposition, etc.) to watershed or global issues [1].We are specifically motivated by contaminant removal during artificial recharge [2,3], which is driven by biochemical reactions.
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 (Figure 2): 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 [4] and Seawater Intrusion [5].The third one refers to non-Fickian transport [6].The fourth one describes a specific surface reaction [7].The four are relevant but not really RT papers by themselves.It is clear that this kind of search Energies 2022, 15, 925 2 of 47 does not portray an accurate picture of the field.In fact, the list does not include the early works of Rubin and James [8], Grove and Wood [9] or Lichtner [10], who set the basis for the current paradigm of RT.However, 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. [11]).
Energies 2022, 15, x FOR PEER REVIEW 2 of 48 to non-Fickian transport [6].The fourth one describes a specific surface reaction [7].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 [8], Grove and Wood [9] or Lichtner [10], who set the basis for the current paradigm of RT.However, 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. [11]).
Figure 1.Reactive transport (RT) describes the transport of solutes (displacement by advection, spreading by dispersion, and dilution by mixing) and the reactions (red arrow) among solutes and with immobile phases.The mathematical nature of RT depends on the type of reactions and species involved [12].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 [13].Some are centered on problems related to CO storage [14,15], which generates acidity, so it is somewhat related to acid mine drainage [16], more formally revised by Amos et al. [17].Nuclear waste disposal has motivated reviews on uranium mobility [18], validity of existing models [19], cement-clay interaction [20] or oxygen transport at depth, which motivated the thorough review of MacQuarrie and Mayer [21] on RT in fractured rock.Global issues motivate the reviews of Dignac et al. [22] and Li et al. [23].Arndt at al. [24] and Regnier et al. [25] concentrate on degradation of organic carbon in marine sediments.Reactive transport (RT) describes the transport of solutes (displacement by advection, spreading by dispersion, and dilution by mixing) and the reactions (red arrow) among solutes and with immobile phases.The mathematical nature of RT depends on the type of reactions and species involved [12].
Energies 2022, 15, x FOR PEER REVIEW 2 of 48 to non-Fickian transport [6].The fourth one describes a specific surface reaction [7].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 [8], Grove and Wood [9] or Lichtner [10], who set the basis for the current paradigm of RT.However, 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. [11]).
Figure 1.Reactive transport (RT) describes the transport of solutes (displacement by advection, spreading by dispersion, and dilution by mixing) and the reactions (red arrow) among solutes and with immobile phases.The mathematical nature of RT depends on the type of reactions and species involved [12].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 [13].Some are centered on problems related to CO storage [14,15], which generates acidity, so it is somewhat related to acid mine drainage [16], more formally revised by Amos et al. [17].Nuclear waste disposal has motivated reviews on uranium mobility [18], validity of existing models [19], cement-clay interaction [20] or oxygen transport at depth, which motivated the thorough review of MacQuarrie and Mayer [21] on RT in fractured rock.Global issues motivate the reviews of Dignac et al. [22] and Li et al. [23].Arndt at al. [24] and Regnier et al. [25] concentrate on degradation of organic carbon in marine sediments.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 [13].Some are centered on problems related to CO 2 storage [14,15], which generates acidity, so it is somewhat related to acid mine drainage [16], more formally revised by Amos et al. [17].Nuclear waste disposal has motivated reviews on uranium mobility [18], validity of existing models [19], cement-clay interaction [20] or oxygen transport at depth, which motivated the thorough review of MacQuarrie and Mayer [21] on RT in fractured rock.Global issues motivate the reviews of Dignac et al. [22] and Li et al. [23].Arndt at al. [24] and Regnier et al. [25] concentrate on degradation of organic carbon in marine sediments.Wetlands are reviewed by Langergraber [26].Energy related subsurface problems are starting to consider RT [27][28][29][30].
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 includes quantitatively oriented hydrologists, engineers and physicists.The latter might sound as though it is 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 working together is not easy.Multidisciplinarity is often hallowed, but interdisciplinary work requires time and a determined effort in developing shared language and understanding.Worse, in practice it is penalized by evaluation panels.Methods cannot be shared by the various specialists, but concepts can.
The goal of this paper is to frame and review basic RT conceptual issues, which we highlight, with special emphasis on the fast-growing biochemical field, both because of our own interest and because we feel the need for shared concepts is acute.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 [41][42][43] nor on the numerous applications 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 2 may be a bit too mathematical for some, but we request the reader to try and concentrate on the concepts).These equations allow us to frame simple reactive transport paradigms (Section 3) that illustrate what we consider the main conceptual issues of RT.In turn, the importance of properly representing mixing motivates Section 4, where we analyze transport equations other than the RT. Biochemical concepts are reviewed in Section 5 and solution methods in Section 6.

Conceptual Issue
No. 1 is that RT is multidisciplinary and requires interdisciplinarity.It is clear that many important questions can only be overcome by overlapping knowledge from different science branches, which requires working together and understanding each other.Thus, we start by reviewing the basic equations and concepts for transport and reactions.

Transport
Transport requires knowing the time and space variability of fluid fluxes.These are obtained from the flow equation, usually solved prior to transport, which requires establishing conservation of fluid mass and momentum.Fluid mass conservation is expressed as (see notation table at the end of the paper): where ρ M w /L 3 w is fluid density, φ [-] is volumetric water content (porosity in saturated media, and φ = 1 in pore scale simulations or when flow is solved in a fluid continuum, e.g., ocean), q L 3 w /L 2 m /T is the volumetric fluid flux (usually called Darcy velocity in porous media, but only a true fluid velocity when φ = 1), f w L 3  w /L 3 m /T represents volumetric fluid sink and sources of inflowing external water with density ρ e M w /L 3  w , which is equal to the resident fluid density ρ for sinks (i.e., when f w < 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 w for water, actually fluid in general, or m 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., [36]), or further simplified to Darcy's Law when the laminar flow is slow (i.e., inertial effects can be neglected) and the porous medium is treated as a continuum [44].
Transport is usually solved with the advection-dispersion equation (ADE, sometimes ADRE, with the R standing for reaction), a partial differential equation (PDE) that expresses the mass balance of each aqueous species when transport is driven by advection and dispersion.Beyond the mathematical complexities, Figure 1 implies that concentration changes due to the following processes: 1.
Advection, that is, the displacement of solutes dragged by fluid flow.Note that, unless flow is solved at the pore scale, advection is represented with a mean fluid flux (mean over a volume that depends on the level of detail with which flow has been solved).

2.
Dispersion includes two processes: molecular diffusion (i.e., Brownian motion causes a solute mass flux away from high concentration regions, given by Fick's Law) and velocity fluctuations with respect to the mean fluid flux, which causes the solute to spread (also written by a Fick's-like Law).We will revisit this term in Section 4.

3.
External sources, solute mass that leaves/enters the medium with fluid sink and sources, f w .
The ADE results from explicitly expressing mass conservation (per unit volume of medium and unit time) for the above four processes, which yields: where c i [M sol /M w ] is the concentration of species i (aqueous or immobile depending on the subscript, a or s, respectively), D[L 2 /T] is the hydrodynamic dispersion tensor (D pq = φD m δ pq + α T qδ pq + (α L − α T )q p q q /q, where q is the modulus of flux, q = |q|, and q p its p-th component; D mp = φD m is the effective molecular diffusion, c ex.i is the concentration of species i in fluid sources (again, c ex,i = c i when f w < 0) and R i M sol /L 3 m /T includes the contributions of chemical reactions to the mass balance of species i, expressed per unit volume of porous medium.This equation is somewhat cumbersome (in part because of the units we have chosen).In fact, it can be written in various forms depending on the context.For example, it is frequent to write D in terms of velocities instead of fluxes for historical reasons, but this is inconvenient in porous media, where velocity is poorly defined (it depends on the model, whereas water flux, given by Darcy's Law, is a robust concept).In addition, we have written concentration as mass fraction (e.g., ppm) to ensure that dispersion does not cause advection and that adding Equation (2) for all aqueous species yields the fluid mass balance (Equation (1), e.g., [45]).However, reactions usually express the solute mass [M sol ] in moles, so that c i can express moles of species i per unit mass of solution (still a mass fraction, for practical purposes).When fluid density is taken as constant then it is valid to write transport in terms of ρc i (molarity).However, chemical databases usually require working with molality (moles solute per kg of H 2 O, or m i = c i /ω w , where ω w is the mass fraction of H 2 O).
The conservative form of the ADE (Equation ( 2)) 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 c i , from Equation (2), and dividing by ρ, which yields: This advective form is conceptually rough because (1) it represents the "flow-equation[removed" [46] balance, but divided by density, and (2) the units of the mass balance are somewhat tricky M sol L 3 w L −3 m M −1 w T −1 , 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 ( 4) by φ, while neglecting its spatial variations for the dispersive term, which yields: where R w,i = R i /φρ is the contribution of reactions, but expressed per unit mass of fluid, D v is similar to D, but expressed in terms of velocity instead of fluid fluxes, and v = q/φ 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 Lagrangian formulations of transport (e.g., Random Walk methods): where Dc i /Dt = ∂c i /∂t + v∇c i is the material derivative (i.e., the rate of concentration change in a flowing water parcel).It is important to point that all these advective forms of the ADE (Equations ( 4)-( 6)) are equivalent to Equation (2).We have spelled them out to illustrate that the ADE can be found in numerous forms.In practice, however, advective forms require being extremely careful in the application of numerical methods.For example, very few authors use the reverse time weighting [47] for transient flow cases, which may lead to severe mass balance errors.Still, for simplicity, we write: where 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 Figure 1.The figure shows that the basic behavior of the ADE is defined by the dispersion length L D = √ Dt, which characterizes the spread over which concentration tends to homogenize, and the advection length L A = qt/φ = vt, distance traveled by conservative solutes.Figure 3 displays the time evolution of these two lengths for different values of D, q, 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., t c = φD/q 2 ), which leads to a characteristic length L c = D/q and the dimensionless ADE: where is a characteristic concentration that depends on the problem, often a boundary or initial concentration, here c c = c ex,i ), f D = f w L c /q (i.e., the fraction of water flux, q, that enters the medium from internal sources, f w , over L c ), and R Di = R i L c /qc c (i.e., the fraction of advective solute flux, qc c , produced by reactions over L c ).This dimensionless form is unusual.The most frequent choice is to select L c as a length of interest and then define the Peclet number as Pe = L int /L c , which appears in front of the advection term.This way, Pe > 1 is the same as saying L int > L c .We prefer using Equation ( 8) because it makes apparent that the Peclet number is not an attribute of the ADE, but of the length of interest, which makes Pe ambiguous 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.where  =   ⁄ , ∇ ≡   ⁄ ,  =   ⁄ ,  =   ⁄ ( is a characteristic concentration that depends on the problem, often a boundary or initial concentration, here  =  , ),  =   / (i.e., the fraction of water flux,  , that enters the medium from internal sources,  , over  ), and  =   / (i.e., the fraction of advective solute flux,  , produced by reactions over  ).This dimensionless form is unusual.The most frequent choice is 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 ( 8) because it makes apparent that the Peclet number is not an attribute of the ADE, but of the length of interest, which makes  ambiguous 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; the Peclet number is ambiguous.Other dimensional forms are possible and, especially, many other transport equations can be chosen.We analyze them in Section 4.

Reactions
The term  represents 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 it (reactions include all biogeochemical processes from mineral dissolution to hormones production).However, from an 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: where  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: Other dimensional forms are possible and, especially, many other transport equations can be chosen.We analyze them in Section 4.

Reactions
The term R i represents the contributions of all reactions to the mass balance of species i, expressed per unit volume of porous medium.This section is devoted to the evaluation of R i .Stated this way, it may sound irreverent, given the complexity of the concepts behind it (reactions include all biogeochemical processes from mineral dissolution to hormones production).However, from an RT perspective, evaluating R i 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: where Nm i is the chemical formula of the i-th species, N sp is the number of species, N r is the number of reactions, and S ij is the stoichiometric coefficient of species i in reaction j.
Then, R i is given by: Energies 2022, 15, 925 7 of 47 where r j is the reaction rate (i.e., moles of reactants that evolve into products, per unit time and per unit volume of medium).S•Nm = 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: R = S t r (12) where S is the N r × N sp stoichiometric matrix, whose rows contain the stoichiometric coefficients in each reaction, and r a vector of size N r , containing rates of all reactions, and Nm a vector of size N sp , containing the chemical formula of all species.To illustrate these concepts, let us consider a simple calcite dissolution problem, which one might define as consisting of three species Nm t = CaCO 3(s) , CO 2− 3 , Ca 2+ and one reaction (CaCO 3(s) − CO 2− 3 − Ca 2+ = 0).Then, the stoichiometric coefficients would be S = (1, −1, −1), so that R 1 = r 1 , R 2 = −r 1 , and R 3 = −r 1 .Note, however, that in most groundwater problems, this definition would be insufficient because most of the carbonate evolves to bicarbonate (HCO − 3 ) under neutral pH conditions, so that both HCO − 3 and H + would have to be included in the set of species and a new reaction (HCO − 3 − CO 2− 3 − H + = 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 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 N sp , and reactions as vectors of the changes in concentrations due to chemical reactions.This view is extensively discussed by Rubin [12] and Friedly and Rubin [48], who attribute it to Prigogine and Defay [49].
For the purpose of computing reaction rates, r, it is convenient to distinguish between kinetic ("slow" in Figure 1) and equilibrium ("fast" in Figure 1) 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 r k = r k (c), where we have introduced the subindex k to indicate that, in general, only a subset of all reactions are kinetic.r k (c) are empirical functions.Lasaga [50] covers them extensively for inorganic reactions.Biochemical reactions are reviewed in Section 5.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: where S e consists of the N e rows of S corresponding to the N e equilibrium reactions, K e is the vector of their equilibrium constants and a 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 databases 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 a i = γ i m i , where γ i is its "activity coefficient".For moderate salinities, the molality is very close to the molarity (i.e., a i ≈ γ i c i ).
In general, however, one would need to write a i = γ i c i /ω w .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 2.1), the activity of adsorbed species becomes a s,i = c si ρφ/(ρ s θ s ), where ρ s and θ s are the density and volume fraction of the solid.For cation exchange the activity is an equivalent fraction, that is, a i = z i c i ρφ/(ρ s θ s CEC), where z i is the charge of the solute and CEC 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 [40].As we see in Section 3.1, model representation of sorption can range from very simple to very complex (Appelo and Postma, [51], Chapter 6).For minerals, one has to understand that the solid comprises 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 numerous species (pure minerals, constant pressure gases, ideal H 2 O) display constant activity.Therefore, from the chemical point of view, they are not unknowns.The number of unknowns in Equation ( 13) is N sp − N c , where N c 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 databases.The selection of the database is a topic in itself.The most widely used databases for geochemical modeling are EQ3/EQ6 [52] and PHREEQE [53].Specific databases are often available for specific purposes.For example, [54] review data for CO2 storage.There are databases in the framework of cement/clay interactions [55] and the NEA (Nuclear Energy Agency) has produced thermochemical databases intended for the radioactive waste management community [56].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 performed by means of a good preprocessor).
Databases are also available for kinetic rate laws, especially for dissolution of minerals [57,58].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 [59].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(accessed on 15 November 2021)).Things get more complicated for biochemical reactions.PubChem (https://pubchem.ncbi.nlm.nih.gov) is a "popular chemical information resource" [60] (accessed on 15 November 2021).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.However, 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 and complexity of RT problems.The issue is illustrated in Figure 4 for a simple gypsum equilibrium example.The chemical system consists of three species (gypsum, Ca 2+ , and SO 2− 4 ).The stoichiometric matrix contains only one row describing the single reaction of the system (Figure 4a).Concentrations (Figure 4b) and activities (Figure 4c) 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 (a 1 = a Ca , in Figure 4c) allows easy calculation of the other (a 2 = a SO 4 ).
Energies 2022, 15, x FOR PEER REVIEW 9 of 48 equilibrium reactions reduce the size and complexity of RT problems.The issue is illustrated in Figure 4 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 (Figure 4a).Concentrations (Figure 4b) and activities (Figure 4c) 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 Figure 4c) allows easy calculation of the other ( =  ).The concepts of Figure 4 are general.Each equilibrium reaction 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 the 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  The concepts of Figure 4 are general.Each equilibrium reaction reduces the dimension of the concentrations space by 1.Therefore, the number of degrees of freedom is N s − N e and one can choose a subset of primary species as the basis of the subspace of "valid" concentrations.The remaining N e 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 S e and vector a into two parts (S e = (S e1 , S e2 ); a t = a t 1 , a t 2 ), where a 1 contains the activities of the (N sp − N e ) primary species and a 2 the remaining (N e ), secondary species, and where S e1 contains the stoichiometric coefficients for the primary and S e2 those for the secondary species.Substituting these definitions in Equation ( 13) and multiplying by S −1 e2 yields: where S * e1 = −(S e2 ) −1 S e1 and logK * = (S e2 ) −1 logK.In this way we obtain a new de- scription of the same chemical system, which is now defined by equilibrium constants K * instead of K and of stoichiometric matrix for the primary species S * e1 instead of S e1 .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 N sp − N e × N s kernel matrix (or null space) of S e .That is, US t e = 0.In vector space nomenclature, this matrix defines the orthogonal complement to the equilibrium subspace of valid concentrations.Matrix U is not unique.In fact, the set of possible component matrices is in itself a vector space of dimension N sp − N e .Hence, there are numerous methods to obtain it.Molins et al. (2004) describe some of them [61].In particular, they define approaches to ensure that U facilitates reactive transport computations.Still, the most widely used method is based on Gauss-Jordan elimination [62], which consists of defining the component matrix as: where I is the (N s − N e ) identity matrix.It is easy to check that, indeed, US t e = 0. Components are then defined as linear combinations of species: 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, c t = c a c s t , components may also contain an immobile part.This can be formalized by splitting U = U a U s , so that u = Uc = u a + u s , where u a = U a c, and u s = U s c.We illustrate these matrices by another simple example involving calcite (Cc) and gypsum (Gyp) dissolution precipitation in equilibrium, as shown in Figure 5.The chemical system consists of two reactions (N r = N e = 2) and five species (Cc, Gyp, Ca 2+ , SO 2− 4 , and CO 2− 3 , i.e., N sp = 5).Therefore, the chemical state can be defined in terms of three components or three primary species (the remaining species, secondary, can be easily obtained from Equation ( 14)).We choose Cc, Gyp, and Ca 2+ to be primary species, and SO 2− 4 and CO 2− 3 to be secondary.The corresponding stoichiometric matrix, S e , is shown in the left column of Figure 5.The choice of primary and secondary species is somewhat arbitrary.The only requirement is that S e2 must be non-singular.In this case, S e2 , which encompasses the last two columns of S e in Figure 5 (blue square), equals the identity matrix (S e1 is the green square).To obtain the components matrix, we first compute S * e1 , defined below Equation ( 14): Application of Equation (15) gives the following component matrix: Total concentrations of components, u = Uc (Equation ( 16)), can be read directly from Equation (18).The first one is total carbonate in the system (it will not change due to reactions, because carbonate must be either dissolved or in calcite).The same can be said about the second component (total sulfate).These two components contain aqueous (u a ) and solid immobile (u s ) parts.The third component (u 3 ) is conservative because gypsum dissolution and/or precipitation will change Ca 2+ and SO 2− 4 by the same amount (and similarly with calcite), so that their difference remains unchanged by these reactions.In short, the three components are indeed conservative.Actually, the problem can be further simplified (conceptually and numerically).As mentioned before, constant activity species are not really unknowns during chemical calculations and it is profitable to eliminate them.This is precisely what was done by Saaltink et al. [63].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.Another possibility is to choose the constant activity species to be primary, which is what we have done in the example of Equations ( 17) and (18), where calcite and gypsum are both constant activity and primary species.As the MALs do not depend on the activity of constant active species, we can decouple the calculation of the concentrations of these primary species and the transport equation of their corresponding components.In practice, the number of unknowns is reduced from  to  −  .
The use of constant activity species is discussed further in Figure 5.If Cc and Gyp are assumed constant activity (as discussed above, this is the standard assumption for pure minerals), they are not chemical unknowns, so that the stoichiometric matrix is reduced to a 2 × 3 matrix, and the components matrix becomes a 1 × 3 matrix (red square in Figure 5, and also in Figure 4).The component space is 1D, so that any species may serve as primary.In the full concentration space, the single component (Ca − +O − −O ) is represented as a plane orthogonal to the equilibrium (and, thus, to both the calcite and gypsum equilibrium subspaces) reactions.Reactions must stay in this plane (red and purple arrows in Figure 5) because they do not change components.
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  = ().However, 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 ( 16)).The solution of these equations can be difficult and may require special solutions methods [64][65][66].A number of codes are available to perform speciation calculations, such as PHREEQC [67], CHEPROO [68] and GEMSK [69].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 [70], which is possibly what the future holds for RT.Actually, the problem can be further simplified (conceptually and numerically).As mentioned before, constant activity species are not really unknowns during chemical calculations and it is profitable to eliminate them.This is precisely what was done by Saaltink et al. [63].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.Another possibility is to choose the constant activity species to be primary, which is what we have done in the example of Equations ( 17) and ( 18), where calcite and gypsum are both constant activity and primary species.As the MALs do not depend on the activity of constant active species, we can decouple the calculation of the concentrations of these primary species and the transport equation of their corresponding components.In practice, the number of unknowns is reduced from N sp to N sp − N c .

Coupled Transport and Reactions
The use of constant activity species is discussed further in Figure 5.If Cc and Gyp are assumed constant activity (as discussed above, this is the standard assumption for pure minerals), they are not chemical unknowns, so that the stoichiometric matrix is reduced to a 2 × 3 matrix, and the components matrix becomes a 1 × 3 matrix (red square in Figure 5, and also in Figure 4).The component space is 1D, so that any species may serve as primary.In the full concentration space, the single component 3 ) is represented as a plane orthogonal to the equilibrium (and, thus, to both the calcite and gypsum equilibrium subspaces) reactions.Reactions must stay in this plane (red and purple arrows in Figure 5) because they do not change components.
Speciation calculations are needed to obtain the vector of concentrations of all variable activity species from a 1 or c 1 .In the first case, Equation ( 14) yields directly a 2 .Concentrations are then obtained from the activity coefficients as c = a/γ.Iterations may still be needed because γ = γ(c).However, this dependence is weak.A more problematic case arises when speciation needs to be performed from components.In this case, the N s concentrations are obtained from the MAL of the N e equilibrium reactions (Equation ( 12)) and the N s − N e definitions of components (Equation ( 16)).The solution of these equations can be difficult and may require special solutions methods [64][65][66].A number of codes are available to perform speciation calculations, such as PHREEQC [67], CHEPROO [68] and GEMSK [69].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 [70], 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, explained in Section 2.1, 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: S e loga = logK e (20) This formulation of reactive transport contains N sp + N r unknowns.These are matched by the N sp mass balances of all species (Equation ( 19)), plus the N e MAL (Equation ( 20)), and the N k definitions of kinetic reaction rates.The latter can be eliminated by substitution of Equation ( 21) into (19).Still, the large number of coupled non-linear equations makes the problem intractable.Therefore, various strategies have been proposed to simplify the problem [71,72].However, most of them are based on eliminating the equilibrium reaction rates by multiplying Equation ( 19) 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, u a = U a c a , is transported.In fact, the transport operator does not apply to the solid (immobile) fraction, u s = U s c s .
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 [12].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 (Figure 6).
The apparent discontinuity in Figure 6 simply reflects that the set of equations in each mineral zone will be different, thus yielding different solutions.A 3 does not participate in Min. 1 and behaves conservatively, up to x = 50, where Min. 2 leads to a new equilibrium.As mentioned above, what happens is that the RT equations are different in the two mineral zones, and a sharp reaction rate (mathematically a Dirac delta) occurs at their interface.While this is known to occur in reality, it hinders numerical solutions.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 Figure 6).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 3.3.Ways to overcome the problem numerically are described in Section 6.This leads to Conceptual Issue No. 5: mineral zone variations lead to sharp fronts in aqueous chemistry.

Basic Reactive Transport Paradigms
In the previous section, we analyzed the ADE for transport, the formulation of general chemical reactions, and how to couple them to write an 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 No. 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 strange (Vw/Ms) and unrelated to its meaning, as they reflect that  is mass per unit volume of water, whereas  is expressed per unit mass of solid.This is why sometimes the reported value is the dimensionless partition coefficient  =    ⁄ , which relates the sorbed to the aqueous solute mass in a volume of aquifer [74].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 ( 19)), while using Equation (23) 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

Basic Reactive Transport Paradigms
In the previous section, we analyzed the ADE for transport, the formulation of general chemical reactions, and how to couple them to write an 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, R i , can produce the rich biochemical diversity of Nature.

Sorption. Conceptual Issue No. 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 a XA = Ka X a A .However, equilibrium is usually obtained experimentally from batch tests where S A (moles of A sorbed per unit mass of solid) is measured for different values of c A .Results (S A vs. c A plots) often display a straight line through the origin (its slope K d is called distribution coefficient), that is: Note that the units of K d are strange (Vw/Ms) and unrelated to its meaning, as they reflect that c A is mass per unit volume of water, whereas S A is expressed per unit mass of solid.This is why sometimes the reported value is the dimensionless partition coefficient K p = ρ b K d /φ, which relates the sorbed to the aqueous solute mass in a volume of aquifer [74].Including sorption into the transport equation simply implies realizing that every mol of A that sorbs onto the solid disappears from solution.Therefore, sorption can be viewed as a sink term equal to R A = −ρ b ∂S A /∂t, where ρ d is the bulk density, mass of solid per unit volume of aquifer.Including this term into the solute transport equation (Equation ( 19)), while using Equation (23) for S A leads to: where R is the "retardation factor", given by R = 1 + ρ b K d /φ = 1 + K p .The impact of R can be deduced from the velocity of c A (coefficient of the gradient divided by the coefficient of the time derivative), which equals to q/φR, compared to the velocity of water, q/φ (i.e., c A is indeed retarded by a factor R with respect to water).Another way to look at it relies on noticing that the lefthand side of Equation ( 24) can be rewritten as φ∂c A /∂(t/R), which causes Equation ( 24) to be identical to that of a conservative solute, if t is substituted by t/R (i.e., the time scale is slowed down by a factor R, hence the name retardation).This is illustrated in Figure 7 for a case with constant R = 2.
Energies 2022, 15, x FOR PEER REVIEW 14 of 48 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 lefthand side of Equation ( 24) can be rewritten as   (/) ⁄ , which causes Equation ( 24) 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 Figure 7 for a case with constant  = 2.A second solute (dashed) enters with the same concentration as resident water, but salinity changes from 2 to 1 its retardation, producing a peak at the salinity front.Note that numerical solution of Sal (lines) with a Lagrangian method is visually identical to the analytical solution (symbols).The same numerical method produces some numerical dispersion for a third solute with constant retardation.
While linear adsorption is a nice simple concept, several remarks can be made.First, recall that, despite its empirical nature, Equation ( 23) is a rewrite of the mass action law.That is,  =   , where  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 [75].Ra tends to sorb strongly and its retardation depends on salinity (Michael et al. [74] 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 ( 24) is not valid because  is not constant.Thus, the lefthand side must be changed to ( +   )  ⁄ = ( )  ⁄ = ( )  ⁄ +  ()  ⁄ .Since  depends on salinity, the last term becomes an Ra source term when salinity increases, which explains the sharp increase in  at the salinity front (dashed lines in Figure 7).While the phenomenon is well-known in the Ra literature, Equation (24) and retardation are so deeply ingrained in the hydrological literature that Equation ( 24) remains the standard.Sorbed ions may be displaced by other solutes [76].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. [77], who explain the chromatographic effect (separation of cations depending on their selectivity).Details on the relationship between adsorption and ion exchange are discussed by Charbeneau [78] and Appelo [79,80].In essence, CE should be considered always, but it is not easy to handle because selectivity coefficients depend on numerous  A second solute (dashed) enters with the same concentration as resident water, but salinity changes from 2 to 1 its retardation, producing a peak at the salinity front.Note that numerical solution of Sal (lines) with a Lagrangian method is visually identical to the analytical solution (symbols).The same numerical method produces some numerical dispersion for a third solute with constant retardation.
While linear adsorption is a nice simple concept, several remarks can be made.First, recall that, despite its empirical nature, Equation ( 23) is a rewrite of the mass action law.That is, K d = FKa X γ A , where F includes the effect of expressing the concentration of X − A as S A per unit mass of solid (i.e., F = a XA /S A ). Therefore, a constant K d 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, K d varies in response to changes in salinity (which changes activity coefficients and competes for sorption sites) or temperature (changes K) or pH (protons also compete).A simple example illustrates the effect of varying K d .Assume that salinity, which reduces K d , 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 [75].Ra tends to sorb strongly and its retardation depends on salinity (Michael et al. [74] report values of K d ranging between 10 and 10 5 for fresh water and between 1.5 and 10 4 for salt water).Under these conditions, Equation ( 24) is not valid because K d is not constant.Thus, the lefthand side must be changed to Since R depends on salinity, the last term becomes an Ra source term when salinity increases, which explains the sharp increase in c A at the salinity front (dashed lines in Figure 7).While the phenomenon is well-known in the Ra literature, Equation (24) and retardation are so deeply ingrained in the hydrological literature that Equation ( 24) remains the standard.
Sorbed ions may be displaced by other solutes [76].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. [77], who explain the chromatographic effect (separation of cations depending on their selectivity).Details on the relationship between adsorption and ion exchange are discussed by Charbeneau [78] and Appelo [79,80].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 [81][82][83].This complexity explains why numerous empirical models have been developed to estimate CE parameters (e.g., [84]), including artificial neural networks [85].However, 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.That being said, 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., [86][87][88][89]) 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 No. 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 buildup or decay reaction (r = kc), 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 n-th order reactions, whose rate is given by: In the absence of transport, the integration of Equation ( 25) yields the evolution of The characteristic time for these reactions is t r = 1/kc n−1 0 .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 Damköhler number, defined as: where τ is the characteristic transport time.The Damköhler 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 c 0 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 is addressed in Section 5).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 Conceptual Issue No.2: diffusion controls microscales).The resulting number would be the product of Da (defined with the advective time) and the Peclet number.The Damköhler number is used for assessing whether a given kinetic reaction needs to be taken into account.When Da is small (much smaller than 1, say Da < 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 Da is large (much larger than 1, say Da > 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 (Figure 8).The analysis is more complex, and the usefulness of the Damköhler 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 Figure 8.It is clear that Da close to 1 implies a sizeable reaction in all cases.However, it is not true that a large Da implies that the solution is close to equilibrium.
Energies 2022, 15, x FOR PEER REVIEW 16 of 48 all cases.However, it is not true that a large  implies that the solution is close to equilibrium.
Figure 8. Dimensionless concentration versus dimensionless distance after a pulse injection of a tracer suffering first order decay (a).The reactions barely affect the solute, if Da < 0.1, but make it disappear for Da > 3.This is not necessarily true for non-linear kinetics (b).Da = 1 implies a sizeable reaction but note that significant amounts of solute remain for Da = 100, if n = 4. Also note 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 ( 25) in an attempt to simplify complex models (recall, again, that simplicity helps robustness).However, Reigner et al. [25] 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 [90] 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. [91]).Valhondo et al. [2] 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 [92].These 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.The issue is especially relevant for biochemical problems, where kinetics may evolve in space and time in response to changes in redox state and growth of the microbial community mediating reactions.When several kinetic reactions co-exist, the fastest ones control the accuracy and efficiency of the solution method.The fastest reaction will generally be the one for which the Damköhler number needs to be computed, but this reaction may change along the simulation.
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 1, but make it disappear for Da > 3.This is not necessarily true for non-linear kinetics (b).Da = 1 implies a sizeable reaction but note that significant amounts of solute remain for Da = 100, if n = 4. Also note 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 ( 25) in an attempt to simplify complex models (recall, again, that simplicity helps robustness).However, Reigner et al. [25] 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 [90] 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. [91]).Valhondo et al. [2] 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 [92].These 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.The issue is especially relevant for biochemical problems, where kinetics may evolve in space and time in response to changes in redox state and growth of the microbial community mediating reactions.When several kinetic reactions co-exist, the fastest ones control the accuracy and efficiency of the solution method.The fastest reaction will generally be the one for which the Damköhler number needs to be computed, but this reaction may change along the simulation.
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 [93].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 5.In this context, the Damköhler number should be viewed as a rough preliminary indicator of the relative importance of kinetic reactions, but cannot be used blindly.An easy paradigm for equilibrium reactions is the one of two ions, A 1 and A 2 , that precipitate as an immobile solid mineral S 3s in equilibrium ( A 1 + A 2 S 3s , example of Figure 4, its solution is the first half, x < 50, of Figure 5).The solution to this problem illustrates the methodology for solving reacting transport [73].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 ( 19) and ( 20): where c i is the concentration of A i , and S 3 is the concentration of the mineral.These four equations contain four unknowns: r, c 1 , c 2 , and S 3 .They are coupled, and there is no explicit expression for the reaction rate, r, 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 Figure 5. Multiplying the vector of the first two equations times the component matrix is equivalent to subtracting Equation ( 28) from (27), which yields: where u = c 1 − c 2 , which is conservative indeed a component, and easy to solve.We reduced the problem from four coupled non-linear PDEs to one 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 (31) for u.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 (31) can be solved for u.
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 2).But in this case, speciation is easy, because it consists of solving the system of equations given by Equation (30) and Step 4: Evaluation of the reaction rate.Substitution of the concentration of the secondary species, A 2 , into its transport equation leads to [73] 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 D. Beyond its basic simplicity, this equation illustrates that under our ideal conditions the reaction rate is always positive (i.e., precipitation occurs because ∂ 2 c 2 /∂u 2 is positive), as discussed by Rubin [12].This observation could be deduced from Figure 4, 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 [94], which affect ∂ 2 c 2 /∂u 2 .As a result, dissolution may occur.For example, changes in ion strength produce dissolution of calcite when calcite is saturated with a fresh and saltwater mix, which explains the development of coastal karst [95,96].
Beyond the limitations of the example, Equation ( 32) is of universal application.This equation has numerous implications.The rate of reaction depends on chemistry, which controls ∂ 2 c 2 /∂u 2 but also on transport processes, controlling the gradient of u.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 come back to this is issue when discussing the "Water Mixing Approach" to transport in Section 6.Another paradoxical feature of Equation ( 32) is that reaction does not necessarily take place where concentrations attain their maximum values.In fact, the reaction rate equals zero when u is maximum or minimum.
It is also interesting to notice that, in the absence of the first contribution on the righthand side of Equation (32), all terms are proportional to the dispersion tensor, D, thus strengthening the relevance of mixing processes to the development of fast reactions.In particular, the term ∇ T uD∇u is a measurement of the mixing rate, which is consistent with the concept of dilution index, as defined by Kitanidis [97] 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 [98], but also the geochemical processes occurring in the system.In fact, De Simoni et al. [73] 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 S 3 of mineral is obtained using r in Equation (30), 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, are described in Section 6.These equations were used in the early part of Figure 8 to compute mineral precipitation in response to a pulse injection, for which the transport Equation (31) can be solved analytically.
The summary of this section is that mixing is fundamental for RT modeling.From a conceptual point of view, this statement is trivial (of course, reactants must come into contact for reactions to occur).However, it may be non-trivial in porous media when reactions may be restricted to specific reaction sites.This is usually ignored in practice because porous media are usually treated as a continuum (an exception are pore scale models), so that access to reaction sites may be well-represented by mixing, but it may not be sufficient in real systems (e.g., biochemical reactions in biofilms or mineral reactions in low permeability locations).As we see in the next section, a problem with the ADE is that it does not simulate mixing or access to reaction sites correctly, which has led to numerous efforts to seek alternative transport formulations.

Representing the Impact of Unknown Heterogeneity on Transport
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 [99,100].Calibrated dispersivity grows with the experiment scale [101,102], and calibrated advective porosities often depend on the residence time [103] or on the flow direction [104,105].This type of behavior is called anomalous, which we find ironic because it is the one always observed in nature.Since an adequate transport equation is critical for RT, we devote this section to discussing 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 the broad range of scales involved in RT (Figure 9), 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, gaining insight into RT, and preventing errors as the ones we have made in the recent past.
Efforts to quantify anomalous transport started during the last part of the 20th 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 [106][107][108].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.), which motivated the search for alternative equations.However, the summary of this introduction is that transport difficulties arise from the fact that velocity varies in uncertain manners.Mean water fluxes can be reproduced, but uncertain fluctuations cannot.The goal of effective equations is not so much to describe these fluctuations as their impact on RT.

Representing the Impact of Unknown Heterogeneity on Transport
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 [99,100].Calibrated dispersivity grows with the experiment scale [101,102], and calibrated advective porosities often depend on the residence time [103] or on the flow direction [104,105].This type of behavior is called anomalous, which we find ironic because it is the one always observed in nature.Since an adequate transport equation is critical for RT, we devote this section to discussing 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 the broad range of scales involved in RT (Figure 9), 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, gaining insight into RT, and preventing errors as the ones we have made in the recent past.
Efforts to quantify anomalous transport started during the last part of the 20th 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 [106][107][108].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.), which motivated the search for alternative equations.However, the summary of this introduction is that transport difficulties arise from the fact that velocity varies in uncertain manners.Mean water fluxes can be reproduced, but uncertain fluctuations cannot.The goal of effective equations is not so much to describe these fluctuations as their impact on RT.

The Search for an Effective Transport Equation
After establishing that the ADE is a poor representation of transport, which is our Conceptual Issue No. 9, researchers focused on finding "effective equations" in the hope

The Search for an Effective Transport Equation
After establishing that the ADE is a poor representation of transport, which is 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 Energies 2022, 15, 925 20 of 47 at different scales to address anomalous transport [109].We make here a short detour into methods to solve the ADE and so as to properly frame the numerous equations that were proposed to overcome the limitations of the ADE.They are schematically described in Figure 10.Solution methods can be divided into Eulerian and Lagrangian (top row of Figure 10).Eulerian methods discretize the equation in space and time over a fixed grid.Accurate Eulerian methods suffer numerical instability.Overcoming this instability is the main motivation for Lagrangian methods, which solve the equation over a grid that moves with the fluid velocity and 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 tracking methods and, specifically, Random Walk (RW) methods, based on Einstein's [110] concept that Fick's Law results from the mean behavior of particles subject to Brownian motion [111,112].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 (Figure 10d) or, far better, by associating a kernel to each particle (Bells method, Figure 10d) [113,114].Note that bells methods implicitly suggest mass exchange between particles, which has not yet been explored.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 (Lagrangian methods are not well-suited for microscales, which are diffusion dominated, Figure 3).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 smoothed particle hydrodynamics (SPH) and 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 [115,116].The Lattize Botzman method is similar in treating particle interactions, but on a fixed grid [117].A description of methods to approximate molecular behavior is presented by Meakin and Tartakovsky [36].The summary of this long detour into solution methods is that the choice should be solely based on numerical convenience, but all methods attempt to solve the same ADE, despite their different looks.
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 an RW in which the time steps for velocity transitions are random [118].Variable time steps allow the modeler to represent non-Fickian behavior (both waiting in immobile zones and fast travelling through preferential flow paths).By appropriately selecting 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 nontrivial observations, such as tailing in field and laboratory experiments [6], proper scaling of variable time push-pull tests [119] and especially pore network model outputs [120].A second option is to substitute the time and/or space derivatives in the ADE by fractional order derivatives [121,122].While using fractional order derivatives sounds awkward, the concept dates back to Leibnitz himself and has been the subject of much mathematical work [123].Its use in hydrology emerges from an 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 [124].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 [125].As it turns out, by appropriately choosing the memory function, the three methods are identical [126,127].From an RW perspective, the memory function can be viewed as the residence time distribution in immobile water.Therefore, these methods are directly related to the Multi-Rate Mass Transfer (MRMT) method [128].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 [129][130][131].The summary of this section is that, as summarized in Figure 11, all these effective formulations (1) attempt to reproduce the effect of unmodeled, uncertain velocity fluctuations, (2) they can all be viewed in PDE or RW forms, and (3) despite their mathematical difficulty, they are not that different (in fact, they are virtually identical [126,127]).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, g, that represents the delay of slowly flowing portions of the medium and the diffusive exchange between them [125].As it turns out, by appropriately choosing the memory function, the three methods are identical [126,127].From an RW perspective, the memory function can be viewed as the residence time distribution in immobile water.Therefore, these methods are directly related to the Multi-Rate Mass Transfer (MRMT) method [128].This method considers a discrete number of well-mixed immobile zones that exchange mass with transfer coefficients, inverse of mean residence times, α j (j = 1, • • • , N im , number of immobile zones).The probability of each zone, β j , 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 [129][130][131].The summary of this section is that, as summarized in Figure 11, all these effective formulations (1) attempt to reproduce the effect of unmodeled, uncertain velocity fluctuations, (2) they can all be viewed in PDE or RW forms, and (3) despite their mathematical difficulty, they are not that different (in fact, they are virtually identical [126,127]).A variation of MRMT is 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 [132][133][134][135][136][137][138].Unlike the "effective equation" approaches, hybrid models do not attempt to simulate "the effect of" velocity variability, but the variability itself.As such, they should be viewed as potentially more accurate, provided an accurate representation of heterogeneity.Methods for this coupling are reviewed by Molins and Knaber [32].

Lessons Learned, Formal Upscaling, and on-Going Efforts
A problem with 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 through natural media, they fail to adequately distinguish mixing from spreading (Figure 12).The distinction is not as trivial as the figure might suggest.Mixing and dispersion are intricately linked: (1) velocity fluctuations,  in Figure 12, increase the interface between different water bodies, thus increasing diffusive mixing, and (2) diffusion feeds back into dispersion by causing particles to jump between flow tubes [139].
Another problem with the above effective equations is that, with the exception of MRMT, they do not provide a value of concentration in its geological setting, which is required for RT.Worse, in our view, reactivity is often largest where velocity is smallest (clay particles, biofilms, etc.).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 [140] or geochemical localization, which stands for the occurrence of reactions that would not be possible in a locally equilibrated medium (e.g., acid injection should lead to carbonate dissolution, but may produce precipitation in isolated zones where other sources of alkalinity are present) [141,142].Still, much of the tailing occurs due to very slow advection and does not necessarily imply mixing [143,144].The need to distinguish between mixing and spreading has prompted Soler-Sagarra et al. [145] to rewrite the dispersion term as ∇( ∇ ), where  is the dispersion flux (Figure 12), 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: A variation of MRMT is 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 [132][133][134][135][136][137][138].Unlike the "effective equation" approaches, hybrid models do not attempt to simulate "the effect of" velocity variability, but the variability itself.As such, they should be viewed as potentially more accurate, provided an accurate representation of heterogeneity.Methods for this coupling are reviewed by Molins and Knaber [32].

Lessons Learned, Formal Upscaling, and On-Going Efforts
A problem with 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 through natural media, they fail to adequately distinguish mixing from spreading (Figure 12).The distinction is not as trivial as the figure might suggest.Mixing and dispersion are intricately linked: (1) velocity fluctuations, J D in Figure 12, increase the interface between different water bodies, thus increasing diffusive mixing, and (2) diffusion feeds back into dispersion by causing particles to jump between flow tubes [139].
Another problem with the above effective equations is that, with the exception of MRMT, they do not provide a value of concentration in its geological setting, which is required for RT.Worse, in our view, reactivity is often largest where velocity is smallest (clay particles, biofilms, etc.).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 [140] or geochemical localization, which stands for the occurrence of reactions that would not be possible in a locally equilibrated medium (e.g., acid injection should lead to carbonate dissolution, but may produce precipitation in isolated zones where other sources of alkalinity are present) [141,142].Still, much of the tailing occurs due to very slow advection and does not necessarily imply mixing [143,144].The need to distinguish between mixing and spreading has prompted Soler-Sagarra et al. [145] to rewrite the dispersion term as ∇ J D ∇ c , where J D is the dispersion flux (Figure 12), instead of ∇(D∇ c ), 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 [146,147]   By 2010, many started to believe that effective equations were failing, which prompted Neuman and Tartakovsky [148] 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 [149,150] 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 along two avenues.The first one is to follow the traditional stochastic upscaling methodology that consists of (1) generating realizations of heterogeneous fields of hydraulic conductivity (and/or chemical properties of immobile phases, pore geometries, etc.), (2) simulating explicitly RT through each realization, and (3) extracting relevant statistics and/or conceptual lessons from the simulation results.The three steps display uncertainties.Proper generation of geological variability remains a challenge [149,151,152].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 [153,154].Even if the procedure is properly followed, formal upscaling remains difficult, because physical and chemical heterogeneity cannot be upscaled separately [155]: 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 [13,142,156].By 2010, many started to believe that effective equations were failing, which prompted Neuman and Tartakovsky [148] 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 [149,150] 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 along two avenues.The first one is to follow the traditional stochastic upscaling methodology that consists of (1) generating realizations of heterogeneous fields of hydraulic conductivity (and/or chemical properties of immobile phases, pore geometries, etc.), (2) simulating explicitly RT through each realization, and (3) extracting relevant statistics and/or conceptual lessons from the simulation results.The three steps display uncertainties.Proper generation of geological variability remains a challenge [149,151,152].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 [153,154].Even if the procedure is properly followed, formal upscaling remains difficult, because physical and chemical heterogeneity cannot be upscaled separately [155]: 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 [13,142,156].
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 [157,158].This brings renewed interest on the role of transverse dispersion, which has not received sufficient attention in the past [159,160].However, 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 [161].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 (J D in Figure 12).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. [162] 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.We believe that two new steps are required: (i) the velocity distribution of the domain has to be discretized [163], and (ii) velocity transitions have to be properly reproduced [164][165][166].This view has led to a new phase space formulation [167], where velocity is considered as another independent variable of solute (much like space or 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 [145].The relevance of velocity fluctuations at fixed spatial steps is reflected by the recent review of spatial Markov approaches [150].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 [168], " . . .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, any C with valence less than +4) decay and redox concepts, prior to discussing the concepts related to kinetics and biofilms.

Organic Carbon Decay and Redox Sequence
Decay of OC involves transforming it into another organic species with higher valence or full mineralization to inorganic carbon (valence +4).This is a redox reaction where the original OC is the (main) electron donor, as shown by the redox half reaction of CH 2 O, which we use as example of OC: R1: 0.25CH 2 O + 0.5H 2 O → 0.25HCO − 3 + 1.25H + + e − Actually, microorganisms catalyze this reaction and grow on the energy it releases.This can be acknowledged by including biomass as a (possibly many) species.A simple way to express cell decay and synthesis is [169] − + 0.125H + When multiple TEAs are present, microorganisms tend to use the thermodynamically most favorable one (i.e., the one yielding most energy).When exhausted, other TEAs are used sequentially in order of their thermodynamic preference.The order of preference, or redox sequence, is O 2 [51,170].Figure 13 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 [171].

Traditional Kinetic Models (Monod and Michaelis-Menten)
A very simple enzyme kinetic rate law was developed by Michaelis and Menten [174].They assumed that an enzyme, E, catalyzes the transformation of a substrate, S, to a product, P, according to two reactions ( E + S ↔ ES → E + P , where ES is an enzyme substrate complex of E and S).Assuming that the first reaction occurs in equilibrium (i.e., the MAL is valid), that the total c E is fixed and much larger than c S , and that the second is kinetic with its rate proportional to c ES , it is easy to show that the rate of production of P is where µ m is the kinetic rate constant of the second reaction and K M,s 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 [175].This equation expresses the dependence of r on limiting factors in a simple way: r grows linearly on c S , with slope r max /K M , when c S is small (approximately up to r = 0.5r max , when c S = K M , which is why it is called half-saturation constant), but tends to the maximum rate, r max , when c S is much larger than K M .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 [169,176]: where c B is the biomass concentration, which can be considered proportional to the term c E + c ES (total enzyme plus enzyme complex).The rate of cell synthesis may be limited not only by the substrate, but also by TEAs, such as O 2 or NO − 3 , which is qualitatively reproduced by a multiplicative Monod equation [169,177]: where A is the TEA.For the rate of endogenous decay, one often assumes a linear relation with the biomass concentration.
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 (I Aj = K inh M,j / K inh M,j + c Aj , where K inh M,j is the inhibition constant and subscript j refers to inhibiting electron acceptor in the redox sequence, (j = 1 for O 2 , j = 2 for NO − 3 , etc.).These factors display a behavior opposite to the Monod factors (constant and equal to 1 when c Aj is too small to inhibit the process and tends to 0, thus inhibiting it, when c Aj becomes larger than K inh M,j ).Several multiplicative factors can be added.The OC decay law becomes (e.g., [178]): where subscript i refers to the TEAP in the redox sequence.This equation replicates the preferential use of TEAs in the redox sequence.Numerous variants of Equation ( 37) have been used.Instead of a Monod factor for c CH 2 O , a linear relation can be used [179].In addition, instead of the inhibition factors, one can use switch functions that are either one or zero, depending on the concentration of an electron acceptor [180].As Equation (37) does not consider Gibbs free energy, it could give rates that are thermodynamically unfavorable [181].Therefore, factors may have to be included to correct this (e.g., [182]).This equation can also be modified to allow for negative concentrations that may arise due to instabilities during transport solution [183].
An alternative approach for organic carbon decay in the redox sequence is the Partial Equilibrium Approach [184].Rather than formulating rate laws for each TEAP, i in Equation (37), 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 N 2 by O 2 (Saaltink and Rodríguez-Escales, 2021, submitted).

Microbial Based Kinetic Models
Models that use Equation (37) and, though to a less extent, the Partial Equilibrium Approach, require many parameters (for each TEAP i one needs a value for r max,i , K M,i and several K inh M,i ).Usually, they are obtained from the literature or by fitting the model to measured concentrations.If we compare Equations ( 35) and (37), we can see that the maximum rate of a particular TEAP, r max,i , is related to the biomass, B i , capable of using the particular electron acceptor i.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 [171,176], 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 r max,i 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 [198].
The basic mass balance for a specific microbe, assumed immobile, simply results from Equation (19), which reads: where superscript j refers to the different substrates that may sustain the growth of the specific microbe i.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 (33) 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" [199].
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 c Bmax .
This term is needed to prevent a boundless increase in c when neither substrate nor TEAs are limited.In fact, terms with this form were proposed by Kindred and Celia [200] and Schäfer et al. [183].They limited biomass growth by the available pore space.However, as we see in Section 5.5, we view microorganisms as species belonging to a specific phase (the biofilm), where c B cannot be limited by the porosity.We conjecture that c Bmax 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 such as where c BTmax is the maximum total microbial concentration.Equation (40) also contains factors to limit growth by the availability of substrate and electron acceptors, similar to Equation (37), 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 [201,202].Xu [203] and Tang [199] studied Monod (r = µc B M OC ), logistic (r = µc B L B ) and hybrid rate laws (r = µc B M OC L B ) and argue in favor of the hybrid rate laws.Many other factors must be acknowledged when simulating microbial growth.Stolpovsky et al. [204] account for the possibility of biomass entering into a "dormant" state during starvation periods, and then reactivating 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 obtain 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 [205].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 Equations ( 38) and (40), are mostly applied to constructed wetlands, which are an alternative for treating waste water.Reviews exist on this type of modeling [206,207].Arguably, they contain the most complex and complete biochemical systems, used for porous media up to now.An example is CMW1 (Constructed Wetland Model No. 1) [201], 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 ( 37), but they do not contain logistic factors, such as in Equation (40).CMW1 has been coupled to flow and transport codes, such as Retraso [208,209] and Hydrus [210,211].There are a few modeling studies in marine and lake sediments that also consider biomass (e.g., [202]), although most of them use the traditional rate laws of Section 5.2.They 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 5.2, but we also take it as indicative of the need for further research.

Genomics Based Metabolic Models
The fact that the 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. [212] and Louca et al. [213] went further and applied a complex 'gene-centric' approach using marker genes as experimental references, the production of which was linked to biogeochemical reactions in the ocean.Pagel et al. [214] 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 multi-objective calibration.
Instead, Scheibe and coworkers propose an alternative approach: constraint-based modeling [215][216][217][218]. 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 [219] and is represented, similar to Equation (19), by a stoichiometric matrix and metabolic and exchange rates, but neglecting the transport terms and assuming steady state within the cells, yielding S t r = 0. Given concentrations of chemical species of the cell's environment, the rates, r, 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., S t r = 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 these kinds of models to RT has been restricted to single microbial species metal redox problems [219], which again suggests that much remains to be done.

Biofilm Modeling
Microorganisms are mostly concentrated in biofilms, a biological phase that covers 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 [220].EPS is a fibrous gel-type matrix composed of polysaccharides, lipids, proteins, and DNA materials excreted by the microbes [221,222].The importance of the biofilm lays in the fact that almost all biochemical reactions, such as organic carbon decay, described in Section 5.1, take place in the biofilm.As the hydraulic conductivity of biofilm is very low [223], exchange of solutes between pore water and biofilm is controlled by diffusion.As a result, overall rates of biochemical reactions become affected by the size and thickness of the biofilm [224].
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 (Multi-Rate Mass Transfer), explained in Section 4. 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 [183].In fact, its inclusion in continuum formulations has been formalized mathematically [225].The biofilm-porewater exchange is assumed proportional to the difference between the concentration in the biofilm and porewater (α(c a − c b ), where α is an exchange coefficient and c a and c b are the concentra- tions in the porewater and biofilm), which is identical to the exchange rate between mobile and immobile zone for MRMT in Figure 11a.
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., [226][227][228][229]).In addition, its effect on the retention curve has been studied and modeled [230][231][232].We contend that the most relevant impact for RT is the effect of biofilm growth on transport properties such as dispersivity, and especially parameters for the porewater-biofilm exchange, which control biochemical reactions and biofilm growth [233].Experimental results have demonstrated that, indeed, biofilm growth induces heterogeneities that affect the transport in porous media [234,235].Moreover, it forms preferential flow paths and stagnation zones [232].Therefore, to understand biochemical RT, we need models that are capable of simulating dynamic biofilm growth and the increase of the heterogeneity of porous media simultaneously.

Solution Methods and Tools
Beyond the conceptual difficulties introduced in Sections 4 and 5, RT modeling is challenging because of the large number of non-linear coupled PDEs (Equations ( 19)-( 21)) 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 4. Since analytical solutions are restricted to very simple conditions (we discussed a few of them in Section 3), 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 (Figure 14).All these methods accept numerous variants.Here we only describe the main features.Specific features are discussed in detail in methodological reviews [38,236].As we see, the choice of solution method may affect the processes that are modeled and the way they are represented.Therefore, surprising as it may sound, our final Conceptual Issue, No. 11, is that one should "select an appropriate RT solution method".
The SIA consists of two steps at each iteration: (1) solve the transport equations of the aqueous component (u a ) 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 Figure 14a).For the first step we rewrite the transport Equations ( 19)-( 21) as: where superscript n refers to iteration number.By using the source-sink term (f) 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 (c n ) from the aqueous components (u a ) 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, Equation ( 12) or (13).This leads to a standard speciation calculation from components, briefly discussed in Section 2.2.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 Equation (41) (Figure 14a).These calculations are performed for each node separately, which facilitates parallelization.After calculation of all concentrations, the sink-source term (f n ) can be updated easily using Equation (42).
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 [67] or be used as a speciation library to solve the chemical calculations [237].Some examples of transport codes that have been coupled to PHREEQC are MODFLOW/MT3D [238], COMSOL [239], or PH-WAT [240].HYDRUS-1D [241] also solves 1D reactive transport (RT) problems along with heat transport and density dependent flow cases.HYDROGEOCHEM [242] solves coupled Thermo-Hydro-Mechanical-Chemical (THMC) 3D models.TOUGHREACT is another THMC and multiphase flow code [243], largely based on CORE-2D [244].Finally, Van Der Lee [245] 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 [246], and OpenGeoSys [247], which is also a THMC code.
The DSA consists of performed one single step (see Figure 14b) per iteration by substituting chemical equations into the transport equations and applying the method of Newton-Raphson.Concentrations of secondary species can be written as a function of primary concentrations by writing the mass action law (Equation ( 13)) as log c 2 = S * e1 log a 1 + log K * + log γ 2 (44) which can be substituted into the transport Equations ( 19)- (21).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 [63].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: the 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 c 2 is an explicit function of c 1 in Equation ( 44)).Therefore, comparison of methods is not straightforward.SIA is simpler to code and leads to smaller systems of equations than DSA, but DSA should converge faster.Initial comparisons favored SIA [38], but more recent comparisons favor DSA in terms of performance [61,236,248].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 [39].Some examples of DSA codes are MIN3P [249], PFLOTRAN [250] and RETRASO [47].CRUNCHFLOW [251] is a code that can use both SIA and DSA.Some codes have been adapted to high performance computing, such as MIN3P-HPC [252] or ParCRUNCH-FLOW [253].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 [65,254,255].In particular, the Method of Lines [256], a general method for solving coupled PDEs, has been used with notable success to solve RT problems [257,258].
The fact that transport equations are virtually identical for all mobile species prompted Soler-Sagarra [259] to propose an alternative concept, the Water Mixing Approach (WMA).The basic idea is to understand transport as water advection and mixing (concentrations being attributes of this water) instead of transport of each individual species.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 [259].The concept can be extended to reactive transport by simply adding the reaction term to the mixing equation, which leads to rewrite Equation (19) as: where λ lm 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.) m to l during the time step, expressed as a fraction of the water volume (or mass) in cell l; N l is the number of cells connected to l, and R l 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. [260] to compute end-member mixing ratios for the interpretation of concentration data.Here, the goal is the opposite; given the mixing ratios (easy to obtain from transport solvers), compute reaction amounts R l ∆t and the updated concentrations.This calculation is relatively easy to implement in existing speciation codes.Soler-Sagarra et al. [259] details the procedure and implemented it on CHEPRO [68].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 Figure 14c).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, as discussed in Section 4, Lagrangian formulations (Equation ( 45)) are more adequate for advection dominated problems.This can be carried out with mixed Eulerian-Lagrangian methods using variants of the method of characteristics [261][262][263][264], Localized Adjoint Methods [265,266] or stream-based grids [267].Most of these methods risk suffering fluctuations or numerical dispersion, which motivates purely Lagrangian methods (Figure 10b,c).Figures 6 and 7 were obtained with a Lagrangian method (isochronal grids) coupled to the WMA. Figure 7 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 3.1 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 (Figure 7).
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 5, 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 [268][269][270][271][272]. The probability of particle collisions, interactions, and transformations may be made to depend on the scale [273].This implies treating particles as individual species, and facilitates assigning each particle its specific transport properties, thus getting around the above problem with sorption [274].The approach has been successful in reproducing bi-molecular reactions in anomalous transport settings [275][276][277], including biomass growth [273].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 [278] and assigns each particle a fraction of the total mass of a given species.Still, the Kernels introduced in Section 4 as part of the method to recover concentration fields from particles (Figure 10d) tend to be used to define a support volume to facilitate the evaluation of particles interaction [114,279,280].A final twist on the method is to extend the concept to total concentrations [281], in the hope that it allows describing transport accurately (as discussed in Section 2.2 and 3.2, 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 [145] and suggests that the RW community may eventually blend with the rest of the RT community.
We conclude this section by adding a summarized discussion on the selection of an "appropriate solution method for efficient and accurate RT".First of all, all consistent methods will converge to the exact solution as the spatial and temporal discretization become small (and all the methods discussed here are consistent).Therefore, the decision must be made on the basis of efficiency and which methods are most tolerant to coarse discretization.Most methods work with aqueous components to handle homogeneous reactions in equilibrium.Therefore, the issue is how to handle kinetic and equilibrium heterogeneous reactions.Equilibrium reactions are effectively treated as kinetic in SIA, which leads to underestimation of reaction rates and to erroneous spatial distribution of reactions [96].Worse, conventional SIA formulations disregard "mineral zonation" (recall Conceptual Issue No. 5, that the set of equilibrium reactions determines the chemical system).For this purpose, it is important to bear in mind that fast kinetics is equivalent to equilibrium.This is especially important in biochemical problems, where kinetics may be controlled by redox state, so that redox zonation will have a similar effect to mineral zones (recall Figure 5) in controlling the effective definition of the chemical system.In short, for fast kinetics, DSA or WMA will generally be more appropriate, unless small time steps are feasible.In cases where several kinetic reactions occur, the time step is controlled by the fastest one.On the other hand, when all heterogeneous reactions are slow, SIA (or even SNIA) may suffice.
A similar comment can be made regarding numerical dispersion.The actual rate of fast reactions is controlled by mixing (recall Conceptual Issue No. 8).Therefore, numerical dispersion causes errors in reaction rates and the spatial distribution of reactions (see example in [259]).This implies that care must be taken in discretization of Eulerian methods and may favor the use of Lagrangian methods. 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 [145] and suggests that the RW community may eventually blend with the rest of the RT community.
We conclude this section by adding a summarized discussion on the selection of an "appropriate solution method for efficient and accurate RT".First of all, all consistent methods will converge to the exact solution as the spatial and temporal discretization become small (and all the methods discussed here are consistent).Therefore, the decision must be made on the basis of efficiency and which methods are most tolerant to coarse discretization.Most methods work with aqueous components to handle homogeneous reactions in equilibrium.Therefore, the issue is how to handle kinetic and equilibrium heterogeneous reactions.Equilibrium reactions are effectively treated as kinetic in SIA, which leads to underestimation of reaction rates and to erroneous spatial distribution of reactions [96].Worse, conventional SIA formulations disregard "mineral zonation" (recall Conceptual Issue No. 5, that the set of equilibrium reactions determines the chemical system).For this purpose, it is important to bear in mind that fast kinetics is equivalent to equilibrium.This is especially important in biochemical problems, where kinetics may be controlled by redox state, so that redox zonation will have a similar effect to mineral zones (recall Figure 5) in controlling the effective definition of the chemical system.In short, for fast kinetics, DSA or WMA will generally be more appropriate, unless small time steps are feasible.In cases where several kinetic reactions occur, the time step is controlled by the fastest one.On the other hand, when all heterogeneous reactions are slow, SIA (or even SNIA) may suffice.
A similar comment can be made regarding numerical dispersion.The actual rate of fast reactions is controlled by mixing (recall Conceptual Issue No. 8).Therefore, numerical dispersion causes errors in reaction rates and the spatial distribution of reactions (see example in [259]).This implies that care must be taken in discretization of Eulerian methods and may favor the use of Lagrangian methods.

Discussion and Conclusions
This review consists of two parts.In the first part (Sections 2 and 3) we revised RT concepts of traditional RT, largely "geochemical", models.In the second, we 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 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 facing 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 (Conceptual Issue No. 1) or that appropriate reactions and species must be included (Conceptual Issue No. 2).We have listed them to point that they deserve further thought.It is not sufficient to bring together different specialists; they need to work together (interdisciplinarity) and to understand each other.For example, when, why, and which aqueous complexes need to be included in an RT model may not be evident to a solute transport specialist.As a result, some published papers appear to ignore Conceptual Issue No. 2. 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 (or, even less, well-modeled).We 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 for RT.Numerous alternatives were proposed to improve on this.They improve on mixing and overcome many of the limitations of the transport equation, such as tailing and scale dependence of dispersion or porosity.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, and because it allows reproducing some of the subtle RT effects, such as chemical localization.We found 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 a large number of non-linear equations.The two traditional methods are DSA (or Newton-Raphson) and SIA (or Picard, 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, which live in biofilms and make them grow.The impact of biofilm growth on transport beyond clogging needs to be properly addressed.Genome based formulations sound like the right approach for biochemical RT.However, their complexity and the difficulties inherent to assigning functionality to genomic data 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 not only expertise on flow, hydrogeochemistry, biology, transport or numerics, but also close interaction between experts to understand interactions among processes.

Diffusion controls microscales; advection controls macroscales; the Peclet number is ambiguous
The ADE characteristic length is D/q.Pe is not a property of the equation, but of one model scale.Several scales may be relevant in practice.Better compare these to D/q, which leads to the above statement for diffusion coefficients and GW flux (Section 3.2). 3 It is important to select the appropriate set of species and reactions Selection of species and reactions requires understanding geochemistry, where it can be helped by speciation codes (Section 2.2), and biochemistry, where selection of biological reaction is more difficult, though at least they should be thermodynamically favorable (Section 5.1). 4

Equilibrium reactions reduce the size and complexity 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 (Sections 2.2 and 2.3). 5

Mineral zones variations lead to sharp fronts in aqueous chemistry
The set of equilibrium reactions in any subdomain (loosely, mineral zone) determines the chemical system and the equations that need to be solved (set of mass action laws, components, etc.) and, thus, water chemistry (Section 3.3).Note that varying microbial communities may have a similar effect (redox zonation). 6 Each species travels at a different velocity Solid species do not move, whereas aqueous species do.Exchanging species are retarded with respect to wa-ter while mobilizing others (Section 3.1).This hinders the accuracy of Lagrangian methods and leads to non-linear dependence between total concentration of a component (u) and its aqueous part (u a ) in the transport equation (Equation ( 22)). 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 (Sections 2.2 and 5).Since kinetic rates are often variable, beware of the Damköhler 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 (32), which implies that mixing rates need to be accurately modeled for accurate RT. 9

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 it 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.Yet, accurate biochemistry requires acknowledging both microbial communities and biofilms.
11 Select an appropriate solution method for efficient and accurate RT Solution methods may control the way mixing, and thus fast reactions, are represented, whether chemical localization can be reproduced, etc.Therefore, not only for efficiency, but also for conceptual consistence, an appropriate solution method needs to be adopted.

Notation
Selected variables in scalar form (bold face for vectors and matrices).

Figure 2 .
Figure 2. The number of papers identified as dealing with RT in the SCIE has increased linearly since 1993.

Figure 1 .
Figure 1.Reactive transport (RT) describes the transport of solutes (displacement by advection, spreading by dispersion, and dilution by mixing) and the reactions (red arrow) among solutes and with immobile phases.The mathematical nature of RT depends on the type of reactions and species involved[12].

Figure 2 .
Figure 2. The number of papers identified as dealing with RT in the SCIE has increased linearly since 1993.

Figure 2 .
Figure 2. The number of papers identified as dealing with RT in the SCIE has increased linearly since 1993.

Figure 3 .
Figure 3.Time evolution of diffusion (red,  for effective porous diffusion), dispersion (brown) and advection (blue) lengths vs. time in arithmetic (a), and (b) and log-log scales (c).Advection dominates sizable (meters and days) field scales (left), diffusion dominates pore scale (center).The frequent practical rule ( = 0.1 L, plotted here for  = 1 m/d) implies a macrodispersion that is unrealistic for mixing.Note the broad range of time and length characteristics scales that can be obtained varying  or .

Figure 3 .
Figure 3.Time evolution of diffusion (red, D mp for effective porous diffusion), dispersion (brown) and advection (blue) lengths vs. time in arithmetic (a), (b) and log-log scales (c).Advection dominates sizable (meters and days) field scales (left), diffusion dominates pore scale (center).The frequent practical rule (α L = 0.1 L, plotted here for v = 1 m/d) implies a macrodispersion that is unrealistic for mixing.Note the broad range of time and length characteristics scales that can be obtained varying α L or q.

Figure 4 .
Figure 4. Example of a gypsum equilibrium chemical system with three species (gypsum,  , and  ). (a) consists of only one row.If gypsum activity is assumed constant, the concentration space,  , and , are reduced to 2 dimensions (red square).Equilibrium restricts  to 1 dimension, both in the  (b),  and log  (c), in log-log scale to highlight linearity in log activity spaces.

Figure 4 .
Figure 4. Example of a gypsum equilibrium chemical system with three species (gypsum, Ca 2+ , and SO 2− 4 ).S e (a) consists of only one row.If gypsum activity is assumed constant, the concentration space, S e , and U, are reduced to 2 dimensions (red square).Equilibrium restricts c to 1 dimension, both in the c (b), u and log a (c), in log-log scale to highlight linearity in log activity spaces.

Figure 5 .
Figure 5. Simplified chemical system for calcite (Cc) and gypsum.The full stoichiometric matrix (blue and green squares) can be simplified (red) if mineral activities are constant (a).Equilibrium imposes constraints on concentrations, which must lie on equilibrium surfaces (b) that become planes in log activities (c).Reactions (red and purple arrows, b) lie on constant  planes (green plane, b).

Figure 5 .
Figure 5. Simplified chemical system for calcite (Cc) and gypsum.The full stoichiometric matrix (blue and green squares) can be simplified (red) if mineral activities are constant (a).Equilibrium imposes constraints on concentrations, which must lie on equilibrium surfaces (b) that become planes in log activities (c).Reactions (red and purple arrows, (b)) lie on constant u planes (green plane, (b)).

Figure 6 .
Figure 6.Transport of pulse injected at x = − 8 m and t = −18 d in a medium consisting of two minerals.A front occurs for A3 at x > 50 m, where it equilibrates with Min. 2, causing its sharp precipitation and an apparent discontinuity in concentrations.The analytical solution (A_r1) [73] is displayed for the precipitation rate of Min. 1. Reaction rates maxima coincide with gradients maxima (black arrows).

Figure 6 .
Figure 6.Transport of pulse injected at x = −8 m and t = −18 d in a medium consisting of two minerals.A front occurs for A 3 at x > 50 m, where it equilibrates with Min. 2, causing its sharp precipitation and an apparent discontinuity in concentrations.The analytical solution (A_r1) [73] is displayed for the precipitation rate of Min. 1. Reaction rates maxima coincide with gradients maxima (black arrows).

Figure 7 .
Figure7.Distribution (at t = 40, red, and 80 d, blue) of salinity (Sal) that enters (Sal = 8) into a medium with low salinity.A second solute (dashed) enters with the same concentration as resident water, but salinity changes from 2 to 1 its retardation, producing a peak at the salinity front.Note that numerical solution of Sal (lines) with a Lagrangian method is visually identical to the analytical solution (symbols).The same numerical method produces some numerical dispersion for a third solute with constant retardation.

Figure 7 .
Figure7.Distribution (at t = 40, red, and 80 d, blue) of salinity (Sal) that enters (Sal = 8) into a medium with low salinity.A second solute (dashed) enters with the same concentration as resident water, but salinity changes from 2 to 1 its retardation, producing a peak at the salinity front.Note that numerical solution of Sal (lines) with a Lagrangian method is visually identical to the analytical solution (symbols).The same numerical method produces some numerical dispersion for a third solute with constant retardation.

Figure 8 .
Figure 8. Dimensionless concentration versus dimensionless distance after a pulse injection of a tracer suffering first order decay (a).The reactions barely affect the solute, if Da < 0.1, but make it disappear for Da > 3.This is not necessarily true for non-linear kinetics (b).Da = 1 implies a sizeable reaction but note that significant amounts of solute remain for Da = 100, if n = 4. Also note the abrupt approach to zero for n < 1, which may cause convergence problems.

3. 3 .
Equilibrium.Conceptual Issue No. 8: The Rate of Fast Reactions Is Controlled by Mixing and Access to Reaction Sites

Figure 9 .
Figure 9.A broad range of scales are involved in RT.Continuum (i.e., Darcy) representations are standard above 1 m scales.Pore scales are becoming frequent below 10 cm scales.Recall that diffusion (double green arrows) dominates sub mm scales, while advection (blue arrows, always spatially variable) dominates large scales.

Figure 9 .
Figure 9.A broad range of scales are involved in RT.Continuum (i.e., Darcy) representations are standard above 1 m scales.Pore scales are becoming frequent below 10 cm scales.Recall that diffusion (double green arrows) dominates sub mm scales, while advection (blue arrows, always spatially variable) dominates large scales.

Figure 10 .
Figure 10.Schematic picture of ADE solution methods.Eulerian methods (a) work on fixed coordinates (accurate but possibly unstable).Lagrangian methods avoid most instabilities by working on moving coordinates, through (b) moving grids or (c) tracking particles (e.g., RW, where dispersion is reproduced by random displacements).Concentration can be obtained from particles locations (d) by binning, or averaging ("Bells" method).

Figure 10 .
Figure 10.Schematic picture of ADE solution methods.Eulerian methods (a) work on fixed coordinates (accurate but possibly unstable).Lagrangian methods avoid most instabilities by working on moving coordinates, through (b) moving grids or (c) tracking particles (e.g., RW, where dispersion is reproduced by random displacements).Concentration can be obtained from particles locations (d) by binning, or averaging ("Bells" method).

Figure 11 .
Figure 11.Effective transport equations can be written in (a) PDE form (multicontinuum with an exchange term, often implicit, expressed in terms of a memory function , or a distribution of immobile regions with residence time 1 α ⁄ ) or (b) RW form (with "anomalous" noise).

Figure 11 .
Figure 11.Effective transport equations can be written in (a) PDE form (multicontinuum with an exchange term, often implicit, expressed in terms of a memory function g, or a distribution of immobile regions with residence time 1/α) or (b) RW form (with "anomalous" noise).
as shown in Figure 12.To us, this is a relevant Conceptual Issue, No. 9: solutes are displaced by advection, spread by dispersion, diluted by mixing.Energies 2022, 15, x FOR PEER REVIEW 23 of 48 proper reproduction of advection, dispersion and mixing, hopefully linked to geological understanding [146,147] as shown in Figure 12.To us, this is a relevant Conceptual Issue, No. 9: solutes are displaced by advection, spread by dispersion, diluted by mixing.

Figure 12 .
Figure 12.Transport through heterogeneous media (a) is characterized by mean advection, dispersion (black dashed arrows, J , in a), which spreads fronts, and molecular diffusion (double yellow arrows), which causes mixing.The ADE (b) fails to distinguish between spreading and mixing, thus overestimating the rates of mixing and fast reactions, even if the solute spatial distribution (vertically integrated during sampling, d) and BTCs (c) look fine.

Figure 12 .
Figure 12.Transport through heterogeneous media (a) is characterized by mean advection, dispersion (black dashed arrows, J D , in (a), which spreads fronts, and molecular diffusion (double yellow arrows), which causes mixing.The ADE (b)) fails to distinguish between spreading and mixing, thus overestimating the rates of mixing and fast reactions, even if the solute spatial distribution (vertically integrated during sampling, (d)) and BTCs (c) look fine.

Figure 13 .
Figure13.Example of the effect of organic carbon oxidation on groundwater chemistry (a and mineralogy (c).Organic carbon is added stepwise and is oxidized assuming equilibrium reactions.Details of the model are given by Appelo and Postma ([51], example 9.6).

Figure 13 .
Figure 13.Example of the effect of organic carbon oxidation on groundwater chemistry (a,b) and mineralogy (c).Organic carbon is added stepwise and is oxidized assuming equilibrium for all reactions.Details of the model are given by Appelo and Postma ([51], example 9.6).

Figure 14 .
Figure 14.Algorithm of the reactive transport methods: (a) Sequential Iteration Approach, (b) Direct Substitution Approach and (c) Water Mixing Approach.

Figure 14 .
Figure 14.Algorithm of the reactive transport methods: (a) Sequential Iteration Approach, (b) Direct Substitution Approach and (c) Water Mixing Approach.

2: diffusion controls microscales; advection controls macroscales; the Peclet number is ambiguous.
This discussion leads to Conceptual Issue No.

Table 1 .
Proposed list of basic RT conceptual issues.