Perple_X Electrolyte Speciation at Elevated Pressure and Temperature

Oct 19, 2017: results presented here were updated with Perple_X 6.8.0. 

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


Aqueous Speciation During Mantle Wedge Decompression Melting (Closed System Model)


Subduction Zone Mass Transfer during Devolatilization and Melting (Two-dimensional Phase Fractionation)

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

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

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) Back-calculation implicitly assumes the solute-bearing phase is in equilibrium with the pure solvent (an assumption that is only rigorous in the infinitely dilute limit). This has the consequence that in the limit of small, but finite, solute concentrations back-calculation always underestimates solute-content. 

2) Back-calculation implicitly assumes the mass of the solute-bearing phase is infinitesimal and therefore cannot affect the solid/melt phase proportions. However, the real mass of the solvent may be such that, if the solvent phase were to have the back-calculated solute load, one or more chemical components would be completely depleted from the solid/melt phase assemblage. E.g., if Gibbs energy minimization predicts stability of 1 mol of H2O (solvent), and back-calculation predicts that this fluid should contain 0.1 mol of Na2O, but the solid assemblage contains only 0.01 mol of Na2O, then the results are inconsistent and may be in substantial error. 

3) 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.

4) 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

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

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.


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. 


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

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) Lagged speciation calculations become invalid if the composition of the impure solvent phase moves outside of the composition space bounded by phase compositions identified in the previous iteration. This is inescapable if a solute component (e.g., K2O, CaO, or Na2O) is dissolved entirely in the solvent phase; it may also occur at high solute concentrations (~5-10 m). When this occurs Perple_X reports that the optimization failed although technically the optimization is possible, though numerically unstable. The results of such optimizations can be obtained by specifying the aq_bad_results option. Characteristically, such optimizations identify both the pure and impure solvent as stable coexisting phases. 

3) 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.

aq_bad_results - allows output of results once a solute component is completely depleted from the solid/melt assemblage (not recommended, see limitations).

aq_vapor_epsilon - allows specification of a threshold for the dielectric constant of the solvent, below this threshold the solvent is assumed to be vapor and incapable of dissolving solutes.

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

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.

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.

Aqueous Speciation During Mantle Wedge Decompression Melting (Closed System Model)


Explanatory text...