UPDATE, Oct 19, 2017: The calculations presented here were made with Perple_X 6.8.0. Prior to the 6.8.0 version, the lagged speciation algorithm suffered from numerical instability if a component was present only in the fluid phase, but close to saturation in the solid/melt phases. In such cases, the algorithm oscillated between the correct solution and a supersaturated solution. In 6.8.0, this source of instability has been eliminated by increasing the number of metastable phases monitored during adaptive optimization.
Forward and Back-Calculated Aqueous Electrolyte Speciation
Thermodynamic Details and Limitations
Thermodynamic Data Files for Electrolytic Fluids
Back-Calculated Electrolyte Speciation (Figure 1)
Lagged Forward-Calculated Electrolyte Speciation (Figure 2)
Brute-Force Forward-Calculated Electrolyte Speciation (Figure 3)
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, 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. 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.
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  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), 
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) 
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.79865·10-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 ). 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.
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 32.10 53.18 87.87 2.90 0.96560 0.00000 0.00000 0.03456 0.51719 en 23.43 14.85 6.04 0.199 0.00000 2.00000 2.00000 0.00000 3.00000 ta 44.46 31.35 6.07 0.200 1.00000 3.00000 4.00000 0.00000 6.00000 mag 0.02 0.01 0.01 0.308E-03 0.00000 1.00000 0.00000 1.00000 1.50000 Phase speciation (molar proportions): COH-Fluid CO2: 0.03448, CH4: 0.00008, H2: 0.00002, CO: 0.00000, H2O: 0.96543 ...output abridged... Chemical Potentials (J/mol): H2 Mg Si C O2 -51274.5 -395061. -423263. 1178.67 -428498. ...output abridged... Back-calculated solute speciation in the solvent: pH = 3.579+/-0.000; neutral_pH = 3.368; ionic_strength = 0.7118E-02; gamma/q^2 = 0.8562 dielectric constant = 16.81; solvent molar mass, g/mol = 18.9107; total solute molality = 0.4542 Solute endmember properties: charge molality mol_fraction g,J/mol* g0,J/mol** SiO2,aq 0 0.360240 0.675435E-02 -851760 -843271 Si2O4,aq 0 0.809606E-01 0.151798E-02 -1703520 -1682620 Mg(HCO3) 1 0.358204E-02 0.671618E-04 -1105143 -1057027 HCO3- -1 0.356513E-02 0.668447E-04 -624326 -576171 HSiO3- -1 0.249123E-02 0.467095E-04 -1048768 -997632 MgOH+ 1 0.166054E-02 0.311345E-04 -677824 -623316 OH- -1 0.814214E-03 0.152662E-04 -197007 -136574 Mg+2 2 0.401311E-03 0.752442E-05 -480816 -410627 H+ 1 0.308045E-03 0.577571E-05 -68515 0 Mg(HSiO3 1 0.209692E-03 0.393165E-05 -1529585 -1457872 CO3-2 -2 0.208688E-06 0.391282E-08 -555811 -422751 MgCO3,aq 0 0.512487E-10 0.960892E-12 -1036628 -839623 HO2- -1 0.298513E-19 0.559700E-21 -411256 -36165 Normalized solvent endmember properties: molality mol_fraction vol_fraction v,cm3/mol* g,J/mol* g0,J/mol*** CO2 1.8228 0.03447 0.06465 32.729 -427318 -401095 CH4 0.0042 0.00008 0.00015 33.019 -101370 -35518 H2 0.0010 0.00002 0.00002 18.858 -51274 26445 CO 0.0002 0.00000 0.00001 30.729 -213070 -124463 H2O 51.0519 0.96543 0.93517 16.904 -265523 -265235 *partial molar, **molal ref. state, ***molar ref. state. Back-calculated fluid bulk composition: mol % wt % H2 62.7985 9.97159 Mg 0.719821E-02 0.137830E-01 Si 0.645428 1.42809 C 2.25586 2.13464 O2 34.2930 86.4519
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.
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).
aq_output - turns back-calculation on/off.
aq_oxide_components - allows speciation calculations with oxide components (not recommended).
aq_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) H2 C Si Al Fe Mg Ca Na K O2 S2 0 -362015.8596 -32514 2507 -424371 -485557 -69476 -406406 -491645 -232901 -266837 -425792 -87860 1 -362018.9189 -32618 2507 -424508 -485556 -69657 -406268 -491849 -233010 -267306 -425656 -87679 1 -362019.0833 -32572 2507 -424416 -485533 -69538 -406317 -491712 -233013 -267008 -425747 -87798 1 -362019.0860 -32572 2507 -424416 -485534 -69538 -406317 -491712 -233013 -267008 -425747 -87798 2 -362019.0904 -32572 2507 -424417 -485534 -69538 -406317 -491713 -233014 -267009 -425747 -87798 2 -362019.1026 -32591 2507 -424454 -485543 -69554 -406334 -491769 -233008 -267118 -425709 -87782 2 -362019.1045 -32591 2507 -424454 -485543 -69554 -406334 -491769 -233008 -267118 -425709 -87782 3 -362019.1045 -32591 2507 -424454 -485543 -69554 -406334 -491769 -233008 -267118 -425709 -87782 3 -362019.1053 -32588 2507 -424449 -485542 -69556 -406327 -491761 -233010 -267104 -425714 -87780 3 -362019.1051 -32588 2507 -424448 -485541 -69552 -406331 -491760 -233009 -267100 -425715 -87784 4 -362019.1050 -32588 2507 -424449 -485542 -69556 -406327 -491761 -233010 -267104 -425714 -87780 4 -362019.1052 -32587 2507 -424447 -485541 -69555 -406327 -491758 -233011 -267097 -425717 -87781 4 -362019.1050 -32587 2507 -424447 -485541 -69555 -406326 -491758 -233011 -267098 -425716 -87781 5 -362019.1051 -32587 2507 -424447 -485541 -69555 -406327 -491758 -233011 -267097 -425717 -87781 5 -362019.1311 -32587 2507 -424445 -485540 -69553 -406326 -491755 -233011 -267093 -425718 -87783 5 -362019.1310 -32587 2507 -424446 -485541 -69554 -406327 -491757 -233011 -267096 -425717 -87783 ...output abridged... Stable phases at: T(K) = 573.000 P(bar) = 10000.0 Phase speciation (molar proportions): wt % vol % mol % mol H2 C Si Al Fe Mg Ca Na K O2 S2 Mica(CF) 34.75 32.19 8.45 0.882E-01 1.00000 0.00000 3.47771 2.04457 0.11375 0.36396 0.00000 0.52229 0.47771 6.00000 0.00000 Stlp 3.67 3.62 0.35 0.366E-02 6.25000 0.00000 8.00000 2.00000 2.54231 2.45769 0.00000 0.00000 0.50000 15.25000 0.00000 Omph(GHP) 13.10 10.42 5.79 0.604E-01 0.00000 0.00000 2.00000 0.42678 0.22675 0.34647 0.49074 0.50926 0.00000 3.00000 0.00000 COH-Fluid 5.05 12.04 26.39 0.275 0.99293 0.00280 0.00206 0.00000 0.00000 0.00000 0.00001 0.00436 0.00022 0.50230 0.00008 law 3.26 2.83 0.99 0.103E-01 2.00000 0.00000 2.00000 2.00000 0.00000 0.00000 1.00000 0.00000 0.00000 5.00000 0.00000 q 29.98 30.15 47.66 0.497 0.00000 0.00000 1.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 1.00000 0.00000 arag 6.63 6.11 6.33 0.660E-01 0.00000 1.00000 0.00000 0.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.50000 0.00000 pyr 3.53 2.60 3.83 0.400E-01 0.00000 0.00000 0.00000 0.00000 1.00000 0.00000 0.00000 0.00000 0.00000 0.00000 1.00000 gph 0.03 0.03 0.21 0.221E-02 0.00000 1.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 Phase speciation (molar proportions): Mica(CF) cel: 0.36396, fcel: 0.11375, mu: 0.00000, pa: 0.52229 Stlp mstp: 0.49154, fstp: 0.50846 Omph(GHP) di: -0.25421, jd: -0.09857, acm: 0.02057, hed: 0.00704, om: 0.92689, cfm: 0.27446, jac: 0.12382 COH-Fluid CO2: 0.67945E-04, CH4: 0.67945E-04, H2S: 0.0000, SO2: 0.0000, H2: 0.0000 CO: 0.0000, H2O: 0.99051, FeC4H6O4: 0.26642E-14, FeC2H3O2: 0.33461E-11, NaCO3-: 0.13683E-02 CO3-2: 0.40603E-02, Al+3: 0.76718E-19, AlO2-: 0.11707E-03, Ca(HCO3): 0.17284E-03, Ca(OH)+: 0.15467E-04 Ca+2: 0.56463E-04, CaCO3,aq: 0.12473E-03, CaSO4,aq: 0.11431E-07, Fe+2: 0.29007E-08, Fe+3: 0.14522E-21 H+: 0.50570E-06, HAlO2,aq: 0.87156E-06, HCO3-: 0.13631, HO2-: 0.21373E-29, HS-: 0.91125E-02 HSiO3-: 0.86903E-01, HSO3-: 0.40874E-09, HSO4-: 0.18557E-08, HSO5-: 0.36435E-37, K+: 0.12139E-01 KOH,aq: 0.51452E-04, KSO4-: 0.82598E-08, Mg(HCO3): 0.13712E-06, Mg(HSiO3: 0.43248E-07, Mg+2: 0.13340E-06 MgCO3,aq: 0.11367E-08, MgOH+: 0.14878E-06, MgSO4,aq: 0.98166E-11, Na+: 0.22949, NaHCO3,a: 0.73048E-02 NaHSiO3,: 0.52260E-02, NaOH,aq: 0.87682E-03, OH-: 0.98122E-02, S2-2: 0.20238E-07, S2O3-2: 0.22912E-10 S2O4-2: 0.13834E-22, S2O5-2: 0.15291E-24, S2O6-2: 0.92015E-28, S2O8-2: 0.72588E-53, S3-2: 0.13133E-10 S3O6-2: 0.19209E-32, S4-2: 0.65246E-14, S4O6-2: 0.47754E-28, S5-2: 0.24770E-17, S5O6-2: 0.31495E-40 Si2O4,aq: 0.20936E-02, SiO2,aq: 0.18943E-01, SO3-2: 0.22319E-10, SO4-2: 0.18533E-05, ionic_st: 0.25096 tot_mola: 56.028, solv_mas: 0.18017E-01, err_lgKw: -0.32330E-02* *molar solvent, molal solute concentrations. ...output abridged... Back-calculated solute speciation in the solvent: pH = 6.481+/-0.000; neutral_pH = 4.337; ionic_strength = 0.2510; gamma/q^2 = 0.6538 dielectric constant = 36.54; solvent molar mass, g/mol = 18.0166; total solute molality = 0.5242 Solute endmember properties: optimized charge molality mol_fraction mol_fraction g,J/mol* g0,J/mol** Na+ 1 0.229522 0.409652E-02 0.409605E-02 -287809 -278773 HCO3- -1 0.136333 0.243329E-02 0.243289E-02 -597564 -586046 HSiO3- -1 0.869006E-01 0.155101E-02 0.155105E-02 -1024519 -1010855 SiO2,aq 0 0.189432E-01 0.338099E-03 0.338100E-03 -850164 -831268 K+ 1 0.121333E-01 0.216557E-03 0.216657E-03 -321894 -298851 OH- -1 0.981200E-02 0.175125E-03 0.175130E-03 -174354 -150300 HS- -1 0.911190E-02 0.162630E-03 0.162641E-03 -5387 19019 NaHCO3,a 0 0.730669E-02 0.130410E-03 0.130377E-03 -885374 -861939 NaHSiO3, 0 0.522635E-02 0.932802E-04 0.932735E-04 -1312328 -1287297 CO3-2 -2 0.406097E-02 0.724804E-04 0.724692E-04 -526472 -492141 Si2O4,aq 0 0.209357E-02 0.373662E-04 0.373662E-04 -1700328 -1670939 NaCO3- -1 0.136864E-02 0.244277E-04 0.244219E-04 -814282 -780843 NaOH,aq 0 0.876885E-03 0.156507E-04 0.156496E-04 -462164 -428628 Ca(HCO3) 1 0.172845E-03 0.308495E-05 0.308482E-05 -1198917 -1155620 CaCO3,aq 0 0.124729E-03 0.222617E-05 0.222617E-05 -1127825 -1084999 AlO2- -1 0.117082E-03 0.208969E-05 0.208954E-05 -856460 -811307 Ca+2 2 0.564585E-04 0.100768E-05 0.100776E-05 -601353 -546651 KOH,aq 0 0.514262E-04 0.917859E-06 0.918325E-06 -496248 -449201 Ca(OH)+ 1 0.154649E-04 0.276019E-06 0.276059E-06 -775708 -720911 HSO3- -1 0.408823E-09 0.729671E-11 0.729517E-11 -643964 -538949 Solvent endmember properties: molality mol_fraction opt_mol_frac vol_fraction# v,cm3/mol* g,J/mol* g0,J/mol*** CO2 0.0038 0.00007 0.00007 0.00014 35.526 -423209 -388964 CH4 0.0038 0.00007 0.00007 0.00014 35.898 -62667 -35054 H2S 0.0000 0.00000 0.00000 0.00000 35.237 -76479 -27157 SO2 0.0000 0.00000 0.00000 0.00000 43.347 -469609 -300964 H2 0.0000 0.00000 0.00000 0.00000 21.086 -32587 20791 CO 0.0000 0.00000 0.00000 0.00000 33.627 -210350 -124737 H2O 55.4966 0.99051 0.99051 0.99972 17.233 -245446 -245411 #normalized, *partial molar, **molal ref. state, ***molar ref. state. Back-calculated vs optimized fluid bulk composition: optimized optimized mol % mol % wt % wt % H2 65.9861 65.9861 10.9493 10.9493 C 0.186198 0.186168 0.184121 0.184092 Si 0.136708 0.136711 0.316097 0.316102 Al 0.139907E-03 0.139897E-03 0.310778E-03 0.310757E-03 Fe 0.344490E-08 0.344457E-08 0.158383E-07 0.158368E-07 Mg 0.550085E-06 0.550002E-06 0.110070E-05 0.110053E-05 Ca 0.438279E-03 0.438279E-03 0.144611E-02 0.144611E-02 Na 0.289768 0.289734 0.548449 0.548384 K 0.144525E-01 0.144593E-01 0.465214E-01 0.465431E-01 O2 33.3808 33.3808 87.9395 87.9395 S2 0.540502E-02 0.540539E-02 0.142685E-01 0.142695E-01
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.
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.
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.
species_output - controls whether phase speciation is output.
aq_output - controls whether the forward result is compared to the back-calculated speciation.
aq_oxide_components - allows speciation calculations with oxide components (not recommended).
aq_solvent_composition - controls solvent concentration units on output.
aq_solute_composition - controls solute concentration units on output.
aq_lagged_iterations - controls number of lagged speciation iterations made for each adaptive minimization iteration.
aq_lagged_speciation - controls whether lagged speciation is computed.
output_iteration_G - controls whether the free energy and chemical potentials of the system are output during iteration.
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.10522E-06, Mg+2: 0.44455E-09, MgOH+: 0.27080E-06, OH-: 0.37691E-06 ...output abridged... Chemical Potentials (J/mol): H2O MgO O2 -340564. -621905. NaN Back-calculated solute speciation in the solvent: pH = 5.300+/-0.000; neutral_pH = 4.969; ionic_strength = 0.2089E-04; gamma/q^2 = 0.9806 dielectric constant = 7.060; solvent molar mass, g/mol = 18.0150; total solute molality = 0.4173E-04 Solute endmember properties: optimized charge molality mol_fraction mol_fraction g,J/mol* g0,J/mol** OH- -1 0.234248E-04 0.421998E-06 0.376910E-06 -218814 -112244 MgOH+ 1 0.131595E-04 0.237069E-06 0.270805E-06 -743653 -631330 H+ 1 0.511649E-05 0.921735E-07 0.105216E-06 -121749 0 Mg+2 2 0.196915E-07 0.354742E-09 0.444549E-09 -524838 -347030 HO2- -1 0.706131E-08 0.127209E-09 -218814 -31359 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 -340563 -340563 *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.237423E-04 0.271250E-04 0.531174E-04 0.606852E-04 O2 0.231070E-05 0.194976E-14 0.410436E-05 0.346324E-14
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).
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.
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.
aq_output - controls whether the forward result is compared to the back-calculated speciation.
aq_oxide_components - allows back- and lagged-forward calculated speciation with oxide components (not recommended, necessary in this example only to obtain the back-calculated result for comparison).
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.