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 CO_{2} 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 CO_{2} 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.

Code | PDE discretization^{1} | Transport and reaction coupling^{2} | Phase conservation and transport coupling^{3} |
---|---|---|---|

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 discretization^{1} | Transport and reaction coupling^{2} | Phase conservation and transport coupling^{3} |
---|---|---|---|

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 |

^{1}DGM, 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.

^{2}DS, direct substitution; OS, operator split.

^{3}COU, 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

### Species and component conservation equation

*i*belonging to phase , which is a particular case of Equation (1) has the following expression: 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: Mobile phase fluxes are calculated according to Darcy's law: where , and are the conductivity tensor, pressure and density of the phase , respectively. Diffusive-dispersive fluxes are calculated according to Fick's law: where and are the diffusion and dispersion tensor for phase , respectively, and is the tortuosity.

*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: 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.

*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.

*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.

*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: Thus the component conservation equation is written only in terms of mobile component conservation: 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

## 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.

H_{2}O_{(l)}, H_{2}O_{(g)}, air_{(l)}, air_{(g)}, Ca, SO_{4}, K, Cl, gypsum, anhydrite |
---|

H_{2}O_{(l)}, H_{2}O_{(g)}, air_{(l)}, air_{(g)}, Ca, SO_{4}, K, Cl, gypsum, anhydrite |
---|

### Initial component definition

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.

*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):

*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):

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

*et al.*2011). Combining the mass actions law for anhydrite and gypsum (Equation (13)), the fixed water activity value can be obtained: Under this scenario, gypsum dissolves and anhydrite precipitates at a rate that ensures this fixed water activity: This implies a singular component definition (see Table 3) which results in the following component conservation equation: 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: , , , .

^{2+}and and no H

_{2}O, 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 Ca

^{2+}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)): 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.

## 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.