Perple_X Electrolyte Speciation at Elevated Pressure and Temperature

Corrections/modifications relevant to this page:


March 15, 2020, Deep Earth Water (DEW) data update: The Huang et al. (2019, Geochem Persp Lett 12:16) data for aqueous Cr-bearing species has been added to DEW19HP622ver_elements.dat. Courtesy of Vincent van Hinsberg (McGill).


February 18, 2020, DEW_2_VER parsing error for aqueous Glutamine (C5H10N2O3) and Isobutane (C4H10): a parsing error caused the DEW_2_VER data conversion program to reject species that included a 2-digit stoichiometric factor with zero as the second digit. For the 20 elements specified in the DEW17 and DEW19 data files, this error caused the omission of aqueous Glutamine and Isobutane. These species have been added to DEW17HP622ver_elements.dat and DEW19HP622ver_elements.dat. Correction courtesy of Andrea Maffeis (Torino).


December 11, 2019, DEW model update, CaCl2(aq) stoichiometric error: The Huang & Sverjensky (2019, GCA 254:192-230) update to the DEW database has been formatted as DEW19HP622ver_elements.dat. The DEW_2_VER program, which converts the DEW spreadsheet format to Perple_X format, parses the spreadsheet species "symbols" for stoichiometry. As of this date, the symbol provided for CaCl2,aq has been CaCl(0); consequently the stoichiometry of this species was incorrect in all previous Perple_X versions of the data. Additionally, no symbol has been provided for Glycinate,aq with the consequence that this species was, and is, not included in previous versions of the data. Refer to Section 1e of README_DEW_2_VER for additional information. Update courtesy of Andrea Maffeis (Torino/ETH).


November 18, 2019, ACETIC-ACID stoichiometry in DEW17HP622ver_elements.dat was not corrected as indicated in the Nov 6 update. Correction courtesy of Mohit Melwani Deswani (JPL, Cal Tech).


November 6, 2019: the DEW_2_VER program, which converts the DEW ( spreadsheet formatted data to Perple_X thermodynamic data file format, incorrectly parsed species formulae if: 1) the formula was terminated by two successive elements without stoichiometric coefficients or a charge or group indicator (e.g., CH3COOH); 2) a stoichiometric coefficient was > 9 (e.g. C6H14); or 3) a final stoichiometric coefficient was followed by a + or - sign (e.g., C5H7O4-). The affected species are:


Species: incorrect -> correct

ACETIC-ACID,AQ C(2)H2(1.5)O2(1) -> C(2)H2(2)O2(1) [H2O(2)CO2(2)O2(-2)]

ETHANOL,AQ C(2)H2(2.5)O2(.5) -> C(2)H2(3)O2(.5) [H2O(3)CO2(2)O2(-3)]

FORMIC-ACID,AQ C(1)H2(.5)O2(1) -> C(1)H2(1)O2(1) [H2O(1)CO2(1)O2(-.5)]

GLUTARAT C(5)H2(3.5)O2(.5) -> C(5)H2(3.5)O2(2) [H2O(3.5)CO2(5)O2(-4.75)]

GLYCOLAT C(2)H2(1.5)O2(.5) -> C(2)H2(1.5)O2(1.5) [H2O(1.5)CO2(2)O2(-1.25)]

HEXANE,AQ C(6)H2(.5) -> C(6)H2(7) [H2O(7)CO2(6)O2(-9.5)]

LACTATE, C(3)H2(2.5)O2(.5) -> C(3)H2(2.5)O2(1.5) [H2O(2.5)CO2(3)O2(-2.75)]

METHANOL,AQ C(1)H2(1.5)O2(.5) -> C(1)H2(2)O2(.5) [H2O(2)CO2(1)O2(-1.5)]

PROPANOL,AQ C(3)H2(3.5)O2(.5) -> C(3)H2(4)O2(.5) [H2O(4)CO2(3)O2(-4.5)]


These species have been corrected in DEW17HP622ver_elements.dat, DEW17HP622ver_oxides.dat, DEW13HP622ver_elements.dat, DEW13HP622ver_oxides.dat, DEW13HP02ver_elements.dat, and DEW13HP02ver_oxides.dat; the DEW_2_VER program has only been corrected in the 6.8.8 beta version of Perple_X. It is probable that this error was responsible for erratic behavior in calculations with organic carbon solute species. I apologize both for the error and for falsely attributing the error to the primary DEW data. Correction courtesy of Mohit Melwani Deswani (JPL, Cal Tech).


October 25, 2019: the lagged speciation error trap:


**warning ver103** pure and impure solvent phases coexist. To output result set aq_bad_result.


was triggered incorrectly during calculations when pure solvent was stable. Correction of this error extends the range of conditions at which the lagged speciation algorithm is feasible.


Jan 27, 2019: The DEW13 aqueous solute data-base was not converted to Perple_X format using the default options specified in the DEW spreadsheet, it is used here only for purposes of demonstration, the DEW17 aqueous solute data-base should be used for quantitative applications.


HKF g-function: as currently configured Perple_X will abort calculations with an H2O-solvent if the physical conditions are outside the range of validity for the HKF g-function indicated by Shock et al (1988); however, for multi-component solvents, Perple_X will zero the g-function at out-of-range conditions and continue calculations irrespective of the actual solvent composition.


Nov 21, 2018: The DEW17HP622ver_elements.dat and DEW17HP622ver_oxides.dat files generated on Nov 19 did not include the c2 heat capacity coefficient. The data files and DEW_2_ver have been corrected.


Nov 19, 2018: The May 19, 2017 revision of the DEW data ( is implemented in Perple_X as DEW17HP622ver_elements.dat and DEW17HP622ver_oxides.dat. The program for converting the DEW spreadsheet data format to Perple_X format, DEW_2_ver, is included with the Perple_X programs and documented at:


Sep 22, 2018: See Errata comment in Connolly & Galvez (2018), the comment also applies to earlier examples reproduced here. 


Aug 30, 2018: Updated to provide the files for Connolly & Galvez (2018) including examples computed with Perple_X 6.8.4.


Oct 19, 2017: results presented here (Figures 1-3) updated with Perple_X 6.8.0.



Forward and Back-Calculated Aqueous Electrolyte Speciation


Thermodynamic Details and Limitations


Molecular Volatiles: Solute or Solvent?


Thermodynamic Data Files for Electrolytic Fluids


Simple back-Calculated Electrolyte Speciation (Figure 1)


Relevant Options

Lagged Forward-Calculated Electrolyte Speciation (Figure 2)


Relevant Options

Files for calculations in Connolly & Galvez (2018)

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, 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 (aka lagged speciation, Connolly & Galvez 2018) is the best method. Simple 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

Both the extended HKF (Helgeson-Kirkham-Flowers) formulation (Perple_X EoS 16; Shock et al. 1992, Sverjensky et al. 2014) and the  Anderson density model (Perple_X EoS 15; Anderson et al. 1991, Holland & Powell 1998) for solute species are implemented in Perple_X; however, at present, only the HKF is linked to speciation calculations. The implementation of HKF/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, the spreadsheet uses, though other options are offered, 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 the Born coefficient (ω) of a solute is only a function of pressure (P) and temperature (T). 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 + (ω - ω0)(1/ε_solvent - 1/ε0_solvent) - (T - T0)ΔS0_solute [verify this equation]

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 standard state solute properties are influenced by solvent chemistry only through the dielectric constant, the conditions used to discriminate the different regions of applicability for the empirical expressions for the HKF correlation function used to compute Born coefficients (Shock et al., 1992) are H2O specific. As the correlation 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 empirical expressions (Harvey & Lemmon 2005; Harvey & Mountain 2017), 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. This expression yields a dielectric constant of 81.97 at 298.15 K and 1 bar for a density of ~0.997 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. Wolery 1990 notes that this approach violates the Gibbs-Duhem relation and 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. The inconsistencies may become significant at high solute concentrations.

Molecular volatiles: Solvent or Solute?

The thermodynamic properties of molecular volatiles are described more accurately by solvent species molecular EoS than they are by solute species thermodynamic extrapolations. However, the resolution on the molar fractions of solvent species is ~ final_resolution/speciation_factor (~ 1e-5 by default), whereas the resolution on solute species is essentially unlimited for both lagged and back-calculated speciation. Consequently, molecular volatiles should only be included as solvent species if their concentrations are likely to be significantly above final_resolution/speciation_factor. If the concentration of a solvent species is at this limit, then the concentrations of all species at or below the limit may become limited by the resolution of the solvent species.

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), hp622ver.dat (Holland & Powell 2011, rev. 2013), and the ba96ver.dat (Berman & Aranovich 1996) thermodynamic data bases. The merged files are DEW13HP02ver_elements.dat, DEW13HP02ver_oxides.dat, DEW13HP622ver_elements.dat, and DEW13HP622ver_oxides.dat, and DEW13BA96ver_oxides.dat, where the suffixes _elements and _oxides indicate whether the data is formulated in terms of element or oxide components. The Berman & Aranovich (1996) data is nominally consistent with the DEW data, but provides a relatively limited description of condensed phases. 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.

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


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.09 53.79 87.88 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.03447, CH4: 0.00008, H2: 0.00002, CO: 0.00000, H2O: 0.96543
...output abridged...
Chemical Potentials (J/mol):
 H2 Mg Si C O2
 -51276.9 -395063. -423268. 1173.80 -428493.
Back-calculated solute speciation in the solvent:
pH = 3.569+/-0.000; Delta_pH = 0.202; ionic_strength = 0.7146E-02; gamma/q^2 = 0.8559
dielectric constant = 16.81; solvent molar mass, g/mol = 18.9107; 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.366108E-02 0.686437E-04 -1104964 -1057027
HCO3- -1 0.349006E-02 0.654373E-04 -624506 -576171
HSiO3- -1 0.243878E-02 0.457261E-04 -1048947 -997632
MgOH+ 1 0.169718E-02 0.318215E-04 -677645 -623316
OH- -1 0.797071E-03 0.149448E-04 -197187 -136574
Mg+2 2 0.419445E-03 0.786443E-05 -480458 -410627
H+ 1 0.314842E-03 0.590316E-05 -68336 0
Mg(HSiO3 1 0.214319E-03 0.401840E-05 -1529405 -1457872
CO3-2 -2 0.200102E-06 0.375184E-08 -556169 -422751
MgCO3,aq 0 0.512487E-10 0.960892E-12 -1036628 -839623
HO2- -1 0.292313E-19 0.548076E-21 -411433 -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 -101380 -35518
H2 0.0010 0.00002 0.00002 18.858 -51276 26445
CO 0.0002 0.00000 0.00001 30.729 -213072 -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.736845E-02 0.141090E-01
 Si 0.645369 1.42796
 C 2.25587 2.13464
 O2 34.2929 86.4517


Simple 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 simple back-calculated speciation in Perple_X follows Galvez et al. (2015). The term "simple back-calculated" speciation is used here to distinguish the back-calculated speciation obtained as in Galvez et al. (2015) from the back-calculated speciation output during lagged speciation calculations (Connolly & Galvez 2018).

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 GLUTARAT (glutarate, C5H7O4-) GLYCOLAT (glycolate, C2H3O3-) and LACTATE (C3H5O3-) which cause convergence to an implausible result if included in calculations for carbonic systems. 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 (trailing commas or otherwise icomplete names result because the Perple_X data conversion program DEW_2_VER truncates the SUPCRT names to 8 characters). Molecular volatiles such as CO2, CH4, and H2S may be included as either solute or solvent species, but should not be present as both. Including molecular volatiles as solute species increases numerical stability and reduces computational cost, but prevents Perple_X from predicting solvent-phase immiscibility (e.g., the coexistence of a water-rich phase with CH4- or CO2-rich phases in carbonic systems).

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

Relevant Perple_X options

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. This example was done with aq_lagged_iterations = 2; thus, for each major iteration, the lagged speciation and chemical potentials of the system are recomputed 2 times. Use of aq_lagged_iterations generally does not result in a significant improvement to the solution, but may be helpful in low resolution calculations. 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), the error in the base 10 log of the water dissociation constant (err_lgKw), the pH, the difference between the pH and neutral pH (Delta_pH), the total solute molality (solute_m), and the solvent dielectric constant (epsilon) 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. Note: the back-calculated result reported for lagged calculations is calculated from the chemical potentials obtained from the impure solvent phase and therefore is not subject to the limitations mentioned for simple back-calculation.


Iteration G(J/mol) H2 C Si Al Fe Mg Ca Na K O2 S2
 0 -362014.3184 -32455 2507 -424252 -485527 -69861 -405900 -491465 -232928 -266660 -425912 -87475
 1 -362019.4133 -32569 2507 -424392 -485518 -69383 -406448 -491675 -232977 -266906 -425772 -87953
 1 -362019.8024 -32554 2507 -424361 -485510 -69508 -406281 -491630 -233041 -266859 -425802 -87828
 2 -362019.8032 -32553 2507 -424361 -485510 -69501 -406289 -491628 -233039 -266855 -425803 -87836
 2 -362019.8189 -32608 2507 -424470 -485537 -69564 -406318 -491792 -233027 -267177 -425694 -87772
 3 -362019.8222 -32600 2507 -424453 -485533 -69557 -406310 -491767 -233030 -267129 -425710 -87779
 3 -362019.8293 -32599 2507 -424451 -485532 -69540 -406327 -491763 -233026 -267116 -425713 -87796
 4 -362019.8293 -32600 2507 -424452 -485533 -69546 -406321 -491766 -233027 -267123 -425711 -87790
 4 -362019.8301 -32597 2507 -424447 -485532 -69545 -406318 -491758 -233029 -267107 -425716 -87791
 5 -362019.8299 -32597 2507 -424447 -485532 -69545 -406318 -491758 -233029 -267107 -425716 -87791
 5 -362019.8356 -32597 2507 -424448 -485532 -69546 -406318 -491759 -233029 -267110 -425716 -87790
 6 -362019.8356 -32597 2507 -424448 -485532 -69546 -406317 -491759 -233029 -267110 -425716 -87790
 6 -362019.8601 -32597 2507 -424446 -485531 -69545 -406317 -491757 -233030 -267106 -425717 -87791
...output abridged...
Stable phases at:
 T(K) = 573.000
 P(bar) = 10000.0
Phase Compositions (molar proportions):
 wt % vol % mol % mol H2 C Si Al Fe Mg Ca Na K O2 S2
 COH-Fluid 5.05 12.04 26.38 0.275 0.99290 0.00270 0.00202 0.00000 0.00000 0.00000 0.00001 0.00443 0.00022 0.50217 0.00009
 Mica(CF) 34.75 32.19 8.45 0.882E-01 1.00000 0.00000 3.47763 2.04473 0.11373 0.36390 0.00000 0.52237 0.47763 6.00000 0.00000
 Stlp 3.68 3.63 0.35 0.366E-02 6.25000 0.00000 8.00000 2.00000 2.54252 2.45748 0.00000 0.00000 0.50000 15.25000 0.00000
 Omph(GHP) 13.09 10.41 5.79 0.604E-01 0.00000 0.00000 2.00000 0.42677 0.22667 0.34655 0.49081 0.50919 0.00000 3.00000 0.00000
 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.64 6.12 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):
 COH-Fluid CO2: 0.45298E-04, CH4: 0.45298E-04, H2S: 0.0000, H2O: 0.99058, FeC4H6O4: 0.26531E-14
 FeC2H3O2: 0.34127E-11, NaCO3-: 0.13357E-02, CO3-2: 0.38938E-02, Al+3: 0.81673E-19, AlO2-: 0.11490E-03
 Ca(HCO3): 0.17634E-03, Ca(OH)+: 0.15771E-04, Ca+2: 0.58796E-04, CaCO3,aq: 0.12484E-03, CaSO4,aq: 0.11425E-07
 Fe+2: 0.30274E-08, Fe+3: 0.15452E-21, H+: 0.51578E-06, HAlO2,aq: 0.87282E-06, HCO3-: 0.13340
 HO2-: 0.20913E-29, HS-: 0.89052E-02, HSiO3-: 0.85012E-01, HSO3-: 0.39968E-09, HSO4-: 0.18151E-08
 HSO5-: 0.35642E-37, K+: 0.12351E-01, KOH,aq: 0.51236E-04, KSO4-: 0.80588E-08, Mg(HCO3): 0.14026E-06
 Mg(HSiO3: 0.44222E-07, Mg+2: 0.13930E-06, MgCO3,aq: 0.11397E-08, MgOH+: 0.15213E-06, MgSO4,aq: 0.98381E-11
 Na+: 0.23349, NaHCO3,a: 0.72751E-02, NaHSiO3,: 0.52021E-02, NaOH,aq: 0.87291E-03, OH-: 0.95995E-02
 S2-2: 0.19353E-07, S2O3-2: 0.21926E-10, S2O4-2: 0.13241E-22, S2O5-2: 0.14639E-24, S2O6-2: 0.88111E-28
 S2O8-2: 0.69529E-53, S3-2: 0.12550E-10, S3O6-2: 0.18380E-32, S4-2: 0.62297E-14, S4O6-2: 0.45653E-28
 S5-2: 0.23632E-17, S5O6-2: 0.30092E-40, Si2O4,aq: 0.20937E-02, SiO2,aq: 0.18946E-01, SO3-2: 0.21376E-10
 SO4-2: 0.17753E-05, ionic_st: 0.25011, tot_mola: 56.029, solv_mas: 0.18016E-01, err_lgKw: -0.40811E-02
 pH: 6.4719, Delta_pH: 2.1369, solute_m: 0.52292, epsilon: 36.545*
 Mica(CF) cel: 0.36390, fcel: 0.11373, mu: 0.00000, pa: 0.52237
 Stlp mstp: 0.49150, fstp: 0.50850
 Omph(GHP) di: -0.25412, jd: -0.09859, acm: 0.02052, hed: 0.00705, om: 0.92694, cfm: 0.27440, jac: 0.12379
*molar solvent, molal solute concentrations.
...output abridged...
Back-calculated solute speciation in the solvent:
pH = 6.472+/-0.000; Delta_pH = 2.135; ionic_strength = 0.2501; gamma/q^2 = 0.6541
dielectric constant = 36.55; solvent molar mass, g/mol = 18.0161; solute molality = 0.5228
Solute endmember properties:
 charge molality mol_fraction mol_fraction g,J/mol* g0,J/mol**
Na+ 1 0.233433 0.416632E-02 0.416736E-02 -287727 -278774
HCO3- -1 0.133353 0.238009E-02 0.238085E-02 -597669 -586049
HSiO3- -1 0.850083E-01 0.151723E-02 0.151729E-02 -1024624 -1010859
SiO2,aq 0 0.189461E-01 0.338151E-03 0.338150E-03 -850164 -831269
K+ 1 0.123596E-01 0.220594E-03 0.220447E-03 -321804 -298851
OH- -1 0.959905E-02 0.171324E-03 0.171331E-03 -174460 -150303
HS- -1 0.890483E-02 0.158933E-03 0.158940E-03 -5497 19016
NaHCO3,a 0 0.727120E-02 0.129776E-03 0.129845E-03 -885397 -861939
NaHSiO3, 0 0.520078E-02 0.928234E-04 0.928466E-04 -1312352 -1287297
CO3-2 -2 0.389249E-02 0.694731E-04 0.694968E-04 -526674 -492150
Si2O4,aq 0 0.209366E-02 0.373677E-04 0.373676E-04 -1700328 -1670939
NaCO3- -1 0.133503E-02 0.238276E-04 0.238396E-04 -814401 -780846
NaOH,aq 0 0.872686E-03 0.155757E-04 0.155796E-04 -462187 -428629
Ca(HCO3) 1 0.176324E-03 0.314703E-05 0.314724E-05 -1198821 -1155622
CaCO3,aq 0 0.124837E-03 0.222810E-05 0.222809E-05 -1127825 -1085003
AlO2- -1 0.114892E-03 0.205058E-05 0.205075E-05 -856551 -811311
Ca+2 2 0.588051E-04 0.104955E-05 0.104938E-05 -601151 -546654
KOH,aq 0 0.512700E-04 0.915066E-06 0.914457E-06 -496264 -449202
Ca(OH)+ 1 0.157739E-04 0.281533E-06 0.281473E-06 -775612 -720912
HSO3- -1 0.399499E-09 0.713025E-11 0.713352E-11 -644073 -538951
Solvent endmember properties:
 molality mol_fraction opt_mol_frac vol_fraction# v,cm3/mol* g,J/mol* g0,J/mol***
CO2 0.0025 0.00005 0.00005 0.00009 35.526 -423209 -388964
CH4 0.0025 0.00005 0.00005 0.00010 35.898 -62687 -35054
H2S 0.0000 0.00000 0.00000 0.00000 35.237 -76493 -27157
H2O 55.5008 0.99058 0.99058 0.99981 17.233 -245456 -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.9941 65.9939 10.9518 10.9518
 C 0.179400 0.179457 0.177419 0.177475
 Si 0.134456 0.134462 0.310923 0.310935
 Al 0.137329E-03 0.137340E-03 0.305085E-03 0.305109E-03
 Fe 0.359553E-08 0.359543E-08 0.165326E-07 0.165321E-07
 Mg 0.565796E-06 0.565921E-06 0.113226E-05 0.113250E-05
 Ca 0.445747E-03 0.445744E-03 0.147091E-02 0.147090E-02
 Na 0.294331 0.294408 0.557147 0.557291
 K 0.147228E-01 0.147130E-01 0.473964E-01 0.473648E-01
 O2 33.3772 33.3772 87.9396 87.9394

S2 0.528289E-02 0.528312E-02 0.139477E-01 0.139482E-01

Lagged Forward-Calculated Electrolyte Speciation


The lagged forward-calculated speciation algorithm (aka, lagged speciation, Connolly & Galvez 2018) 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. Note: the "back-calculated" result reported for lagged-speciation calculations corresponds to the speciation obtained from Eq 12 of Connolly and Galvez (2018) at h = n, i.e., once iteration has ceased after n iterations; this differs from the "simple back-calculated" result reported if the aq_lagged_speciation is set to false. The latter corresponds to the speciation obtained from Eq 12 at h = 0.


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 eight mock-species in the computational results (Figure 2): ionic_st (molal ionic strength), tot_mola (total molality), solv_mas (solvent molar mass, kg/mol), err_lgKw (the error in the base 10 log of the water dissociation constant), pH, Delta_pH (the difference between the pH and neutral pH), solute_m (total solute molality), and epsilon (solvent dielectric constant). The err_lgKw parameter is diagnostic of the precision of the 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

Relevant Perple_X options

Files for the calculations presented in Connolly & Galvez (2018)

The 6.8.1 version of Perple_X used for the calculations shown in Connolly & Galvez (2018) had an error that caused the solvent dielectric constant to be computed by the MRK EoS regardless of the molecular EoS specified by the user. The error was not corrected because its consequences are minor, but the error has the consequence that the current (6.8.4) version does not exactly reproduce all the published figures (Figure 2 and A.1 are exceptions because these figures were calculated for the revised version of the paper). The COHS-solvent model used in the paper, while perhaps useful for purposes of illustration, is not the optimal practical configuration. In most of the examples below this model has been replaced by a CO2-H2O-solvent model. The distinction between these models is clarified below and has no significant impact on computational results. The files listed below include some output for testing purposes and/or comparison. The collected files are in

Notes on the calculations and Matlab plotting*:

*Similar plots can be generated with PYWERAMI or the Perple_X program PSTABLE.

Files common to multiple calculations:

Figures 2 and A.1. K/Si-solubility, Closed system, H2O-Solvent Model, Lagged Speciation (Figures 2, A.1):

Lagged Speciation (Figures 4a, 5, 6ab), GLOSS Closed-System, CO2-H2O-Solvent Model:

Lagged Speciation (Figures 5c, 6c), GLOSS Open-System, CO2-H2O-Solvent Model:

Lagged Speciation (Figure 7), GLOSS Open-System + infiltration, CO2-H2O-Solvent Model:

Simple Back-Calculated Speciation (Figures 4b, 5a), GLOSS Closed-System, CO2-H2O-Solvent Model:

Lagged Speciation (Figures 4a, 5b, 6b), GLOSS Closed-System, H2O-Solvent Model:

Lagged Speciation (Figures 4a, 5, 6ab), GLOSS Closed-System, COHS-Solvent Model:

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


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

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