Reactive transport modeling involves solving several nonlinear coupled phenomena, among them, the flow of fluid phases, the transport of chemical species and energy, and chemical reactions. There are different ways to consider this coupling that might be more or less suitable depending on the nature of the problem to be solved. In this paper we acknowledge the importance of flexibility on reactive transport codes and how object-oriented programming can facilitate this feature. We present PROOST, an object-oriented code that allows solving reactive transport problems considering different coupling approaches. The code main classes and their interactions are presented. PROOST performance is illustrated by the resolution of a multiphase reactive transport problem where geochemistry affects hydrodynamic processes.

INTRODUCTION

Reactive transport models are tools that help to understand the hydraulic and chemical behavior of natural and artificial porous media. They have been used to solve a broad range of problems, such as groundwater remediation (Loomer et al. 2010), nuclear waste disposal (MacQuarrie & Mayer 2005) and CO2 sequestration (Zhang et al. 2012) among others, from micro-scale (Trebotich et al. 2014) to field-scale problems (Sassen et al. 2012).

Modeling reactive transport in porous media involves simulating several coupled phenomena: phase flow, solute transport, and reactions. It may also involve multiphase flow, heat transport, and porous media deformation (Steefel et al. 2015). These phenomena may be complex to model individually, and modeling together brings new difficulties associated with coupled effects (Lichtner 1996). Which coupled effects have to be considered and the optimal solution strategy for the coupled equations depend on the nature of the problem to be solved and may vary significantly from case to case (Zhang et al. 2012).

The ideal reactive transport code would have to use an accurate, robust, and efficient numerical approach. However, it is difficult to obtain these goals with a single numerical approach. Therefore, concessions have to be made and different coupling alternatives have to be chosen at different levels. Numerical accuracy is generally preferred above other issues when solving modeling research applications. On the other hand, when solving field-scale problems, efficiency and robustness have priority while accuracy remains within the bounds of the uncertainty associated with model parameters (Yeh et al. 2012).

Two big families of methods were addressed to account for the coupling between solute transport and chemical reaction processes: (1) the operator splitting (or sequential iterative (or non-iterative) approach and (2) the global implicit or direct substitution approach (Saaltink et al. 2001). As regards the first one (i.e., the sequential methods), whether iterative (SIA) or not (SNIA), operator splitting techniques should be adopted that effectively decouple component transport equations. As regards the last one, direct substitution approaches (DSA) solve both transport and chemical reactions simultaneously. A number of authors have studied the numerical performance of these methods (Steefel & MacQuarrie 1996), and they conclude that in spite of the fact that the DSA is more accurate and robust, there are cases where the SIA is more convenient from an efficiency–accuracy point of view. In addition, SNIA may be appropriated for scenarios with Courant number smaller than 1 (Xu et al. 2012). Some reactive transport codes are able to work with both of these approaches (CRUNCHFLOW, Steefel 2009; DUMUX, Flemisch et al. 2011; HYDROGEOCHEM, Yeh et al. 2010; PFLOTRAN, Lichtner et al. 2013; RETRASOCODEBRIGHT, Saaltink et al. 2004), while others use the fully implicit approach (NUFT, Hao et al. 2012; MIN3P, Mayer et al. 2012), or different variants of operator splitting techniques (CORE, Samper et al. 2009; HYDRUS-PHREEQC (HP1), Jacques et al. 2011; HYTEC, Lagneau & Van Der Lee 2010; IPARS, Wheeler et al. 2012; OPENGEOSYS, Li et al. 2014; ORCHESTRA, Meeussen 2003; PHAST, Parkhurst et al. 2004; PHREEQC, Parkhurst & Appelo 2013; PHT3D, Prommer & Post 2010; RT3D, Johnson & Truex 2006; STOMP, White & McGrail 2005; TOUGHREACT, Xu et al. 2011).

On a more complex level is the coupling between phase conservation and reactive solute transport. Most reactive transport codes decouple phase conservation (i.e., flow equation) from reactive transport calculations (RT3D, MIN3P, PFLOTRAN, PHAST, RETRASOCODEBRIGHT, HYTEC, TOUGHREACT).

This approach is convenient in most cases, but a numerically coupled solution will generally be more suitable when the phenomena involved are highly physically coupled. One example of this could be found in problems related to the CO2 sequestration in brine aquifers, which has prompted the development of codes that solve coupled multiphase flow and reactive transport (Fan et al. 2012) and even mechanical deformation (Zhang et al. 2012). Likewise, Wissmeier & Barry (2008) showed that the consumption of water due to hydrated mineral precipitation can have impacts on flow and solute transport for unsaturated flow problems. These impacts can be even more important if gas transport is also considered because water activity, which controls vapor pressure, is affected by capillary and osmotic effects. Moreover, certain mineral paragenesis can fix water activity (producing an invariant point), causing the geochemistry to control vapor pressure, which is the key variable for vapor flow (Risacher & Clement 2001). In such cases, decoupling is not appropriated. Formulations that are able to represent these effects are complex to implement since they should consider all coupled phenomena and a variable number of components in space and time.

While most reactive transport codes consider a single technique for the resolution of the partial differential equation some codes can adopt more than one. In Table 1 the supported discretization method and coupling strategies for different reactive transport codes are detailed.

Table 1

Supported discretization method and coupling strategies for different reactive transport codes

Code PDE discretization1 Transport and reaction coupling2 Phase conservation and transport coupling3 
CORE FEM OS SEQ 
CRUNCHFLOW FVM OS, DS COU 
DUMUX FEM, FVM, MFDM OS, DS COU, SEQ 
HYDROGEOCHEM FEM, MMC OS, DS ITER 
HYDRUS-PHREEQC (HP1) FEM OS SEQ 
HYTEC FVM OS SEQ 
IPARS MFEM, DGM OS INDP 
NUFT FVM DS SEQ 
MIN3P FVM DS SEQ 
OPENGEOSYS FEM OS SEQ 
ORCHESTRA MC OS INDP 
PFLOTRAN FVM OS, DS SEQ 
PHAST FDM OS INDP 
PHREEQC MC OS INDP 
PHT3D FVM, MMC OS INDP 
RETRASOCODEBRIGHT FEM OS, DS SEQ, ITER 
RT3D FDM OS INDP 
STOMP FEM OS SEQ 
TOUGHREACT FVM OS SEQ 
Code PDE discretization1 Transport and reaction coupling2 Phase conservation and transport coupling3 
CORE FEM OS SEQ 
CRUNCHFLOW FVM OS, DS COU 
DUMUX FEM, FVM, MFDM OS, DS COU, SEQ 
HYDROGEOCHEM FEM, MMC OS, DS ITER 
HYDRUS-PHREEQC (HP1) FEM OS SEQ 
HYTEC FVM OS SEQ 
IPARS MFEM, DGM OS INDP 
NUFT FVM DS SEQ 
MIN3P FVM DS SEQ 
OPENGEOSYS FEM OS SEQ 
ORCHESTRA MC OS INDP 
PFLOTRAN FVM OS, DS SEQ 
PHAST FDM OS INDP 
PHREEQC MC OS INDP 
PHT3D FVM, MMC OS INDP 
RETRASOCODEBRIGHT FEM OS, DS SEQ, ITER 
RT3D FDM OS INDP 
STOMP FEM OS SEQ 
TOUGHREACT FVM OS SEQ 

1DGM, discontinuous Galerkin method; FDM, finite difference method; FEM, finite elements method, FVM, finite volume method; MC, mixing cell; MFDM, mimetic finite difference method; MFEM, mixed finite element method; MMC, modified method of characteristic.

2DS, direct substitution; OS, operator split.

3COU, coupled; INDP, independent; ITER, iterative; SEQ, sequentially.

Reactive transport modeling in fractured media might also require flexibility regarding the way the medium is considered. Important changes in fluid pressures and solute concentrations will propagate rapidly through the fracture system, while exchanges with the matrix blocks will occur slowly. To account for this, some reactive transport codes have included multiple interacting continua modeling (TOUGHREACT, PFLOTRAN).

In short, for reactive transport modeling the adopted coupling techniques, the partial differential equation discretization method, and the way the domain is considered, may be problem dependent. Therefore, a reactive transport code should include several solution approaches to be used in a broad range of problems. Moreover, in order to ensure its use for present and future problems, it must have an extensible design. A number of authors have pointed out that object-oriented (OO) programming facilitates the implementation of these features (Filho & Devloo 1991; Commend & Zimmermann 2001).

The scientific community has been adopting OO techniques for problem solving since the end of the last century (Forde et al. 1990; Wang & Kolditz 2007; Slooten et al. 2010). However, only in the last decade have OO codes been developed for reactive transport modeling. Meysman et al. (2003) developed an OO reactive transport code for a single fluid phase. Gandy & Younger (2007) developed an OO multiphase reactive transport code for pyrite oxidation and pollutant transport in tailing ponds. Shao et al. (2009) included reactive transport calculations into a thermo-hydro-mechanic OO framework adopting a sequential noniterative approach (SNIA). Bea et al. (2009) developed an OO module capable of solving reactive transport for a single phase considering the SNIA, SIA, or DSA approach. However, all of these codes, and most of the procedural reactive transport codes, have a predefined strategy for dealing with coupling effects. Particularly, they do not allow for changing number and definitions of chemical components when solving flow and reactive transport in a coupled way.

The objective of this paper is to present an OO structure for reactive transport that can accommodate different levels of physical and chemical processes coupling. The structure presented here is capable of modeling from single-phase SIA problems to fully coupled multiphase reactive transport problems. In addition, to the best of our knowledge, it is the first OO tool capable of considering the occurrence of invariant points (e.g., for reference see Risacher & Clement 2001) in a reactive transport problem. This is an extreme case where geochemical processes significantly affect fluid flow and the number and definitions of chemical components may vary significantly in space and time. This structure has been implemented in PROOST which was programmed in FORTRAN 95 following the OO paradigm, and until now could solve single phase reactive transport by the SIA method and a fully coupled multi-phase reactive transport by the DSA method.

EQUATIONS TO SOLVE

Reactive transport modeling implies establishing several conservation principles, like mass or energy conservation expressed as partial differential equations (PDEs), and several constitutive and thermodynamic laws (such as retention curve or mass actions laws) expressed as algebraic equations (AEs). Darcy's law is used to represent momentum conservation. In this section we present a generic conservation equation to represent conservation principles in reactive transport problems. We consider in detail the species and component conservation and we briefly present the constitutive and thermodynamic laws.

General conservation equation

Conservation of a physical entity can be expressed as: 
formula
1
where is the amount of per unit volume of medium, is the flux of due to the driving force (e.g., advection or diffusion), and is a sink source term. Since time and spatial derivatives are involved, conservation equations usually take the form of a PDE.

Species and component conservation equation

The conservation of a species i belonging to phase , which is a particular case of Equation (1) has the following expression: 
formula
2
where is the volumetric content of phase , is the species i concentration in phase, is the stoichiometric coefficient of the equilibrium reaction j for the species i, is the reaction rate of the equilibrium reaction j, and Ne is the number of equilibrium reactions. , and Nk are analogous to , and Ne but for kinetic reactions. is an external sink-source term, and is the linear transport operator for the mobile phase involving advective and diffusive-dispersive processes: 
formula
3
Mobile phase fluxes are calculated according to Darcy's law: 
formula
4
where , and are the conductivity tensor, pressure and density of the phase , respectively. Diffusive-dispersive fluxes are calculated according to Fick's law: 
formula
5
where and are the diffusion and dispersion tensor for phase , respectively, and is the tortuosity.
Note that the general sink source term of Equation (1), , involves several different terms in Equation (2): 
formula
6
There is no explicit expression for the equilibrium reaction rates , their value has to be such that the corresponding mass action law is satisfied. Therefore, values can be written as a function of both transport and chemical processes (DeSimoni et al. 2005). A common approach to avoid dealing with these terms is to formulate the conservation of components as a linear combination of species that remain unaffected by equilibrium reactions. As such, equilibrium reactive rates disappear from the conservation equations of components (Steefel & MacQuarrie 1996). However, components may involve species belonging to different phases, therefore conservation equation for components have to be written: 
formula
7
where and are the i component concentration in mobile phases and immobile phases , respectively, and is a linear combination of the kinetic terms that affect the species composing the component. We consider as immobile phases minerals and fluid–solid interface, despite the fact an interphase is not a phase from a thermodynamic point of view. Note that the component conservation Equation (7) has the same structure as Equation (2). The main difference is that a component may be present in more than one phase, while a species belongs to a single phase. There are several ways of defining components and therefore some freedom in the choice of components. This has led to formulations that try to define components that do not affect each other, such as those proposed by Molins et al. (2004), Kräutle & Knabner (2005) and Hoffmann et al. (2010). Saaltink et al. (1998) introduced a definition that eliminates species whose activities are known and constant. That is the case of minerals, that are considered as pure phases, so that their activity equals unity. Also, the activity of water can be assumed unity for the case of diluted solutions. Minerals, often considered as constant activity species, might appear or disappear from portions of the domain due to precipitation–dissolution processes. Therefore, under equilibrium assumption, the dimension of the component vector, the number of components, may be different at each discrete point in space and vary in time. This increases the difficulty of solving Equation (7) since the matrix system to be solved has a dynamic size, which significantly affects the code.

Once all component conservation and geochemical equations have been solved, all species concentrations are known. Equilibrium reaction rates are then calculated from species conservation Equation (2). If constant activity species have been eliminated from the component definition, their concentration must also be calculated from Equation (2).

Constitutive and thermodynamic laws

The literature provides several models for density, viscosity, and diffusion coefficients of mobile phases. These parameters are usually expressed as an explicit function of phase composition, pressure, and temperature. Several models express saturation and relative permeability as an explicit function of capillary pressure and surface tension. All these relations lead to a local system of equations, which is valid at every point of the domain.

Thermodynamic relations also form part of this local system of equations. The most important of these are the chemical equilibrium reactions, which may be expressed by means of mass action laws, as often done in reactive transport. Also required are models for the calculation of activity, such as Debye-Hückel (1923) or Pitzer (1973), and expressions for kinetic rate laws (such as Monod or Lasaga) (Mayer et al. 2002).

Minor changes on the solid matrix, like porosity changes due to mineral dissolution–precipitation or clogging, may also be expressed as algebraic equations (Soleimani et al. 2009). More complex mechanical processes, like deformation or consolidation, involve momentum conservation equation and have to be solved as a PDE (Villar et al. 2008).

Constitutive and thermodynamic relationships define a set of AEs that have to be solved together with the conservation equations (PDEs).

Numerical solution of the equations

Methods such as finite element or finite differences, among others, are normally used to approximate time or space derivative terms in PDEs. Application of such methods leads to a set of equations that represent the conservation principle for discrete portions of the domain (representing nodes or cells). The current version of PROOST supports two methods: the finite elements and the mixed finite elements. Contrary to constitutive or thermodynamic laws, these equations are not local, that is, equations at a discrete point are a function of variables at other discrete points. As constitutive and thermodynamic models (AEs) involve variables that appear in the PDE, both AEs and PDEs may have to be solved simultaneously. Generally, the resulting set of equations is nonlinear, which makes their solution more difficult. As mentioned in the Introduction, different approaches can be adopted for solving these coupled sets of equation: independently, sequentially, iteratively, or coupled.

OO ANALYSIS OF REACTIVE TRANSPORT MODELING AND PROOST CLASS ORGANIZATION

According to the OO philosophy, the numerical solution of reactive transport can be represented by a group of interacting objects. These objects belong to classes which define common types of data and functionality. According to Filho & Devloo (1991), defining suitable classes is the first and perhaps the most important step in software design under OO.

Our analysis was based on the following abstraction: reactive transport modeling is considered as a set of equations (PDEs and AEs), representing the conservation of chemical species, that need to be solved in a certain domain. These equations involve several variables or fields (such as concentrations, density, or porosity) which are also defined over portions of the same domain. The domain is discretized and fields are defined over the discretized space (nodes or cells). Using discretization techniques (such as finite element or finite differences methods), PDEs are converted into a set of algebraic equations which represent a discrete version of the PDE. For each discretized time interval, this set of equations can simultaneously be solved with the AE or using an operator splitting technique.

The above description points to a natural class structure for our problems. The PDEs share attributes such as terms in the equation, state variables, or domain definitions, and also share functionalities such as computing the balance or the matrices for the discretized PDE. Therefore, we find it natural to define a class, termed Phenomenon, to identify PDEs. In the same fashion, we define Process as the class whose instances will be specific terms in the PDE (e.g., advection, dispersion, etc.). The class Meshfields defines objects representing various properties defined over space (and time). To deal with the geochemical processes we use the class CHEPROO (CHEmical PRocesses Object Oriented, Bea et al. 2009). All these objects produce the terms for the (nonlinear) discretized PDEs, which are solved with the functions of the class Solver. The class organization described above is shown in Figure 1 and its detailed description is given below.
Figure 1

Organization of the main classes in PROOST. Each box represents a class with its attributes and methods. A paradigm is shown below each class.

Figure 1

Organization of the main classes in PROOST. Each box represents a class with its attributes and methods. A paradigm is shown below each class.

Phenomenon class

PDEs are a central ingredient of reactive transport modeling. All PDEs represent a conservation principle. All of them consist of different terms, like storage, flux divergence, or source terms and are subject to initial and boundary conditions. Therefore, we define a class for representing PDEs. We term this class Phenomenon. Note that a number of authors have also defined similar classes in their analysis (Meysman et al. 2003; Boivin & Ollivier-Gooch 2004; Kolditz & Bauer 2004). But the main difference here, is that in our case, the Phenomenon object is composed of several objects of the class Process which represent the different terms of the PDE. This is a key aspect that facilitates code reuse, as will be shown below in the Process class description.

Beside the Processes that define the PDE, the initial and boundary conditions are also the main attributes of the Phenomenon class. Methods include the computations of balances or the contribution to matrices comprising the discrete version of the PDEs. The values of the solution variables or unknowns will be obtained from the solution of this matrix system.

Initial conditions and Dirichlet boundary conditions are defined as a Meshfields and are handled by the Phenomenon class. The rest of the boundary conditions, that can be expressed as different terms of the PDE, are represented by instances of the Process class. The Dirichlet boundary conditions have the particularity of imposing the state variable value over different parts of the domain. For this reason they are handled directly by the Phenomenon.

A Phenomenon object can be used to represent a single conservation principle (such as species mass or energy) or several conservation principles with similar equations, like components concentrations. For this latter case, the Phenomenon class makes use of the fact that the same conservation equation applies to all components, and therefore only one PDE has to be defined which applies to all components.

Process class

The terms that compose the PDE (e.g., storage or advection) and the boundary conditions that constrain it (e.g., leakage) are represented by the Process class. The actual nature of this term is defined via inheritance by specialization classes (Figure 1).

The main attributes of this class are the time and space where the Process is applied (e.g., the location of a pumping well for a sink-source Process) and the fields it involves (the pumping rate in this example). Methods include the computation of the process contribution to the system matrix or to the global balance. All these are performed by using methods of the class Mesh, where all discretization-integration information and methods are encapsulated.

The Processes objects are the terms that constitute the conservation equations. A Phenomenon can be formulated by combining different Processes. This class facilitates the extensibility of the code because only the new terms (new specialization of the class Process) have to be programmed to extend the set of equations that can be solved. It also allows reusing code, since the same type of Processes can be used for different conservation equations. For example, a diffusive process for a mass conservation equation and an energy conservation equation are different instances of the same class. Another example is the extension of the component conservation equation from single phase to multiphase (Equation (7)). For this case, all Processes have to be replicated for each mobile phase. This can be easily done by considering new instances of the same Process objects.

There are certain limitations regarding the kind of Processes that can be added to a Phenomenon, and are related to the numerical method chosen for solving it. The nature of the considered Process has to be supported by the numerical method. For example, in its current implementation, advective terms cannot be considered when solving a PDE with the mixed finite element method.

Most boundary conditions are represented with objects of the class Process. Imposed fluxes and variable dependent fluxes are considered through Sink-Source objects, which are specialization of the class Process. As mentioned before, Dirichlet boundary conditions are handled by the Phenomenon class.

Mesh class

There are different techniques to solve PDEs numerically. All these techniques share an approach for discretizing the spatial domain (such as nodes, elements, or cells) and methods to integrate (or differentiate) the terms (Process) of a PDE (Phenomenon) to produce a matrix system from which the discrete solution of the PDE can be obtained.

Thus, all the data and functionality regarding spatial discretization and the discretization–integration methods for solving PDE (such as finite element or finite differences) define a class that we term Mesh. A number of authors have defined similar classes in their analysis. However, most of them separate the domain discretization from the integration methods in different classes (Commend & Zimmermann 2001; Wang & Kolditz 2007). This integration was made because, despite the fact that both methods can share a mesh (elements and nodes), the mesh topological data required might be different. For example, mixed finite elements and finite elements can both use the same mesh, but mixed finite elements need extra information about edges for 2D problems or sides for 3D. Another difference between these two methods is that while finite elements give a continuum scalar field for the solution over the mesh, mixed finite elements give a vector field. Therefore, some aspects of the spatial discretization are related to the integration method, and that is why both are considered a single object in PROOST.

The main attributes of the Mesh class are the domain discretization information (such as nodes or cell coordinates and connectivity between these discrete elements). Methods include yielding information of space discretization (such as the number of discrete elements and their geometrical information), integrating the different terms of the conservation equation (Processes) over the domain, and evaluating spatial properties of variables such as gradients.

The Mesh class allows incorporating new discretization–integration numerical methods by adding new specializations of the class. Two specializations of the class Mesh are currently implemented in PROOST: the finite elements and the mixed finite elements.

Meshfield class

Another important element of reactive transport modeling is the AEs that represent constitutive and thermodynamic laws. Constitutive laws express one field as a function of others. Thus a class termed Meshfield is defined to represent the projection of different scalar, vector, or tensor fields (such as pressure, flux, or conductivity) in the discrete domain. The main attributes of this class are the values and derivatives of a field for the discrete entities (nodes, elements, or cells) and the parameters of the function or constitutive laws they represent. The main methods of the class are to calculate its values and derivatives, and to interpolate its values over any point of the domain. Among others, Meshfield is used to represent retention curves, relative permeability curves, and dispersion coefficients.

For example, a flux Meshfield object defined as , can calculate its values and its derivatives to transmissivity T and head h fields. When a Meshfield represents one of the solution variables of the problems, like head in the previous example, its values are set by the solver class.

The Meshfield class facilitates code extension since new constitutive laws can be easily added to the code by creating new specializations.

CHEPROO class

Geochemical calculations for the component concentrations and kinetic rate laws of Equation (4) are, in fact, constitutive laws. Hence, we treat them as a specialization of a Meshfield, which we term Chemical Meshfield. Many geochemical variables affect the evolution of the system but do not appear explicitly in any PDE (e.g., the activity of aqueous species). For this reason and also because of the complexity of some geochemical calculations, all geochemical models and computations are encapsulated into a single object of a class termed CHEPROO. Only the chemical variables that appear in PDE (such as component concentration or density) are stored in a Chemical Meshfield.

The CHEPROO class uses a module with the same name, with an internal class hierarchy including classes like species, phase, and reaction (Bea et al. 2009). CHEPROO attributes include the geochemical models, such as those for activity coefficients, density or kinetic rates laws, and the chemical data associated with each discrete point of the problem, such as concentrations or components definition. CHEPROO includes methods for calculating the values and derivatives of chemical variables (like component concentration) with respect to the solution variables of the PDE, and to dump them into Chemical Meshfield objects.

CHEPROO also controls the number of chemical components. For some formulations, like that of Saaltink et al. (1998), the number of components may change in time and space. Thus, CHEPROO has to provide information about the components in order to establish the dimension of the final matrix system to be solved.

Solver class

A coupling strategy (coupled or decoupled) needs to be chosen when solving several PDEs. A solution technique for nonlinear systems (Newton-Raphson or Picard) is also needed. An object of the Solver class will be in charge of solving a number of PDEs with a chosen solution strategy:

  • Independently, there are no crossed influences between Phenomena (for example, changes on porosity due to chemical changes are not considered when solving fluid phase conservation).

  • Sequentially, influences between Phenomena are considered lagged in time (for the porosity example, changes due to chemistry in time t are considered for flow in time t+dt).

  • Iteratively, all Phenomena are alternately solved until no significant changes on linking variables occurs (for the porosity example, flow and transport are solved alternately until no significant changes in porosity occurs).

Coupled, all Phenomena are solved at once. Solver attributes include the set of Phenomena, the coupling strategy, the time discretization parameters or the convergence criteria. Methods are required for assembling and solving the discretized PDE system, for time integration. To address these, Solver uses other classes. For instance, matrix systems are handled by a class termed Matrix that encapsulates matrix data and solution techniques for linear systems.

Solver is the class that contributes most to the flexibility of the code since it can be used to solve several conservation equations following different strategies. For example, it might be used to solve first a steady state phase conservation equation (for phase flow calculation) and then a transient component conservation. Or it can be used to solve simultaneously the component and energy conservation.

Component conservation Phenomenon for the SIA and DSA approach

Despite the fact that the SIA and DSA are two approaches for solving the same Phenomenon (the component conservation equation), the way this Phenomenon is formulated in PROOST depends on the chosen approach.

When solving component conservation equations with the DSA approach, the input Phenomenon for PROOST should be the same as in Equation (7). However, for the SIA approach, immobile species storage and kinetic reactions are treated as a sink source term: 
formula
8
Thus the component conservation equation is written only in terms of mobile component conservation: 
formula
9
The PROOST class organization allowed implementing the SIA method without many modifications. The SIA sink source term was represented with the pre-existing sink source Process class. This process evaluates the values of the sink source term, which are given by a Meshfield, and calculates its contribution to the discretized PDE system. By doing this, all the complexity of this term is encapsulated in the class CHEPROO, which sets the values of the SIA source term in a Chemical Meshfield.

SOLUTION PROCEDURE SCHEME FOR A TIME STEP

The interaction between PROOST objects can be illustrated by the solution of a time interval for a reactive transport problem considering the DSA method. The flow diagram is shown in Figure 2, from which 15 relevant points have been identified.
Figure 2

Flow diagram of a time interval resolution for a reactive transport problem in PROOST. 1. Solver establishes the size of the matrix system to be solved. This size depends on the number of coupled phenomena and the dimension of each state variable. Recall that component conservation dimension can be different for each discrete point and may change among the iterative process. 2. Solver assembles the matrix system to be solved. To this end, Solver requests each Phenomenon for its contribution. 3. Phenomenon requests the contribution of all its Processes. 4. Each Processes requests the values of all the Meshfields to which it is related. 5. Meshfield computes its values. 6. CHEPROO calculates Chemical Meshfield values. 7. Mesh computes the contribution of the Process to the matrix system. 8. Matrix solves the matrix system. 9. Solver updates the calculated solution variables (concentrations, temperature, or pressures) in CHEPROO. 10. CHEPROO calculates the concentration of all species from these values (speciation). If there are significant changes on chemical composition, like complete dissolution of minerals in equilibrium or precipitation of new ones, geochemical calculation might not converge. If that is the case, the length of time interval is reduced and the resolution procedure is restarted. The user sets the ideal time step, but if the resolution of the matrix system (which goes from step 2 to 11) exceeds a certain number of iterations, also set by the user, the time step is reduced. 11. Solver controls the convergence of the PDEs linearization and resolution process. When convergence is reached all variables involved in the Phenomena are known, except equilibrium reactions rates that were eliminated when solving component conservation (Equation (7)). These rates can be calculated from the species conservation equations (Equation (2)). In order to avoid the formality of formulating both components and species Phenomenon, this is done by considering an alternative component definition; each mobile species is considered a component. Therefore, the result of the balance of the new component conservation will be the product of the stoichiometric coefficient and the equilibrium reaction rates. These aspects are illustrated with an example in the next section. 12. CHEPROO changes component definition (each mobile species is considered a component). This step allows solving species conservation equations with the same structure used for component conservation equations. This is one of the advantages of PROOST class organization. More details on this particular aspect will be given in the section ‘Application’ below. 13. Phenomenon computes balance (similar to step 3). 14. CHEPROO calculates equilibrium reaction rates from Phenomenon balance. Some reactive transport formulations, like the one of Saaltink et al. (1998), eliminate constant activity species, like minerals, from component composition. These species concentration can be calculated once the equilibrium reaction rates are known. 15. If the formulation considered eliminates constant activity species, the number of components is affected by the disappearance or appearance of minerals. Therefore, component definition has to be controlled after the eliminated species are calculated. If component definition changes the resolution procedure has to be started for the new definition, if not the resolution procedure for the time step is finished.

Figure 2

Flow diagram of a time interval resolution for a reactive transport problem in PROOST. 1. Solver establishes the size of the matrix system to be solved. This size depends on the number of coupled phenomena and the dimension of each state variable. Recall that component conservation dimension can be different for each discrete point and may change among the iterative process. 2. Solver assembles the matrix system to be solved. To this end, Solver requests each Phenomenon for its contribution. 3. Phenomenon requests the contribution of all its Processes. 4. Each Processes requests the values of all the Meshfields to which it is related. 5. Meshfield computes its values. 6. CHEPROO calculates Chemical Meshfield values. 7. Mesh computes the contribution of the Process to the matrix system. 8. Matrix solves the matrix system. 9. Solver updates the calculated solution variables (concentrations, temperature, or pressures) in CHEPROO. 10. CHEPROO calculates the concentration of all species from these values (speciation). If there are significant changes on chemical composition, like complete dissolution of minerals in equilibrium or precipitation of new ones, geochemical calculation might not converge. If that is the case, the length of time interval is reduced and the resolution procedure is restarted. The user sets the ideal time step, but if the resolution of the matrix system (which goes from step 2 to 11) exceeds a certain number of iterations, also set by the user, the time step is reduced. 11. Solver controls the convergence of the PDEs linearization and resolution process. When convergence is reached all variables involved in the Phenomena are known, except equilibrium reactions rates that were eliminated when solving component conservation (Equation (7)). These rates can be calculated from the species conservation equations (Equation (2)). In order to avoid the formality of formulating both components and species Phenomenon, this is done by considering an alternative component definition; each mobile species is considered a component. Therefore, the result of the balance of the new component conservation will be the product of the stoichiometric coefficient and the equilibrium reaction rates. These aspects are illustrated with an example in the next section. 12. CHEPROO changes component definition (each mobile species is considered a component). This step allows solving species conservation equations with the same structure used for component conservation equations. This is one of the advantages of PROOST class organization. More details on this particular aspect will be given in the section ‘Application’ below. 13. Phenomenon computes balance (similar to step 3). 14. CHEPROO calculates equilibrium reaction rates from Phenomenon balance. Some reactive transport formulations, like the one of Saaltink et al. (1998), eliminate constant activity species, like minerals, from component composition. These species concentration can be calculated once the equilibrium reaction rates are known. 15. If the formulation considered eliminates constant activity species, the number of components is affected by the disappearance or appearance of minerals. Therefore, component definition has to be controlled after the eliminated species are calculated. If component definition changes the resolution procedure has to be started for the new definition, if not the resolution procedure for the time step is finished.

CODE IMPLEMENTATION

The code presented the results from merging and expanding two existing codes: PROOST and CHEPROO. The original design of PROOST was already capable of solving different phenomenon, in a coupled or decoupled way, by considering different techniques for the resolution of nonlinear systems (such as Newton-Raphson or Picard). However, such a design only allowed solving Phenomenon objects that had one scalar field as unknown. Also Phenomenon Processes had to be written explicitly as a function of the unknown variable. These featured clashed with the resolution of component conservation, especially when the DSA approach is considered.

The solution of component conservation equations involves considering the conservation equation of several components. As the number of components and its definition might change in time and space (because of complete dissolution or appearance of new mineral species), the number of Phenomenon considered would also have to vary. In order to avoid this difficulty, and as the same Processes affect all component concentrations, only one Phenomenon is considered which applies to a vector variable: the component concentration vector. Therefore, Phenomenon and Process classes were expanded to handle a vector variable whose size may change in time and space.

Processes were originally designed to represent terms of PDEs that directly involve the unknowns of the problem (i.e., main state variables of the Phenomenon: pressure for flow equation, concentration for transport equation). For example, all Processes in a conservative transport problem involve the solute concentration variable, which is also the unknown of this problem. When solving reactive transport by the DSA method, Processes are formulated in terms of component concentrations, but the unknowns of the problem are the primary species concentrations. Therefore, Phenomenon and Process classes were expanded so they can be formulated in terms of any variable and not necessarily the unknown.

Originally, CHEPROO was capable of solving single phase reactive transport problems. CHEPROO uses a matrix system calculated by another conservative transport code to formulate and solve the reactive transport problem (Bea et al. 2009). In order to take advantage of PROOST's flexibility we choose to formulate and solve the multiphase reactive transport equations in PROOST instead of CHEPROO. Therefore, CHEPROO was added to the PROOST structure with the only purpose of performing the chemical calculations (speciation) and provide geochemical variables values and derivatives.

Besides adding new services to make chemical variables available outside its module, several improvements were made in CHEPROO. Phase properties like density, viscosity, and enthalpy, and capillary effect on water activity were added. The PROOST class organization allowed representing all these chemical variables in the class Chemical meshfield. By doing this, all the work related to the evaluation, update, and dependency of these fields to others (like pressure or temperature) is done by pre-existing methods.

Also a new speciation algorithm that uses the Newton-Raphson method had to be programmed in CHEPROO due to the high nonlinearity of concentrated solutions. CHEPROO and PROOST were programmed in FORTRAN 95 following the OO paradigm. This language was chosen for its high popularity among hydrogeologists and its excellent performance reputation. Even though FORTRAN is not a full object-oriented language it can directly support many of the important concepts of OO programming. Details about OOP concepts in FORTRAN can be found in Maley et al. (1996), Decyk et al. (1998), Norton et al. (1998), Akin (1999), Carr (1999), and Gorelik (2004).

APPLICATION

Application description

In order to illustrate the classes introduced before, some aspects of the solution procedure scheme for a time step (generically described in the section ‘Solution procedure scheme for a time step’) are shown for a concrete application. We present the modeling of a 24 cm column of porous gypsum subjected to a constant source of heat, in which significant evaporation occurs. We will focus on the component conservation equation. This synthetic example was designed for illustrating the interaction between hydrodynamic and geochemical processes and it is described in detail by Gamazo et al. (2012). Owing to this interaction, a compositional formulation was adopted and therefore no phase conservation equations are explicitly solved. The finite element method was used for the spatial discretization. One of the most interesting aspects of the application is how the equilibrium reaction rates were calculated. This implies solving a different conservation equation, species conservation instead of components. The PROOST structure allowed calculation of the equilibrium reaction rates by using pre-existing methods.

This application example includes gypsum, liquid water, and vapor, dissolved and gaseous air, calcium and sulfate (main components of gypsum besides water) and two conservative species, potassium and chloride (see Table 2). It also considers the occurrence of anhydrite, which may precipitate as a result of gypsum dehydration. Note that the coexistence in equilibrium of anhydrite and gypsum can fix water activity and therefore produce invariant points (Risacher & Clement 2001). As was shown above, to our knowledge, PROOST is the first multiphase reactive transport capable of modeling this scenario.

Table 2

Chemical species and reactions considered

H2O(l), H2O(g), air(l), air(g), Ca, SO4, K, Cl, gypsum, anhydrite 
 
 
 
 
H2O(l), H2O(g), air(l), air(g), Ca, SO4, K, Cl, gypsum, anhydrite 
 
 
 
 

Initial component definition

Pore solution was initially considered in equilibrium with gypsum with a mineral volumetric content of 0.6. The incoming energy heats the column, which increases evaporation and reduces saturation degree at the top. This induces an ascending non statured flow of liquid water. At a certain point a descending evaporation front appears followed by an also descending gypsum dehydration front in which anhydrite precipitates (see Figure 3). Note that this second front has a significant effect on vapor flow.
Figure 3

Liquid saturation, mineral mass, evaporation rate, and vapor flux for the upper 8 cm of the column. Note that besides the typical evaporation front associated with the drying front there is a second evaporation front associated with hydrated mineral dissolution. This second front has a significant effect on vapor flow.

Figure 3

Liquid saturation, mineral mass, evaporation rate, and vapor flux for the upper 8 cm of the column. Note that besides the typical evaporation front associated with the drying front there is a second evaporation front associated with hydrated mineral dissolution. This second front has a significant effect on vapor flow.

When the simulation starts, the whole domain has the same mineral composition and therefore the component conservation equations for all nodes are the same (see Table 3 for component definition): 
formula
10
Table 3

Component definition for different mineral combinations (from top to bottom: only gypsum, gypsum and anhydrite, only anhydrite) and the ‘one component per mobile species’ component definition U1

 
 
 
 
 
 
 
 

This implies that the number of components is the same for the entire domain. This aspect is controlled by a single object of the CHEPROO class, and affects almost all classes: from the Solver, in charge of calculating the dimension of the system to be solved, to the Meshfield, in charge of storing field values and their derivatives to state variables.

Despite having several components, each with its own conservation equation or Phenomenon, PROOST treats components as entities pertaining to one Phenomenon. This simplifies the code's internal operability and problem definition, since it allows advantage to be taken of the fact that several Processes affect species in the same way. For example, the storage, advection, and diffusion–dispersion Processes in Equation (10) affect all species from a phase in the same way. For these Processes the contributions to the system matrix are calculated for all components together. Encapsulation allows confining to the Process class all the complexity associated with the fact that Processes can be part of one or a set of partial differential equations. Currently, the only Process that acts differently over each species is the ‘sink/source’ Process.

The CHEPROO object also defines which species and variables will be considered as state variables for the Newton-Raphson system. When only gypsum is present in the system, the states variables associated with component conservation equations are: , . The rest of species (, , ) are secondary and values are calculated by CHEPROO by considering mass actions laws. Reaction rates and nonmobile species concentrations are calculated in a subsequent step.

In order to understand the physical meaning of component conservation equations, it is helpful to associate state variables with specific components. For example, each of the species chosen as state variables (, , and ) can be considered as the constituents of the four first components of the conservation, Equation (10). The association of the liquid pressure () with a specific component is not straightforward. Liquid pressure is related to liquid saturation which affects all components. Since the last component in Equation (10) contains all water species and only involves secondary species, liquid pressure can then be associated with its main variable. However, variables like activity coefficients, density, viscosity, gas pressure, and liquid saturation, depend on all state variables and make the system fully coupled. Nevertheless, the exercise of defining a main variable for every component provides a more profound knowledge about variables dependency, which may be relevant for some cases, as is shown in the section ‘Anhydrite precipitation’. For that case water species is eliminated from the component equation and both calcium and sulfate are defined as secondary variables.

As can be seen in Figure 2, matrix system assembling is the core of a time interval resolution. It involves all the classes shown in Figure 1.

Once the system is solved, there are still unknown variables to be calculated: the eliminated species concentration and the equilibrium reaction rates.

These variables can be calculated by considering the species conservation equation. In order to avoid formulating a different Phenomenon, the PROOST class organization allows using the same structure as used for calculating component conservation for species conservation. This is one of the advantages of the PROOST class organization. The same Phenomenon is considered and only the component definition is changed. The new component definition considers every mobile species as a component (see Table 3): 
formula
11
Note that all the Processes in Equation (11) are analogous to Equation (10), except the last one. This is the only term in Equation (11) that has unknown variables (, , ); the other terms involve known variables. In order to calculate these unknown variables, the component definition is considered and a general method of the Process class, balance, is used to calculate all terms at the right-hand side of Equation (12): 
formula
12

The result is used by CHEPROO to calculate the reaction rates (evaporation, volatilization of dissolved air, and gypsum precipitation). Note that the number of equations exceeds the number of unknowns (eight and three, respectively). In theory, solution of all equations should give the same reaction rates. For simplicity we used the least square method for the solution of Equation (9). Once the reaction rates are calculated, the mole variations of mineral species can be computed (gypsum for this case). If a mineral is completely depleted or if the solution has become saturated for a new mineral, components should be redefined and calculations for the time step recalculated.

Anhydrite precipitation

As the system evolves over time, water activity decreases at the top of the column due to osmotic and capillary effect, and anhydrite starts to precipitate. When anhydrite and gypsum coexist a singularity, known as ‘invariant point’, occurs and water activity remains constant (Risacher & Clement 2001; Gamazo et al. 2011). Combining the mass actions law for anhydrite and gypsum (Equation (13)), the fixed water activity value can be obtained: 
formula
13
Under this scenario, gypsum dissolves and anhydrite precipitates at a rate that ensures this fixed water activity: 
formula
14
This implies a singular component definition (see Table 3) which results in the following component conservation equation: 
formula
15
Note that all forms of water have been eliminated from the components conservation equation. For all nodes where anhydrite and gypsum coexist these new components have to be considered, and the state variables associated with conservation equation will be: , , , .
This conservation equation, and the corresponding state variables, may make it difficult to associate state variables with specific components. As in the previous system, the state variables , , can be associated with the first, second, and fourth component, respectively. Hence, the remaining variable, liquid pressure, must be associated with the third component of Equation (15). Although this third component only contains Ca2+ and and no H2O, it still depends on liquid pressure through volumetric content, θaq, and retention curve. Hence, there is no problem in using it for the calculation of liquid pressure. However, since neither or are state variables, it is not straightforward to understand how the system can manage balance of these two species, especially when advective and diffusive fluxes are considered. In nodes where gypsum and anhydrite coexist water activity is fixed. The main mechanism for this is the sink source term of water produced by the interaction of these two minerals as shown in Equation (14). However, this interaction can also affect calcium and sulfate concentration through differences in precipitation rates. For instance, if the rate of gypsum dissolution is the same as anhydrite precipitation then only water is released, but any differences between these rates can release or consume dissolved Ca2+ and . Hence, gypsum–anhydrite interaction will release the necessary amount of water, calcium, and sulfate to keep water activity constant and to conform these two species conservation equations. This can be seen by considering the species conservation equations of calcium and sulfate (Equations (16) and (17)), and the sum of liquid and gaseous water (Equation (18)): 
formula
16
 
formula
17
 
formula
18
The most interesting aspect of this is that the precipitation dissolution rates were eliminated and therefore were not involved in the component conservation equations. This situation continues as long as the two minerals coexist. If that is not the case, and one mineral is depleted, the component definition is changed and the time step is recalculated.

When a mineral content is exhausted during the resolution of a time step, the new component definition considers that this mineral is no longer present. However, it was still present at the beginning of the time step. In order to keep track of this remaining amount of mineral, this amount is distributed among the species that form it and treated as a fixed source term. In the present application, when gypsum disappears, a new component definition is considered () and a source term is added for calcium, sulfate, and water equal to the amount of remaining gypsum times the corresponding stoichiometric constant, 1 for calcium and sulfate and 2 for liquid water.

The resulting component conservation equations are similar to the ones shown in Equation (10). 
formula
19
Note that the last component only involves liquid and gaseous water.

SUMMARY AND CONCLUSIONS

An object-oriented multiphase reactive transport class organization has been presented. It was designed to ensure extensibility and flexibility. Its main classes are: Mesh (contains all the discretization information and integration methods, such as finite elements or finite differences); Meshfield (represents spatial fields and the constitutive laws that relate them, like saturation or concentrations); Phenomenon (represents the conservation of a physical magnitude expressed as a PDE, such as mass or energy conservation); Process (represents a term of a Phenomenon PDE, like advection or storage, and non-Dirichlet boundary conditions); Solver (controls the coupling strategy for solving all Phenomena and assembles the matrix system to solve); and CHEPROO (encapsulates all thermodynamic data and performs geochemical calculations).

The flexibility and extensibility of PROOST come from the following particularities of its design. Several Phenomenon can be formulated by combining the available Processes. In order to solve a new kind of Phenomenon, only new Processes have to be programmed. The Solver class can be set to solve all Phenomena independently, sequentially, or coupled. New constitutive laws can be easily added to the code by creating new specialization of the Meshfield class, and new numerical methods for discretization–integration of PDE can be added by implementing new specializations of the Mesh class. The main challenging task, for solving reactive transport problems, was to implement the ability of using changing definitions of components. This could be achieved by considering the components as entities pertaining to one Phenomenon. The flexibility of the structure allowed the implementation of the SIA method by mainly creating a new specialization of the Meshfield class. The performance of PROOST is illustrated by describing the solution procedure for a concrete application: the modeling of a column of porous gypsum subjected to a constant source of heat. The problem involves important interaction between hydrodynamic and geochemical processes like the occurrence of invariant points. The flexibility of the structure is shown in the example. In this regard, it highlights the fact that a single Phenomenon object is considered for representing both component and species conservation in two different steps of the resolution procedure. This allows the calculation of the equilibrium reaction rates using pre-existing methods.

REFERENCES

REFERENCES
Akin
J. E.
1999
Object oriented programming via Fortran 90
.
Engineering Computations
16
,
26
48
.
Bea
S. A.
Carrera
J.
Ayora
C.
Batlle
F.
Saaltink
M. W.
2009
CHEPROO: a Fortran 90 object-oriented module to solve chemical processes in earth science models
.
Computers & Geosciences
35
(
6
),
1098
1112
.
Boivin
C.
Ollivier-Gooch
C.
2004
A toolkit for numerical simulation of PDEs. II. Solving generic multiphysics problems
.
Computer Methods in Applied Mechanics and Engineering
193
(
36–38
),
3891
3918
.
Commend
S.
Zimmermann
T.
2001
Object-oriented nonlinear finite element programming: a primer
.
Advances in Engineering Software
32
(
8
),
611
628
.
Debye
P.
Hückel
E.
1923
The theory of electrolytes. I. Lowering of freezing point and related phenomena
.
Physikalische Zeitschrift
24
,
185
206
.
Decyk
V. K.
Norton
C. D.
Szymanski
B. K.
1998
How to support inheritance and run-time polymorphism in Fortran90
.
Computer Physics Communications
115
,
9
17
.
DeSimoni
M.
Carrera
J.
Sanchez-Vila
X.
Guadagnini
A.
2005
A procedure for the solution of multicomponent reactive transport problems
.
Water Resources Research
41
(
11
),
W11410
.
Flemisch
B.
Darcis
M.
Erbertseder
K.
Faigle
B.
Lauser
A.
Mosthaf
K.
Müthing
S.
Nuske
P.
Tatomir
A.
Wolff
M.
Helmig
R.
2011
Dumux: DUNE for multi-{phase, component, scale, physics, …} flow and transport in porous media
.
Advances in Water Resources
34
(
9
),
1102
1112
.
Forde
B.
Foschi
R. O.
Stiemer
S. F.
1990
Object-oriented finite element analysis
.
Computers and Structures
34
(
3
),
355
374
.
Gamazo
P.
Bea
S. A.
Saaltink
M. W.
Carrera
J.
Ayora
C.
2011
Modeling the interaction between evaporation and chemical composition in a natural saline system
.
Journal of Hydrology
401
,
154
164
.
Gamazo
P.
Saaltink
M. W.
Carrera
J.
Slooten
L.
Bea
S. A.
2012
A consistent compositional formulation for multiphase reactive transport where chemistry affects hydrodynamics
.
Advances in Water Resources
35
,
83
93
.
Gorelik
A. M.
2004
Object-oriented programming in modern Fortran
.
Programming and Computer Software
30
,
173
179
.
Hao
Y.
Sun
Y.
Nitao
J. J.
2012
Overview of NUFT: A versatile numerical model for simulating flow
. In:
Groundwater Reactive Transport Models
(
Zhang
F.
Yeh
G. T.
Parker
J. C.
, eds).
Bentham e-Books, Bentham Science Publishers, Sharjah, UAE
. .
Jacques
D.
Šimůnek
J.
Mallants
D.
van Genuchten
M. T.
Yu
L.
2011
A coupled reactive transport model for contaminant leaching from cementitious waste matrices accounting for solid phase alterations
. In:
Sardinia 2011 Proceedings – Thirteenth International Waste Management and Landfill Symposium
.
Johnson
C. D.
Truex
M. J.
2006
RT3D Reaction Modules for Natural and Enhanced Attenuation of Chloroethanes, Chloroethenes, Chloromethanes, and Daughter Products, PNNL-15938
.
Pacific Northwest National Laboratory
,
Richland, Washington
,
USA
.
Kolditz
O.
Bauer
S.
2004
A process-oriented approach to computing multi-field problems in porous media
.
Journal of Hydroinformatics
6
,
225
244
.
Lagneau
V.
Van Der Lee
J.
2010
HYTEC results of the MoMas reactive transport benchmark
.
Computational Geosciences
14
,
435
449
.
Lichtner
P. C.
1996
Continuum formulation of multicomponent–multiphase reactive transport
. In:
Reactive Transport in Porous Media, Reviews in Mineralogy
(
Lichtner
P. C.
Steefel
C. I.
Oelkers
E. H.
, eds).
Vol. 34
,
Mineralogical Society of America
,
Washington, DC
,
USA
, pp.
1
81
.
Lichtner
P. C.
Hammond
G. E.
Lu
C.
Karra
S.
Bisht
G.
Andre
B.
Mills
R. T.
Kumar
J.
2013
PFLOTRAN User Manual: A Massively Parallel Reactive Flow and Transport Model for Describing Surface and Subsurface Processes. http://www.pflotran.org/docs/user_manual.pdf
.
Loomer
D. B.
Al
T. A.
Banks
V. J.
Parker
B. L.
Mayer
K. U.
2010
Manganese valence in oxides formed from in situ chemical oxidation of TCE by KMnO4
.
Environmental Science & Technology
44
(
15
),
5934
5939
.
MacQuarrie
K.
Mayer
K. U.
2005
Reactive transport modeling in fractured rock: a state-of-the-science review
.
Earth-Science Reviews
72
(
3/4
),
189
227
.
Maley
D.
Kilpatrick
P. L.
Schreiner
E. W.
Scott
N. S.
Diercksen
G. H. F.
1996
The formal specification of abstract data types and their implementation in Fortran 90: implementation issues concerning the use of pointers
.
Computer Physics Communications
98
(
1–2
),
167
180
.
Mayer
K. U.
Amos
R. T.
Molins
S.
Gérard
F.
2012
Reactive transport modeling in variably saturated media with MIN3P: Basic model formulation and model enhancements
. In:
Groundwater Reactive Transport Models
(
Zhang
F.
Yeh
G. T.
Parker
J. C.
, eds).
Bentham e-Books, Bentham Science Publishers, Sharjah, UAE
. .
Meysman
F. J. R.
Middelburg
J. J.
Herman
P. M. J.
Heip
C. H. R.
2003
Reactive transport in surface sediments. I. model complexity and software quality
.
Computers & Geosciences
29
(
3
),
291
300
.
Molins
S.
Carrera
J.
Ayora
C.
Saaltink
M. W.
2004
A formulation for decoupling components in reactive transport problems
.
Water Resources Research
40
(
10
),
W10301
.
Norton
C. D.
Decyk
V.
Slottow
J.
1998
Applying Fortran 90 and object-oriented techniques to scientific applications
.
Object-Oriented Technology
1543
,
462
463
.
Parkhurst
D. L.
Kipp
K. L.
Engesgaard
P.
Charlton
S. R.
2004
PHAST – A program for simulating ground-water flow, solute transport, and multicomponent geochemical reactions: US Geological Survey Techniques and Methods 6–A8, 154 pp.
Parkhurst
D. L.
Appelo
C. A. J.
2013
Description of input and examples for PHREEQC version 3 – a computer program for speciation, batch-reaction, one-dimensional transport, and inverse geochemical calculations
.
US Geological Survey Techniques and Methods
,
book 6, chap. A43
,
Denver, CO
,
USA
.
Prommer
H.
Post
V.
2010
PHT3D, A Reactive Multicomponent Transport Model for Saturated Porous Media. User's Manual v2.10
.
Saaltink
M. W.
Ayora
C.
Carrera
J.
1998
A mathematical formulation for reactive transport that eliminates mineral concentrations
.
Water Resources Research
34
(
7
),
1649
1656
.
Saaltink
M. W.
Carrera
J.
Ayora
C.
2001
On the behavior of approaches to simulate reactive transport
.
Journal of Contaminant Hydrology
48
(
3–4
),
213
235
.
Saaltink
M. W.
Batlle
F.
Ayora
C.
Carrera
J.
Olivella
S.
2004
Retraso, a code for modeling reactive transport in saturated and unsaturated porous media
.
Geologicaacta
2
(
3
),
235
251
.
Sassen
D. S.
Hubbard
S. S.
Bea
S. A.
Chen
J.
Spycher
N.
Denham
M. E.
2012
Reactive facies: an approach for parameterizing field-scale reactive transport models using geophysical methods
.
Water Resources Research
48
,
W10526
.
Shao
H.
Dmytrieva
S. V.
Kolditz
O.
Kulik
D. A.
Pfingsten
W.
Kosakowski
G.
2009
Modeling reactive transport in non-ideal aqueous-solid solution system
.
Applied Geochemistry
24
(
7
),
1287
1300
.
Slooten
L. J.
Batlle
F.
Carrera
J.
2010
An XML based problem solving environment for hydrological problems
. In:
XVIII Conference on Computational Methods in Water Resources (CMWR)
. .
Soleimani
S.
Van Geel
P. J.
Isgor
O. B.
Mostafa
M. B.
2009
Modeling of biological clogging in unsaturated porous media
.
Journal of Contaminant Hydrology
106
(
1–2
),
39
50
.
Steefel
C. I.
2009
Crunch Flow Software for Modeling Multicomponent Reactive Flow and Transport User's Guide
.
Earth Sciences Division, Lawrence Berkeley National Laboratory
,
Berkeley, CA
,
USA
.
Steefel
C. I.
MacQuarrie
K. T. B.
1996
Approaches to modeling of reactive transport in porous media
.
Reviews in Mineralogy and Geochemistry (Reactive Transport in Porous Media)
34
,
85
129
.
Steefel
C. I.
Appelo
C. A. J.
Arora
B.
Jacques
D.
Kalbacher
T.
Kolditz
O.
Lagneau
V.
Lichtner
P. C.
Mayer
K. U.
Meeussen
J. C. L.
Molins
S.
Moulton
D.
Shao
H.
Šimůnek
J.
Spycher
N.
Yabusaki
S. B.
Yeh
G. T.
2015
Reactive transport codes for subsurface environmental simulation
.
Computational Geosciences
19
(
3
),
445
478
.
Trebotich
D.
Adams
M. F.
Molins
S.
Steefel
C. I.
Shen
C.
2014
High-resolution simulation of pore-scale reactive transport processes associated with carbon sequestration
.
Computing in Science & Engineering
16
(
6
),
22
31
.
Villar
M. V.
Sánchez
M.
Gens
A.
2008
Behaviour of a bentonite barrier in the laboratory: experimental results up to 8 years and numerical simulation
.
Physics and Chemistry of the Earth, Parts A/B/C
33
,
S476
S485
.
Wang
W.
Kolditz
O.
2007
Object-oriented finite element analysis of thermo-hydro-mechanical (THM) problems in porous media
.
International Journal of Numerical Methods in Engineering
69
,
162
201
.
Wheeler
M. F.
Sun
S.
Thomas
S. G.
2012
Modeling of flow and reactive transport in IPARS
. In:
Groundwater Reactive Transport Models
(
Zhang
F.
Yeh
G. T.
Parker
J. C.
, eds).
Bentham e-Books, Bentham Science Publishers, Sharjah, UAE
. .
White
M. D.
McGrail
B. P.
2005
STOMP Subsurface Transport Over Multiple Phases, PNNL-15482
.
Pacific Northwest National Laboratory
,
Richland, Washington
,
USA
.
Xu
T.
Spycher
N.
Sonnenthal
E.
Zhang
G.
Zheng
L.
Pruess
K.
2011
TOUGHREACT Version 2.0: a simulator for subsurface reactive transport under non-isothermal multiphase flow conditions
.
Computers & Geosciences
37
,
763
774
.
Xu
T.
Spycher
N.
Sonnenthal
E.
Zhang
G.
Zheng
L.
Pruess
K.
2012
TOUGHREACT: A simulator for subsurface reactive transport under non-isothermal multiphase flow conditions
. In:
Groundwater Reactive Transport Models
(
Zhang
F.
Yeh
G. T.
Parker
J. C.
, eds).
Bentham e-Books, Bentham Science Publishers
,
Sharjah, UAE
. .
Yeh
G. T.
Fang
Y. L.
Zhang
F.
Sun
J. T.
Li
Y.
Li
M. H.
Siegel
M. D.
2010
Numerical modeling of coupled fluid flow and thermal and reactive biogeochemical transport in porous and fractured media
.
Computational Geoscience
14
,
149
170
.
Yeh
G. T.
Tripathi
V. J.
Gwo
J. P.
Cheng
H. P.
Cheng
R. J.
Salvage
K. M.
Li
M. H.
Fang
Y. L.
Li
Y.
Sun
J. T.
Zhang
F.
Siegel
M. D.
2012
HYDROGEOCGEM: A coupled model of variably saturated flow, thermal transport, and reactive biogeochemical transport
. In:
Groundwater Reactive Transport Models
(
Zhang
F.
Yeh
G. T.
Parker
J. C.
, eds).
Bentham e-Books. Bentham Science Publishers, Sharjah, UAE
. .
Zhang
R.
Yin
X.
Winterfeld
P. H.
Wu
Y.
2012
A fully coupled model of nonisothermal multiphase flow, geomechanics and chemistry during CO2 sequestration in brine aquifers
. In:
Proceedings, TOUGH Symposium 2012
,
Lawrence Berkeley National Laboratory
,
Berkeley, CA
,
USA
,
September 17–19
.