•  Introduction
•  Consistency of a Reaction System
•  Rate Equations and Parameter Determination
•  Reactive Transport and Transport Velocities
•  Numerical Methods Specific to Reactive Transport
•  Conclusion

•  Assess the Consistency of a Reaction System
•  If a system is not consistent, modeling exercise is meaningless.
•  Formulation of Rate Equations - Determination of Parameters.
•  Values of parameter should not contain state variables (p, T, C's, V).
•  Formulation of Transport Equations - Transport Velocities.
•  True transport velocities depend mainly on heterogeneous reactions.
•  Numerical Methods: The most overlooked importance.
•  Methods should be robust - if not, no solution.
•  Methods should be accurate - if not, distort physics.
•  Methods should be efficient - if not, no application to field sites.

("cheap" algorithm on a fine grid vs accurate "expensive" algorithm on a coarse grid)

•  Basic Formulation of Concentration-vs-time curves


•  Diagonalization of Reaction Matrix


•  If the kinetic data (concentration-vs-time curves of all M species) indicates that any of the NC mass conservation equations is violated, then the system is inconsistent: either the number of species identified is too many or the number of linearly independent reactions is too few. The diagonalization is not unique and we should take advantage of the non-uniqueness to explore the consistency of the system.

•  Conventional Approach


•  Diagonal Approach No. 1


•  Diagonal Approach No. 2

• Ad Hoc Approaches

A first- or second-order or any other exotic rate equation is postulated and reaction parameters including rate constants are determined from the experimental analysis of concentration-vs-time plots of a target species. The major problem is that the reaction parameters are not "true" constants.
 
• Mechanistic Approaches

Rate laws are derived from proposed mechanisms (pathways). A mechanism (pathway) consists of a sequence of elementary steps (an elementary step is a reaction that can be described by the forward and backward rates with the order of the rate given by the molecularity of the reaction) that describes how the final products are formed from the initial reactants. A mechanism determines the overall reaction. The major problem is the determination of the pathway is extremely difficult. The major advantage is that the reaction parameters are "true" constants.
 

• Direct simulation of the proposed mechanism - all elementary reactions in the pathway are included: The advantage is that either mass action equations for fast reactions or only one type of rate law for slow reactions needs to be considered and the dis-advantage is that all intermediate species have to be included.

• Indirect simulation of an overall rate law for the proposed mechanism: The advantage is that all intermediate species can be eliminated and the dis-advantage is that one needs to derive a rate law for every reaction.

• Ad Hoc Approach

NOTE: As shall be seen that the first order rate constant in fact contains the effect of the concentration of the electron acceptor [O].
 
• Proposed Mechanism


E+S <--> ES                                                              (r1)
E+O <--> EO                                                             (r2)
ES+O <--> ESO                                                        (r3)
EO+S <--> ESO                                                        (r4)
ESO <--> E+P                                                            (r5)
 

• Indirect simulation of the overall rate - The Mono kinetics


• Conventional approach


• Diagonalized approach


NOTE: The indirect simulation needs to calculate species, [S], [O], and [P] only. This approach depends on our ability to derive the over-all rate law and the validity of simplifying assumptions: Reactions (r1) through (r4) are fast reactions and reaction (r5) is ir-reversible slow reaction with the first order rate.
 
• Direct simulation of the proposed mechanism.

There are 7 species and 4 chemical reactions all of which are linearly independent. The conventional approach is not a good way to simulate a reactive system when mixed kinetic and equilibrium reactions are involved. The diagonalized approach should be preferred even when all the reactions are slow. For the direct simulation, one of the diagonalized approaches would be set up as:


NOTE: The direct approach would require the simulation of all species [S], [O],[E], [ES], [EO], [ESO] and [P]. It is clear that Eqs. (16) and (17) are equivalent.

• Equilibrium constants for all linearly independent equilibrium reactions can bedetermined one reaction by one reaction.
• Reaction rate constants (parameters) for those kinetic reactions that are linearly dependent on equilibrium reactions are irrelevant to the problem.
• Rate constants (parameters) for kinetic reactions that are linearly independentto other kinetic reactions can be determined based on only one curve of concentration rate-vs-time.
• Rate constants (parameters) for those kinetic reactions that are linearly dependent on each other must be determined simultaneously based on a number of curvesof concentration rate-vs-time that are be exactly equal to the number of those reactions.

• Type 0: Nonlinear Algebraic Equations

From Ne Equilibrium Reactions
 

• Type I: Ordinary Differential Equations

From (NI - Ne) transport-reaction equations with those rows that reduced [B] is zero.
 
• Type II: Linear Transport equations

From NC component equations with those rows that the reduced [A] is the same as the reduced [B]. Vt = Vf
 

• Type III: Nonlinear Transport Equations

From NC component equations with those rows that the reduced [A] is different from the reduced [B]. Vt = (C/T)Vf

• Type IV: Linear Transport-Reaction Equations

From (NI - Ne) transport-reaction equations with those rows that the reduced [A] is the same as the reduced [B]. Vt is not easy to define.

• Type V: Nonlinear Linear Transport-Reaction Equations

From (NI - Ne) transport-reaction equations with those rows that the reduced [A] is different from the reduced [B]. Vt is not easy to define.

NOTE: If the system of Eq. (18) is reducible to the system of Type 0, Type I, and Type II, then the transport can be completely decoupled from reactions. Under such circumstances, "true" transport velocities for chemical components are just the fluid velocity and the "true" transport velocities for chemical species are unnecessary to know. We would like to know what kind of properties the reaction matrix [] must have in order to render the system of Eq. (18) reducible. If Type V equations are not created, then the diagonalized approach is always superior to the conventional approach.

NOTE: The conventional approach contains Type I (for surface species) and Type IV (for aqueous species) equations only.

First Decomposition

• A transport velocity can be defined for every equation.

• Not only a complete decoupling has not been accomplished. Furthermore, it creates one undesirable type IV equation.

• Specific Requirements
• No negative concentration is allowed, however small. Preclude oscillatory numerical schemes.
• Solutions should be accurate for large as well as small values.
• High and low concentrations created due to reactions should be captured.
 
• The LEZOOMPC Scheme
• Backward node tracking: to obtain concentrations at global nodes to within the error tolerance.
• Forward node tracking: to determine rough elements for adaptive local grid refinement and to capture peak/valleys for shape preservation.
• Adaptive Local grid refinement: to refine rough elements, one by one, with regularly spanned subelements based on "true error" estimates.
• Peak/Valley Capture: to capture peak/valleys within each subelement so that the shape can be preserved and the loss of peak and valleys can be prevented.
• Recording of forward-node-travel-through elements: to handle boundary sources.

This benchmark problem compares 1st order upwinding with the TVD scheme for nonreactive flow at 45o angle to the grid. Parameters used in the calculation are: a fluid flow velocity of 15.768 m y-1 along the positive x- and y-axes, diffusion coefficient equal to 10-5 cm2 s-1, and 100% porosity. The grid consists of 50 nodes in the x- and y-directions with an equal spacing of 0.2 m. This gives a grid Peclet number of 100. The initial condition consists of a concentration mound of unit concentration with a width of 2 m occupying the square region 1-3 m by 1-3 m corresponding to nodes 5-15 in each direction. The concentration is zero everywhere else. Zero gradient boundary conditions are imposed around the periphery of the domain. A Courant number of 0.1 is used in the TVD calculation. The calculation is carried out for 5 x 106 s.

The result for upwinding are shown in Figure 1 along with the initial condition. In Figure 2 is shown the result of the TVD calculation. Both calculations conserve mass essentially exactly. Clearly the TVD scheme does a much better job compared to the upwinding scheme whcih reapidly disperses the concentration mound with time. The TVD scheme is able to very accurately track the mound without spreading it out as in the upwinding scheme, although some distortion can be seen perpendicular to the direction of flow.

Peter C. Lichtner's Problem Solved with LEZOOMPC

Cases 1A and 1E

Peter C. Lichtner's Problem Solved with LEZOOMPC

Cases 1D and 1H

• Nonuniform Flows

The equations describing the advection-reaction problem are fairly well known [see for example, Bear, 1972]. Given N_M mobile components {c_i}, i = 1, N_M,

Here, V (LT^-1) is the velocity vector; K (LT^-1) is the hydraulic conductivity tensor; P (L) is the potential head; q (T^-1) is the strength of the fluid source/sink, positive for a sink and negative for a source; and c_i (ML^-3) is the concentration of the i-th mobile component. c_i* (ML^-3) is the concentration of the i-th mobile component in a source/sink. It is equal to the specified injection concentration for a source and is equal to c_i for a sink. r_i (ML^-3T^-1) is the rate of production by reaction for the i-th component. is the porosity taken here to be invariant in time.

Parabolic Flow

Consider again a rectangular domain (0 <= x <= Lx, 0 <= y <= Ly), wih an incompressible flow given by,
 
 

V_x(x,y) = m - gx + gy,  V_y(x,y) = -m + gx - gy; m, g > 0
 
• Reactions

Second order, first order, and half order reactions are considered for advection. The reaction expressions are simple so that the rate expressions can be analytically integrated. The following are only a sample. Any number
of other reactions can also be considered. Fell free to pick your favorite reaction system (but need all mobile species).

 
 

A Two-D Nonlinear Reaction Problem under Nonuniform Flows

 
• Parabolic Flow Fields
• Case A Boundary Value Problem
• Case B Initial Value Problem
 

Solution of B. V. Problem for Non-Reactive Tracer by Chilakapati

• Breakthrough between x = 0.5 and 1.0 for injection at y = 0.5 and 1.0

Soluiton of B. V. Problem for Reactive Solute by Chilakapati

• Breakthrough between x = 0.5 and 1.0 for injection at y = 0.5 and 1.0

 

Solution of Initial Value Problem by Chilakapati

Solution of IC and BV Problems for Non-Reactive Tracer with LEZOOMPC

Solution of B.V. Problem for Reactive Solute with LEZOOMPC - 1

Solution of B.V. Problem for Reactive Solute with LEZOOMPC - 2

Fraction of Reaction Mass for I. V. Problem with LEZOOMPC

Fraction of Reaction Mass for B. V. with LEZOOMPC

• Consistency of Reactive System

• Once a reactive system is specified, its consistency should be evaluated with experimental data before rate laws are proposed.

• Formulation of Rate Law and Parameter Determination

• A rate law corresponding to each biogeochemical reaction must be formulated based on a rigorous theory so that the corresponding reaction parameters are free of state variables. The determination of reaction parameters can be achieved reaction by reaction when all the reactions are linearly independent.
 
• Reactive Transport Formulation and Transport Velocities.

• The transport velocity for each of the original set of reactive transport equations is not clearly defined. Thus, it is preferable to reduce the original set of equations to a new set of equations, for which a transport velocity of an entity can be clearly defined for the corresponding equation.

• The original set of transport equations can be decomposed to a new set of transport equations of five types. The decomposition is not unique. The guiding principle is try to generate as many Type I and II equations as possible.
 
• Numerical Methods.

• There are specific requirements of numerical methods in simulating reactive chemical transport.

• The LEZOOMPC algorithm can, accurately and efficiently, simulate both non-reactive and reactive chemical transport. Its robustness has been demonstrated.