Perple_X Electrolyte Speciation at Elevated Pressure and Temperature


Forward and Back-Calculated Aqueous Electrolyte Speciation


Thermodynamic Details and Limitations


Thermodynamic Data Files for Electrolytic Fluids


Back-Calculated Electrolyte Speciation (Figure 1)


Relevant Options

Lagged Forward-Calculated Electrolyte Speciation (Figure 2)


Relevant Options

Brute-Force Forward-Calculated Electrolyte Speciation (Figure 3)


Relevant Options

Forward- and Back-Calculated Aqueous Electrolyte Speciation

Electrolyte speciation can be calculated in Perple_X 6.7.9 in three different ways: in lagged (Figure 2) and brute-force forward-calculations (Figure 3) electrolytic fluid speciation is evaluated during optimization (as in programs such as EQ3 and GEMS), whereas in back-calculation (Figure 1) the chemical potentials obtained from a forward calculation or from optimization without considering the effects of electrolytes are used to compute electrolyte speciation. For most purposes lagged forward-calculation is the best method, but it cannot be used to treat certain types of phase separation. Back-calculated speciation based on a pure-solvent forward model (Galvez et al. 2015) often provides an excellent approximation for solution chemistry in rock-dominated hydrothermal systems and is useful for verification of lagged and brute-force forward-calculations. Brute-force forward-calculation is the normal mode of treating solutions in Perple_X and is robust with regard to any kind of phase separation, but has limited precision and the configuration of electrolytic solution models requires user expertise.

Thermodynamic Details and Limitations

Perple_X's implementation is generic, but the extended HKF (Helgeson-Kirkham-Flowers) formulation (Perple_X EoS 16; Shock et al. 1992, Sverjensky et al. 2014) is the only formulation that can account for the change in dielectric constant due to non-aqueous solvent species (Galvez et al., 2015). The implementation of SUPCRT data was benchmarked against Harrison & Sverjensky's (DEW, 2013) spreadsheet calculator, minor discrepancies arise because of the choice of the equation of state for water. Specifically, Harrison & Sverjensky (2013) use, though they offer other options, Zhang & Duan (2005) to compute the density of water. This density is then used in the Sverjensky et al. (2014) expression for the dielectric constant and the Gibbs energy of water is taken from Delaney & Helgeson (1978). In Perple_X all properties of water are obtained by whatever equation of state is specified (the initial implementation of the HKF in Perple_X 6.7.3 did use Zhang & Duan's [2005] EoS for density, this EoS will be restored once the routine has been modified to compute water fugacity). Of the available water EoS in Perple_X: the CORK and HSMRK EoS should not be used at ambient conditions because they do not reproduce liquid water properties; the Haar EoS provides the most accurate representation of low pressure (< 1 GPa) water properties; and the PSEoS and MRK EoS give reasonable liquid densities at 298.15 K and 1 bar. I have not checked the boiling curves predicted by the latter two EoS.

The extension of the HKF to non-aqueous solvents follows if it is assumed that the Born coefficient (ω) of a solute is only a function of pressure (P) and temperature (T), this is formally an assumption of the HKF. The HKF then gives the change in the reference Gibbs energy of the solute due to a change in the dielectric constant (ε) of the solvent as

ΔG0_solute = ω0(1/ε0_solvent - 1/ε0_H2O),               [1]

where 0 denotes the reference P-T condition. The solvent contribution to the solute Gibbs energy at elevated pressure and temperature is:

ΔG_solute = ΔG0_solute + ω(1/ε_solvent - 1) -  ω0(1/ε0_solvent - 1) - (T - T0)ΔS0_solute

which simplifies to

ΔG_solute = ω(1/ε_solvent - 1) - ω0(1/ε0_H2O - 1) - (T - T0)ΔS0_solute, 


ΔS0_solute = ω0([∂ln(ε0_H2O)/∂T]/ε0_H2O - [∂ln(ε0_solvent)/∂T]/ε0_solvent)               [2]

is the change in the entropy of an HKF solute due to a change in solvent chemistry. This entropic effect is relatively minor and is not currently evaluated in Perple_X to avoid the computational overhead associated with evaluating [∂ln(ε0_solvent)/∂T]/ε0_solvent for every solvent composition considered. Instead the latter term is approximated by the value for [∂ln(ε0_H2O)/∂T]/ε0_H2O commonly cited in SUPCRT applications (-5.7986510-5/K). This practice should be changed for applications where the solvent is not water-rich. The S0 values in SUPCRT data do not include the solvent term -ω0([∂ln(ε0_H2O)/∂T]/ε0_H2O), hence a small discrepancy between the S0 values indicated in HKF thermodynamic data files and the solute entropy obtained by Perple_X's FRENDLY thermodynamic calculator. 

The HKF includes several additional parameters (e.g., the characteristic temperature, θ, and pressure, ψ) that may mask dependencies on solvent chemistry. Even if Perple_X's treatment, which assumes that standard state solute properties are influenced by solvent chemistry only through the dielectric constant, is justified, some of the conditions used to discriminate the different regions of applicability for the empirical expressions for the HKF correlation function "g" used to compute Born coefficients (Shock et al., 1992) are clearly H2O specific. As the g function is non-zero only for water densities < 1 g/cm3, this particular limitation is likely to be unimportant for calculations at high pressure. For dilute non-aqueous solvents, the influence of solvent chemistry on the bulk dielectric constant and solute thermodynamic properties (as predicted by the HKF) is minor, thus these effects may well be insignificant in comparison to other sources of uncertainty.

Solvent dielectric constants for pure non-aqueous volatile species are computed from the empirical expressions given by Harvey & Lemmon (2005), but using the partial molar volume of the species in the impure solvent rather than the volume of the pure species (Galvez et al. 2015). The dielectric constant of pure water is using the expression given by Sverjenksy et al. (2014), which is indicated to be valid to 6 GPa and 1473 K. The expression of Sverjensky et al. (2013) yields a dielectric constant of 81.39 at 298.15 K and 1 bar for a density of ~1 g/cm3, this value differs from the value of 78.47 usually used in SUPCRT and leads to a change in the standard state Gibbs energy of solutes equal to ω0(1/81.39 - 1/78.47) (Eq [1]). The dielectric constant of the impure solvent is computed using the Looyenga mixing rule (Mountain & Harvey 2015) with volume fractions computed from the partial molar volumes of the solvent species.

A minor contradiction between the HKF and Debye-Hueckel theory, is that the latter identifies charged species and solvent, whereas the former differentiates the solvent (water), neutral species, and charged species. In the context of Debye-Hueckel theory, neutral solute species are properly considered solvent species. In Perple_X: solvent species are considered to be any neutral species that is assigned a thermodynamic reference state in which the species is pure and neutral solute species are any neutral species described by a solute reference state. Activity coefficients for electrolytes are computed by the Davies extension of the Debye-Hueckel equation, other formulations can easily be incorporated. Neutral solutes are assumed to be ideal, though, if available, interaction terms can be accommodated. Solvent species activities are computed from the molecular EoS selected by the user. As has been noted (e.g., Wolery 1990) this approach violates the Gibbs-Duhem relation. Wolery (1990) uses the Gibbs-Duhem relation to derive an expression for the activity of water in a Debye-Hueckel fluid. This approach cannot be extended to a multicomponent solvent. These inconsistencies could become significant at high solute concentrations.

Thermodynamic Data Files for Electrolytic Fluids

To permit calculations with silicate solution models, a file containing Harrison & Sverjensky's (DEW, 2013) revised SUPCRT/HKF data (DEW13ver.dat) has been merged with the hp02ver.dat (Holland & Powell 1998, rev. 2002) and hp622ver.dat (Holland & Powell 2011, rev. 2013) thermodynamic data bases. The merged files are DEW13HP02ver_elements.dat, DEW13HP02ver_oxides.dat, DEW13HP622ver_elements.dat, and DEW13HP622ver_oxides.dat, where the suffixes _elements and _oxides indicate whether the data is formulated in terms of element or oxide components. The Holland & Powell data bases include electrolyte data, but it appears the data is not maintained. To avoid conflicts with the revised SUPCRT data, the Holland & Powell electrolyte data has been moved to the file hpAQver.dat. The Holland & Powell aqueous data is for an equation of state that assumes solute thermodynamic properties at elevated pressure can be predicted from ambient pressure data as a function of temperature and solvent density (Manning 1998, Anderson et al. 2002). The characterization of the solvent by only two state variables implies the solvent is either pure (i.e., H2O) or that its properties are independent of chemistry, which is demonstrably false in the case of carbonic solvents (Harvey & Lemmon 2005). The density extrapolation method is supported in Perple_X (EoS 15). For additional equation of state and data file documentation refer to Perple_X_Thermodynamic_Data_Files.

Stable phases at:
                             T(K)     =  1000.00
                             P(bar)   =  30000.0
Phase Compositions (molar  proportions):

             wt %      vol %     mol %     mol        H2       Mg       Si       C        O2
 COH-Fluid   30.83     52.21     87.62     2.83      0.97704  0.00000  0.00000  0.02313  0.51148
 en          15.27      9.71      4.02    0.130      0.00000  2.00000  2.00000  0.00000  3.00000
 tan         52.17     36.90      7.28    0.235      1.00000  3.00000  4.00000  0.00000  6.00000
 mag          1.73      1.18      1.08    0.350E-01  0.00000  1.00000  0.00000  1.00000  1.50000

Phase speciation (molar proportions):

 COH-Fluid         CO2: 0.02305, CH4: 0.00009, H2O: 0.97687

...output abridged...

Chemical Potentials (J/mol):

      H2            Mg            Si            C             O2
    -51299.1      -391935.      -426749.      -2307.90      -428257.

...output abridged...

Back-calculated solute speciation in the solvent:

pH = 4.791+/-0.000; neutral_pH = 3.327; ionic_strength = 0.8067E-01; gamma/q^2 = 0.6667
dielectric constant = 17.41; solvent molar mass, g/mol = 18.6139; 
total solute molality = 0.4518

Solute endmember properties:

          charge   molality     mol_fraction      g,J/mol*     g0,J/mol**
SiO2,aq      0   0.253024       0.307103E-02      -855006       -843579
Mg(HCO3)     1   0.801347E-01   0.972622E-03     -1128344      -1103988
HCO3-       -1   0.582164E-01   0.706591E-03      -604275       -577261
Si2O4,aq     0   0.374689E-01   0.454773E-03     -1710012      -1682705
OH-         -1   0.205700E-01   0.249665E-03      -173710       -138047
HSiO3-      -1   0.186328E-02   0.226152E-04     -1028717       -973086
Mg(HSiO3     1   0.473607E-03   0.574833E-05     -1552786      -1485767
MgOH+        1   0.405550E-04   0.492230E-06      -697780       -610327
H+           1   0.242859E-04   0.294766E-06       -91716             0
CO3-2       -2   0.122701E-04   0.148926E-06      -534197       -426692
Mg+2         2   0.521524E-06   0.632990E-08      -545708       -411944
MgCO3,aq     0   0.386801E-06   0.469473E-08     -1036628       -913863
HO2-        -1   0.106986E-17   0.129852E-19      -387839        -40426

Solvent endmember properties:

         molality   mol_fraction  vol_fraction  v,cm3/mol*  g,J/mol*   g0,J/mol***
CO2        1.2381     0.02305        0.04369      32.742    -430564     -401095
CH4        0.0047     0.00009        0.00017      33.040    -104906      -35518
H2O       52.4805     0.97687        0.95614      16.903    -265427     -265235

*partial molar, **molal ref. state, ***molar ref. state.

Back-calculated fluid bulk composition:

             mol %           wt %
    H2      63.8068        10.2806
    Mg     0.978873E-01   0.190188
    Si     0.400894       0.900067
    C       1.67638        1.60962
    O2      34.0181        87.0195

Figure 1. Back-calculated solute speciation at 1000 K and 3 GPa obtained by running MEEMUM with the aq3_el.dat problem definition file.


Back-Calculated Electrolyte Speciation

Given the chemical potentials of the thermodynamic components of a system, Perple_X formulates a closed system of equations for the molal concentrations of all solute species by using the charge balance constraint. Perple_X can solve the speciation in terms of any arbitrarily chosen ion. The algorithm currently solves for speciation in terms of H+, should numerical error become significant the code can easily be modified for a more rational basis ion (e.g., OH-). A peculiarity of the contradiction between the HKF and Debye-Hueckel theory, is that the concentrations of electrolytic solute species are independent of the concentrations of neutral solute species. This independence reflects that the neutral solutes have no interactions with the electrolyte species and that the partial molar Gibbs energies of the solvent species are assumed to be identical to those in the solute-free solvent. 

Except for some minor technical details, the implementation of back-calculated speciation in Perple_X follows Galvez et al. (2015). The calculations are surprisingly robust (Figure 1), except that they are critically dependent on the availability of chemical potentials from a forward calculation. C-O-H-S-N solvents have a strong tendency to order for certain bulk compositions (pure H2O and the H2O-CO2 join) at low temperature. These degenerate compositions may prevent Perple_X from determining elemental chemical potentials, when this occurs the problem should be formulated in terms of oxide components. Three examples illustrate how the choice of components controls the feasibility of back-calculated speciation: aq1_ox.dat is the input file for H2O solvent; aq2_ox.dat is the input file for a binary H2O-CO2 solvent; and aq3_el.dat is the input file for a ternary C-O-H solvent. 

NOTE: for calculations with oxide components it is necessary to include O2 as a component so that Perple_X will include nominally reduced aqueous species, e.g., H+ = 1/2 H2O - 1/4 O2 - 1 electron. The amount of O2 may be specified as zero (i.e., to preclude fluid-rock redox processes), in which case its chemical potential is undefined (and irrelevant). An electron/proton component is not necessary in Perple_X because charge balance is implicit. 

Perple_X automatically does back-calculated speciation if: 1) solute species data is present in thermodynamic data file; 2) a recognized solvent is stable; and 3) the Perple_X option keyword aqueous_output is set to true (it is true by default). Recognized solvents are either pure H2O, a generic hybrid mixed-volatile solution model (solution model type 39), or an aqueous electrolyte solution model (solution model type 20, e.g., Aq_Solvent). At present, saturated phases, and saturated or mobile components cannot, or at least should not, be specified in back-calculated speciation problems.

Back-calculated speciation can be obtained directly with the program MEEMUM or by running VERTEX and then using WERAMI (Mode 1) to recover the speciation. In the latter case, speciation and some solvent properties can be recovered as a function of one or two variables with WERAMI (Mode 2-4, Property choice 4). The resulting spread-sheet format tab file can be analyzed with Perple_X_plot (MatLab), PSTABLE, or PYWERAMI. By default the concentrations of the 20 most abundant aqueous species in the chemical system of interest are output. The number of species output can be changed by setting the aqueous_species keyword to the desired number. Concentration units can be specified via the aq_solvent_composition and aq_solute_composition keywords.

In tests I found that the HKF formulation may substantially overestimate the concentrations of molecular volatiles if they are entered as neutral solute species (e.g., ETHANE,Aq; CO2,Aq; CO,Aq; H2,Aq etc.) rather than as solvent species. This problem is egregious for some high molecular weight organic compounds for which the HKF predicts implausible concentrations in water-rich carbonic fluids at 1200 K and 0.5 GPA. As I lack the expertise necessary to identify plausible predictions, I excluded all volatile neutral solute species, i.e.: CO2,aq; H2,aq; O2,aq; CO,aq; METHANE,; ACETIC-A; PROPANE,; HEXANE,A; ETHANE,A; ETHANOL,; ETHYLENE; FORMIC-A; BENZENE; METHANOL; LACTIC-A; LACTATE,; GLUTARAT; GLUTARIC; GLYCOLAT; GLYCOLIC; PROPANOL; PROPANOI; and TOLUENE, (the trailing comma in some names results because the Perple_X data conversion program DEW_2_VER truncates the SUPCRT names to 8 characters). Because, as noted above, neutral solute species have no influence on electrolyte speciation, these exclusions are cosmetic, i.e., I made them to avoid being distracted by the presence of these species in the output. Organic ions are a different matter because they affect electrolyte speciation. As in the case neutral organic species, I lack the expertise to identify plausible concentrations of the organic ions in the SUPCRT data base and therefore excluded all organic ions in my test calculations, i.e.: LACTATE,; GLUTARAT; GLYCOLAT; FORMATE,; PROPANOA; and ACETATE (glutarate and glycolate are particularly problematic). To facilitate identification of potentially anomalous species both MEEMUM and VERTEX print the aqueous species considered in a calculation to the user console at the beginning of a calculation.

The error on the back-calculated pH is half the mismatch between log(K_w) and ΔG_w/R/T. If this error is non-zero it is diagnostic of computational problems. 

Limitations of Back-Calculation

1) The phase proportions obtained by obtained optimization do not account for the solute-content of the solvent phase. Consequently, the method cannot be used for reactive transport calculations without introducing ad-hoc corrections.

2) If no forward model for solute speciation is included in the calculation, then an implicit assumption of back-calculation is that solute species have no significant influence on the partial molar free energies of the solvent species. While much is made of activity-composition phase relations and the limitations of the Debye-Huckel equation, the error introduced by this assumption is probably smaller than the error in the data if the total mole fraction of the solutes exceeds 10% (~5 molal solute concentration).

Relevant Perple_X options

aqueous_output - turns back-calculation on/off.

aqueous_species - controls the number of solute species output. 

aq_solvent_composition - controls solvent concentration units on output.

aq_solute_composition - controls solute concentration units on output.

Iteration    G(J/mol)           Na       Mg       Al       Si       K        Ca       Fe       O2       H2       C        S2
    0       -357971.7408       -237981  -421260  -513391  -460188  -278037  -514523  -114093  -389976   -50242   -28337   -43243
    1       -357980.2460       -242330  -420923  -511967  -460354  -280552  -512574  -109530  -389810   -50325   -30536   -47807
    2       -357996.6729       -243717  -421177  -511637  -460667  -281276  -512116  -107361  -389496   -50642   -31463   -49975
    3       -357997.6793       -243573  -420984  -511359  -460280  -281155  -511874  -108426  -389884   -50461   -31124   -48910
    4       -357997.7436       -243265  -420583  -510812  -459513  -280998  -511433  -110611  -390650   -50078   -30415   -46725
    5       -357997.7882       -243504  -420904  -511242  -460122  -281162  -511793  -108894  -390041   -50382   -30969   -48442
    6       -357997.7987       -243521  -420927  -511272  -460166  -281173  -511819  -108772  -389998   -50403   -31009   -48564
    7       -357997.8176       -243513  -420916  -511258  -460145  -281168  -511806  -108831  -390019   -50393   -30990   -48505
    8       -357997.8176       -243516  -420919  -511261  -460150  -281170  -511810  -108816  -390013   -50396   -30995   -48521
    9       -357997.8197       -243522  -420927  -511272  -460165  -281172  -511818  -108773  -389998   -50403   -31008   -48563
   10       -357997.8197       -243522  -420927  -511273  -460167  -281172  -511819  -108769  -389997   -50404   -31010   -48567

...output abridged...

Stable phases at:
                             T(K)     =  573.000
                             P(bar)   =  10000.0

Phase Compositions (molar  proportions):
                   wt %      vol %     mol %     mol        Na       Mg       Al       Si       K        Ca       Fe       O2       H2       C        S2
 Chl(HP)            0.96      0.90      0.13    0.175E-02  0.00000  4.98522  2.00181  2.99910  0.00000  0.00000  0.01388  9.00000  4.00000  0.00000  0.00000
 Mica(CF)          37.19     31.88      7.27    0.974E-01  0.59566  0.40362  2.19132  3.40434  0.40434  0.00000  0.00072  6.00000  1.00000  0.00000  0.00000
 Omph(HP)           5.91      4.41      2.13    0.286E-01  0.50110  0.48673  0.49517  1.98806  0.00000  0.49890  0.03004  3.00000  0.00000  0.00000  0.00000
 COH-Fluid          6.02     18.84     35.15    0.471      0.01194  0.00004  0.00000  0.00039  0.00976  0.00027  0.00000  0.33891  0.63165  0.00194  0.00509
 law                0.44      0.35      0.11    0.141E-02  0.00000  0.00000  2.00000  2.00000  0.00000  1.00000  0.00000  5.00000  2.00000  0.00000  0.00000
 q                 34.16     31.58     43.12    0.578      0.00000  0.00000  0.00000  1.00000  0.00000  0.00000  0.00000  1.00000  0.00000  0.00000  0.00000
 arag               7.78      6.60      5.90    0.791E-01  0.00000  0.00000  0.00000  0.00000  0.00000  1.00000  0.00000  1.50000  0.00000  1.00000  0.00000
 pyr                6.23      4.21      5.37    0.720E-01  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000  1.00000  0.00000  0.00000  0.00000  1.00000
 any                1.31      1.24      0.83    0.111E-01  0.00000  0.00000  0.00000  0.00000  0.00000  1.00000  0.00000  2.00000  0.00000  0.00000  0.50000

Phase speciation (molar proportions):

 Chl(HP)           daph: 0.00278, ames: 0.00090, clin: 0.99632
 Mica(CF)          cel: 0.40362, fcel: 0.00072, mu: 0.00000, pa: 0.59566
 Omph(HP)          di: 0.01543, jd: 0.00000, cats: 0.01194, acm: 0.02980, hed: 0.00023, omph: 0.94259
 COH-Fluid         CO2: 0.50632E-04, CH4: 0.0000, H2S: 0.48998E-05, H2O: 0.96438, FeC4H6O4: 0.11083E-27
                   FeC2H3O2: 0.18554E-17, NaCO3-: 0.32332E-03, CO3-2: 0.18677E-04, Al+3: 0.52715E-19, AlO2-: 0.31079E-04
                   Ca(HCO3): 0.97386E-03, Ca(OH)+: 0.13435E-02, Ca+2: 0.71635E-04, CaCO3,aq: 0.12492E-03, CaSO4,aq: 0.21668E-01
                   Fe+2: 0.61665E-10, Fe+3: 0.79961E-24, H+: 0.28478E-05, HAlO2,aq: 0.10937E-05, HCO3-: 0.46208E-01
                   HO2-: 0.22615E-28, HS-: 0.43784E-02, HSiO3-: 0.42669E-02, HSO3-: 0.96490E-05, HSO4-: 0.50232E-02
                   HSO5-: 0.64519E-29, K+: 0.23094E-01, KOH,aq: 0.85410E-05, KSO4-: 0.83638, Mg(HCO3): 0.68585E-04
                   Mg(HSiO3: 0.41143E-05, Mg+2: 0.53286E-06, MgCO3,aq: 0.12343E-05, MgOH+: 0.60757E-06, MgSO4,aq: 0.37023E-02
                   Na+: 0.92351, NaHCO3,a: 0.11980, NaHSiO3,: 0.67421E-02, NaOH,aq: 0.63221E-03, OH-: 0.20955E-02
                   S2-2: 0.54974E-07, S2O3-2: 0.54178E-05, S2O4-2: 0.26290E-15, S2O5-2: 0.14360E-15, S2O6-2: 0.46964E-17
                   S2O8-2: 0.40904E-38, S3-2: 0.34413E-08, S3O6-2: 0.70260E-20, S4-2: 0.16400E-09, S4O6-2: 0.44942E-13
                   S5-2: 0.59742E-11, S5O6-2: 0.61571E-24, Si2O4,aq: 0.20937E-02, SiO2,aq: 0.18948E-01, SO3-2: 0.38529E-07
                   SO4-2: 0.25192E-01, ionic_st: 0.97443, tot_mola: 57.552, solv_mas: 0.18016E-01, err_lgKw: -0.15752E-01*

 *molar solvent, molal solute concentrations.

...output abridged...

Back-calculated solute speciation in the solvent:

pH = 5.768+/-0.000; neutral_pH = 4.334; ionic_strength = 0.9744; gamma/q^2 = 0.5993
dielectric constant = 36.55; solvent molar mass, g/mol = 18.0164; total solute molality = 2.047

Solute endmember properties:
          charge   molality     mol_fraction   mol_fraction      g,J/mol*     g0,J/mol**
Na+          1   0.923511       0.108573E-01   0.160466E-01      -281592       -278774
KSO4-       -1   0.836379       0.983290E-02   0.145326E-01     -1047381      -1044090
NaHCO3,a     0   0.119795       0.140838E-02   0.208152E-02      -884731       -874622
HCO3-       -1   0.462075E-01   0.543240E-03   0.802886E-03      -603138       -586051
SO4-2       -2   0.251918E-01   0.296168E-03   0.437724E-03      -742060       -714765
K+           1   0.230945E-01   0.271511E-03   0.401282E-03      -319243       -298851
CaSO4,aq     0   0.216677E-01   0.254736E-03   0.376490E-03     -1316098      -1297842
SiO2,aq      0   0.189485E-01   0.222768E-03   0.329242E-03      -850164       -831269
NaHSiO3,     0   0.674213E-02   0.792640E-04   0.117149E-03     -1313888      -1290070
HSO4-       -1   0.502316E-02   0.590548E-04   0.872806E-04      -791410       -763751
HS-         -1   0.437839E-02   0.514746E-04   0.760774E-04       -11416         16897
HSiO3-      -1   0.426687E-02   0.501636E-04   0.741397E-04     -1032295      -1003858
MgSO4,aq     0   0.370228E-02   0.435260E-04   0.643296E-04     -1225206      -1198532
OH-         -1   0.209545E-02   0.246352E-04   0.364098E-04      -182130       -150306
Si2O4,aq     0   0.209373E-02   0.246150E-04   0.363800E-04     -1700328      -1670939
Ca(OH)+      1   0.134351E-02   0.157950E-04   0.233443E-04      -770090       -736148
Ca(HCO3)     1   0.973856E-03   0.114492E-04   0.169214E-04     -1191098      -1155623
NaOH,aq      0   0.632210E-03   0.743259E-05   0.109851E-04      -463723       -428629
NaCO3-      -1   0.323318E-03   0.380109E-05   0.561786E-05      -821458       -780730
HSO3-       -1   0.964897E-05   0.113438E-06   0.167657E-06      -596412       -538953

Solvent endmember properties:

         molality   mol_fraction  opt_mol_frac  vol_fraction# v,cm3/mol*  g,J/mol*   g0,J/mol***
CO2        0.0028     0.00005        0.00005        0.00011      35.526    -421007     -388964
CH4        0.0000     0.00000        0.00000        0.00000      35.898    -131819      -35138
H2S        0.0003     0.00000        0.00000        0.00001      35.237     -74688      -27157
H2O       53.5281     0.96438        0.96438        0.99988      17.233    -245403     -245230

#normalized, *partial molar, **molal ref. state, ***molar ref. state.

Back-calculated vs optimized fluid bulk composition:

                          optimized                     optimized
             mol %          mol %          wt %           wt %
    Na      1.23561        1.23561        2.18218        2.18218
    Mg     0.444086E-02   0.444086E-02   0.829134E-02   0.829134E-02
    Al     0.378240E-04   0.378240E-04   0.783970E-04   0.783970E-04
    Si     0.401474E-01   0.401474E-01   0.866170E-01   0.866170E-01
    K       1.01045        1.01045        3.03490        3.03490
    Ca     0.284291E-01   0.284291E-01   0.875257E-01   0.875257E-01
    Fe     0.724968E-10   0.724968E-10   0.311007E-09   0.311007E-09
    O2      33.9103        33.9103        83.3564        83.3564
    H2      63.0433        63.0433        9.76097        9.76097
    C      0.200242       0.200242       0.184759       0.184759
    S2     0.527065       0.527065        1.29828        1.29828

Figure 2. Lagged forward-calculated solute speciation at 573 K and 1 GPa obtained by running MEEMUM with the gloss0.dat problem definition file. The output of the Gibbs energy and chemical potentials after each iteration indicates well resolved solution. In contrast to the simple back-calculated model results in Figure 1, the composition of the solvent phase (COH-Fluid) accounts for solute species. In addition to the complete speciation in the "Phase speciation" section of the output, the lagged ionic strength (ionic_st), total molality (tot_mola), solvent mass (solv_mas), and the error in the base 10 log of the water dissociation constant (err_lgKw) are output at the end of the list of species. The latter parameter is indicative of the quality of the solution. The comparison of the back-calculated and lagged (optimized) speciation at the end of the output also demonstrates convergence.


Lagged Forward-Calculated Electrolyte Speciation


The lagged forward-calculated speciation algorithm exploits the iterative nature of adaptive optimization calculations in Perple_X. After each iteration, the chemical potentials are used to formulate expressions for the solute concentrations in terms of the solvent concentrations and partial molar Gibbs energies. These expressions together with an equation of state for the solvent species and the charge balance and closure constraints comprise a closed system of equations which is solved iteratively for the concentrations of both solute and solvent species. The solution becomes exact in the limit that the chemical potentials of the system do not change significantly during successive iterations. The degree to which this limit has been realized can be assessed by setting the output_iteration_G keyword to True and/or by simultaneously evaluating the back-calculated speciation by setting the aqueous_output keyword to True. The difference between the back- and lagged-speciation is attributable to the use of lagged chemical potentials to compute the latter, i.e., the back-calculated composition is a more precise solution, but the composition from the lagged solution is used to compute phase proportions and is therefore to be preferred for reactive transport modelling. 


Perple_X automatically does lagged speciation if: 1) solute species data is present in the thermodynamic data file; 2) a recognized solvent solution model is present; and 3) the Perple_X option keyword lagged_aq_speciation is set to true (it is true by default). Currently the only recognized solvent solution model is the generic hybrid mixed-volatile solution model (solution model type 39). In contrast to other solution models, generic hybrid solution models may have only one endmember (e.g., the WADDAH solution model) so that lagged speciation calculations can be done for a pure water solvent. At present, saturated phases, and saturated or mobile components cannot, or at least should not, be specified in lagged speciation problems.


Lagged speciation calculations will not converge if the equilibrium solute concentrations are extraordinarily high. Thus, at least with the DEW version of the SUPCRT data, the organic species mentioned re back-calculation must be excluded. If lagged speciation calculations do not converge, back-calculation with lagged speciation turned off (i.e., lagged_aq_speciation  set to False), will identify the problematic species. For the gloss0.dat example (Figure 2) it emerges that the lagged speciation calculations fail at high temperature because the equilibrium concentrations of alkali sulfates become excessive (I cannot judge whether these concentrations are plausible).


In addition to species concentrations, in lagged speciation calculations Perple_X includes four mock-species in the computational results (Figure 2): ionic_st (molal ionic strength), tot_mola (total molality), solv_mas (solvent molar mass, kg/mol), and err_lgKw (the error in the base 10 log of the water dissociation constant). The error in the base 10 log of the water dissociation constant is diagnostic of the precision of the lagged speciation calculation. Results involving a lagged speciation electrolyte solution models are analyzed using WERAMI or MEEMUM in the same manner as a non-electrolytic solution model.

Limitations of Lagged Forward-Calculation

1) Lagged speciation calculations will identify electrolytic phase separation between phases with different solvent compositions. For example, it will identify CH4- and CO2-rich solvent phases coexisting with a water-rich solvent phase, but it will not identify an electrolyte-rich H2O brine coexisting with an electrolyte poor H2O solvent (explicit solution models such as F(salt) should be used in such cases). This limitation is not truly algorithmic, but reflects the limited range of solute concentrations that can be modeled with the HKF formulation, at least as currently implemented in Perple_X.

2) As currently implemented, in lagged speciation calculations there is no accounting for solvent-solute interaction parameters other than Debye-Huckel activity coefficients. This has the consequence that the solvent activity coefficients are not a function of solute content. This behavior is likely to be a reasonable approximation if the total mole fraction of the solutes is less than 10% (~5 molal solute concentration). This limitation is not algorithmic. 

Relevant Perple_X options

species_output - controls whether phase speciation is output.

aqueous_output - controls whether the forward result is compared to the back-calculated speciation.

aq_solvent_composition - controls solvent concentration units on output.

aq_solute_composition - controls solute concentration units on output.

lagged_aq_speciation - controls whether lagged speciation is computed.

Stable phases at:
                             T(K)     =  1200.00
                             P(bar)   =  5000.00
Phase Compositions (molar  proportions):

                   wt %      vol %     mol %     mol        H2O      MgO      O2
 Aq_solvent        57.28     75.00     75.00     3.00      1.00000  0.00000  0.00000
 per               42.72      0.00     25.00     1.00      0.00000  1.00000  0.00000

Phase speciation (molar proportions):

 Aq_solvent        H2O: 1.0000, H+: 0.11286E-06, Mg+2: 0.55849E-09, MgOH+: 0.22506E-06, 
                   OH-: 0.33904E-06

...output abridged...

Chemical Potentials (J/mol):

      H2O           MgO           O2
    -340891.      -621905.        NaN

Back-calculated solute speciation in the solvent:

pH = 5.214+/-0.000; neutral_pH = 4.976; ionic_strength = 0.1550E-04; gamma/q^2 = 0.9833
dielectric constant = 7.060; solvent molar mass, g/mol = 18.0150; 
total solute molality = 0.3718E-04

Solute endmember properties:
          charge   molality     mol_fraction   mol_fraction      g,J/mol*     g0,J/mol**
OH-         -1   0.185885E-04   0.334871E-06   0.339037E-06      -221095       -112244
MgOH+        1   0.123738E-04   0.222913E-06   0.225059E-06      -741700       -628789
H+           1   0.620632E-05   0.111807E-06   0.112861E-06      -119795             0
Mg+2         2   0.707702E-08   0.127492E-09   0.558489E-09      -534941       -347030
HO2-        -1   0.576149E-08   0.103793E-09                     -221095        -31637

Solvent endmember properties:

         molality   mol_fraction  vol_fraction  v,cm3/mol*  g,J/mol*   g0,J/mol***
H2O       55.5093     1.00000        1.00000      28.502    -340890     -340890

*partial molar, **molal ref. state, ***molar ref. state.

Back-calculated vs optimized (forward) fluid bulk composition:

                          optimized                     optimized
             mol %          mol %          wt %           wt %
    H2O    100.0000       100.0000       100.0000        99.9999
    MgO    0.223041E-04   0.225617E-04   0.498997E-04   0.504761E-04
    O2     0.518966E-08   0.106624E-15   0.921809E-08   0.189390E-15

Figure 3. Forward and back-calculated solute speciation at 1200 K and 0.5 GPa obtained by running MEEMUM with the aq1_ox_forward.dat problem definition file.


Brute-Force Forward-Calculated Electrolyte Speciation

In brute-force forward-calculated electrolyte speciation models it is currently assumed that the activity coefficients of the solvent species are the same in the solute-bearing fluid as in a fluid consisting only of the solvent species. The concentrations of neutral and electrolyte solutes influence the partial molar free energies of the solvent species, unlike the back-calculated case, by affecting the concentrations of the solvent species.

Brute-force forward-calculation of electrolyte speciation involves specification of a solution model for the electrolytic fluid. Electrolytic solution models are identified as solution model type 20. Three examples, Aq_solvent, Aq_solven0, and Aq_Solven1 are present in solution_model.dat (input files for the forward solution models that correspond to the three back-calculated examples are included with those examples, the input files were run with MEEMUM at 1200 K and 5000 bar). The conditions for the use of these models are the same as for generic hybrid mixed-volatile solution models. In contrast to other solution model types in Perple_X, the initial resolution of the model is not controlled by initial resolution keyword; rather, the initial resolution is taken as specified in the solution model subdivision scheme. If it is desired to resolve concentrations that span many orders of magnitude the a non-linear subdivision scheme is recommended and the stretch factor should be adjusted to obtain the desired density of compositions. The lowest concentration (as a mole fraction) that Perple_X can resolve is final_resolution * speciation_factor. It is futile to include species present in concentrations below this limit.

Results involving brute-force forward-calculated electrolyte solution models are analyzed using WERAMI or MEEMUM in the same manner as a non-electrolytic solution model.

The accuracy of forward speciation (i.e., "optimized") calculations can be assessed by comparison with the back-calculated speciation and bulk composition (Figure 3). 

Limitations of Brute-Force Forward-Calcuation

1) Brute-force calculation is computationally intensive and is impractical for >> 10 species. Selecting the species to be included in the model requires expertise (which may be gleaned by doing a preliminary calculation with back-calculated speciation). Additionally, it may be necessary to adjust the default resolution and subdivision parameters to obtain precise solutions. In the future, an adaptive algorithm will be introduced to identify the dominant species as a function of physico-chemical conditions. 

Relevant Perple_X options

species_output - controls whether phase speciation is output. NOTE: the species concentrations output in brute-force calculation is in molar molar units regardless of the aq_solvent_composition and aq_solute_composition keywords.

aqueous_output - controls whether the forward result is compared to the back-calculated speciation.

stretch_factor - controls the assymetry of non-linear subdivision models, the subdivision is quasi-logarithmic in the limit that the stretch factor goes to zero.