Element balance formulation in reactive compositional flow and transport with parameterization technique

We present a novel nonlinear formulation for modeling reactive-compositional flow and transport in the presence of complex phase behavior due to a combination of thermodynamic and chemical equilibria in multi-phase systems. We apply this formulation to model precipitation/dissolution of minerals in reactive multiphase flow in subsurface reservoirs. The proposed formulation is based on the consistent element balance reduction of the molar (overall composition) formulation. To predict the complex phase behavior in such systems, we include the chemical equilibrium constraints to the multiphase multicomponent negative flash calculations and solve the thermodynamic phase and chemical phase equilibria simultaneously. In this solution, the phase equilibrium is represented by the partition coefficients, whereas the chemical equilibrium reaction is represented by the activity coefficients model. This provides a generic treatment of chemical and thermodynamic equilibrium within the successive substitution loop of mulmultiphase flash to accommodate chemical equilibrium reactions (precipitation and dissolution reactions). Equilibrium Rate Annihilation matrix allows us to reduce the governing component conservation equations to element conservation equations, while the coupling between chemical and thermodynamic equilibrium is captured by a simultaneous solution of modified multiphase flash equations. The element balance equation written in terms of overall component mole fractions is modified and defined in terms of element mole fractions. Therefore, the primary set of governing equations are the element balance equations and the kinetic equations. This element composition of the mixture serves as an input to the modified multiphase flash computations, whereas the output is fractions of components in each phase, including solids. The nonlinear element–based governing equations are solved with the modified version of the operator-based linearization (OBL) approach where the governing equations are formulated in terms of space- and state-dependent parameters constrained by the solution of the extended multiphase flash. The element balance molar formulation along with the modified multiphase flash has been tested in a simple transport model with dissolution and precipitation reactions. The simulation of multidimensional problems of practical interest has been performed using the adaptive OBL technique. This is the first time when a robust multiphase multicomponent flash based on element fractions is coupled with an element balance–based compositional formulation and tested for multidimensional problems of practical interest. The proposed technique improves both robustness and performance of complex chemical models.


Introduction
information into the simulator to make better production forecast.
Fluids in subsurface reservoirs are a complex mixture of different types of hydrocarbons which can be present as either gas, liquid, or even solid phase (asphaltenes, coke), depending on the pressure, temperature, and composition of the reservoir fluids. Along with the hydrocarbon vapor and liquid phases, there is always a presence of an aqueous phase which contains mainly water but can also have dissolved hydrocarbon components and minerals. In addition to the in situ reservoir components, other components are also introduced into the reservoir mixture during injection operations to enhance the production and recovery from the reservoir. Examples include low salinity brine, surfactants, polymers, and CO 2 . The distribution of components in different phases are usually controlled by thermodynamic equilibrium, chemical equilibrium, and chemical kinetics. Therefore, for accurate modeling of such systems, we need to couple these effects with flow and transport which requires solving multiple nonlinear governing equations simultaneously.
The modeling of flow and transport in the subsurface is divided into two different research categories. The first direction is being supported by the reservoir engineering community which focuses on complex equation of state (EOS) models to resolve the multiphase multicomponent flow and transport. The second category of models is being developed by the hydrological community which almost ignores multiphase phenomena due to the limited presence of other phases, but incorporate a wide range of chemical equilibrium and kinetic reactions. One of the first efforts to rigorously couple chemical reactions with the flow has been done in [14]. He reduced the species into primary and secondary components and decoupled the local chemical equilibrium reactions from the solution of single-phase flow [13].
From another side, compositional reservoir simulators put a major emphasis on phase resolution due to the presence of multicomponent fluid. Traditional reservoir simulators determine the phase state of the fluid using the Gibbs energy minimization technique proposed in [16], which is based on the equality of chemical potentials of species in different phases at equilibrium. Using this to determine the stability of the fluid system, flash calculations can be run if the system has multiple phases as suggested in [17] to determine the phase fractions and compositions.
The phase split problem at equilibrium assumptions includes the solution of the nonlinear system of equations. In [22], they reformulate this system in one nonlinear equation which yields phase fraction of phase. Later, in [12], authors allowed the phase fractions to become negative which indicates single-phase mixture. Next, the range of vapor fraction values was bounded in [28]. Since any point in sub-critical single-phase region can be parameterized by the tie-line, the negative flash procedure can be used as an effective phase-state identification method which does not require phase stability test [7].
With the use of modern production technologies like enhanced oil recovery, well acidization, and CO 2 sequestration, various phenomena involve dissolution and precipitation in the reservoir, mostly in the near-well region. This introduces a need to effectively couple the chemical reactions with multiphase flow and transport using reactive compositional simulation. The first attempt to introduce this coupling was performed in [5] where he applied the element-based formulation for multiphase reactivecompositional flow and transport in reservoir simulation. Later, in [6], the model was extended to multiple equilibria and kinetic reactions using both natural-and molar-based formulations in the Automatic Differentiation General Purpose Research Simulation (ADGPRS). In [18], an element and species balance formulation for chemical reactions was introduced. The Gibbs-Helmholtz constrained (GHC) equation of state solver for the simultaneous solution of thermodynamic and chemical equilibria [15], was incorporated in ADGPRS framework to model multiple salt precipitation/dissolution reactions [26].
There are different techniques which are used to couple the non-linear compositional transport equation with the chemical reactions. In [29], authors studied these different approaches and then suggested the sequential iterative approach (SIA) to model the reactive transport, which is faster and has lower numerical dispersion but requires smaller time steps for stability. Another way to solve the system is by using the sequential non-iterative technique where first the flow and transport is solved and later, using the transported concentration, the reaction terms are solved. This method assumes that the transport is fast and the reactions can be applied after the transport solution is complete, which is not an accurate assumption for, e.g., in kinetic reactions which has varied time scales. In [5] and [6], authors solved reactive flow and transport problems using fully implicit formalism also known as the global implicit technique in the hydrological community, where the transport and the reactions are solved simultaneously. The fully implicit technique gives the liberty to perform large time steps without stability issues. However, the phase behavior evaluation was not tightly coupled with chemical equilibrium computations in these approaches which affect the robustness and efficiency of the simulation process.
In this work, we utilize a recently developed fully implicit operator-based linearization (OBL) technique [25]. The OBL method controls the nonlinearity of the problem with the multi-linear representation of state-dependent operators in governing equations. During the course of the simulation, these operator values are linearly interpolated on the mesh with a predefined accuracy. The values of the state operators are calculated adaptively using conventional property estimators based on correlations or solution of equation of state [9]. Thus, this method additionally improves simulation time by skipping the routine evaluation of computationally expensive phase behavior calculations performed by the conventional technique. The OBL approach has been implemented and successfully tested for the solution of complex geothermal [8] and multiphase multicomponent flow and transport problems with buoyancy [10].
In this paper, we develop an algebraic framework for the simultaneous solution of thermodynamic and chemical equilibria and apply it to reactive multiphase multicomponent flow and transport problem. The first section describes governing equations for multiphase compositional simulation in chemical and thermodynamic equilibrium assumptions. In the second section, we explain the reduction of governing equations to element formulation. The third section describes the simultaneous solution of thermodynamic and chemical equilibrium using an extended negative flash technique and includes an illustrative example. This approach is implemented in Delft Advanced Reservoir Terra Simulation (DARTS) [4] to adaptively parameterize the governing equations for simulation which is described in the fourth section. The last section presents several numerical simulations performed in DARTS and used to test and validate the developed framework.

Conventional governing equations
In this section, we present the conventional governing equations which define the mass balance of species in reactive flow and transport at thermodynamic and chemical equilibrium assumptions.

Conservation of components
We start with the basic mass balance equations including the effect of chemical reactions as source/sink term [11]: where n c is the overall mass of component, l c is the total flux associated with that component, q c is the total well flow rate associated with that component, v ck is the stoichiometric coefficient associated with kinetic reaction k for the component c and v cq is the stoichiometric coefficient associated with equilibrium reaction q for component c, r k is the rate for kinetic reaction, and r q is the equilibrium reaction rate.
In this study, we only consider equilibrium reactions. The overall mass of components is defined as Here, P stands for the total number of fluid phases and M stands for total number of solid phases. Here, the first term indicates total mass of component c in all the fluid phases; whereas the second term is the mass of component c in the solid phases. The term l c defines the flux of component c and is given as: where the term d cj corresponds to the dispersion of component c in phase j. For simplicity, this term is neglected in our study. The term u j is the velocity of the phase j and is defined by Darcy's law: The rest of the properties are split into the statedependent (ω) and space-dependent (ξ ) relations.
• Space based (properties altered in space): In our work, we assume that the solid phase is not moving; therefore, we take the velocity of solid phase zero and it is not considered in the flux term. With all above-mentioned assumptions, Eq. 1 can be written in a vector form: where n = (n 1 ,. . . ,n C ) T , l = (l 1 ,. . . ,l C ) T , q = (q 1 ,. . . ,q C ) T is the well flow rate, V is the stoichiometric matrix in a reaction, and r = (r 1 ,. . . ,r Q ) T is the reaction rate vector.

Closing assumptions
Equilibrium reactions are defined as reactions with the forward reaction rate the same as the rate of backward reaction. There are no observable changes in the properties of the system since the reactants and the products do not change with time but, the system is in dynamic equilibrium.
Since the reaction proceeds very fast and is reversible in nature, the instantaneous local equilibrium assumption can be made. The equilibrium relations are defined by the law of mass action and are given as: Here, α c is the activity of component c, Q q is the reaction quotient, and K q is the equilibrium reaction quotient or equilibrium solubility limit in the case of dissolution/precipitation of minerals. The next closing relation is the fugacity constraint which is based on the equality of fugacity of a particular component across different phases. A component is at thermodynamic equilibrium if the chemical potentials of the component in both phases are equal. Numerically, it can be written as The fugacity of a component in a particular phase is given by where φ ij is fugacity coefficient. These relations (7) can also be written in terms of partition coefficients (K) There are two methods to resolve the thermodynamic phase behavior: EOS-based approach and constant K value approach, see [19] for details. The EOS-based approach starts with an initial guess of K values and they are updated in every iteration using the cubic EOS until the fugacity constraint is satisfied. In the second approach, K values are considered only a function of pressure and temperature but not the composition. This representation is valid for hydrocarbon mixtures at low-pressure conditions. However, it can be extended for higher pressure as well with additional composition dependency [23]. For simplicity, we are using the constant K value assumption in this work.
The final set of closing relations are given below where v j is defined as and z i is the overall mole fraction of component i. Equation 10 is the definition of overall composition of component c, Eq. 11 is the phase composition, and Eq. 12 is the overall phase fraction constraint.

Nonlinear formulation
There are different nonlinear formulations which can be used for the solution of the compositional flow problem. These formulations depend on the type of primary equations and unknowns selected for a fully implicit system. There are two major types of formulations which are being used in reservoir simulation community: (1) natural formulation [2] and (2) molar formulation [1,3]. An extended analysis of the different types of formulations has been covered in [26]. In this study, we are using the molar formulation as an alternative to the natural formulation suggested in [5]. Comparison of reactive simulation using both natural and molar formulation has been performed in [6] and can be referred for further details. The governing Eq. 1 in molar formulation are solved for pressure (p) and overall composition (z) where the overall molar mass of component in (1) can be written as and ρ t is the total density of the fluid.
In a fully implicit model, when more that one phase exists, we need to determine the secondary variables using the given primary variable p and z. The set of secondary Eqs. 9 to 12 are solved using multiphase flash procedure. The governing equations are linearized using the Newton-Raphson approach given as J y k y k+1 − y k = −r y k .
Here, J stands for the Jacobian matrix with respect to primary unknowns, k stands for the nonlinear iteration, and y stands for a vector of nonlinear unknown which are pressure and overall compositions. Finally, r is the residual of the mass balance Eq. 5 at k-th iteration. The linear Eq. 15 is solved on every nonlinear iteration to obtain an update of primary unknowns. If the residual is below the pre-defined tolerance, then the system is converged, and the Newton loop is terminated. In this approach, the derivatives of secondary variables with respect to the primary nonlinear unknowns are calculated using the inverse theorem [26].

Element balance reduction
In this section, we describe the reduction of component mass balance equations to element mass balance equations with thermodynamic and chemical equilibria constraints used as secondary equations.

Stoichiometric and annihilation matrices
The stoichiometric matrix takes into account the mass balance of components in chemical reactions. Stoichiometric matrix S, when written in canonical form, consists of component (primary species) and non-components (secondary species). Secondary species are components which are uniquely defined by a chemical reaction and can be written in terms of primary species. The general form of the stoichiometric matrix as suggested in [6] is given below (16) Here, R stands for the total number of reactions, Q is the number of equilibrium reactions, K is the number of kinetic reactions, and C is the total number of components. The rows represent the components involved in chemical reactions, whereas the columns represent the chemical reactions itself. In this study we are using the equilibrium reactions only; therefore, R is equal to Q and K = 0. Equilibrium Rate Annihilation matrix E removes the equilibrium reaction rates from the governing mass balance equations when they are pre-multiplied by it. In addition, this multiplication reduces the C component mass balance equations to E element mass balance equations. The E matrix can be visualized as the distribution of elements in different components and can be written in the matrix form as The multiplication by E matrix also linearly combines the kinetic reaction components so that they appear only in the corresponding governing equations. Therefore, C component balance equations are reduced to E element mass balance equations, K differential kinetic reaction relations, and Q algebraic chemical equilibrium relations which sum to C (C = E + K + Q). The formulation of E matrix is dependent on the stoichiometric matrix S. The general form of E is given in [6] as The rows of the E matrix are chemical elements which combine to form the components. These elements are the smallest chemical species which do not disassociate into the smaller entities. The column of matrix E represents the components involved in the system. To illustrate this approach, we consider a simple system consisting H 2 O, CO 2 and one dissolution-precipitation reaction of CaCO 3 . Since H 2 O and CO 2 do not disintegrate into smaller species, they can be treated as elements for this example. Whereas CaCO 3 disassociates into Ca 2+ and CO 2− 3 ions, therefore these ions are considered as elements. The equilibrium rate annihilation matrix for this system is given as The stoichiometric matrix for this system is given as

Reduction to element formulation
The addition of chemical reactions to compositional formulation posses serious challenges. Firstly, the instantaneous chemical equilibrium reaction rates make the transport problem very stiff, which requires very small time steps to resolve the nonlinearity. Therefore, we use the local equilibrium assumption and decouple the chemical equilibrium reactions. Secondly, equilibrium reactions and kinetic reactions are both functions of concentrations of reactants and products, which are in turn functions of transport. As a result, we need a robust mechanism to capture this coupling. Next, we formulate the element balance approach as suggested in [24] and write the element-based governing equations in terms of mole fractions of elements. To write the governing equation in terms of elements, we pre-multiply the vector Eq. 5 with E which yields the element balance equation The above equation is now reduced to E mass balance equations written as and K kinetic reaction relations The above relations are still written in terms of overall mole fractions and need to be converted to the governing equation in terms of element mole fractions (z E ) using the relations described below. Element to composition transformation shows how the elements combine to form the components and are written as Total molar component density is given by and total molar element density is given as total moles of the element to the total volume Introducing ρ E T in Eq. 20 transforms the mass balance equation in terms of z E and is written as Using the above relations, we have additional E unknowns apart from the C already present. This E unknowns are supplemented by the E element to component transformation Eqs. 22 which can also be written as The flux terms are still written in terms of component compositions in each phase which can be determined from the multiphase flash coupled with the chemical solver. As there are a fewer number of elements compared to the total number of components, the above system becomes under-defined. Therefore, we also include the chemical equilibrium as the closing relations along with the phase equilibrium relations in a Newton loop with the element to component transformation in Eq. 26. The unknown variable set for the element formulation is given below: . . . , E, c = 1, . . . , C, There are in total E + C + P + CP + 1 unknowns in this formulation which has additional E equations in comparison to the overall formulation. The system is closed by E element conservation equation, K kinetic reaction relations, and Q equilibrium relations which together adds up to C. Then there are E element to component transformation relations (26), C(P − 1) fugacity relations, C overall mole fraction relation, P − 1 phase composition relations, and one phase constraint and overall mole fraction relations. So adding all of them up, the total number of equations becomes C + E + CP + P + 1 which is equal to the total number of unknowns.

Thermodynamic and chemical equilibria
The phase constraint relations, when solved with chemical equilibrium relations and element to component transformation, suggests how elements partition in different phases and how these elements are combined to form components. The chemical equilibrium relations, which follow the law of mass action (6), are solved with thermodynamic relations (9) to (12) for phase behavior prediction. The chemistry of precipitation and dissolution reactions are simplified by several adjustments which can be later relaxed.
If the value of Q q > K q , then the reactants (minerals) will form instantaneously in order to reduce Q q and the opposite occurs when Q q < K q . Any change in the system at equilibrium will cause the composition to move in the direction which neglects the change according to Le Chatelier's principle. Therefore, if Q q is greater than K q , precipitation occurs in cases of mineral reactions as more reactants are formed to reduce Q q . Since the equilibrium reactions occur instantaneously (independent of time), and they are the only function of the concentration of components in the grid blocks, we can decouple them from the mass balance and consider them as secondary equations.
As shown in [5], the equilibrium relations can be written in terms of mole fractions instead of the activities. The activity of the pure mineral phases and water is considered as one. The activity of components can be written in terms of molality as given below: Here, standard molality of solute species is taken as one mol/kg, m cw stands for molality of component c, and γ cw is the activity coefficient of component c in water. Activity coefficient for very dilute solutions can be taken as one. The molality of a component is related to the mole fractions of the components as where M wm = 55.508 is the moles of H 2 O per kilogram of the aqueous phase, x cw stands for mole fraction of component c in the aqueous phase, and x ww stands for mole fraction of water in the aqueous phase. The initial guess for phase compositions is determined by solving the RR equation for thermodynamic equilibrium as shown in Appendix A.

Rigorous equilibrium validation
In this section, we resolve the phase behavior at local thermodynamic and chemical equilibrium assumptions using the method described earlier. We first start with a simple one equilibrium reaction system with two fluid phases and one solid mineral phase present.
Using the law of mass action, we can algebraically write the above reaction in terms of activities of present components as shown in (6). The molality relations can be used to convert the relations in terms of component mole fractions For very dilute solution, we can assume activity coefficient to be one and reduce the above equation to To resolve the combined thermodynamic and chemical equilibria simultaneously, Eqs. 9 to (12) are solved along with (26) and (6) in a nonlinear loop. If the solution gives all the three-phase fraction values as positive, then the system is in a three-phase region. If any of the phase fractions are negative, then that phase is missing, and we reduce the system of equations to two-phase case; if both two-phase fractions are negative, then there is only one phase present. Figure 1 shows the phase distribution generated using the element flash and chemical equilibrium solver for different values of K sp from 10 −08 to 5800. A range of K sp values is taken to check the robustness of the coupled thermodynamic and chemical equilibrium solver and to understand how the phase regions change with changing K sp values. Therefore, starting with the element mole fractions, partition coefficients, and equilibrium constant values, we can resolve the phase behavior of chemical and thermodynamic equilibrium system in a coupled manner. The yellow region in Fig. 1 is the three-phase region where all the three phases (aqueous, gas, and solid) are present. As the value of K sp increases, it implies that more ions can dissolve in the aqueous phase, hence the solid phase region begins to shrink. This can be clearly seen in the figure, as we increase the K sp values.

Operator-based linearization
OBL framework has been used before to solve compositional and geothermal problems with buoyancy [8,10], but has never been tested for flow with chemical reactions. The OBL approach provides a flexible solution for nonlinear formulations with complex derivatives in a fully implicit manner. In element balance formulation, all derivatives need to be evaluated with respect to element concentration z E which is a non-trivial procedure.
Equation 31 gives the description of conventional OBL approach for non-reactive compositional problem. The finite-volume discretization of compositional governing equations is given as This equation can be rewritten in parameterized form with space-and state-dependent operators as where ω is a state-dependent parameter and ξ is a space dependent parameter. Equation 31 is translated into Eq. 32 using the following operators: Here, α and β are state-dependent operators. The value of these operators can be determined for parametrized set of p and z points, by evaluating different properties at these points. These points are defined as base nodes or nodal values. During the course of the simulation, the operators are evaluated using multi-linear interpolation in parameter space. The larger the number of base node points we use, the higher is the accuracy of the interpolation.
To improve the performance of OBL approach for a large number of species, an adaptive extension has been proposed in [10]. This is a minor extension of OBL approach in which the tables are generated during the course of a simulation. The grid is uniform but the values of the operators at the nodes are calculated during the course of simulation when it is required. Adaptive OBL is useful for simulating a multicomponent system where only a limited number of compositions defines hyperbolic transport. Therefore, the phase behavior is resolved only at those nodes which are used for interpolation. The chemical reaction system was

Adding chemical reaction in OBL
Until now, we discussed the OBL approach for compositional transport, i.e., where the operators do not take chemical reactions into account. Here, we modify the operators to include equilibrium reactions. We start with a simple CaCO 3 dissolution-precipitation reaction example where the reaction is considered as an equilibrium one. This system has been defined in the previous section. Since the governing equations in reactive transport are written in terms of elements, the operators have slightly different form as compared to the compositional OBL formulation. The element-based governing equations can be described as Here, φ T is the total porosity of the rock which includes the reactive (mineral) and the fluid part, and it is always constant throughout the course of the simulation as the nonreactive part of the rock does not take part in any reaction or flow. The fluid porosity keeps on changing with the change in reactive porosity, hence affecting the permeability of the system. Since the elements are also part of solids, the total mole fractions of elements also include the solid component. Therefore, the accumulation term should also include the reactive part of the reservoir. Hence, the total porosity will remain constant in the accumulation term. The alpha operator for a chemical reaction is given as Using the definition of ρ E T from Eq. 24, α i operator can be can written as Space-dependent parameter a is modified to include the total porosity instead of fluid porosity. The total porosity is constant throughout the course of the simulation, unlike the fluid porosity which changes as the reaction proceeds. As a result of considering total porosity in the accumulation term, φ T o becomes a space-dependent parameter instead of statedependent parameter reducing the nonlinearity in the alpha operator Here, φ T o is the total initial porosity which does not vary with time but only in space. The concept of total porosity and how it depends on the mineral mole fraction is described in Appendix B.
The alpha operators generated for the current system are shown in Fig. 2. The last operator for the current system represents the combined accumulation of both the ion species Ca 2+ and CO 2− 3 since they behave exactly the same when there is no source for these ions. Hence, we are able to represent a four-element system with just three elements. This is not the case when there are multiple reactions or there is some source of individual ions. In real cases, the number of operators would be equal to the number of elements as the individual ion concentration will be different in a particular grid block for some cases.
Here, z E is an E-element vector, whereas ρ E T and compressibility are scalar quantities. ρ E T depends on the overall mole fraction distribution of the components in a grid block. For the current system, we assume the density of all phases to be equal to one so we can reduce the alpha operator to The β i operators are similar to the conventional compositional case with the only difference that it should be pre-multiplied by the Equilibrium Rate Annihilation matrix to get the element-based flux terms. The E matrix takes into account the flux of all components formed by a particular element. Therefore, it is a sum of overall fractional flow of all the components which contain the individual elements. The β i operators for the current system is given below where E is the equilibrium annihilation matrix and l is the flux term. The above equation can be written in terms of phase compositions as where e ic is the amount of element i in component c. In total, there should be E number of β i operators but for this example, we just need three similar to α i operators, see Fig. 3. The above-described operators linearize the reactive flow and transport governing equations. This problem is solved using the fully implicit scheme. A uniform grid is used for parameterizing and the operators are generated and stored adaptively.

Numerical simulation results
In this section, we present the results generated using the newly developed element-based formulation and adaptive OBL for different cases with single and multiple reactions. The K sp values for the equilibrium reactions are highly scaled up to make the reaction more pronounced in the simulation for the testing purposes.

Reactive flow and transport in 1D
The results are generated using the adaptive OBL simulation framework in the 1D test case for one calcite dissolution and precipitation reaction system shown below: . The rest of the parameters for the simulation can be found in Appendix C. The simulation was performed for 2000 days. In the first simulation run, we model just one mineral reaction; therefore, only the first reaction affects the porosity. Figure 4 shows the element mole fraction    The right part of Fig. 4 shows the overall composition in which the green line represents the calcite mole fraction. It is clear, that near the injection well, there is vaporization of water which leads to the calcite deposition. As a result, the fluid porosity of the system decreases and the reactive porosity increases while keeping the total porosity of the system constant. Due to the variation in fluid porosity, the permeability of the reservoir will also be affected which in turn will affect the transport.
Next, we consider two dissolution/precipitation reactions in the system by adding the MgCO 3 equilibrium reaction. The simulation results are shown in Fig. 5. Here, the element and component mole fractions at the end of the simulation are shown. In the composition plot, we can also see the fluid porosity variation due to dissolution and precipitation reactions, which depends on the amount of both minerals present.

Reactive flow and transport in 3D
Next, we use the developed framework for a modification of the 3D Brugge field model [21]. In our models, we look into CO 2 sequestration to the depleted gas field with only one CO 2 injection wells and two gas production wells. We start with the one-layer model and then extend it to the 3D model.
The simulation was performed for one layer (layer 5) of the Brugge field only. Figure 6 shows the dynamic at various times.
The element-based governing equations are solved using the fully implicit adaptive OBL method. We can see that as the water mole fraction decreases, the CaCO 3 deposition occurs based on the equilibrium quotient of the reaction. Since the K sp for this reaction is taken unrealistically high, the initial mole fraction of calcite is zero at the initial conditions. As the CO 2 front proceeds, water is vaporized and CaCO 3 precipitates.
Next, the rest of the layers were activated in the Brugge model. The pressure profile along with mole fraction of CO 2 and calcite are shown for three layers in Fig. 7.
From Fig. 7, it is clear that the CO 2 front has mostly progressed in the layer with the higher permeability. The change in the fluid porosity with precipitation is also determined as described in Appendix B but its effect on permeability has not been considered in the performed simulations and will be investigated in future research.

Conclusions
We present a simulation framework which includes chemical equilibrium reactions fully coupled with flow and compositional transport. The newly developed nonlinear formulation reduces the component-based governing equations to the element-based mass balance equations and provides an effective coupling between chemical equilibrium reactions, thermodynamic equilibrium, and compositional flow and transport. For the computation of phase behavior, we expand the negative flash technique and include chemical reaction into the nonlinear loop. This approach was rigorously validated for a single equilibrium reaction coupled with the two-phase thermodynamics.
For reactive flow and transport, we use an element molar formulation where the secondary unknowns are fully resolved by the extended multiphase flash combined with thermodynamic and chemical equilibrium relations. Recently proposed adaptive operator-based linearization (OBL) technique is employed to solve the nonlinear mass balance equations in a fully implicit manner. The developed multiphase flash, together with other nonlinear relations, provides an effective parametrization of operators in governing equations based on the element composition. The reactive-compositional simulations were performed in the Delft Advanced Research Terra Simulation (DARTS) framework based on the pressure and element compositions. This choice can significantly improve the run time of simulation due to a reduction in the size of the algebraic system in realistic settings. An additional improvement in simulation time can be achieved due to an adaptive parametrization of the reactive multiphase flash result in the OBL approach.
Simple mineral dissolution and precipitation system with two equilibrium reactions are modeled using DARTS framework. To demonstrate the applicability of the developed framework to realistic problems, we modify the Brugge field and transform it to the problem of CO 2 storage in the depleted gas field with one CO 2 injector and two producers. We successfully model the process of calcium carbonate deposition due to the changes in the composition of the dynamic system. The approach will be later extended to include chemical kinetics and permeability alteration to model practical field applications. Furthermore, benchmark studies will be performed using the newly developed framework by comparing the experimental and simulation results. convergence rate can be slow. Figure 8 above shows the solution for a three-phase three-component system. This technique provides a robust mechanism to resolve the thermodynamic equilibrium by solving the RR equation. To improve the convergence rate, we include a Newton loop inside the bisection. In the case when the solution based on Newton iterations fails to converge or moves out from the physical bounds, the solver switches from Newton method back to bisection.
The comparison between the original bisection method and the combination of bisection and Newton methods is given in Table 1. We parameterized the complete compositional space for a three-phase three-component system. The first column in Table 1 shows the total number of iterations required for 18997 flash computations with a pure bisection strategy. The second column uses a combination of Newton and bisection approaches for only one of the two RR equation. The last column shows the result when applying Newton for both the RR equations. Here, the residual tolerance was set to ε = 10 −12 .
From the table, we can see that there are some waste Newton iterations in the case when Newton solver fails to converge and, hence, the system switches back to bisection. Even though some of the Newton iterations were wasted, the total number of iterations for the system with Newton and bisection is almost an order of magnitude less than the total iterations for the pure bisection method. This shows that including Newton method along with bisection for phase computations can significantly improve the computational efficiency of the thermodynamic flash solver and can be extended further to find more robust and effective methods to solve the RR equation. The compositional diagram generated using the flash calculations is shown as the third figure, where the yellow tie-triangle represents a threephase region, the blue regions are the single-phase regions, and the green regions are the two-phase regions.

Appendix B: Solid phase treatment
This appendix describes the concept of fluid and reactive porosity which combine to give the total porosity of the system. We treat the volume, occupied by the mineral component, as a part of the pore volume. The classic porosity, which represents the volume occupied by fluids, is what we call a fluid porosity. For a system without chemical reactions, the porosity only varies with changes in pressure due to the compressibility of the rock. But in the case of chemical reactions when mineral precipitation and dissolution are present, we have continuous changes in the pore space depending on the concentration of minerals. Therefore, the reactive porosity varies with mineral mole fraction. The bulk volume of the model is defined here by three parameters: non-reactive volume (V nr ), reactive volume (V r ), and the pore volume (V φ ). The non-reactive volume is the part of the rock which is not involved in any of the chemical reaction hence its volume is always constant. The reactive volume is the mineral part of the rock, and the pore volume is the volume occupied by the fluids in the rock; both of these volumes are changing depends on the amount of mineral present. Therefore, we can define the total volume as the sum of all the three components as shown below From the above equation, it can be seen that the total porosity of the system is always constant irrespective of the concentration of mineral. If there is less mineral deposited, the pore volume will be higher; otherwise, the reactive volume will be higher. The mineral saturation is defined as Using the definition of total porosity, the above equation can be written as The relation for fluid porosity can be written as where M stands for the number of mineral species. Therefore, knowing the initial total porosity of a control volume, we can calculate the porosity based on the saturation value of the mineral. As the reaction progresses, the dissolution/precipitation process occurs which alters the fluid porosity φ p . This porosity value can be used to determine the permeability values using empirical relations and update the velocities in the governing equations.

Appendix C: Simulation parameters
In this section, we describe the fluid and rock properties as well as numerical parameters of OBL used in the simulation. Table 2 shows the rock and fluid properties. The rock is considered low compressible and the total mineral and fluid porosity are taken as 0.3.  Table 3 shows the initial pressure condition of the reservoir and the OBL parameters which are used for adaptive simulation.  Table 4 shows the initial and injection reservoir mole fractions.  Table 5 shows the thermodynamic partition coefficient and the chemical equilibrium constants for two equilibrium reactions. These values were slightly elevated in comparison with real physical values.