Modeling Phase Behavior in Complex Systems Using the Gibbs-Helmholtz Constrained Equation of State

Modeling the phase and thermo-physical behavior of multi-component fluid systems using a cubic equation of state (EOS) is important for many industrially relevant applications. For example, simulations involving hydrocarbon recovery or carbon dioxide sequestration in geological formations rely on thermodynamic models to determine the phase behavior and density of the relevant reservoir fluid mixture. Accurate models for the density phase behavior of these fluids is required to make reliable predictions of the hydrocarbon production capability, or the carbon sequestration capacity of a given formation. Cubic EOS models have remained the industry standard in thermodynamic modeling for fluid phase behavior for the last 50 years, probably due to the relative ease with which these models can be implemented and generally acceptable accuracy for many systems. However, the current models include empirical parameters that are regressed to pure component saturated liquid density and saturation pressure data, as well as binary (and higher order) mixture composition data. Conducting experiments to collect pure component data for the regression of EOS parameters can be expensive, especially at elevated temperatures and pressures. Further, collecting mixture data over the entire composition space at varying temperatures and pressure quickly becomes intractable (especially for mixture of three or more components). Clearly, a cubic EOS model which requires less data for the regression of parameters is desirable. The multi-scale Gibbs-Helmholtz constrained (GHC) equation of state (EOS) is an innovative approach to EOS modeling which uses molecular scale information about the component(s) of interest to calculate bulk scale EOS parameters. More specifically, the molecular attraction parameter (‘a’) in the two parameter Soave-Redlich-Kwong (SRK) EOS is calculated as a function of temperature using an expression derived from the Gibbs-Helmholtz equation (a classical thermodynamic relationship). Further, the GHC expression for the attraction parameter incorporates molecular level information using results of isobaric-isothermal Monte Carlo (NPT-MC) molecular simulations. In this work, some aspects of the GHC EOS performance and thermodynamic consistency are investigated. Further, novel modeling frameworks are developed for the application of the GHC EOS to systems capable of forming simple structure I (sI) gas hydrates and molecular salts. Finally, the GHC EOS is incorporated into a fully compositional and thermal reservoir simulator. The GHC EOS is then used as the thermodynamic model for the underlying reservoir fluid in novel reservoir simulations relevant to enhanced oil recovery, carbon sequestration, and groundwater contamination modeling.


Thermodynamic Consistency: Theory and Verification
In this section, the basic theory behind the multi-scale GHC equation and verification of thermodynamic consistency are presented.

Theory
Starting (1. 4) where is the critical temperature, is the critical pressure, ( , ) = 0.42748 2 2 / , is the molecular co-volume, is the gas constant, and is a molecular length scale internal energy of departure for the liquid phase given by = ( , ) = ( , ) − ( ), where ( ) is the ideal gas internal energy. Note that serves as a natural bridge between the molecular and bulk phase length scales.
When eq. 1.4 is used to determine the energy parameter in Eq. 1.1, the EOS is called the multi-scale Gibbs-Helmholtz Constrained (GHC) equation of state.
The derivation for mixtures, follows essentially the same steps as that for pure components. Starting with the fundamental expressions where denotes a mixture property and is the pure component fugacity coefficient In summary, the multi-scale GHC equation of state starts with the exact same expression for ln( ) or ̂, as the SRK EOS, calculates the temperature derivative of , and then uses the Gibbs-Helmholtz equation to constrain or to derive the novel expressions for the energy parameter given by eq. 1.4 and 1.8. The question central to this article is whether or not the a posteriori development of the multi-scale GHC expressions for and preserve thermodynamic consistency. As shown in section 1.2.3, and illustrated throughout this article, that depends on the methodology for estimating and .

Internal Energies of Departure from NPT Monte Carlo Simulations
Before turning to the issue of thermodynamic consistency, the procedure for constructing look-up tables of pure component internal energies of departure using NPT Monte Carlo simulations is briefly described. Pure component look-up tables contain discrete sets of as a function of and , at varying intervals of temperatures from 250 -600 K (or 750 K in the case of water) and pressures from 1 to 600 bar. In generating one data point for a given and , we typically use a small number of particles (e.g., N = 32 particles) and run 4 sets of 50,000 equilibration cycles + 400,000 sampling cycles. Results for all 4 sets are then averaged and entered as a single data point in a look-up table. This procedure is described in considerable detail in   The key points in generating and using molecular length scale information are (1) look-up tables contain discrete sets of data with uncertainty, (2) defining for points between values of and in look-up tables is open and can be done in any number of ways (e.g., by averaging, linear interpolation, etc.), and (3) each way of defining impacts thermodynamic consistency differently.

Verification
To verify the thermodynamic consistency of the GHC equation, the following for mixtures. It is important for the reader to understand that eqs. 1.11 and 1. 12-13 only hold locally (at a given state of the system) and therefore can only be verified at a given thermodynamic state. This means that while the derivation of requires that = ( ) , the condition of thermodynamic consistency only holds locally. This means that = ( ) can be relaxed locally and we are free to define the functionality of differently as long as we do not violate thermodynamic consistency.
The procedure for numerically verifying thermodynamic consistency is as follows: 1) Fix for pure components or for mixtures. Set to some small number (e.g., 10 −4 ).
2) Choose . The critical step with regard to thermodynamic consistency is the methodology used to determine ( , ), which can be conveniently separated into two parts -a temperature part and a pressure part. Some of the ways in which can be determined are as follows: 1) For any temperature, , linear interpolation between the appropriate isothermal data in a look-up table can be used to determine ( , ).
2) There are various ways to include the pressure effect once the temperature effect has been calculated. a) Average over the pressure range of interest, say ( ). In this case, ( , ) = ( ) for all pressures at the given temperature, . There will generally be a different average value at each temperature.
c) For the given pressure , use linear interpolation in pressure.
We recommend using linear interpolation in temperature, as in step 1, followed by setting ( , ) = ( , ), which is the equivalent of constructing isothermal staircase functions for ( ).
Clearly 2a) yields thermodynamic consistency since the use of average ( ) at all pressure for any given implies depend only on temperature. Also, 2b) yields thermodynamic consistency because the condition for thermodynamic consistency, = ( ) , involves only point functions, and thus holds locally.
However, using linear interpolation in pressure is not recommended.

Pure Components
In this sub-section three separate cases for determining for pure components are described.
Case 1: Direct Use of in Look-Up

Mixtures
In this sub-section, verification of the thermodynamic consistency of the multiscale GHC equation for mixtures is presented.

Methane-Water
Mixtures of light gas and water are usually challenging so we have selected the methane-water system as a first example to illustrate the thermodynamic consistency of the multi-scale GHC equation for mixtures. were computed by forward finite differences using ∆ = 10 −3 . Note that in all cases the condition of thermodynamic consistency is satisfied to a tolerance of at least = 10 −4 . This next section of the paper provides a brief summary of the underlying mathematical framework for the pressure dependence of . Use of the methodology ( , ) = ( , ) is equivalent to defining as a staircase function in pressure, where the staircase function is simply a mathematical construct that permits any quantity, in this case internal energy of departure, to vary with respect to any independent variable (e.g., pressure) and has the following mathematical properties: 1) It is not continuous everywhere. In fact, the smaller the step width (run) and step height (rise) of the staircase, the finer the granularity and number of points of discontinuity for a given range.
2) It has either a derivative or one-sided derivative of zero everywhere, except for at any discontinuity.   Note that by adopting an isothermal staircase representation of the pressure dependence of the internal energy of departure the following statements are true: 1) clearly depends on pressure but the derivative (or one-sided derivative) of ( , ) with respect to at constant temperature, ( / ) = 0.
2) The staircase can have any desired level of granularity.
3) Finer granularity can be adding more NPT Monte Carlo simulations data to the look-up table. Each step of the staircase would then have a smaller height (rise) and smaller run (width). Fig. 1.1 can also be interpreted as quasi-linearization of with respect to , as described by Bellman (1973). Quasi-linearization in this context is simply a Taylor series of isothermal ( , ) expanded about a discrete data point with all pressure derivatives set to zero. 5) Since ( / ) = 0 is everywhere zero, clearly implies that ( / ) = 0

4) The staircase function in
and therefore the multi-scale GHC equation is thermodynamically consistent.
The proposed staircase procedure for isothermal values of molecular with respect to pressure can be contrasted to other possible approaches where 1) A least-squares, chi-squared, or maximum likelihood function is used to fit internal energies of departure so that bulk phase can be treated as a smooth function of pressure.
2) Multiple NPT Monte Carlo simulations are run at different pressures, averaged, and then numerical differentiation is used to approximate bulk phase ( / ) .
In contrast and regardless of granularity, using a staircase representation of , gives molecular values of ( / ) that are always zero for practical uses (eg. when evaluated at the pressure and temperature of interest).
Recall that molecular internal energies of departure for mixtures, , are calculated using the simple linear mixing rule = ∑ (1.14) where and are the mole fraction and internal energy of departure for the i th component in the mixture. If each in eq. 1.14 is given by an isothermal staircase function in pressure, then is also an isothermal staircase function. Moreover, there is no restriction that the data points for each component need to coincide or that the height (rise) or width (run) of each staircase function for be the same. These facts are easily proved. However, in our opinion, it is more informative and instructive to give a clear geometric interpretation of these facts. Thus Fig. 1

Density & Phase Equilibrium Results
A given equation of state must also provide accurate densities and phase equilibrium to be useful in practice. Table 1.14 gives a comparison of liquid density predictions for the SRK equation  with Peneloux volume translation (SRK+) (Peneloux et al., 1982), the volume translated Peng-Robinson (VTPR) equation , and the GHC equation with experimental data for a number of different compounds.  Kell and Whalley (1965). Note again that the average AAD% error in molar volume is quite small -less than 1%.

Pure Components
Extrapolating the results for average in Table 1

Mixtures
In this sub-section, density and phase equilibrium results for a number of mixtures are presented and, where possible, calculated results are compared with experimental data.

Density
In this section, density predictions for CO2-water and NaCl-water are presented to illustrate that the GHC equation provides accurate mixture densities in addition to being thermodynamically consistent. Density results for the GHC equations are also compared with predictions using the SRK+, PSRK, VTPR, and, where applicable, the ePSRK equations of state, and with experimental density data.  Chen et al. (1977) Nonetheless, Table 1.17 shows that the multi-scale GHC equation predicts the experimental density of the mixtures investigated here with < 1 % AAD.

Phase Equilibrium
The phase equilibrium calculations of mixtures exhibiting vapor-liquid and vapor-liquid-liquid behavior are described in the next section.

Discussion of Results
Theoretical analysis and numerical verification using finite difference values of ( ∆ ∆ ) ( ∆∆ ) , were used to clearly show that the multi-scale GHC equation is thermodynamically consistent for pure components and mixtures. The novel idea of using an isothermal staircase function representation of molecular scale internal energies of departure, , as a function of pressure was proposed and illustrated both geometrically and by reporting tabulated values of NPT Monte Carlo simulation data.
We have shown that the GHC approach relaxes the assumption that the energy parameter is pressure independent (after the fact), and with an appropriate methodology for adjusts the energy parameter in a way that preserves thermodynamic consistency and is supported by the fact that the condition of thermodynamic consistency, = ( ) , involves only point functions, and holds locally. the densities of hexagonal water ice and structure I gas hydrates as well as phase equilibrium for methane-water and CO2-water mixtures.

Introduction
The general warming of land and oceans around the world has spawned interest in understanding long term environmental impacts associated with the melting of ice sheets and thawing of permafrost regions. Permafrost, which accounts for approximately 25% of the land mass in the Northern Hemisphere, is the general term used to describe land that remains below 0 C for two or more years at depths that can range from less than a meter to around 1000 meters. While temperatures and pressures in permafrost regions vary, typical ranges are -20 C to 5 C and 0.1 to ~10 MPa respectively.
Phase behavior in permafrost regions is quite complicated and can involve light gases (e.g., methane and carbon dioxide), liquid water, brines, hexagonal ice, and gas hydratesdepending on conditions. Understanding physical propertys and phase behavior in permafrost regions is important for several competing reasons: 1) The large amount of methane sequestered in permafrost is a potential source of low carbon fuel if it can be captured. According to Kvenvolden (1998) there is an estimated 2.1 x 10 16 standard cubic meters (SCM) of methane contained in hydrates on the ocean bottom, which is twice the energy sum of all other fossil fuels on Earth. There is also an estimated 7.4 x 10 14 SCM in permafrost regions (see, MacDonald, 1990).
2) Thawing of permafrost increases microbial activity, which in turn has the potential to release large amounts of greenhouse gases (CO2 and methane).
3) Deeper locations in permafrost have been suggested as potential sites for carbon storage.
The focus of this work is the modeling of physical properties, specifically density, and phase behavior of mixtures of light gas and water at conditions typically found in permafrost regions (i.e., temperatures ranging from 250 to 280 K and pressures from 0.1 to ~10 MPa) using an equation of state. It is demonstrated that the multi-scale Gibbs-Helmholtz Constrained (GHC) equation , can be used to predict thermo-physical properties and phase equilibrium in systems involving gas, liquid, hexagonal (1h) ice and gas hydrate. Accordingly, this paper is organized in the following way. Section 2.

Ice and Gas Hydrates
Most people have some intuitive feel for or understanding of ice, and hexagonal ice is the only ice that occurs naturally on Earth. As its name implies 1h ice has a hexagonal lattice structure. Gas hydrates, on the other hand, are less well known and understood, despite the fact that they were discovered in 1810 (see, Ghiasi, 2012) and later identified as the cause of flow assurance problems in gas lines by Hammerschmidt (1934 cages, two 4 3 5 6 6 3 cages, and one 5 12 6 8 cage and have 34 water molecules. Here the general notation e n describes a polygon such that has e edges and n faces.

A Brief Review of Relevant Literature
Many traditional equations of state (EOS) have been used to model various thermo-physical properties and phase equilibrium involving mixtures of light gases and water including the Soave-Redlich-Kwong (SRK) equation , the Peneloux volume translation (Peneloux et al., 1982) of the SRK equation (SRK+), the Peng-Robinson (PR) equation , the volume translated PR or VTPR equation , and various forms of the Statistical Associating Fluid Theory (SAFT) equation (Huang & Radosz, 1991). Other equations of state such as the cubic plus association (CPA) equation Voutsas et al., 2000) and the Elliott-Suresh-Donahue (ESD) equation (Elliott et al., 1990) have also been used. There are also specialized equations of state for water such as the The IAPWS-95 model for water has been extended to 1h ice by Feistel & Wagner (2006) and is also a Gibbs potential equation of state. Like the model for liquid water (Wagner & Pruss, 2002), the extended IAPWS-95 equation for 1h ice has been fit to a large amount of experimental data (i.e., 14 independent parameters fit to 522 data points from 32 separate categories of experimental data). This experimental data includes specific Gibbs free energy data, specific entropy data, (∂p/∂T) data along the ice melting curve, heat capacity data, volumetric and volume-temperature derivative data, isentropic compressibility data, and isentropic compressibility-temperature derivative data. See use a more traditional solid-liquid equilibrium formulation in which the fugacity of ice is expressed in terms of the fugacity of sub-cooled water and an exponential correction based on enthalpy and volumetric differences due to fusion.
The most commonly used models for predicting the properties of gas hydrates are based on the cell theory model developed by van der Waals and Platteeuw (1959) and given the acronym vdWP. A detailed derivation of the model from statistical mechanics as well as a very good survey of the relevant hydrate literature can be found in textbook by Sloan and Koh (2007). Here a brief overview of the modifications of the vdWP model is given starting with the work of Parrish and Prausnitz (1972) who extended the model to allow for the practical calculation of hydrate dissociation pressures for mixtures of guest molecules. Klauda and Sandler (2000) presented a fugacity-based model that removed the need for reference energy parameters for the empty hydrate lattice and relaxed the incorrect vdWP assumption that the volume of the crystal lattice is independent of guest molecule type. A similar model was developed by Ballard and Sloan (2002) and Jager and coworkers (2003). More recently, Bandyopadhyay and Klauda (2011) presented an updated fugacity model that uses the PSRK equation to represent the fluid phases in equilibrium with the hydrate. Other modifications include models that account for multiple cage occupancy (see, Klauda and Sandler, 2003;Martin, 2010), those that correct for the effects of guest-guest interactions on the Helmholtz energy (Zhdanov et al., 2012), and others (e.g., Jäger et al., 2013) that adapt the model originally developed by Ballard and Sloan (2002) so highly accurate reference EOS can be used. This last approach has also been tested against a large data set for CO2 hydrate.
The vdWP model with all of its recent adaptations has enjoyed great success and is used heavily in industry because of its accuracy and relative simplicity. However, it is a semi-empirical model that requires adjustable parameters to be regressed to experimental data. For this reason, it is difficult to make predictions about the phase behavior and properties of potentially hydrate-forming systems outside the ranges of experimental data.

The Multi-Scale Gibbs-Helmholtz Constrained (GHC) Equation
The multi-scale Gibbs-Helmholtz Constrained (GHC) equation, which is a modification of the Soave form of the Redlich-Kwong (RK) equation , is given by (2.1) where is pressure, is temperature in Kelvins, is the molar volume, is the universal gas constant and and are the energy and molecular co-volume parameters respectively.

Pure Components
The multi-scale GHC equation is a radically different approach to EOS modeling that is based on three simple ideas: 1) The molecular co-volume for the liquid, , is set equal to the molar volume of the pure solid or high density (glassy) liquid.
2) The liquid phase energy parameter, ( , ), is constrained to satisfy the Gibbs- is what makes the GHC equation a multi-scale equation of state.

Non-Electrolyte Mixtures
For  Lucia et al., 2012a). Mixture critical properties are calculated using Kay's rules where is either critical temperature or critical pressure.
More detailed information associated with the GHC equation, its derivation, expressions for pure component and partial fugacity coefficients, and density and phase equilibrium predictions for a number of systems including pure liquid water, water-light gas mixtures, n-alkanes of varying chain length, and electrolyte solutions can be found in references , Lucia et al, 2012aLucia et al., 2012b). More recently, Lucia and Henley (2013) have shown that the GHC equation is thermodynamically consistent by demonstrating that it satisfies the relationship given by ( ) = .

Modeling Hexagonal Ice Using the GHC Equation
To model hexagonal ice in the GHC framework, values for and , are needed along with a means of calculating the pure component fugacity of ice.

The Molecular Co-Volume and Internal Energy of Departure for Ice
Two straightforward but important physical interpretations are required to apply the multi-scale GHC equation (i.e., eqs. (2.1) and (2.2) with appropriate values of and ) to 1h ice: 1) In our opinion, the correct value for the molecular co-volume for 1h ice is the molar volume of liquid water (i.e., = 18.015 3 / ) because ice melts under pressure.
2) The internal energy of departure for 1h ice can be calculated in two different ways: a) Using the expression These two methods for determining , are described and compared in considerable detail in Appendix A.

A Reference State for Ice
The fugacity coefficient of 1h ice must be different from that of liquid water at phase equilibrium. Appendix B shows that the natural log of the fugacity coefficient for hexagonal ice can be rigorously expressed in the form where and are the fugacity coefficients of 1h ice and liquid water respectively and ∆ can be interpreted as a measure of the impact of long-range structure on the fugacity coefficient of ice and is assumed to be constant. This leads to the following expression for the standard state fugacity of ice

Modeling Gas Hydrates Using the GHC Equation
Although the approach for modeling gas hydrates using the GHC equation is similar to that of ice, it involves considerably more detail. This is because 1) Gas hydrates are unusual heterogeneous structures so it is not at all clear that the usual mixing rules apply.
2) Empty hydrate cages are not stable.
3) Standard states only apply to pure components.
The approach presented here starts with the idea of estimating the bulk density of an empty hydrate. Why? The reason is that some physically sensible way of estimating the fugacity (or Gibbs free energy) of pure or empty hydrate is needed in order for this framework for addressing phase equilibrium to be consistent with classical theory. In order to predict empty hydrate density using the GHC equation values of the internal energy of departure for an empty hydrate, , , and molecular co-volume, ℎ , are required.

The Molecular Co-Volume and Internal Energy of Departure for Hydrate
The value of the empty hydrate molecular co-volume is taken to be the same as that of a completely filled hydrate and given by ℎ = 0.148148 4 + 0.851852/ (2.9) where 4 = 29.614 3 / and where the subscript denotes water and the superscript ℎ represents hydrate. While it may appear that eq. (2.9) is the same as the usual linear mixing rule for , it is not. Eq. (2.9) is interpreted as the high pressure limit of separable gas and solid phases in a hydrate. Additionally, within the GHC approach to hydrates, eq. (2.9) is used for both empty and filled gas hydrates regardless of the gas occupancy in the cages. Thus, there is no mole fraction weighted average of molecular co-volume for hydrate phases; is held fixed at the value given in eq. (2.9).
Because empty hydrates are not stable, considerable care must be taken to get estimates of internal energies of departure , . To do this, the following procedure was followed: 1) Placed guest molecules in the hydrate cages.
2) Equilibrated the simulations with the guest molecules present.
3) Removed the guest molecules.
4) Ran production cycles with restricted volume moves so that the empty cages would not collapse.
NPT Monte Carlo simulations were performed for empty hydrate using the aforementioned procedure over relevant ranges of temperatures (250 to 280 K) and pressures (1 to 100 bar). Values of , ranged from −4.65 10 5 to −4.58 10 5 3 / .

Gas Hydrate Density
Next, a theory for predicting density for physically meaningful gas hydrates, which have gas molecules occupying the cavities or cages was developed. Because gas hydrates are heterogeneous structures, the density of an S1 gas hydrate is calculated using the expression ℎ = (5.75 + ) /5.75 (2.10) where is the fractional occupancy of gas in the hydrate cages and ranges from 0 to 1 and where is the empty hydrate density computed from the GHC equation [i.e., eq.
(2.2)]. Note that eq. (2.10) is still an application of the multi-scale GHC equation to gas hydrates and, in our opinion, provides a straightforward, yet physically meaningful, way of computing hydrate density. Moreover, eq. (2.10) covers the full range of gas occupancy; thus, the more gas there is in the cages, the denser the hydrate. At the same time, it avoids many of the complications associated with composition dependence in determining hydrate properties like fugacity coefficients.

A Reference State for Pure Water in Gas Hydrate
To calculate a standard state for pure water in a hydrate phase we let ℎ = 1 and use a procedure identical to that in Appendix B for ice. This leads to the expression 0,ℎ = (2.11) where = exp (∆ ) and = exp(∆ ). The quantity ∆ = − represents the difference in long range structure between 1h ice and empty hydrate and can be interpreted as a second structural contribution. Appendix C gives all of the details of the derivation of eq. (2.11).

Numerical Results
In this section, numerical results are presented for density and phase equilibrium for light gas-water mixtures in regions where ice and gas hydrate can form. All numerical results presented in this section use the GHC equation of state (except for comparisons). Critical property and other relevant data are summarized in Appendix D.

Methane-Water Mixtures
Methane and water can exhibit a number of different phase equilibrium -gasliquid, gas-ice, gas-ice-hydrate, and gas-liquid-hydrate at conditions typical of permafrost.

Methane-Water in Gas-Liquid Equilibrium
Numerical results for methane-water have been reported by Lucia and Henley (2013), compared to experimental data of Servio & Englezos (2002), and are reproduced here for the reader's convenience.

Hexagonal Ice Density and the Ice-Water Melting Curve
Before studying gas-ice equilibrium and to validate the theory developed in sections 2.4.1 and 2.4.2 and Appendix B, results for the density and ice-water melting curve as predicted by the GHC framework are presented.    predicting the 1h ice melting curve.

Figure 2.1: Melting Curve for Ice Predicted by the IAPWS-95 and GHC Equations
The results in Table 2.2 and 2.3 and Figure 2.1 show that the expression for the standard state fugacity for water in a hexagonal ice phase derived in Appendix B leads accurate predictions of ice/water phase equilibrium and density over the region investigated.

Methane Gas-Ice Equilibrium
At low enough temperature and pressure, methane and water will exhibit gassolid equilibrium with essentially pure phases. However, it is important to understand that other 'false' equilibrium (e.g., gas-liquid equilibrium) can be found at the same conditions. Table 2.4 presents results for the phase equilibrium of a mixture of 10 mol% methane and 90 mol% water at 272 K and 1 bar. At the given temperature, the pressure is too low for methane hydrate to form so this example represents a reasonable test of the capability of the multi-scale GHC EOS to find gas-solid equilibrium. Also, note that the correct global minimum Gibbs free energy solution is methane gas with a very small amount of water vapor in equilibrium with 1h ice (GSE). However, there is a false gas-liquid equilibrium (GLE) solution with a value of G/RT that is only slightly higher than the gas-ice equilibrium solution. It is also possible to find a false gas-liquid-ice (GLSE) equilibrium at these conditions with a value of G/RT between that of the GS and GL equilibria and a very small amount of solid (~3.5 mol%/mol feed).

Methane Hydrate Density
The next example gives numerical results for GHC-predicted density of empty and filled hydrate and, in doing so, tests the validity of the theoretical material presented in sections 2.5.2 and 2.5.3 and Appendix C. Empty and filled methane hydrate densities are shown in Table 2.5. Table 2.5 also shows a comparison of the GHC-predicted filled methane hydrate density with the analytical expression given in Sloan and Koh (2007) and the Sloan and Koh expression with the temperature and pressure dependence of empty cell volume given by Klauda and Sandler (2000) included. Note that there is reasonably good agreement among all three methods.  Klauda & Sandler (2000) temperature & pressure dependent empty cell volume.

The Quadruple Point for Methane-Water Equilibrium
Below the melting point of ice, the phase equilibrium that exists is either methane gas-ice or methane hydrate-gas/ice equilibrium. That is, at low pressure, below pressures for which hydrate can form, methane gas is in equilibrium with 1h ice. See Table 2.4 for an example of low pressure GSE. As one raises the pressure at fixed temperature, methane hydrate will form and gas hydrate will be in equilibrium with gas and/or ice -depending on the overall amounts of methane and water in the system.
Above the melting point of ice the equilibrium changes from gas-liquid to methane hydrate-gas/liquid as a function of increasing pressure. One of the more challenging tasks in the neighborhood of the ice melting point is predicting the quadruple point -the temperature and pressure at which four phases (ice, methane hydrate, liquid water, and gas) co-exist. Anderson (2004) reports a quadruple point of 272.9 K and 25.63 bar.  Because the quadruple point is a singular point, it is very difficult to compute accurately.
However, note that the GHC prediction of the quadruple point gives all three two-phase equilibrium that have dimensionless Gibbs free energies that are quite close (i.e., differing by < 10 -3 ). Note also the predicted fractional occupancy of the gas in the hydrate phase is 0.9437, which is quite close to the value of occupancy of 0.9461 predicted by the correlation in Parrish and Prausnitz (1972) but less close to the value of occupancy of 0.90 ± 0.01 predicted by Gibbs ensemble Monte Carlo simulation (GEMC).

Sub-Cooled Methane Hydrate-Ice Equilibrium
This last methane-water example in this article is used to illustrate that phase equilibrium and gas occupancy for conditions away from phase boundaries can be determined using the theoretical framework described in sections 2.3 -2.5 and Appendices A, B, and C. Consider permafrost conditions for a mixture of 13 mol% methane and 87 mol% water at 272 K and 27 bar. For the given problem (1) there is an excess amount of water, (2) methane hydrate should be in equilibrium with pure ice, and (3) the resulting phase equilibrium solution is away from the phase boundary.
Numerical results for this example are shown in Table 2.7. It is important to describe the phase equilibrium computations for this ice-gas hydrate equilibrium flash solution in a bit more detail. The water-to-methane ratio in the feed is equal to 6.6923, well in excess of the value of 5.75 for a completely filled methane hydrate. Therefore, if methane hydrate is one of the equilibrium phases, there should also be excess water and this implies that either liquid water or ice will also exist. In addition, from Table 2.7 it is clear that ice should be the phase in equilibrium with methane hydrate at the given conditions. This is true because the false GSE solution has a lower value of G/RT than that for the false VLE solution. Moreover, in any icemethane hydrate equilibrium, methane is confined to the hydrate phase. Thus, defining equilibrium becomes a bit more challenging since, in theory, there is only a single chemical potential to use to define conditions of phase equilibrium. This is where the use of structural contributions to the ice and empty hydrate standard state fugacities for water become important. To define equilibrium the following condition is used It is also important to understand that eq. (2.12) is singular unless trace amounts of methane are permitted in the ice and empty hydrate phases. However, even if trace amounts of methane are present in these phases to avoid singularity, the problem is still so nearly singular that quadratic acceleration (see, Eq. (A2), p. 2562 in Lucia and Feng, 2003) is required to get convergence to a reasonable tolerance.
To give the reader some appreciation for this last point, we have shown the calculated equilibrium solution to eq. (2.12) in the presence of trace amounts of methane at 272 K and 27 bar in Table 2.8.  (7.23x10 -5 , 0.999928) Using quadratic acceleration, the computations converge to the solution shown in Table   2.8 in 20 iterations. However, note how very close in composition both phases in Table   2.8 are; both contain less than 0.01 mol% methane. The two critical pieces of information that come from solving eq. (2.12) in this manner are the ice and hydrate phase fractions. From these phase fractions, it is straightforward to determine the methane occupancy, gas hydrate density [i.e., using Eq. (2.10)], and corresponding value of G/RT for the hydrate phase -as well as G/RT for the ice-methane hydrate equilibrium solution.
For this particular example, the fractional occupancy of methane in the hydrate phase using Gibbs-ensemble Monte Carlo simulations was also calculated for the purpose of comparison. The value obtained from GEMC was 0.89 ± 0.01. Table 2.7, on the other hand, shows that the calculated fractional occupancy using the GHC-based framework developed in the work is 0.8759, whereas the occupancy determined using the correlation in Parrrish and Prausnitz (1972) is 0.9497. Thus, it is clear that the proposed GHC-based framework for gas hydrates gives reasonable values of gas occupancy directly from phase equilibrium.

Carbon Dioxide-Water Mixture
Like methane-water mixtures, carbon dioxide-water mixtures can exhibit a number of different types of phase equilibrium -gas-liquid, liquid-liquid, gas-ice, hydrate-ice, and so on.

Gas-Liquid Equilibrium
In previous work, Lucia and co-workers (2012a)

Carbon Dioxide Gas-Ice Equilibrium
CO2-ice equilibrium is straightforward to compute using the GHC-based framework developed in this work. However, depending on conditions, there can be many other false equilibrium with values of G/RT that are close to the global minimum value of G/RT, as illustrated in Table 2.9. The presence of many equilibrium makes the prediction of the correct solution quite challenging.

Carbon Dioxide Hydrate Density
The density of CO2 hydrate is considerably higher than that of methane hydrate due to the fact that CO2 is much heavier than methane. However, because CO2 hydrate is also a structure 1 hydrate, eq. (2.10) can be used to calculate hydrate density. Table   2.10 compares GHC-predicted CO2 hydrate densities with densities calculated using the analytical expression given in Sloan and Koh (2007) with and without the Klauda and Sandler (2000) correction for empty hydrate cell volume.  al. (1999), the molecular weight of CO2 hydrate is 21.59 g/gmol. From Eq. (2.10), the actual hydrate density for a fractional CO2 occupancy of 0.9583 is 0.051679 mol/cm3, which when multiplied by the molecular weight of the CO2 hydrate gives a mass density of 1115.57 kg/m 3 . This GHC-predicted value of CO2 hydrate density represents an error of 1.43%. For the exact same conditions, the analytical expression in Sloan and Koh (2007) gives a density of 1101.02 kg/m 3 while the expression in Sloan and Koh with the Klauda and Sandler (2000) correction for empty hydrate cell volume results in a mass density of 1162.25 kg/m 3 .

Phase Equilibrium in Regions Where CO2 Hydrates Can Exist
One of the interesting differences between methane hydrates and CO2 hydrates is the presence of liquid CO2 at high pressure so we consider an example that would potentially result in storage of CO2 in a hydrate phase. Let the temperature and pressure Finally, for hydrate in equilibrium with fluid phases (and not 1h ice), the conditions defining equilibrium are the equality of chemical potentials for CO2 and water in all phases. Numerical results are shown in Table 2.11.  Table 2.11 is consistent with the hydrate phase sinking in water.
Finally, the GHC-predicted density and fractional occupancy agree quite well with the density of 1100 kg/m 3 and occupancy of 0.9583 reported in Brewer et al. (1999).
The numerical aspects associated with computing the HLLE flash solution in Table 2.11 are somewhat complicated. Because there are two components and three phases, the solution is near singular and quadratic acceleration is required for reliable convergence. In this example, the HLLE flash computations converge to an accuracy of 10 -6 in the 2-norm of the equality of chemical potentials in 9 iterations.

Conclusions
The multi-scale Gibbs-Helmholtz Constrained (GHC) equation of state-based framework has been extended for determining densities and phase equilibrium in light gas-water mixtures at conditions where ice and/or hydrates can exist. This extended framework is built around the use of the multi-scale GHC EOS, and the novel ideas that (1) make use of the GHC equation to determine densities of ice and empty hydrate, (2) derive and employ structural contributions to the standard state fugacities of ice and empty hydrate in phase equilibrium, (3)

Appendix A: Internal Energies of Departure for Ice
Internal energies of departure for hexagonal ice can be determined either from for liquid water plus energies of fusion or by direct Monte Carlo Simulation.

Internal Energies of Departure for Hexagonal Ice from Energies of Fusion
Along the melting curve, internal energies of departure associated with fusion,

Simulation
We have also calculated internal energies of departure for 1h ice using the TIP4P-Ew force field model initialized from a hexagonal structure. The optimized force field parameters for liquid water used in the MC simulations of ice were taken from Horn et al. (2004) and are given in Appendix D. Figure A1 gives a snapshot of the crystal structure predicted by Monte Carlo simulation for N = 96 water molecules, = 272 and = 10 during the sampling (or production) phase. For these simulations, 50,000 equilibration cycles and 200,000 production cycles were used. The minimum energy structure for hexagonal ice at 0 is also shown in Fig It is also important to note that the estimation of for hexagonal ice from direct NPT Monte Carlo simulation is not restricted to the melting curve but can be used for any temperature and pressure. However, it is unclear if this distinction matters in practice since Table A1, which gives a comparison of The results in Table A2 show that there is strong agreement between the two methods for determining internal energies of departure for hexagonal ice. Specifically, the % , which is defined as % = 100| , ( ) − ( + , )/ , ( )|, is less than 3% and the standard deviation, , is small -despite the fact that estimations of , using the Clausius-Clapeyron equation are really not strictly applicable at conditions that are not on the melting curve. Note that Table A1 also shows that the pressure effect on , is small.

Appendix B: A Reference State for Water in a Hexagonal Ice Phase
The expression for ln in the multi-scale GHC framework for a pure component is where = / 2 and = / . Application of Eq. (B1) to 1h ice and liquid water gives Subtracting Eq. (B2) from Eq. (B3) and some algebra gives where ∆ denotes the difference between the natural log of the fugacity coefficients for water and 1h ice due to long-range structural differences between water and ice.
Note that Eq. (B4) is a rigorous relationship between and and permits the calculation of the difference between the two fugacity coefficients from information that is readily available from density calculations for 1h ice and water.
At conditions of ice-water phase equilibrium, where denotes the Gibbs free energy, the subscript denotes water and the superscripts and represent hexagonal ice and liquid respectively. Equation (B5) can be expanded in the form where is the universal gas constant. Because = = 1, Eq. (B6) reduces to 0, All that remains is to evaluate at some reference state, say 0 = 273.15 and 0 = 0.1013 , which is easily accomplished using information from density calculations at reference state conditions. Table B1 shows the values of the quantities needed in Eq.

Remark
It turns out that the value, = 3.18531, is very close to the second peak distance in the O-H partial radial distribution function for 1h ice (see Fig.1, p. 514 in Chau and Hardwick, 1998). However, it is unclear if there is a rigorous relationship between ∆ and the O-H partial radial distribution function or if this is just a coincidence.

Appendix C: A Reference State for Pure Water in a Gas Hydrate Phase
At conditions of ice-hydrate phase equilibrium, To calculate a standard state of pure water in the hydrate phase we let ℎ = 1 in Eq.
(C5), which means the hydrate cage is empty. This gives where the quantity ∆ = − represents a second structural component that captures the differences in long range structure between 1h ice and empty hydrate.
To calculate this second structural component contribution to the standard state for pure water in the hydrate phase we use the same strategy that we used for determining the structural component for 1h ice. That is, we apply Eq. (B1) to determine the natural log of the fugacity coefficients for ice and empty hydrate and simply subtract the latter from the former to get ∆ in terms of compressibility and equation of state parameters, and . The resulting expression for ∆ is which clearly shows that only EOS information is needed to determine ∆ .
Straightforward computations using the GHC equation of state shows that ∆ varies very little with temperature and pressure (i.e., values of between -2.51 and -2.21 over the temperature and pressure ranges shown in Table 1). Table C1 gives a sample calculation of ∆ at 270 and 50 . We also repeated the same derivation of the structural component for empty hydrate by using water as the reference condition. Specifically, we re-wrote Eq. (C5) as where ∆ = − now measures the structural contribution to the hydrate standard state with respect to pure liquid water. This, in turn, led to an expression very similar to Eq. (C7). That is, However, this choice for the expression for ∆ did not yield good results.   (3) and a unique formulation of the combined chemical and multi-phase equilibrium flash problem that accounts for salt deposition but decouples the chemical and phase equilibrium aspects of the flash.

Appendix D: Critical and Other Relevant Physical Property Data
Examples from real EOR and CO2 sequestration applications are presented.
Results clearly show that the proposed numerical approach is reliable, robust, and efficient and can be used to determine salt deposition in multi-phase flash problems.
Several geometric illustrations and numerical details are used to elucidate key points of the proposed approach.

Introduction: Motivation and Background
Groundwater aquifers and seawater generally contain a number of ions such as Na + , K + , Ca 2+ , Mg 2+ , Fe 2+ , Cl -, SO4 2-, CO3 2-, and others. At high enough concentrations, these ions can combine and form salts and precipitate out of solution causing operational difficulties in a wide range of applications. Salt deposited on heat exchange equipment can cause poor heat transfer. A common example of this is scaling on a household water heater from 'hard' water (i.e., water with high levels of calcium, magnesium and carbonate ions). In enhanced oil recovery (EOR) operations, production water containing a wide range of ions is often re-injected into a reservoir in order to avoid the use of municipal water and/or costly water treatment. However, re-use of production water can result in salt deposition, which in turn can plug pores, block injection wells, and/or re-direct or sometimes stop flow through porous media. In CO2 sequestration, deposition of calcium carbonate can be either a blessing or a curse, depending on where deposition occurs. Obviously, it is desirable if deposition in CO2 sequestration occurs far from injection well bores. Salt deposition is also important for reservoir modelers to understand, predict, and quantify so that reservoir simulation results correctly reflect physics. In EOR and CO2 sequestration problems, multiple equilibrium phases often exist within a reservoir due to the presence of light gases and water. Thus, the quantification of salt precipitation when multiple salts can form must be considered in conjunction with multi-component, multi-phase equilibrium flash problem, which is a combined chemical and phase equilibrium problem.

Organization of This Article
In this article, a new equation-solving approach to determining simultaneous chemical and phase equilibrium for mixtures of light gases and aqueous electrolyte solutions with the potential for multiple salt deposition is described. Concluding remarks are given in section 3.6.

Literature Survey
The open literature on multi-phase flash abounds and dates back many years.
Thus citation of all papers is not possible. Some of the more notable methodologies and articles for multiphase flash include the inside-out algorithm (Boston and Britt,1978), tangent plane analysis (Michelsen, 1982a,b), the negative flash , the Mixed Integer Nonlinear Programming (MINLP) approach of McDonald and Floudas (1995), interval Newton-bisection methods (Stadtherr et al.,1995), homotopy-continuation (Sun and Seider, 1995), and the successive quadratic programming approach . Notable methods and articles that address single and homogeneous multi-phase chemical reaction equilibrium include linear programming (White et al., 1958), an equation-tearing approach (Sanderson and Chien, 1973), nonlinear programming (Castillo and Grossmann, 1981), and the work of Smith and Missen (1982).
Elementary information regarding heterogeneous chemical equilibrium, equilibrium solubility products, and common ion effects can be found in many standard textbooks on chemical engineering thermodynamics. The reader is referred to the textbooks by  and Elliott and Lira (2012), which give basic theory for dilute solutions and/or common ion effects involving only a single common ion. In On the other hand, there are empirical correlations for estimating properties of electrolyte solutions (e.g., Lam et al., 2008). There are also more rigorous approaches for determining the solubility of ions in solution and their associated properties such as the use of (1) Debye-Huckel theory, (2) Pitzer expansions, (3) equation of state models (i.e., variants of the electrolyte Predictive Soave-Redlich-Kwong (ePSRK) by Li et al., 2001 and, (4) variants of the SAFT equation (Behzadi et al., 2005), and (5) various activity coefficient models as studied by Raatikainen and Laaksonen (2005) and Mohs and Gmehling (2012). Moreover, there is a large body of experimental data for highly soluble chloride, nitrate and some sulfate salts; salts like NaCl, KCl, Na2SO4, etc. Unfortunately, all of the computational papers cited deal strictly with the determination of solubility. More importantly, none of the cited papers explicitly address the problem of combined chemical and phase equilibrium where multiple fluid phases and multiple solid salts can co-exist nor do they address the associated numerical issues that can arise.

Problem Formulations
The problem of interest in this work can be conveniently divided into two parts -equilibrium flash and salt deposition.

Multi-phase Equilibrium Flash
The approach to the multi-component, multi-phase isothermal, isobaric equilibrium (TP) flash problem used in this work in the absence of chemical equilibrium has been described in detail in the previous work of Lucia et al. (2012aLucia et al. ( , 2012bLucia et al. ( , 2013 and is based on finding the global minimum of the dimensionless Gibbs free energy.

Salt Deposition
Salt deposition is a heterogeneous chemical equilibrium problem and its combination or coupling to multi-phase equilibrium flash has not been addressed to any great extent in the Process Systems Engineering (PSE) literature -particularly for cases where there are potentially a large number of salts that can form.
In general, salt deposition is approached by calculating and comparing equilibrium solubility products, Ksp j , to actual ion solubility products, Qsp j , for all possible salts, j = 1, ...,ns, where ns denotes the number of molecular salts. That is, if the conditions Qsp j > Ksp j (3.2) then salt j will precipitate out of solution. Moreover, from equilibrium solubility limits one can readily determine how much of each salt will precipitate. In these applications, equilibrium solubility products are generally computed from standard Gibbs free energies of formation, ∆Gf 0 , and adjusted for temperature using either a van't Hoff correction using heats of formation data, ∆Hf 0 and, if necessary, differences in heat capacities, ∆Cp. The resulting chemical equilibrium formulation leads to a system of nonlinear equations of the form where the component functions can be a mixture of polynomial, linear, and transcendental equations and the unknown variables are ion concentrations.
The wide variety of possible ions and salts in many applications leads to challenges in automating the formulation and solution of salt deposition problems.
However, there are two key equation-solving issues that need consideration (1) the basic formulation of the equations and variables and (2) an appropriate reduction strategy, when needed, that leads to a deterministic subset of equations that, when solved, gives a correct solution.

Basic Chemical Equilibrium Formulations
There are at least two different ways equations and unknown variables can be defined to determine equilibrium concentrations of ionic species in aqueous solution.
For example, the dissociation of an ideal solution of calcium and magnesium carbonate can be written in the form of 2 quadratic equilibrium solubility equations and 1 linear mass balance constraint (from charge neutrality) given by where

and 3.4 with ion activities.
The equations governing the dissociation of these carbonate salts, eqs. 3.4 and 3.5, can also be written in the form where now the unknown variables are zj = ln[cj], j = 1, 2, 3. This second formulation has 2 linear equations and 1 transcendental equation. Note that irrespective of formulation, the result is a set of nonlinear equations.

Equation Reduction
Since the salt deposition problem must consider the formation of all possible salts, the number unknowns, nC + nA, is always less than or equal to the number of possible salts, nCnA, when both nC and nA >

Chemical Equilibrium and Equilibrium Solubility Products
Equilibrium solubility products, Ksp, can be computed directly from Gibbs free energies of formation and related to the reaction equilibrium constant by assuming that the solid salt is a pure phase. For example In general, given the salt dissociation reaction CvCAvA → vCC + vAA (3.16) where C and A denote cation and anion respectively and vC and vA are the associated stoichiometric numbers, the equilibrium solubility product of an ideal solution at standard conditions is given by Values of ∆Gf 0 for a number of common ions and salts are listed in Table A1 in Appendix A.

Nonideal Solutions
When needed, mean ionic activity coefficients are used to account for non-ideal aqueous phase behavior. That is, where ± is the mean fugacity coefficient that represents non-ideal solution behavior for a composite fluid consisting of all ions dissolved in the solvent. The variable ϕ ± in Eq. (3.18) takes the place of γ ± in the usual representation of the product of the ion activities in electrolyte solutions. See Sandler (Eqs. 9.2-5 and 9.2-6, p. 668-669, 1999).
The actual value of ϕ ± for any ionic species, j, dissolved in water is determined by the simple ratio ϕ ± = (ϕj aq /ϕj w ), where ϕj aq and ϕj w represent the fugacity coefficient of j th ion in an aqueous solution at the given concentration and the fugacity coefficient of the j th ion in an aqueous phase at infinite dilution respectively and the same temperature and pressure. The value of lnϕ ± starts at an upper bound of 0 at infinite dilution and becomes more positive as the molality of the aqueous electrolyte solution increases.
In this work, we use the GHC equation to determine fugacity coefficients.
However, the electrolyte Predictive SRK (ePSRK) equation of , the Cubic Plus Association (CPA) equation of , or variants of the SAFT equation suitable for electrolyte solution (e.g., Behzadi et al., 2005) could alternatively be used to compute fugacity coefficients. Thus, the methodology presented here is not dependent on the EOS model.

Temperature Effects
When T is not too far from T0 = 298.15 K and the heat of formation is constant, temperature effects can be included using heat of formation data, ∆Hf 0 , and the van't Hoff equation A more rigorous approach requires the calculation of the heat of reaction as a function of temperature using heat capacity data and the integrals ∆Hf = ∫[∆Hf 0 + ∆Cp]dT (3.20) and lnKsp T = lnKsp + ∫[∆Hf/(RT 2 )]dT (3.21) where ∆Cp = vCCp(C) + vACp(A) -Cp(CvCAvA). Values of ∆Hf 0 used in this work are given in Table A1 in Appendix A.

Combined Chemical and Phase Equilibrium
When light gases are present and there is the possibility of ions in aqueous solution precipitating as molecular salts, the problem becomes a simultaneous chemical and phase equilibrium problem with interesting features. Of particular interest are the facts that 1) Light gases make it necessary to use an equation of state.
2) Ions and molecular salts are generally considered non-volatile so they do not exist in the vapor phase. b) If, on the other hand, precipitation of one or more molecular salts is indicated, then the amount of each salt that will precipitate must be calculated and the computations must move to step 3.
3) Re-solve the TP flash problem for a new flash solution. Because precipitation has occurred, the multi-phase TP flash results from step 1 will generally no longer satisfy equality of chemical potentials for the fluid phases. Therefore, a new feed must be calculated from the initial feed and the amounts of salt precipitates. The multi-phase TP flash is re-solved to high accuracy (|| .|| < 10 -10 ) to ensure closure of the solution.

Salts
In this section, a novel tearing algorithm for determining the equilibrium solubility limits is presented. Remember, equilibrium solubility limits and equilibrium solubility products are needed in steps 2 and 3 of the overall algorithm for systems that involve multiple salts. This algorithm is preferred over Newton's method because it can be proved to be globally convergent under very mild conditions as shown in Appendices B and C.

General Equation Formulation
The general formulation of equations used to define the equilibrium solubility limit problem is as follows. where Ci denotes i th cation, Aj denotes the j th anion, and vi,j and vi,nc+j are the stoichiometric coefficients for the dissociation of the k th salt, where k = (i-1)*nA+j. The total number of possible molecular salts is ns = nCnA.
3) Adjust all chemical reactions for non-ideal solution behavior and temperature effects.

Algorithmic Considerations
The steps in Section 3.4.1 lead to a deterministic set of (nC+nA) equations in (nC+nA) unknowns, written in the general form F(z) = 0.

A Novel Equation-Solving Algorithm
Here the steps of a novel tearing algorithm are presented.

Attributes of the Proposed Algorithm
The advantages of the proposed tearing algorithm are four-fold. Specifically, the proposed algorithm 1) Only requires the solution of linear subsets of equations in steps 7 and 10.
2) Is globally convergent for very large domains of attraction. See Appendices B and C.
3) Has a domain of attraction that can be quantified.
These facts are important because they show theoretically that it is straightforward to choose a starting point that is guaranteed to converge and converge quickly.

Convergence Results
Convergence results for the tearing algorithm depend on the salt reaction(s) and the initial values chosen for the anion concentrations. General results are not possible in the strictest sense. Appendices B presents convergence results for the proposed tearing algorithm for some simple cases. A number of illustrations supporting the convergence results in Appendix B are given in Appendix C.

Numerical Examples
In this section, two examples are presented that use the proposed numerical strategies given in sections 3.3.3 and 3.4.2 to solve multi-phase equilibrium flash problems in which the deposition of one or more salts is possible. The equation of state that is used in these computations is the multi-scale GHC equation of Lucia and coworkers. Table 3.1 gives a list of the conditions for the two example problems studied in this work. Example 1 is a simple model CO2 sequestration problem while Example 2 represents actual problems from a real enhanced oil recovery process in which production water is re-injected to cool an in situ burner and generate steam.

Example 1
This first example is a problem in understanding phase behavior in CO2  The next step of the overall algorithm is to determine equilibrium solubility limits, molalities, equilibrium solubility products (Ksp), and ion solubility products (Qsp) using the proposed tearing algorithm and to determine which, if any, molecular salts precipitate.    Table 3.2 is supersaturated. The final step of the overall algorithm is to resolve the flash problem following salt precipitation and combine the results of the re-solved flash with the salt precipitation results. Table 3.5 gives the final numerical solution for this multi-phase flash/salt precipitation problem. Note the following: 1) The dimensionless Gibbs free energy for the VLE + salts solution in Table 3.5 is -4.960854 and is lower than that for the supersaturated VLE flash solution shown in Table 3.2.
2) All ion molalities in the VLE + salt flash solution correspond to saturated values of molality.
3) For slightly soluble salts like CaCO3 and MgCO3, there is very little change in the amount of vapor because there is essentially no boiling point elevation due to the presence of these salts. However, this is not the case for highly soluble salts as illustrated in the next example.

Example 2
The second example in this article is taken from a real EOR process for a reservoir located in Saskatchewan, Canada. The EOR process is CO2 plus steam injection, where the steam is generated from production water and the objective of this example is to determine if salts will precipitate in the well bore or near wellbore region. Table 3.6 gives the feed condition for this second example, which is taken from an actual analysis of production water from the reservoir. Note that the molalities of the dominant ions, Na + and Cl -, are 1.09 and 1.16 m respectively and the electro-neutrality of the aqueous mixture is -7.36x10 -7 , which is very close to zero. Moreover, because of the relatively higher molalities of Na + and Clions in the production water, it is not unreasonable to anticipate the precipitation of NaCl. Table 3.7 gives results for the computation of the first VLE flash solution for this problem.  Table 3.8. From Table 3.8, one can easily see that the equilibrium solubility product computations predict that NaCl will precipitate at the given conditions of temperature and pressure. Again, this is because approximately half of the aqueous liquid has been vaporized, leaving the remaining aqueous liquid incapable of accommodating all of the ions. The actual amount of NaCl that precipitates is 7.2830x10 -3 moles/mole of feed. The equilibrium ion solubility limits shown in Table 3.8 converged to an accuracy of 5.575x10 -9 in 29 iterations. Table D1 in Appendix D gives the iterative values of all ions. The next step in the overall algorithm is to re-solve the VLE flash to determine the correct fluid phase equilibrium following precipitation. These results are shown in Table 3.9.  Finally, note that the VLE + salt solution has a lower dimensionless Gibbs free energy than the supersaturated VLE flash solution shown in Table 3.7, as it should, and that it is lower by approximately the dimensionless Gibbs free energy of formation for NaCl.

Conclusions
Combined chemical and phase equilibrium of mixtures involving light gases, aqueous electrolytes and the potential for solid salt precipitation were studied. A new tearing algorithm for determining equilibrium ion solubility limits was proposed and shown to generate a Cauchy sequence and therefore to be convergent under very mild conditions. A number of numerical and geometric illustrations of aqueous electrolyte solutions were presented that validate the robustness and convergence properties of the proposed tearing algorithm for finding equilibrium ion solubility limits. An overall multi-phase equilibrium flash algorithm that includes salt precipitation was also presented. Perhaps the most unique feature of this overall flash algorithm is that it decouples the chemical and phase equilibrium aspects of the flash problem. Two examples -one from CO2 sequestration and the other from EOR -were presented and used to illustrate the reliability and efficiency of the proposed multi-phase flash algorithm with salt precipitation.

Acknowledgement
This work was supported in part by a grant from the 'Petroleum Institute, Abu Dhabi' under Subaward Agreement No. 60371833-49954-B.

for Some Common Ions and Salts
The values of standard Gibbs free energies and enthalpies of formation  for the compounds studied in this article are given in Table A1.
To prove convergence, we use the contraction mapping theorem given in Ortega and Rheinboldt (p, 120, 1970) which states that if G(y): D ⊂ R n → R n is contractive on a closed set D0 ⊂ D, then G(y) converges to a unique fixed point in D0.

Remarks
Note that the results for this case are quite similar to those for the single cation case. In particular, 1. Note that all large values of y1,0 and y2,0 are guaranteed converge to y*.
2. The results given by Eq. B20 and B21 hold for ideal and non-ideal aqueous electrolyte solutions alike since no restrictions have been placed on any Ksp.
3. The results given by Eq. B20 1nd B21 also hold for any relevant temperature.

Appendix C: Illustrations of Convergence
In this appendix, a number of illustrations are presented to support the proof of convergence in Appendix B.
Illustration 1 Table C1 gives all iterative values of yk, G(yk), and xk for NaCl at 298.15 K and 1 atm and clearly corroborates the convergence results given in Appendices B.   It also shows that the condition G'(y*) = 0 is satisfied, which means asymptotic convergence is very fast.

Introduction
Gas hydrates are non-stoichiometric crystalline compounds composed of a network of hydrogen bonded water molecules forming cavities. Small guest molecules can occupy the cages, and in doing so stabilize the structure. The three most common naturally forming hydrate structures have been classified as structure I (sI), structure II (sII), and structure H (sH). This work is only concerned with sI hydrates. The cubic sI hydrate unit cell consists of two small (5 12 ) and six large (5 12 6 2 ) cavities. The notation e f denotes a polygon with f faces consisting of e edges, so that the sI small cages are polyhedra made of twelve pentagons, while the large sI cages are made of twelve pentagons and two hexagons. Since their discovery by Sir Humphrey Davies in 1810, gas hydrates have been of interest to the scientific community. They were later rediscovered by Hammerschmidt (1934) as the root cause of blockages in gas lines.
Estimates made over the last decade concerning the amount of methane currently stored in gas hydrates vary widely (3 -120 x 10 15 m 3 STP)(Sloan and Koh, 2007,

Procedure
Constant pressure Gibbs ensemble Monte Carlo (GEMC) (Panagiotopoulos, 1987) simulations were performed to measure the equilibrium occupancy of pure structure 1 (sI) gas hydrates. The constant pressure Gibbs ensemble was preferred over the grand canonical ensemble (GCMC) because in this approach the pressure rather than the chemical potential is specified as an input parameter. In the grand canonical ensemble, the pressure must be calculated from the specified chemical potential either (1) using an EOS model, (2) determined during the simulation by calculating the molecular virial (Frenkel and Smit, 2001) or (3) by applying a method similar to Widom test particle insertion (Widom, 1963) in which test volume moves are attempted, but no volume move is ever actually applied to the system (Eppenga and Frenkel, 1984).
Method 1 requires a choice of equation of state for the system of interest, and both of the later methods involve ensemble averages, and thus are subject to statistical uncertainty. Simulations in the Gibbs ensemble have been previously applied to study adsorption isotherms (eg. Bai et al., 2014;Demir and Ahunbay, 2014) but as far as we know this is the first application of the method to a gas hydrate system.
In this work, the simulations were initialized using two boxes, one containing an empty sI hydrate lattice and the other containing guest molecules. The box containing guest molecules was initialized to a simple cubic structure. Monte Carlo (MC) moves available to both boxes included translation and rotation of the molecule about the center of mass. Isotropic volume moves were available to the box containing only guest molecules, while anisotropic unit cell displacement moves were applied to the box containing the hydrate lattice. In addition, particle swap moves were available only to the guest molecules in the simulation, so that no water molecules were removed from the hydrate phase or inserted into the gas phase, but guest molecules could be removed from the gas phase and inserted into the hydrate phase and vice versa. In the rigid hydrate model, the rotation and translation moves were only applied to guest molecules, and no volume displacement was attempted on the box containing the hydrate phase.
Volume displacement, unit cell displacement, and particle swap moves were all attempted with a frequency of 0.01, while center of mass translation and rotation were attempted with frequencies 0.57 and 0.4 respectively. The simulations were first run in 'equilibration' mode where the maximum displacements (volume, unit cell, translation and rotation) were updated every 10 MC cycles to achieve an acceptance ratio of 0.5 for each type of move. Properties were sampled in a 'production' mode in which the maximum displacements were updated four times during the simulation. A classical Lennard-Jones (LJ) potential was applied using Lorentz-Berthelot mixing rules for unlike sites. Long range electrostatic interactions were incorporated using Coulomb's law with the Ewald summation method. A LJ cutoff of 10 A was applied and the standard LJ tail correction was applied for distances greater than the LJ cutoff. For the Ewald summation five inverse space vectors were used in each direction and the electrostatic cutoff was allowed to adjust to half of the current box length. The simulation box containing the structure I hydrate lattice studied in this work was comprised of 8 unit cells (2x2x2) (consisting of 368 water molecules). The TIP4P/Ew potential (Horn et al., 2004)

Computations
All computations were performed using the open source Monte Carlo software MCCCS Towhee version 7.0.6 (Martin, 2013). Towhee input files for any of the simulations in this work will be made available by the authors on request. In addition, all relevant simulation details that would be required to reproduce the work are given in section 4.2.

Validation of Method using Phase Equilibrium Computations
Typically in the vdW-P model, hydrate occupancy is determined at a given temperature and pressure by the expression (assuming a simple hydrate system) = 1+ (4.1) where the subscript refers to cavity type, is guest fugacity, and is the cavity Langmuir constant. Normally, a distinction is made between the small and large cavities in structure I hydrates when determining the Langmuir constant, but Lasich et al. (2014) have shown that it is possible to match experimental methane hydrate data without making such a distinction between the two cavity types. This is because methane exhibits one site adsorption behavior in the hydrate phase, filling both cavity types without preference (Glavatskiy et al., 2012). In general, Langmuir constants can be calculated from Kihara potential parameters, but commonly (for example: Parrish and Prausnitz, 1972;Munck et al., 1988) their temperature dependence is described by the following relationship as a function of the parameters A and B that are regressed to various data sources From the work of van der Waals and Platteeuw (1958) and Parrish and Prausnitz (1972) the criterion for gas hydrate equilibrium is defined as where refers to either liquid water or ice, depending on which phase is present at the current conditions. The hypothetical empty hydrate phase is chosen as a reference state, the difference in the chemical potential of water in the filled and empty hydrate at the given temperature and pressure can be expressed as where is the ratio of cavities of type to water molecules in a unit cell. The subscript runs over all cavity types and runs over components. The difference in chemical potential of water in the − ℎ (eg. liquid or ice) and the empty hydrate reference state can be expressed as where Δ and Δ are the enthalpy difference and volume difference between thephase and the empty hydrate reference phase. The volume difference is assumed to be constant over the temperature and pressure range of interest. The temperature dependence of the enthalpy difference is incorporated using the difference in isobaric heat capacity as All of the required reference parameters can be found in the literature (Parrish and Prausnitz, 1972;Munck et al. 1988).
In the normal procedure that is followed using models based on vdW-P theory, the parameters A and B in equation (4.2) for the calculation of Langmuir constants, or Kihara potential parameters would be regressed to the experimental gas hydrate phase equilibrium data. In this work, following the procedure developed by Lasich et al.
(2014), the parameters A and B were fit to gas hydrate occupancy data obtained using GEMC simulations. Then these parameters were used with equations (4.1-6) to calculate hydrate dissociation conditions. The equations were solved iteratively using Newton's method, and converged to a tolerance of 10 −12 .

Simple structure I methane hydrate occupancy results
In this section, results of simulations for pure structure I methane hydrate with both a rigid lattice and flexible lattice are described in detail.  The flexible hydrate model is a more rigorous description of the hydrate phase because in reality the hydrate lattice is distorted by the presence of guest molecules and changes in temperature and pressure (Klauda & Sandler, 2000;Tse, 1987). A comparison of occupancy predicted by both models to experimental data (Uchida et al., 1999) is given in  Overall hydration number of the hydrate structure is simple to calculate from gas hydrate total occupancy information. A comparison between calculated methane gas hydrate hydration number and GEMC simulation results for methane hydrate is shown in Table 4.3. The data in Table 4.3 is taken directly from Anderson (2004) (see Table 6 on page 1124), and was calculated using a wide range of experimental hydrate dissociation data and the Clapeyron equation. The resulting AAD is 6.97% and 2.44% for the flexible and rigid lattice models respectively. The rigid hydrate lattice model again seems to predict values that are in much better agreement both with experimentally measured values (Uchida et al., 1999), and the results of Anderson (2004). A comparison of the simulation results for both lattice models and predictions made by vdw-P theory for methane hydrate occupancy is shown in

Methane Hydrate Phase Equilibrium Predictions
The procedure described in Section 4.2.2 was used to calculate simple sI methane gas hydrate dissociation conditions. The governing equations were solved iteratively using Newton's method and converged to a tolerance of 10 −12 . The results for both the rigid and flexible model are compared to experimental data in  Table 4, page 53) present results for AAD in predicted dissociation temperature using similar potential models and grand canonical simulations of 0.8-3.7%. The smaller deviation from experimental dissociation temperatures in this work may be a direct result specifying pressure explicitly. The heat of dissociation of the hydrate phase can be estimated using the Clausius-Clapeyron equation Assuming that Z (the gas compressibility) is 1 over the range of interest Δ can be easily obtained by plotting versus 1 (see figure 4.4).

Figure 4.4: Ln(p) versus (1/T) and a Linear Fit used with Equation (4.7) to Determine the Heat of Dissociation
The resulting value for Δ is 88.4 kJ/mol, in error of about 63% from the experimentally measured values of 54.19 ± 0.28 kJ/mol reported by Handa and coworkers (1986) and 53.5 ± 1.3 kJ/mol reported by Anderson (2004). One possible cause of the strong deviation is the assumption that Z = 1 in equation (4.7), which is not valid at high pressures and it essentially removes any temperature dependence of the heat of dissociation. Anderson (2004) discusses the weaknesses associated with making such an assumption for methane hydrate systems, and performs calculations in which this assumption is removed. The assumption is retained in this work for simplicity, and to facilitate comparisons between this method and the work of Lasich and coworkers

Conclusions
The following conclusions were drawn from the results presented and discussed in section 4.3: • Simulation of sI gas hydrate occupancy using the constant pressure Gibbs ensemble is an attractive alternative to the grand canonical ensemble because pressure, instead of fugacity or chemical potential, is given as an input parameter. This removes the need to calculate the pressure using an appropriate equation of state model, calculation of the molecular virial, or test volume moves. The main disadvantage of the constant pressure Gibbs ensemble compared to the grand canonical ensemble is the simulations are more computationally expensive because of the need to calculate interactions in both the box containing the hydrate lattice, and the fluid phase. However, this allows the calculation of properties for the fluid phase in equilibrium with the hydrate phase, and will be much more interesting in an extension of the work in which a mixture of guests is introduced. Another common disadvantage occurs when dealing with dense phases (high pressure and/or low temperature) or bigger molecules; in this case the acceptance rate of particle swap moves can be rather low. Biasing techniques have been developed to deal with this situation (Snurr et al., 1993). However, in this work they were not needed as acceptance rates of the particle swap moves, though low, were manageable (0.3-1.1%). It is possible that in order to efficiently sample with a more complicated guest molecule in the simulation biasing techniques will be required, and if so they will be implemented in future work. showed that in a GCMC simulation framework the number of unit cells included in the simulation does not affect the results. However, the simulations from this work were rerun using a single unit cell and the resulting occupancy was systematically higher (by as much as 5% at higher temperatures) than the occupancy predicted using eight cells. Therefore, the increased accuracy of the results is likely due to a combination of both including pressure explicitly and using a larger number of unit cells. Using a single unit cell also presents an issue because the lattice parameter, and thus the box size for the simulation box containing the hydrate phase, is 12.03 Angstroms. This limits the possible length for the Lennard-Jones cutoff to ≤ 6.015 Angstroms because it must be less than half of the box size to be consistent with the minimum image convention (Allen and Tildesley, 1991). This leads to fewer important interactions being considered in the single cell simulation. For example, the average distance between guest molecules in a structure I hydrate is ~6 Angstroms (Sloan and Koh, 2007), which may explain why the occupancy results from single cell simulations are systematically higher. That is, the interactions between guest molecules with guests in neighboring cavities and beyond, as well as with water molecules past the first coordination shell, are not being properly accounted for.
Another attractive aspect of the constant pressure Gibbs ensemble is that it is easy to imagine including additional simulation boxes. For example, a box could be added containing a liquid-like phase in a mixed hydrate system to study hydrate occupancy with a vapor and liquid phase in thermodynamic equilibrium. Future work on this project will involve the study of simple sI hydrates of CO2, mixed structure I hydrates of CO2 + CH4, the addition of another simulation box to simulate mixed hydratevaporliquid equilibrium, and finally structure II hydrates.

Abstract
Fully compositional and thermal reservoir simulation capabilities are important in oil exploration and production. There are significant resources in existing wells and in heavy oil, oil sands, and deep-water reservoirs. This article has two main goals: (1) to clearly identify chemical engineering sub-problems within reservoir simulation that the PSE community can potentially make contributions to and (2)  Numerical results for several chemical engineering sub-problems and reservoir simulations for two EOR applications are presented. Reservoir simulation results clearly shows that the Solvent Thermal Resources Innovation Process (STRIP) outperforms conventional steam injection using two important metrics -sweep efficiency and oil recovery.

Introduction
Modeling and simulation to predict long-term performance of oil recovery methods (i.e., reservoir simulation) is a topic dating back to the 1950s, late 1960s, and the early 1970s (see Douglas Jr. et al., 1959;Price and Coats, 1974;Todd et al., 1972).
Early reservoir models (e.g., black-oil reservoir models) were typically based upon rigorous mass balance equations for key species (oil, water, and gas) but only used approximate phase equilibrium (e.g., no oil dissolved in the water phase) and/or There remains considerable oil in place (OIP) in many reservoirs that are either in current operation or have been shutdown (often with infrastructure remaining in place). There are also large amounts of fossil fuels in heavy oil, oil sands, and deep sea reservoirs -but these hydrocarbons are more challenging and more costly to produce.
An increase in production of just 1% represents a 6. Models consisting of differential-algebraic equations (DAEs).
7. Nonlinear equation-solving using Newton and trust region methods.

Iterative linear equations solving.
Thus it is our opinion that chemical engineers, particularly those in the process systems engineering (PSE) community, are in a unique position to make significant contributions to various aspects of reservoir simulation.
In this paper, we present an advanced reservoir modeling and simulation framework for fully compositional and thermal reservoir simulation and subsequently apply this simulation framework to a comparative study of steam injection and STRIP in EOR applications. We also identify those sub-problems that are, in our opinion, problems that the PSE community can contribute to. Accordingly, this paper is

Literature Survey
The main focus of this article is numerical reservoir simulation, which comprises a vast body of literature and thus it is not possible to survey all relevant scientific papers.
Therefore in this section only a summary of those papers and numerical methods directly relevant to the modeling and simultaneous solution of numerical reservoir models is presented. We refer the reader to the book by Peaceman (2000) for an introduction to the fundamentals of reservoir modeling and simulation and a description of some of the earlier underlying numerical methods that have been used. A secondary focus of this manuscript is to identify sub-problems within a larger reservoir simulation that are clearly within the skill set of process systems engineers.
Some of the earliest work in numerical reservoir simulation dates back to 1959 and the pioneering work of Douglas Jr. et al. (1959), who developed numerical methods for the simultaneous solution of time dependent two-phase flow problems in one and two spatial directions. Governing partial differential equations (PDEs) describing conservation of mass and flow were converted to nonlinear algebraic equations using difference approximations and the resulting nonlinear algebraic model equations were then solved using various numerical methods including alternating-directions implicit, Jacobi iteration, successive over-relaxation, Gauss-Seidel iteration, and other established techniques. We refer the reader to the book by  for a comprehensive description of these numerical methods. The work of Douglas Jr. et al. (1959) was later extended to three spatial dimensions by Coats et al. (1967) and three-phase flow problems by Peery and Herron (1969) and Sheffield (1969). Other journal articles that address additional physics in reservoir simulations and solve model equations simultaneously include those by Snyder (1969), Settari and Aziz (1974), and Trimble and McDonald (1981). Key differences among many of the early approaches to reservoir simulations reside largely in model formulation and the methods used to solve the resulting algebraic model equations. These differences persist today. Note that all of the topics that have just been described are topics that are very familiar to the PSE community -dynamical equations describing conservation of mass and energy, differencing, and nonlinear equation solving.
Current state-of-the-art reservoir simulation has moved to two basic nonlinear formulations -a natural formulation (Coats, 1980) and a molar formulation (Acs, 1985).
Large sets or subsets of nonlinear algebraic equations result from discrete representations of the governing PDEs that describe the spatial and temporal evolution of the system. The most commonly used approaches for discrete representation are finite-difference or finite-volume approximations on structured or unstructured grids.  (Cao, 2002). A good survey of the numerical characteristics of FIM, IMPES, and AIM is given by Marcondes et al. (2009).
Regardless of the formulation, many current solution methods use some form of iterative linear equation solving (e.g., GMRES or other Krylov subspace methods) with pre-conditioning to solve the linear system of equations that determines the Newton correction to the variables at each time step.

Reservoir Model Equations
The equations describing the time evolution of fluid compositions, temperatures, and pressures in a reservoir comprise a set of coupled, nonlinear PDEs that describe conservation of mass, energy, and momentum. In addition, various thermo-physical properties, equilibrium (or non-equilibrium) behavior of fluid phases, properties of porous media, and well configuration specifications are included as algebraic constraints to the governing PDEs. In this article, the governing PDEs are represented in discrete form using finite-volume discretization and when used with additional constraints that form a large set of nonlinear algebraic equations. In this section, the reservoir equations as well as other constitutive equations are described.

Reservoir Model Equations in General Form
The nonlinear time-dependent PDEs that represent conservation of mass and energy in a reservoir are given by where ϕ is the porosity of the porous media, denotes molar density, is composition in mole fraction, is saturation, is volumetric flow, is molar diffusion flux, which is usually ignored for large scale applications, and is a mass source or sink term. In eq. (5.2), denotes internal energy, is enthalpy, is heat conduction flux, is an energy source or sink term, and the summation is also over all phases k = 1, …, P. The subscript i denotes a given component while the superscript k denotes a given phase. C is the total number of components in the mixture and P is the total number of phases.
The subscript M in eq. (5.2) denotes the porous media whereas the symbol ∇ denotes the gradient of a vector.

Phase Equilibrium in General Form
Phase equilibrium in any given finite volume (also referred to as a cell or grid block) is described by the equality of partial fugacities for all components in all phases, clearly a topic on which publications in the chemical engineering literature abound. In particular, 1 = 2 = ⋯ = , = 1, … , ; = 1, … , (5.3) where denotes the partial fugacity of component i in phase k and is given by where denotes the mole fraction, denotes the fugacity coefficient of component i in phase k, and p is pressure.

Conservation of mass within any grid block is represented by a set of component mass balance equations
− ∑ = 0, = 1, … , (5.5) where and are the total density and mole fraction of component i in the cell. Note that there is some overlap in symbols because standard notation in reservoir engineering and chemical engineering thermodynamics use the same symbol to denote different quantities. We caution the reader to pay careful attention to context so the meaning of a symbol is clear.  have a lower computational overhead and provide results that are within acceptable accuracy. As described later in this article, our implementation provides the user with a number of more commonly used cubic EOS.

Other Constitutive Equations
Other constitutive equations are also needed to close the numerical model and allow proper integration of eqs. (5.1) and (5.2). These constitutive equations include Darcy's Law, heat conduction, and when relevant, diffusion equations -again all topics familiar to the PSE community in traditional applications such as heat and mass transfer in catalyst pellets, non-equilibrium models in multi-stage distillation as well as more recent applications in bio-medical modeling and simulation of the brain, and others.

Darcy's Law
Darcy's law describes the volumetric flow of each phase, , through porous media and is given by where is an intrinsic rock or soil permeability, is relative permeability, is viscosity, is the acceleration due to gravity, and is the coordinate in the direction of gravity.

Heat Conduction Equations
The heat conduction equation is where K is the conduction coefficient for phase k and is absolute temperature.

Equation Coupling
The conservation of mass and energy, flow, and conduction through porous media described by eqs. (5.1)-(5.7), and the equations describing the conservation of mass, conservation of energy with heat losses to the surroundings, phase equilibrium, and holdup in injection and production wells form a large system of nonlinear algebraic equations that are strongly coupled. In a hierarchical sense, the EOS lies at the innermost level of the computations and provides the phase densities. Phase densities are used to calculate fugacity coefficients, and fugacities to determine the type and amounts of each phase present in a grid block (i.e., by solving the traditional chemical engineering Tp flash). Calculated phase densities and composition from the flash are then used to determine the unknown variables at the reservoir level (e.g., pressures, saturations, and temperatures) as well as heat conduction fluxes, and the flow of phases through the porous media.

Implementation
As

AD-GPRS
AD-GPRS is an advanced reservoir simulator with wide ranging capabilities that There are, of course, many details associated with AD-GPRS ; here we only summarize its main features.

Formulations
Both natural and molar formulations are available in AD-GPRS . Regardless of formulation, the primary dynamical model equations 2) summation equations for the mole fractions in each phase.
3) saturation summation equations Similarly heat (energy) flux can be expressed as where , 1 is the geometrical part of the transmissibility coefficient assuming Two-Point Flux Approximation (TPFA) for the conduction term .

Solution of Nonlinear Algebraic Equations
The Newton-Raphson method is used where a Jacobian matrix is assembled and a linear system of equations solved for each iteration. The relations, other than the mass conservation equations, are treated as constraints that are local to a grid block (cell). To minimize the size of the global linear system, a Schur-complement procedure is applied to the full Jacobian matrix of each block to express the primary (mass conservation) equations as a function of the primary variables only (Cao, 2002). After the size of the system is reduced, the resulting global linear system of equations is solved for the primary variables. An iterative linear equation solver with pre-conditioning is used to solve the linear system. After the linear system is solved, the computed changes in the primary variables are used with the secondary equations to obtain changes in the secondary variables locally in each grid block. Next, the nonlinear variables can be updated using different strategies and safeguards to ensure that the solution remains within physical boundaries.
Convergence of Newton-Raphson iteration depends on several aspects that include (1) any corrections to updated variables that employ safeguards, (2) various chopping strategies for different unknowns, and (3)

Phase Behavior Computations
In this section, we describe different approaches to phase behavior computations in AD-GPRS, which include the use of intermittent flash solutions and CSAT.

Intermittent Flash Problem Solutions
For phase behavior computations, AD-GPRS uses a two-stage procedure. In the first stage, the number of phases that exist in each cell is determined. This can be obtained using Gibbs energy minimization or phase stability analysis . In the second stage, flash calculations are performed to determine the compositions of the existing phases . At both stages a combination of Successive Substitution Iteration and Newton's method is used.
As an alternative to this two-stage strategy, a generalization of the negative-flash based approach of  can be used .
Here it is assumed that the number of phases present is the maximum possible, and then eqs. (5.3) and (5.5) are solved, allowing for phase fraction to be less than zero, or greater than one. When the phase fractions of a converged negative flash procedure are negative, fewer existing phases are assumed and a similar procedure for this reduced system may be required .

Compositional Space Adaptive Tabulation (CSAT)
Solving flash problems for all grid blocks over all nonlinear iterations and time steps is computational demanding. To improve the performance of phase behavior computation in reservoir simulation, the CSAT approach originally developed by Voskov and Tchelepi (2009a,b) is used. CSAT adaptively stores a discrete set of tielines at different pressures and temperatures to represent phase behavior during reservoir simulation. This collection of tie-lines is interpolated and used to look up the phase state of the mixture at a particular pressure and temperature. In addition, the number of tie lines can be collected adaptively based on the specific attributes of a compositional solution during a reservoir simulation.
CSAT completely replaces the need for phase stability tests and provides good initial guesses for the standard Tp flash computations.

Compositional Space Parameterization (CSP)
The compositional space parameterization (CSP) method (Voskov and Tchelepi, 2009a;Zaydullin et al., 2013) is based on casting the nonlinear governing equations (1) and (2), including thermodynamic phase equilibrium constraints (3), in terms of the tiesimplex ( ) space. During a simulation, the space is adaptively discretized using supporting tie-lines. The coefficients for the governing system of equations, including the phase compositions, densities, and mobilities, are computed using multi-linear interpolation in the discretized space.
Using the CSP methodology, phase behavior computations can be replaced by an iteration-free look-up table procedure during the course of a reservoir simulation, removing the need for standard EOS computations (phase stability and flash). Also, it is important to note that the error associated with multi-linear interpolation is bounded and decreases with grid (or space) refinement  and therefore only a limited number of supporting tie-lines are needed for the accurate representation of phase behavior. That, in turn, leads to significant gains in computational efficiency.
Tie-lines or tie-simplexes needed for CSAT and CSP can be parameterized using the generalized negative flash procedure  or using GFLASH (Section 5.4.2).

GFLASH
GFLASH is a FORTRAN suite that models and solves the traditional chemical engineering multi-phase, multi-component isothermal, isobaric (Tp) flash problem. That is, given an overall composition for a fluid mixture, , a temperature, , and pressure, , GFLASH determines the number of phases that exist at equilibrium and their corresponding compositions, fugacities, densities, and enthalpies. In this section, formulations, overall solution strategies, and methods of solution are described.

Equations of State
A number of the commonly used cubic EOS with and without volume translation are implemented in GFLASH. The EOS available include the 1) Soave-Redlich-Kwong (SRK) equation ,

Formulation and Solution
All EOS are formulated as cubic polynomials in compressibility factor, , in the complex plane in the form ( ) = 1 3 + 2 2 + 3 + 4 = 0 (5.14) The resulting single variable function, ( ), is solved using Newton's method in the complex plane to find any root to an accuracy of | ( )| ≤ 10 −12 . The cubic polynomial is then deflated to a quadratic equation, which is solved using the quadratic formula to determine the other two roots. This approach removes the need to use an accurate initial guess for Newton's method, guarantees that all three roots will always be found, and is actually faster than using the analytical solution to a cubic polynomial.

Root Assignment
Correctly determining which root is liquid and which root is vapor is as important, if not more important, than computing roots to EOS and is particularly challenging under harsh conditions (i.e., high and high ).The current approach used to assign roots in GFLASH is as follows: For a set of roots given by { 1 , 2 , 3 }, where any root has the complex variable form = ± , = 1, 2, 3, we define where | | denotes the complex absolute value function given by | | = √ 2 + 2 and the superscripts L and V denote liquid and vapor respectively. Phase densities are easily computed from the compressibility factors, and , using the expression = (5.16)

Flash Problem Formulations and Method of Solution
The flash problem is really two problems -a phase stability problem and a phase equilibrium problem. In GFLASH, the formulations of the phase stability and phase equilibrium conditions use the dimensionless Gibbs free energy of mixing, ∆ / , and the dimensionless Gibbs free energy, / , respectively.

Phase Stability
Minima in ∆ / often turn out to be inexpensive and good approximations for points of tangency. The necessary conditions for a minimum in ∆ / are formulated in terms of the equality of dimensionless chemical potentials, , = 1, … , . For the phase split (or phase stability) problem, which is always a two-phase determination, the model equations are given by and then projected onto the mass balance constraints in eq. (5.18) to reduce the size of the phase equilibrium problem and to ensure that conservation of mass is satisfied at each iteration.
Other Modeling Capabilities in GFLASH GFLASH also has the capability of solving chemical reaction equilibrium problems and combined chemical and phase equilibrium problems, topics that are both familiar to the PSE community and important in various applications of EOR. For example, in applications of STRIP, partial oxidation of methane is used to generate in situ CO2 and steam. In other EOR applications, where production water is re-used in order to defray the high cost of purchasing municipal water or expensive water treatment, salt precipitation in the presence of multiple fluid phases, which is a combined chemical and phase equilibrium problem, can be a serious operational problem. These sub-problems are also solved by Gibbs free energy minimization within the GFLASH framework.

Method of Solution
GFLASH uses a trust region method to solve both the phase stability and phase equilibrium model equations. This methodology is a simple version of the terrain methodology developed by Lucia and Feng (2003). When applied to phase stability and phase equilibrium, we restrict the terrain method to look for only one stationary point in each of ∆ / and / respectively.
In addition, when solving flash problems, GFLASH alternates between phase stability and phase equilibrium sub-problems, maintaining a monotonically decreasing sequence of values of / until a global minimum identifying the number and type of phases as well as their associated mole numbers is found. Phase stability problems [i.e, eq. (5.17)] are solved to an accuracy of || ( )|| 2 ≤ 10 −6 , where || . ||2 denotes the 2norm or Euclidean norm. In contrast, phase equilibrium problems (eq. 5.20) are solved to an accuracy of || ( )|| 2 ≤ 10 −4 for two-phase equilibria and 10 −5 for three-phase equilibria.

The Connection Between CSAT and GFLASH
In this section, we describe the connection between rigorous phase stability and flash computations using GFLASH and CSAT. We also describe the interface between AD-GPRS and GFLASH.

Conventional phase behavior computations
The number and types of phases (or phase state) of a mixture in a given grid block can vary. For example, for mixtures that exhibit three-phase behavior, there are seven different possible phase states -three single phase states (i.e., water-rich liquid, vapor, or oil-rich liquid), three different two-phase states (i.e., LLE, water-rich VLE, or oilrich VLE), or vapor-liquid-liquid equilibrium (VLLE). Thus, the phase state as well as all corresponding phase compositions need to be determined for every grid block on each Newton iteration. For the natural formulation, a three-step procedure is usually used for these computations: 1. For any grid block, the current phase state is determined using a phase stability test.
2. If the current phase state is different from one on a previous Newton iteration, flash computations are performed to obtain phase compositions.
Because of the complexity of the ADGPRS-GFLASH interface, both a phase stability test and flash computations are performed simultaneously in the current model framework realization.

Phase behavior computations with CSAT
As noted, CSAT can significantly improve the time required for phase behavior computations in fully compositional reservoir simulation (Voskov and Tchelepi, 2009a;Voskov and Tchelepi, 2009b). The general multiphase implementation of CSAT Voskov and Tchelepi, 2009b) is a two-step procedure: 1. Computation of supporting tie-simplexes (i.e., tie-triangles for three-phase systems).
In the original CSAT implementation of , a generalization of the negative-flash idea  for Step 1 and geometrical parameterization (i.e.,. tracking tie-lines from each side of a tie-triangle) for Step 2 was used. While this approach proves to be robust for challenging three-phase systems, it requires some preliminary knowledge of a multiphase mixture under investigation because the geometry of tie-simplex subspace can be quite complicated.
In this work, we have used a different strategy. First, we use GFLASH to provide fugacities for given pressure, temperature, and phase compositions and the generalized negative flash approach  is used to find a supporting tie-simplex for the CSAT procedure. Next, an extension of the tie-simplex is adaptively discretized and GFLASH determines the phase state of a model cell. Finally, the collection of tiesimplexes and their extensions are interpolated for a particular pressure and temperature and used to look up the phase state of the mixture.

AD-GPRS/GFLASH Interface
Because AD-GPRS is written in C++ and GFLASH is a FORTRAN suite, the proposed modeling and simulation framework is necessarily mixed language and therefore an interface is needed to communicate information between the two programs.
The interface program is described in Appendix A.

Thermal EOR Methodologies
In this section, steam injection and STRIP are introduced along with common performance metrics used to evaluate thermal EOR methods. We refer the reader to the work of Aziz et al. (1987), which is the 4 th SPE Comparative Solution Project: Comparison of Steam Injection Simulators, for an introduction to steam injection and associated simulation challenges.

Steam Injection
Steam injection is generally implemented using surface facilities to generate superheated steam, which is injected into a reservoir through a well. The entering steam heats the formation and lowers oil viscosity, which allows the oil to flow more easily to production wells. In all steam injection methods, surface generation of steam suffers from a number of disadvantages, not the least of which is energy losses (up to 50%) to the walls of the injection well.

Solvent Thermal Resource Innovation Process (STRIP)
STRIP, which has been developed by RII North America, is an environmentally friendly approach to EOR, which is typically deployed into existing wells, so there is little or no disruption of land. Unlike other steam injection processes, STRIP generates steam and CO2 by in situ combustion of methane in oxygen, which eliminates energy losses to the injection well and delivers steam directly to the formation. STRIP also provides a co-solvent, CO2, which is well known to enhance oil recovery by swelling oil and lowering viscosity. The STRIP burner can be placed in a number of configurations, but in this work the STRIP burner resides in a vertical section of the injection well. Because the combustion temperature can approach 3,000°C, the STRIP burner is typically cooled using production water, significantly reducing and often removing the need for municipal water. The nominal composition of lumped gases entering a reservoir formation is around 10 mol%, with roughly 6.7 mol% being CO2.

Performance Metrics
Several common metrics are used to evaluate the performance of a thermal EOR methodology. These metrics include (1) sweep and (2) oil recovery, which, of course, is of primary interest.

Numerical Examples
In this section, two reservoir simulation examples are presented to elucidate key points, to compare the performance of steam injection and STRIP, and to demonstrate the reliability and computational efficiency of the numerical tools in GFLASH and AD-GPRS. However, prior to presenting results for reservoir simulation with STRIP, a number of traditional chemical engineering sub-problems needed to be solved, including a chemical equilibrium problem, an adiabatic flame temperature problem, and a salt precipitation problem, to clarify and quantify various aspects of the reservoir simulations.

Example 1: Chemical Equilibrium of STRIP Combustion
As noted, STRIP generates in situ CO2 and steam by partial oxidation of methane. In a typical application, the reactants are fuel-rich and thus there are a number of 'major' syngas by-products such as H2 and CO, and un-reacted methane and O2 in addition to the CO2 and steam. However, the composition of the combustion product stream is a function of both the O2/CH4 ratio and the reaction temperature, the latter of which we do not know. The O2/CH4 ratio, which is denoted by r, is an operational decision based on extensive laboratory experimentation and can vary between at 1.6 and 1.9 depending on the application. Note that the stoichiometric ratio of oxygen/methane for complete combustion is 2 and thus STRIP combustion is fuel-rich and thus will produce some syngas (H2 and CO).

Example 2: Adiabatic Flame Temperature
While Fig. 5.1 gives chemical equilibrium results for a wide range of combustion temperatures, the actual temperature of the combustion products for a given set of conditions should be estimated using an adiabatic flame temperature calculation and is coupled to the previous chemical equilibrium problem because the temperature and composition are interdependent. However, because product gas compositions are relatively weak functions of temperature above about 1000 ⁰C, we can effectively decouple the problems and solve the adiabatic flame temperature problem using the product gas compositions shown in Table 1. The corresponding problem formulation, which is rather simple and can be found in several undergraduate level thermodynamics textbooks, is shown in eq. (5.26).
where ∆ 0 is the standard heat of reaction, 0 is a reference temperature and equal to 25 ⁰C, and , is the heat capacity for the ℎ component. Data for ∆ 0 and , are given in Appendix B. The calculated adiabatic flame temperature for an oxygen-tomethane ratio of 1.6 is 3343.975 ⁰C.

Example 3: Salt Precipitation
Use of production water in EOR processes can reduce the cost of purchasing municipal water or operating an on-site water treatment plant and concomitantly lower environmental impact. The main challenge associated with the use of production water is the presence of ions and the potential for salt precipitation. To lower the potential for precipitation, production water can be mixed with clean water. In the case of STRIP, production water is mixed with in situ generated steam for two purposes -to generate additional steam and to cool the STRIP combustion burner. Table 5.2 gives an illustration of the compositions of production water and the combined feed for a STRIP application to a real reservoir in Saskatchewan, Canada before and after mixing. Note that production water analyses are generally reported as molality, whereas mass or mole fractions are generally used in flash calculations. The real concern regarding precipitation comes from the fact that the combustion products from STRIP are very hot and thus the amount of liquid available to dissolve ions, even after mixing, might be quite small if too much vaporization occurs.
Remember, the main purpose of STRIP is to inject enough steam and CO2 into the reservoir for improved oil recovery. What this means is that the desired fluid stream entering the reservoir (i.e., at the injection well bore) should have a relatively high vapor fraction, say between 0.7 and 0.8. To determine whether or not salt precipitation will occur, we must therefore solve a combined chemical and phase equilibrium flash problem at high temperature. Salt precipitation is a heterogeneous chemical equilibrium problem and must be determined by comparing equilibrium solubility products, , to ion solubility products, , for all possible molecular salts as shown in eq. (5.24  Note that the supersaturated VLE solution has a lower value of G/RT than the single phase solution. where is the molar amount of in the feed, [ ] + is the solubility limit of + in the aqueous liquid, 2 is the number of moles of water in the aqueous liquid, and 2 is the molecular weight of water.  1. The VLE + salt solution has a lower dimensionless Gibbs free energy than either the single phase solution or the supersaturated VLE solution.
2. The STRIP criterion of 0.7 -0.8 vapor fraction has been met in the final solution.
3. Salt precipitation is potentially a serious concern in this application of STRIP unless production water is mixed with clean water. composition points are generated for each temperature and pressure. This is to ensure that phase boundaries are smooth and that changes in V-only, L-only, VLE, LLE, and VLLE regions make physical sense. Table 5.6 gives computational details for the rigorous flash tests.

Example 5: Comparison of Steam Injection and STRIP
This first EOR example compares model results for steam injection and STRIP for a 3D heterogeneous reservoir formation containing light oil. Input data are listed in Table 5.7. In EOR applications it is typical to 'lump' components to reduce computational costs; thus for STRIP all light gases were treated as CO2, which is the solvent of interest.
In this example, the model is based on a fragment of the up-scaled SPE10 porosity and permeability fields (see Christie and Blunt, 2001). Here we used a grid size of 30×60×3 m 3 with a uniform grid block volume of 12×6×1.2 m 3 . The injection and production wells were placed at opposite corners of the reservoir. A single component, n-decane, was used to model the oil and the EOS used was the SRK equation. Steam injection was modeled using an injection stream of 1 mol% CO2 and 99 mol% steam while STRIP, which contains more CO2 from combustion, injected 10 mol% CO2 and 90 mol% water. The heat and water input were the same for steam injection and STRIP so an equitable comparison could be made.

Main Simulation Results
The performance of steam injection and STRIP are compared using the metrics of sweep and oil recovery.

Sweep Efficiency
The sweep efficiency can be deduced from oil saturation at the end of the operating period. Figure 5.3 shows the oil saturation in the reservoir for steam injection and STRIP after 2,000 days of operation.

Figure 5.3: Oil Saturation for Steam Injection (upper) and STRIP (lower) for Example 5
Note that the blue regions are much larger for STRIP than for steam injection, indicating that STRIP removes more oil. One can also make more quantitative measures of sweep efficiency using the following expression = Δ (5.29) where denotes the sweep ratio, Δ is the porous volume for which the oil composition has changed by 1% or more and is the total porous volume available to the oil. Figure   5.4 shows quantitative results for sweep ratio for steam injection and STRIP as a function of time. The sweep ratios for steam injection and STRIP after 2,000 days of operation are 60% and 83%, respectively. Clearly the sweep ratio of STRIP is superior to steam injection.  To compare oil recovery, the Original Oil in Place (OOIP) at surface conditions, which is a common assumption in the petroleum industry, must be computed. OOIP is calculated using formula OOIP = ∑ ϕ oil oil (5.30) where V is the block volume, is oil saturation, and is the surface-to-reservoir formation volume factor. Here again, STRIP outperforms steam injection -recovering 16,764 m 3 (105,442 barrels) more oil and leaving less OIP after 2,000 days. Table 5.8 summarizes the performance of steam injection and STRIP for this first example.

Example 6: Comparisons between Conventional EOS and CSAT
This second example compares a conventional reservoir simulation approach, which uses an EOS, to one that uses CSAT. For this example, pore volumes and permeability fields were taken from the upper layer of the original SPE10 model (Christie and Blunt, 2001). The simulations were performed using an initial reservoir composition of 1 mol% CO2, 49 mol% n-decane, 20 mol% n-hexadecane, and 30 mol% water and the initial reservoir pressure and temperature were 31 bar and 300 K, respectively. One injection and one production well were placed at opposite corners of the reservoir. The injection well operates under constant pressure and constant temperature conditions of 60 bars and 500 K. The STRIP injection fluid consisted of 15 mol% CO2, and 85 mol% water. The production well was set to a constant pressure of 3.45 bars. Input data for this example are shown in Table 5.9.  Christie and Blunt (2001)

Main Simulation Results
The details of the performance of STRIP are discussed along with the features of the simulator.

Sweep Efficiency
Oil and gas saturations provide enough information to quantify sweep efficiency. Figure 5.6 shows the oil and gas saturation in the reservoir for STRIP after 7,000 days of operation, where the x and y axes denote grid blocks and the color bar shows saturations. As expected, both CSAT and GFLASH provide identical results for gas and oil saturation. This is because CSAT only skips phase identification and rigorous flash computations when compositions are far from phase boundaries.
Simulation Statistics Table 5.10 summarizes the simulation statistics for this example and shows that rigorous flash solutions take the bulk of the computer time for the conventional EOS approach. CSAT, on the other hand, significantly decreases the number of EOS solves and, therefore, reduces total flash solution time by almost two orders of magnitude.

Coda
In addition to the chemical engineering sub-problems described in this work, there are also others, including sub-problems that require (1) The characterization of oils with many components.
(2) Development of better methods for determining viscosity and relative permeability in harsh conditions.
(3) Understanding chemical EOR methods and the associated phase equilibrium in the presence of surfactants and other chemical additives.
(5) Reaction kinetics models for gas hydrate formation and CO2 sequestration.
(6) Improved numerical methods for flash and for solving 'stiff' DAE systems.
In our opinion, the PSE community is ideally positioned to make contributions in these and other areas as they relate to reservoir simulation. densities (gprs_rho), densities derivatives (gprs_der_rho), equilibrium phase enthalpies (gprs_enth), and associated phase enthalpy derivatives (gprs_der_enth).
One variable that need some clarification is the integer variable RE_ENTER, which provides re-entry facilities in GFLASH, is defined as follows: RE_ENTER gives AD-GPRS complete control of GFLASH and is an important feature for reducing computational workload. RE_ENTER allows AD-GPRS to determine when it is necessary to solve a rigorous flash problem for a given grid block or when to skip solving the flash for that grid block. More specifically, when the rigorous solution of a flash problem is needed, AD-GPRS sets RE_ENTER = 0 or 1 and GFLASH solves a rigorous flash problem. As a result, the fugacity constraints for the given grid block are satisfied for the conditions of temperature, pressure and equilibrium phase compositions for the grid block and communicated back to AD-GPRS. On the other hand, when the temperature, pressure and equilibrium phase compositions for that grid block change but it is anticipated that the number of equilibrium phases in that grid block is unlikely to change, then AD-GPRS ask GFLASH to simply evaluate fugacities, densities, enthalpies, and their derivatives without solving a flash problem by setting RE_ENTER = 99. This approach works because AD-GPRS includes fugacity and other constraints in the equation set for a given time step, it converges the fugacity conditions defining equilibrium along with the conservation of mass and energy equationswhether or not they are satisfied at the beginning of the time step.

Appendix B: Standard Gibbs Free Energy, Standard Heat of Formation, and Ideal Gas Heat Capacity Data
The standard heat of reaction, ∆ 0 , is calculated using standard heats of formation data, ∆ , 0 , plus the simple equation Pure component ideal gas heat capacities are computed using a polynomial in temperature of the form Table B1 gives the heat of formation data, which was taken from Appendix IV in , and ideal gas heat capacity coefficient data, which has units of J/mol, taken from Reid et al. (1987).  Table B2 shows the standard Gibbs free energy and heat of formation data needed for computing equilibrium solubility products of molecular salts. Some of the data was taken from  while the remaining data was taken from the CRC Handbook of Chemistry and Physics.  3 -435.8 -748.1 -795.8 -1265.2 -1382.8 -1316.4 -1433.7 -1320.3 -1432.7 Nomenclature number of atoms j in molecule i total amount of atom j oil surface-to-reservoir formation volume factor 1 , 2 , 3

Background and motivation
Thermal multiphase flow and compositional reactive transport in porous media is the basis for simulation of almost all energy and environment-related industrial processes. The development of a simulation framework capable of modeling this class of problems on a continuous scale has been an important task in both the reservoir engineering and hydrology communities. Reservoir engineers usually deal with problems involving thermal multiphase flow and multi-component transport tightly coupled with complex phase behavior . These problems include different enhanced oil recovery processes such as steam or gas injection. Usually, chemical reactions are not treated as having a first-order impact on these models.
On the other hand, the hydrology community has been concerned with subsurface modeling of multiple components and multiple chemical reactions. The work by Lichtner (1985) laid the theoretical foundation for continuum models for mass transport and chemical interactions. Current chemical models include a wide variety of different reactions, including dissolution-precipitation and adsorption-desorption (Steefel et al., 2005). However, these models mostly deal with only the aqueous phase in slightly heterogeneous reservoirs. Reactive transport modeling in subsurface hydrology has never been fully coupled with equilibrium phase behavior of complex hydrocarbon mixtures in highly heterogeneous formations, despite some recent attempts (Flemisch et al., 2011). Due to the emerging interest in complicated subsurface dynamic processes like CO2 sequestration, methane hydrate recovery, and geothermal processes, there is a growing need in integrating full chemical reaction modeling capabilities with compositional reservoir simulation (Marchand and Knabner, 2014;Farshidi, 2016). Any heterogeneous structure of subsurface formations and the multiple scales of governing processes requires implicit time approximation for numerical solutions to be unconditionally stable on simulation time steps appropriate for the problem of interest.
The main purpose of this study was to develop, for the first time, capabilities for reactive transport modeling in subsurface hydrology fully coupled to equilibrium phase behavior of complex mixtures in highly heterogeneous formations within a numerical reservoir simulator. We tested our framework on a problem of particular practical importance -CO2 sequestration in aqueous aquifers. One of the major challenges in modeling this class of problems is accurate representation of dissolution trapping (Elenius et al., 2014 andElenius et al., 2015). Macroscopic dissolution rates can be enhanced significantly by gravity-driven currents (up to an order of magnitude). This complex behavior is strongly affected by many factors, including the chemical composition of the brine, different impurities in the injection stream of CO2, changes in pressure and temperature, and simultaneous chemical reactions. In addition, small scale precipitation and dissolution of minerals impacts the dynamics of gravity-driven flows and, in turn, effects the dissolution as well. This work is the first attempt to create a universal tool for predictive reservoir simulation of CO2 sequestration in aqueous aquifers that takes into account all of the complex mechanisms that effect the macroscopic dissolution rate. This dissolution rate can then be used in a realistic, largescale reservoir model using simplified physical models (Gasda et al., 2011 andLagasca, 2014) to predict the dynamics of CO2 trapping for medium time scales (i.e., tens to a hundred thousand years).

EOS modeling
Additional complexity in compositional modeling stems from phase behavior computations. An Equation of State (EoS) model is usually employed to describe the phase behavior of the multi-component system (Lake, 1989). For given temperature, pressure, and overall composition, EoS computations define the phase state and composition of each phase  The GHC EOS ) is a recent modification of the SRK  EOS that constrains the energy parameter, a, to satisfy the Gibbs-Helmholtz equation and uses Monte Carlo simulations to incorporate molecular length scale information. In particular, the energy parameter in the GHC EOS for pure components is given by equation (6.1).
( , ) = (0.42748 * 2 + ln(2) + 2 ( ) ln (2) ) − ln(2) − ( 2 ln (2) ) ( ) (6.1) where and are critical properties, is the molecular co-volume parameter, and is the gas constant, and is the molecular scale internal energy of departure at the given temperature and pressure. Pure component is determined a priori over wide ranges of temperature and pressure using Monte Carlo molecular simulation, stored in look up tables, and readily up-scaled to bulk phase EOS calculations using eq. (6.1).
The novel features of the GHC equation include the use of molecular information in the energy parameter expression and the estimation of b from pure component density data.
The details of the derivation of the GHC EOS, the extension to non-electrolyte and electrolyte mixtures can be found in the literature Lucia et al., 2012a,b;Lucia et al., 2015). It is important to note that the GHC EOS only uses parameters based on pure component properties (e.g., mixture critical properties from Kay's rules) and pure component Monte Carlo molecular simulation (mixture from a linear mixing rule), even in systems containing electrolytes, and is truly predictive.

Monte Carlo simulation details
Application of the multi-scale GHC EOS requires prior knowledge of pure component internal energy of departure for the components in the system at relevant conditions. Monte Carlo molecular simulations in the isobaric-isothermal ensemble are used to generate the necessary information a priori and create lookup tables. The specific details of the simulation change based on the component and molecular model being used, however generally a 6-12 Lennard-Jones potential is used to account for van der Waals interactions, along with the recommended cut off (depending on the model) and tail corrections. The Coulomb potential with an Ewald summation is used to account for electrostatic interactions. A cubic box with periodic boundary conditions is used in all bulk fluid simulations. Standard center of mass translation and rotation moves, as well as isotropic volume moves are applied. Simulations are run in equilibration mode for a number of cycles (depending on the system) in which the maximum displacement and rotation are frequently adjusted to maintain 50% acceptance rates of translation and rotation moves, respectively. The system is then switched to production mode, in which the maximum displacement and rotation are fixed, for sampling. Typically, four parallel sets are run for the same system and the results are averaged. References for the simulation parameters for the components used in this work can be found Table B3.
Simulations are typically run either using the open source Towhee MCCCS software (Martin, 2013) or an in-house FORTRAN program. MC simulation runtimes can range from a few hours to a few days, depending on the complexity of the molecular model, the number of particles included in the simulation, and the potential model(s) used.

GFLASH
The GFLASH library is a multi-component, multiphase, isothermal, isobaric (TP) flash calculation program written in FORTRAN. Given an overall composition, temperature, and pressure of a fluid mixture, GFLASH has the capability to determine the number of existing phases at equilibrium and calculate their compositions, densities, enthalpies, fugacities, and all property derivatives with respect to pressure, temperature, and composition. The main capabilities and details are outlined in previous publications Lucia, et al., 2015). The main reasons for the use of GFLASH in this work are the implementation of (1) a robust stability and flash algorithm and (2) the GHC EOS, and (3) the capability to handle simultaneous phase and homogeneous/heterogeneous chemical reaction equilibrium with mineral deposition/dissolution.

Description of Reaction Equilibrium Model
A full detailed description of the numerical methods used for solving the equilibrium reactions for molecular salt formation included in this work can be found in the paper by Lucia et al. (2014). For the examples studied in this work, the formation/dissolution of solid salts was limited to NaCl, Na2CO3, CaCl2, and CaCO3, described by the following reactions: In addition, the following reaction of dissolved carbon dioxide with water to generate carbonate ion was included: Reaction (6.6) was obtained by summing the reactions in the carbonate series, and equilibrium concentrations of carbonic acid and bicarbonate ion, which were not of interest in the examples studied. Also, the formation of hydronium ion, and therefore changes in pH due to the dissolution of CO2, were neglected. This is because Soong and coworkers (2004) have found that brine with a pH of 11 produces the most CaCO3 when reacted with CO2. Therefore, it is assumed that the original pH of formation brine was high enough to fully support CaCO3 precipitation. For a brine solution with pH of 11, the change in pH due to the amount of CO3 2produced from CO2 in the examples in this study was, in fact, negligible, as shown in Appendix A. The equilibrium constants for the reactions were calculated from tabulated standard Gibbs free energies of formation data and corrected for temperature using tabulated standard enthalpies of formation and the van't Hoff equation. See Appendix B.
The primary goals of this study were (1) to demonstrate that the coupled ADGPRS/GFLASH software system has the capability to accurately model mixtures in which minerals dissolve and/or precipitate and (2) to show that ADGPRS/GFLASH can accurately model CO2 sequestration with residual, dissolution, and mineral trapping.

ADGPRS
The reservoir modeling software used in this work is called Automatic Differentiation General Purpose Research Simulator (ADGPRS) and was developed and maintained by SUPRI-B research group at Stanford University. ADGPRS is written mainly in C++, and widely used throughout the reservoir and petroleum engineering communities because of its wide-ranging capabilities, which include ): 1. Flexible treatment of all nonlinear physics.
2. A fully thermal-compositional formulation for any number of phases.
3. Multi-phase CSAT for efficient and robust computation of phase behavior.  respectively for additional information regarding the implementation of the GHC EOS and the multiphase flash algorithm and handling of salt precipitation within GFLASH.

Natural formulation for reactive systems with precipitation and dissolution
The basic information flow between ADGPRS and GFLASH is given in Fig. 6.1.
Specifically, for each grid block in the reservoir ADGPRS passes the current estimate of the temperature, pressure and overall composition of that block to GFLASH.
GFLASH, in turn, uses that information to determine the number of equilibrium phases, their amounts and compositions, in this case using the GHC equation of state, and returns this information to ADGPRS. Many of these details are described in previous publications . In this section, specific modifications for the treatment of flow and transport in the presence of chemical (equilibrium) reactions and details of the coupling between ADGPRS and GFLASH are described. For problems with precipitation and dissolution with equilibrium reactions, new (solid) phases must be introduced into the general natural variables logic ).

Governing equations
The typical governing equations for CO2 sequestration include conservation of mass for each species and closure assumptions and constraints. The mass conservation equations for isothermal compositional systems can be written as where P is the number of co-existing thermodynamic phases, R is the number of reactions, and C is the number of species. Also, ϕ is the porosity of the porous media, ρ is molar density, x denotes composition in mole fraction, S is phase saturation, V is volumetric (Darcy) flow, J is molar diffusion flux, q denotes a mass source or sink term, r is reaction rate, and v are stoichiometric coefficients.

Rearrangement of equations
Eq. (6.7) can be written in the general matrix form given by where a corresponds to an accumulation vector of length C, I is a C vector of fluxes, q is the well source term vector, also of length C, r is the reaction rate vector of length Q, and V is a CxQ stoichiometric matrix.
Following the logic described in Farshidi et al. (2013), we introduce the ExC matrix, E, which represents the stoichiometry for each element associated with the reactions of each species. In general, this matrix can be determined by solving the equations × = 0 (6.9) Multiplying Eq. (6.8) by E gives E element mass conservation equations of the form + + = 0 (6.10) To close the system, an additional C−E independent equilibrium constraints are needed and take the form where the multi-phase flash procedure described in Section 6.2.2 is used to define the equilibrium solubility product, ( , ), using Gibbs free energies of formation and the identity of the minerals that precipitate. The ion solubility product, ( , , ), on the other hand, is defined using the actual ion concentrations in the brine determined from the GHC equation of state.

Mineral precipitation
When mineral precipitation is controlled by chemical equilibrium, a procedure is needed to account for the appearance and/or disappearance of solids. Since precipitation only occurs when the ion solubility product for any particular salt exceeds its equilibrium solubility product, a more general form of eq. (6.10) is needed.
( , , ) > ( , ); (6.12) ( , , ) > ( , ); Once the potential for precipitation is identified, new conservation equations for those mineral components and mineral phases of the form = (6.13) must be added to the original system given in eq. (6.8). Additional unknown variables must also be added to the set of unknowns corresponding to the concentrations of the solid species. Finally, the matrix E should be modified and include a new rate corresponding to the reaction for the precipitated mineral, see Farshidi (2016) for details.

Illustrative example
In this section, a simple example is presented to illustrate the extension of the compositional AD-GPRS/GFLASH framework to reactive systems with five components and three phases as shown in Table 6.1.
For simplicity, we assume that the brine phase always exists. Treatment of the disappearance of water or brine phases can be found in Farshidi (2016). For the illustrative example, there are only four possible combinations of co-existing phases, as shown in Table 6.2.
The matrix E for this example is Eqs. (6.14) and (6.15), each row corresponds to an element (i.e., carbon, hydrogen, oxygen, etc.) while each column corresponds to a component (i.e., molecular or ionic species). Finally, a correct set of unknowns (and equations) can be easily constructed for each combination of phases shown in Table 6.2.

Numerical results and discussion
In this section, numerical results for four separate CO2 sequestration examples are presented to illustrate the robustness of the proposed methodology in capturing the correct physics of solubility and mineral trapping, carbonate chemistry, and mineral precipitation and dissolution.

Example 1 -Large scale model with solid precipitation and dissolution
This first example is a simple reservoir model with the fluid system of CO2-H2O-Ca 2+ -Na + -Cl -CO3 2-. For the conditions shown in Table 6.3, there are many different equilibrium phase states possible (i.e., VLE, VLLE, LLE, SLE, SLLE, etc.), depending on temperature, pressure and composition throughout the reservoir. The main purpose of this example was to demonstrate the ability of the ADGPRS/GFLASH system to model solid precipitation and dissolution; therefore, carbonate chemistry was not included and represents a system in which a CO2-rich injection stream is introduced into a formation containing a single-phase brine. It is set up in a 50 × 50 × 1 grid to study horizontal propagation of the injection stream. This particular example contained an injection well in the lower left corner, and a production well in the upper right-hand corner. As shown in Table 6.3, the reservoir was initialized with solid NaCl in each block and a fluid with a composition different than the composition of the injected fluid.
As the simulation of this first example progressed, the pressure changed until supercritical CO2 broke into the production well in the upper right-hand corner. Also, the flow of CO2 from the injection block caused an increase in CO2 composition in many blocks of the reservoir, eventually causing a second CO2 rich phase to appear in the system (see Fig. 6.2). Fig. 6.3, on the other hand, shows the equilibrium phase state (i.e., LLE, SLLE, VLE, etc.) of each block in the system as the simulation evolved. Clearly, Fig. 6.3 shows that as fluid flowed out of the feed block into the reservoir, NaCl dissolved before the Ca 2+ in the reservoir could react with CO3 2from the feed block to form CaCO3.

Example 2 -Large scale simulation of CO2 injection
The purpose of this second example was to demonstrate that the coupled ADGPRS/GFLASH system could successfully model CaCO3 precipitation/dissolution. This second example focused on buoyancy driven vertical migration of a CO2 plume, which is a key feature of many carbon sequestration studies, the use of approximate carbonate chemistry, and its impact on the precipitation of CaCO3 as described in Section 6.2.2.1. Here we assume that the CO2 plume is trapped in the geological formation and monitor short time-scale mineralization processes. This example used the same fluid system that was used in Example 1 with conditions shown in Table 6.4. More specifically, the reservoir in this example was homogeneous with a porosity of 0.18 and discretized with a 25 × 1 × 25 grid. The grid spacing was also homogeneous and set to 8 m in each direction, with absolute permeability equal to 150 mD and 220 mD in the x and z directions respectively. This system was first equilibrated to approximate hydrostatic equilibrium so that the pressure in the initial system varied with depth from 250 bar to 270 bar. Pure CO2, which was less dense than the surrounding reservoir fluid, was injected into the middle block in the bottom row of the reservoir. The increased amount of CO2 that dissolved in reservoir brine resulted in an increase in dissolved carbonate ions in the aqueous phase, and under the model conditions defined in this example, also resulted in the precipitation of CaCO3. As this CO2 plume rose, it was observed that the amount of dissolved CO2 in the aqueous phase in the blocks near the plum increased. This, in turn, led to the generation of CO3 2-, which then reacted with dissolved Ca 2+ in the reservoir to form solid CaCO3.   It is clear from Figs. 6.4-6.6 that the CO2 plume migrates upward, causing an increase dissolved CO2 in the brine phase, leading to precipitation of CaCO3. Moreover, in this example the only source of CO3 2− ions was the equilibrium reaction between dissolved carbon dioxide and water and the resulting precipitation of CaCO3 in the reservoir clearly illustrates that the coupled ADGPRS/GFLASH system has the capability to model salt precipitation in the presence of carbonate chemistry.

Example 3 -Small scale model with solubility trapping
In this section, we use a simulation model from Elenius et al. (2015) to (2015), we used the GHC EOS for property evaluation at different thermodynamic conditions to account for the presence of aqueous ions. Table 6.5 gives the initial compositions for the model. We used an initial pressure distribution starting from p = 240 bar at the lower part of the model and constant temperature T = 345 K. No minerals initially precipitated in the model. Due to molecular diffusion, the initial portion of CO2 in the plume starts dissolving in the brine, and the difference in density of the brine with the dissolved CO2 (since it is heavier than the original brine) initiates the formation of unstable fingers of CO2-rich brine. These fingers enhance the dissolution rate of CO2 several fold and significantly increase the trapping capability of the aquifer due to the higher dissolution.
The numerical simulation of this process requires very fine resolution of the simulation grid which makes it prohibitive for a full field simulation (Elenius et al., 2015). Small scale models can predict CO2 dissolution rate quite accurately, which can be used in various up-scaling models (Gasda et al., 2011;Lagasca, 2014) to predict the migration distance of the CO2 plume in a large aquifer over medium time-scale (tens to one hundred thousand years). The left part of Fig. 6.7 shows the composition of CO2 in the brine phase at different times. In the right part of Fig. 6.7 is shown the dissolution rate of CO2 as a function of time (years). The dissolution rate and trend are similar to ones reported in Elenius et al. (2015).

Example 4 -Small scale model with combined solubility/mineral trapping
It is believed that due to the time-scale of chemical reactions in brine-CO2 systems, mineral trapping does not affect the early stages of the CO2 sequestration process. However, small-scale precipitation may affect dissolution trapping due to the change in the dynamics of the instabilities. Precipitation and dissolution can change the porosity, which in turn can affect diffusion and the formation of fingers. In the presence of a capillary transient zone, this process will be coupled to phase behavior and become quite challenging to predict without a reliable simulation tool. Any inaccuracy in the prediction of the small-scale dissolution rate can change the predicted capacity of aquifers used for a large-scale sequestration several fold (see the example in Elenius et al., 2015). Here we present unique simulation results where the physics of all CO2 trapping mechanisms are taken into account. For the simulation of CO2 fingering in the presence of chemical precipitation, the second type of small-scale model from Elenius et al. (2015) was used. In this model, the entire two-phase region was located in the high volume area, which maintains the original CO2 profile. This model represents the leading part of the CO2 plume. For simplicity we used the same configuration and composition of the plume as in the previous simulation, but increased the temperature to T = 380 K, which decreased the CO2 solubility to xCO2 = 0.016. At conditions described in Table 6.5, the higher concentration of CO2 triggers the carbonate reaction (Eq. (6.6)), which increases the concentration of carbonate ion and, in turn, initiates precipitation of CaCO3 based on Eq. (6.5).
The influence of chemical precipitation on the generation of fingers is shown in  Note that the difference in CO2 distribution is insignificant and clearly suggests that the influence of early-time mineralization can be ignored in the accurate estimation of the CO2 dissolution rate for the system studied. However, this does not guarantee that at different thermodynamic or chemical reaction conditions, the same conclusion will hold.
In the next simulation the same system is kept, but amplify the molar volume of minerals to update the porosity 1000 fold. It can be seen in Fig. 6.9 that the porosity variation in this case is more significant due to the larger pore volume occupied by minerals resulting from precipitation. In this case, the dynamics of the fingering process have changed, which affects the subsequent macroscopic dissolution rate. Obviously, the last case is hypothetical since the molar volume of CaCO3 for the porosity update is unrealistically large.
In the next simulation, we assume the presence of NaCl with an initial concentration CNaCl = 1 kmol/m 3 and the same under-saturated brine. The dynamics of the instability completely changes in this case, and no fingers are formed in the reservoir formation. See Fig. 6.10. Here, the increase in porosity due to the dissolution of NaCl in the diffusion zone and the reduction of solubility due to the presence of Na+ and Cl− ions stabilizes the CO2-brine interface and prevents fingers from growing.

Conclusion
A multi-scale framework for modeling and simulating reactive transport in subsurface hydrology with fully coupled equilibrium phase behavior of complex mixtures in homogeneous to highly heterogeneous reservoir formations was developed.
The key attributes of this numerical reservoir simulation framework include the capabilities to model (1) simultaneous chemical/phase equilibrium, (2) homogeneous and heterogeneous chemical reactions, (3) aqueous electrolytes, and (4) mineral precipitation/dissolution. Simultaneous multi-component, multi-phase/chemical equilibrium was modeled using a combination of the multi-scale Gibbs-Helmholtz constrained (GHC) equations of state to describe the behavior of all fluid phases and Gibbs free energies and enthalpies of formation to predict equilibrium solubility products. Precipitation was identified using comparisons of equilibrium and ion solubility products. Governing partial differential conservation and constraint equations were solved using a fully implicit method (FIM). All modeling capabilities were implemented in the coupled software system ADGPRS/GFLASH.
Four different types of processes related to subsurface CO2 sequestration were used to test the efficacy of the proposed framework. Using the multi-component mixture CO2 and brine with dissolved ion Na + , Ca 2+ , Cl − and CO3 2− , the first example demonstrated that the modeling framework successfully predicted the formation and dissolution of solid salts NaCl and CaCO3 throughout an isothermal and homogeneous reservoir formation in the presence of brine and a CO2-rich super-critical fluid. The second example was similar to the first but also included two additional modeling challenges (1) the equilibrium reaction between CO2 and water to form carbonate and hydronium, H3O + , ions and (2) vertical migration of a CO2 plume. Therefore, this example was used to demonstrate combined homogeneous/heterogeneous chemical reactions, salt dissolution, buoyancy driven flow, and simultaneous phase and chemical equilibrium. It was also shown that the carbonate chemistry resulted in insignificant changes in the initial pH of the reservoir, which was assumed to be 11. Accordingly, pH changes were neglected in this second example. Here again, numerical results showed that the proposed framework was capable of successfully modeling all of the relevant physics. The third and fourth examples were used to demonstrate that ADGPRS/GFLASH can capture all of the higher fidelity physics of the interplay between residual, solubility, and mineral trapping in the presence of unstable CO2 fingers, which enhance the macro-scale CO2 solubility in brine, following buoyancydriven flow of a CO2 plume.
Accurate modeling of all physics in this work (i.e., rigorous, complex EOSbased simultaneous phase and chemical equilibrium involving homogeneous and heterogeneous reactions) is essential for developing a high-fidelity model for carbon sequestration in any reservoir-specific environment. To our knowledge there is no software system currently available with all of these capabilities. The use of a rigorous EOS allows us to treat impurities such as O2, Ar, SO2, CH4, N2, H2S, etc. Impurities can have strong effects on the density of the CO2-rich phase as well as the solubility of CO2 in the aqueous phase, which in turn can affect the migration and spread of the injected CO2 plume and CO2 storage capacity (Sin, 2015). Therefore, the next phase of this work will include the development of models that include impact of impurities on CO2 sequestration.

Acknowledgements
We acknowledge financial support from the Petroleum Institute and Abu Dhabi National Oil Company. We also acknowledge the SUPRI-B program at Stanford University for permission to use ADGPRS and Sara Farshidi and Maria Elenius for helpful discussions.

Appendix A -pH Change due to CaCO3 formation
Let the pH of the brine solution be 11. It is demonstrated that the amount of hydronium ion produced is small so that associated changes in pH can be neglected.
The maximum CaCO3 concentration obtained in examples in this work was ∼1.3 × 10 −3 mol/L. From reactions (5) and (6), the concentration of hydronium ion corresponding to this amount of CaCO3 is [H3O + ] = 2.6 × 10 −3 mol/L. Also, the equilibrium constant for the dissociation of water into hydronium and hydroxide ions at 350 K and 250 bar using linear interpolation is Kw = 1 × 10 −12.58 (see, Bandura and Lvov, 2006). At a pH of 11, the initial hydroxide concentration is Therefore, for the examples studied in this article, the amount of hydronium ion produced as carbonate ion is formed does not significantly change the pH of the formation brine. A Corey-type relative permeability model was used to define the relative permeability as a function of saturation (Eqs. B1-B3), the parameters are given in Table B In the third example, the residual water and CO2 saturations were each set to 0.2, and the relative permeability and saturation distribution were determined by specifying a static pressure distribution in each phase (see Elenius et al., 2014, Elenius et al., 2015 to create a capillary transition zone at the interface between the two phases. The viscosities of the brine phase and the CO2-rich phase were assumed to be constant

Appendix B: Fluid and Rock Properties
in all examples in this work. Values of 0.511 and 0.061 cP were chosen for the brine and CO2 phases, respectively, based on values used in a similar study (Elenius et al., 2015).

Introduction
Understanding the flow of pure fluids and their mixtures in porous media is of interest to a wide variety of practical applications including oil and gas exploration/recovery, carbon dioxide sequestration, and ground water flow. Several software tools exist that allow one to solve the set of partial differential equations Recently, AD-GPRS has been adapted for modeling CO2 sequestration application (Voskov et al., 2017) such as CO2 injection into a deep saline aquifer, where mineralization of CO2 can impact the dissolution rate of CO2 and hence affect the storage capacity of the formation.
MODFLOW is a mature open source numerical groundwater simulation program developed by USGS. It has been used to study the interaction between ground and surface water, the long and short-term effects of pumping, and the change in the water table as a function of variations in seasonal inflow and outflow. MODFLOW has a wide user base, and many extensions have been developed to add features to the basic MODFLOW package. Including the ability to adjust for variable density fluid flow (a posteriori), track particles originating in some point source for a given system, and include many different boundary conditions representing surface features such as rivers and streams.
The tools and extensions developed for MODFLOW are very useful for studying groundwater systems. However, the main MODFLOW package makes several major assumptions about the composition of the fluid and its physical properties, which limit its applicability. The specific assumptions include that (1) the fluid is a single component (water), (2) the density of the fluid is independent of pressure, temperature or composition, and (3) the fluid is stable in a single phase (liquid) state. For many cases, these assumptions are completely reasonable. However, in systems where pressure change is expected to be large, or additional phases may be present in the system, these assumptions are not adequate to represent the physics of the system. For example, a carbon sequestration example was recently presented in Voskov et al. (2017) using AD-GPRS in which the variation in density with composition and pressure strongly effects the dissolution rate CO2 from a supercritical carbon dioxide phase into a brine phase. In contrast, the advanced features implemented in AD-GPRS allow the study of groundwater systems in which density variation due to pressure or fluid composition may strongly affect the results, multiple phases may be present, and/or the distribution of the fluid composition at a given time step is of interest to the user.
Due to differences in the formulation of the model equations, implementation, and availability of various boundary conditions, a direct comparison between MODFLOW and AD-GPRS is not feasible. Therefore, the purpose of this study is to demonstrate the efficacy of using AD-GRPS to model groundwater flow in a realistic groundwater system, to highlight the advantages gained from using AD-GPRS in these systems, and to identify the weaknesses and/or challenges that should be overcome to apply AD-GPRS to these systems more efficiently. To accomplish these goals, a simulation model for the AD-GPRS simulator was developed using a previously developed MODFLOW model as a guide. The reservoir model used in this study is the Big River Management area in central Rhode Island, published in the USGS report of Masterson and Granato (2012).

Procedure
The AD-GPRS reservoir simulator  developed by the SUPRI-B group at Stanford University was used in the work. AD-GPRS is a state-ofthe-art, multi-phase, fully thermal and compositional reservoir simulation program. It has been used to simulate many industrial processes including steam injection and steam injection with propane co-injection for enhanced oil recovery , in situ CO2-steam co-injection for heavy oil recovery, and plume migration in carbon dioxide sequestration in the presence of convective dissolution and gravity currents (Elenius et al., 2015;Voskov et al., 2017). In this work, AD-GPRS is used to model and perform fully compositional simulations of groundwater flow.
Treating the groundwater flow problem in this way enables one to study problems in which the composition of the formation water is not constant throughout the reservoir.
Thus, it is possible to monitor the concentration of any one species in the mixture at any location in the model at any given time step.
In this section, we will briefly highlight some of the features of the software used in this work and discuss the development of the simulation model using the published MODFLOW model of the Big River Management area as a guide. Fluid properties were calculated using a program developed by Professor A. Lucia at URI called GFLASH, which has been interfaced with AD-GPRS.

AD-GPRS Overview
The subsurface flow numerical simulation program used in this work is called AD-GPRS. AD-GPRS is a simulation program written primarily in C++ and maintained by the SUPRI-B project at Stanford University. It is used extensively in the reservoir engineering community due its wide-ranging capabilities, which include , discretization schemes are described in ), Voskov (2011 gives a description of the methods for solving non-linear and linear systems of equations, and the many approaches for fluid phase behavior computations are presented in the work of .

GFLASH Overview
The GFLASH software developed by Professor A. Lucia was used for all the fluid properties calculations in this work. GFLASH is a FORTRAN program suite which given temperature, pressure, and overall mixture composition calculates the number of equilibrium phases and their compositions, along with their respective densities, fugacities, enthalpies. Several commonly used cubic equations of state (EOS) models are implemented in GFLASH including Soave-Redlich-Kwong (SRK) equation , the Peng-Robinson (PR) equation (Peng & Robinson, 1979), and the Gibbs-Helmholtz constrained (GHC) equation of . For this work, the GHC EOS was used exclusively because it predicts the density of water more accurately than the traditional cubic equations of state  and because it is applicable to aqueous electrolyte solutions. Neither the SRK nor PR equations can treat aqueous electrolytes. Additional information about the methodologies in GFLASH for the solution of the classical isothermal, isobaric flash problem can be found in the literature , as can more information on the GHC EOS model .

Converting the MODFLOW Model
The AD-GPRS model used in this work was developed from the input files from

Results and Discussion
Simulation results using the model described in the previous section are detailed in this section for three distinct cases. In the first case, the inflows and outflows to/from the active model area represent rivers or streams entering or exiting the model area.
There is no other pumping/injection in the model. The second case is a modification of the first in which pumping is performed at designated locations for a period of 5 days.
The effects of pumping at a specified rate for 5 days is analyzed. This type of analysis is important to determine amount of pumping that is acceptable for an aquifer of interest, and to understand the effects of pumping on the water level in the surrounding areas.
Finally, the third example aim to highlights the utility of running a fully compositional simulation by introducing a contaminant into one or more of the inflows. The composition distribution can then be analyzed over time to understand the spread of contaminant, the potential size of the contaminated area, and possible remediation strategies.

Example 1: Simple case with no pumping
This case was initialized using the same initial pressure distribution as the MODFLOW model. Locations of inflows and outflows across the boundary of the model area were also adapted from the MODFLOW model. The simulation was run for 1000 days in AD-GPRS to allow it to reach steady state and for hydrostatic equilibrium to be reached. The composition of the fluid in the reservoir is given in Table 7.1. The GHC EOS was used to calculate the fluid density as a function of pressure and the necessary derivatives.  reference. Note that in the first case, there was no pumping from the wells, all inflow and outflow sources were located on the top model layer, and the pumping wells were located in the fifth layer as in the MODFLOW model. The results from the 1000 day equilibration simulations with no pumping depicted in Fig. 7.1 were used to initialize a second simulation in which water was pumped out of the system at the two production wells (the gray dots in Fig. 7.1) for 5 days. Note that the results using AD-GPRS/GFLASH shown in Fig. 7.1 are in qualitative agreement with those shown in Fig. 9 of Maserson & Granato (2012).

Example 2: Pumping response in 5 days
The model parameters for this example are exactly the same as in the previous example, with the exception that the two production wells were added at the locations indicated in Fig. 7.1. Figures 7.2 and 7.3 show the change in head after pumping for 5 days in layers 1 and 2 respectively.   Contaminant concentration as function of distance from the source is depicted in Fig.   7.8, for 1000, 3000, and 5000 days. Clearly, the contaminant concentration decreases as distance from the source increases, as expected.
Examples such as the previous two can also be used to predict the concentration of contaminant at any point within the model at any given time. This type of analysis could prove useful in evaluating different contamination scenarios by varying the location and concentration of contaminant introduced into the reservoir. In addition, the AD-GPRS/GFLASH modeling framework makes it possible to introduce contaminant levels greater than the corresponding solubility limits of those contaminant in water, which in turn could result in the formation of a second fluid phase. In addition, the migration of multiple fluid phases throughout the reservoir is easily modeled in the AD-GPRS/GFLASH framework by including a model of the relative permeability as function of saturation. Figure 7.9 shows the saturation of the octane-rich phase in the two-phase example described previously. The octane concentration introduced into the source block is greater than the solubility of octane in water and a second liquid phase forms and propagates through the reservoir model.

Conclusion
The numerical reservoir simulator AD-GPRS was used to investigate (1)