PROOST: object-oriented approach to multiphase reactive transport modeling in porous media

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

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.

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.

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

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

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.

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.

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.

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.

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.

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.

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.

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.

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.

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

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 .

View Large

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.

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.

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.

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.