Abstract

Carbon is subducted to depths where metamorphism liberates water-bearing fluids. The C-bearing fluids facilitate partial melting of the upper mantle, generating magmas that may erupt as arc volcanics. Degassing of the magmas releases CO2 and other volatile species to the atmosphere. Over geological time, this process contributes to the composition of the atmosphere and planetary habitability. Here I summarize the background needed to carry out theoretical geochemical modelling of fluids and fluid–rock interactions from surficial conditions into the upper mantle. A description of the general criteria for predicting equilibrium and non-equilibrium chemical reactions is followed by a summary of how the thermodynamic activities of species are related to measurable concentrations through standard states and activity coefficients. Specific examples at ambient conditions involving dilute water are detailed. The concept of aqueous speciation and how it can be calculated arises from this discussion. Next, I discuss how to calculate standard Gibbs free energies and aqueous activity coefficients at elevated temperatures and pressures. The revised Helgeson–Kirkham–Flowers equations of state are summarized and the revised predictive correlations for the estimation of equation of state coefficients in the Deep Earth Water (DEW) model are presented. Finally, the DEW model is applied to the solubility and speciation of aqueous aluminium.

The motivation for developing an understanding of water–rock interactions in the deep Earth comes from long-standing suggestions that deep fluids link subducting plates and volatile fluxes to the atmosphere (Ulmer & Trommsdorff 1995; Manning 2004; Evans et al. 2012). Indeed, the fluxes of C, N and S to the atmosphere are thought to help maintain long-term planetary habitability and influence the redox state of the atmosphere. As an example, a schematic illustration showing the connections between the near-surface carbon reservoirs and the deep carbon cycle is shown in Figure 1 (Manning 2014). Carbon in a variety of forms is subducted to depths where metamorphism of the subducting materials liberates water-bearing fluids represented by the wiggly arrows in Figure 1. The C-bearing fluids are thought to facilitate partial melting of the upper mantle, generating magmas that may in part erupt as arc volcanics. Degassing of the magmas releases CO2 and other volatile species to the atmosphere. Over geological time, this process contributes to the composition of the atmosphere and planetary habitability.

But what do the slab fluids depicted in Figure 1 really contain? Numerous modelling studies of metamorphic dehydration and decarbonation have relied on a carbon–oxygen–hydrogen, or COH, fluid model that treats the fluid as a mixture of dissolved gas species, typically CO2, CO, CH4, H2 and H2O (Kerrick & Connolly 1998; Kerrick & Connolly 2001a, b; Hacker et al. 2003a,b; Hacker & Abers 2004; Gorman et al. 2006; Hacker 2008; Zhang & Duan 2009; van Keken et al. 2011). However, experimental studies of the solubilities of minerals (Manning 2013), the solubilities of rocks (Kessel et al. 2005a,b, 2015; Dvir et al. 2011; Dvir & Kessel 2017; Tsay et al. 2017; Tumiati et al. 2017; Tiraboschi et al. 2018) and of aqueous speciation (Martinez et al. 2004; Mysen 2010; Louvel et al. 2013; Mysen et al. 2013; Facq et al. 2014, 2016) now clearly show that the fluids leaving subducting slabs must contain all the rock-forming elements and trace elements in addition to the dissolved gases mentioned above. All the above evidence emphasizes the need for comprehensive models of aqueous speciation, solubility and chemical mass transfer under subduction-zone conditions.

Theoretical geochemical modelling of water–rock interactions has long been possible under crustal conditions. The foundations of theoretical aqueous geochemistry were laid by the pioneering work of Harold C. Helgeson in the 1960s (Helgeson 1964, 1969). For the properties of aqueous species, these studies culminated in a monumental, four-part series of papers in the American Journal of Science (Helgeson & Kirkham 1974a,b, 1976; Helgeson et al. 1981) that provided the HKF (Helgeson–Kirkham–Flowers) equations of state for aqueous species. Subsequent research built on this foundation, resulting in a series of papers documenting the revised HKF equations (Shock & Helgeson 1988, 1990; Tanger & Helgeson 1988; Shock et al. 1989, 1997). A major achievement of the revised HKF equations was the development of predictive estimation schemes for the equation of state parameters of aqueous species, which allowed widespread estimation of the characteristic properties of many aqueous species. The standard partial molal properties referring to 25°C and 1.0 bar and the equation of state coefficients of these species were incorporated into the database for the code SUPCRT92 (Johnson et al. 1992). The revised HKF approach has subsequently been used extensively in many geoscience and geotechnical applications. Its upper pressure limit, however, was 5.0 kbar for more than two decades. This limitation resulted from lack of knowledge of the dielectric constant of pure water at pressures >5.0 kbar (Fig. 2a). Consequently, for many years, theoretical geochemical modelling of water–rock interactions was limited to pressures less than or equal to 5.0 kbar, too low for consideration of fluid–rock interactions deep into subduction zones.

The dielectric constant of water is a vital parameter in the HKF approach because it appears in the Born equation for solvation that is part of the equation of state for computing the standard partial molal properties and the non-standard state contributions for aqueous species as functions of temperature and pressure (Helgeson et al. 1981). To remove the limitation of 5.0 kbar in the HKF approach, a characterization of the dielectric constant of water at greater pressures was developed by building on a semi-empirical equation for polar solvents (Franck et al. 1990) and ab initio molecular dynamics results to about 100 kbar (Pan et al. 2013). The new model for the dielectric constant of water (Fig. 2b) led to the development of the Deep Earth Water (DEW) model (Sverjensky et al. 2014a). By calibrating the model with experimental speciation data for ions referring to high pressures (Facq et al. 2014, 2016), and experimental data for the solubility and speciation of neutral species of aqueous silica and aluminium (Manning 1994; Newton & Manning 2002; Tropper & Manning 2007; Mysen 2010) revised predictive correlations were developed (Sverjensky et al. 2014a) allowing calculations of the standard partial molal properties of many aqueous species to 60 kbar.

The DEW model allowed prediction of the equilibrium constants for mineral hydrolysis reactions, aqueous complexing reactions, and many other types of reactions involving hydration and organic species. Quantification of the dielectric constant of water also allowed calculation of the parameters in the long-range electrostatic interaction term of the Debye–Hückel equation for aqueous ionic activity coefficients (see below). Together with equations for mass balance and charge balance, it became possible to quantitatively model aqueous speciation and solubilities and chemical mass transfer under upper mantle conditions. Examples of the application of DEW model predictions include the following: models of the importance of organic species with oxidation states of carbon ranging from −4 to +4 coexisting in subduction-zone fluids at pressures greater than about 30 kbar (Sverjensky et al. 2014b); prediction of the variable speciation of nitrogen in subduction-zone fluids and implications for the origin of the Earth's atmosphere (Mikhail & Sverjensky 2014); a new theory for the formation of diamond by pH drop at constant redox conditions (Sverjensky & Huang 2015); and a model for the breakdown of antigorite in subduction zones leading to a dramatic increase of fO2, a mechanism for generation of highly oxidizing fluids as part of the subduction process (Debret & Sverjensky 2017). Additional applications have focused on predictions of rock solubility and pH changes in metamorphic fluids (Galvez et al. 2015, 2016), the use of theoretical modelling to interpret solubility and speciation experiments (Huang 2017; Tumiati et al. 2017; Tiraboschi et al. 2018), new high-pressure experimental NMR studies of aqueous species (Pautler et al. 2014; Augustine et al. 2017), and the predictions of ab initio molecular dynamics for aqueous complexing reactions (Mei et al. 2018).

The purpose of the present study is to summarize the background needed to be able to carry out theoretical geochemical modelling of the type mentioned in the examples above. To do this, the presentation starts with a summary of the criteria for predicting equilibrium and non-equilibrium chemical reactions. In turn, this leads to the discussion of two important topics. First, I discuss how the thermodynamic activities of species are related to measurable concentrations through the definitions of standard states and activity coefficients. Specific examples at ambient conditions involving dilute water are detailed. The concept of aqueous speciation arises from this discussion and leads to a brief summary of why aqueous speciation models are needed and how they can be calculated. Second, I discuss how to model upper mantle fluids; specifically, how to calculate the standard Gibbs free energies of species and aqueous activity coefficients at elevated temperatures and pressures. The revised HKF equations of state are summarized and the revised predictive correlations for the estimation of equation of state coefficients are presented. Finally, an application involving the solubility and speciation of aqueous aluminium over a wide range of conditions is given.

Prediction of nonequilibrium reactions

For the completely general case of a single chemical reaction, and whether or not it should take place in the direction written, it is necessary to evaluate the sign of the Gibbs free energy of the reaction (ΔGr) according to 
ΔGr=ΔGr0+2.303RTlogQ
(1)
where ΔGr0 is the standard Gibbs free energy of reaction and Q is the thermodynamic activity quotient, which are defined by 
ΔGr0=iviΔG¯f,i0
(2)
and 
Q=i(ai)vi.
(3)
In equation (2), vi represents the stoichiometric reaction coefficient of the ith species in the reaction and ΔG¯f,i0 represents the standard Gibbs free energy of formation of the ith species. The standard Gibbs free energy of formation of each species refer in each case to the irreversible free energy of reaction forming the species in its standard state from its constituent chemical elements in their chosen reference state. Any standard state is a special set of conditions and behaviour for a species with a fixed chemical composition (see below). Therefore values of ΔG¯f,i0 for a given species in a chosen standard state depend only on pressure and temperature. The values are derived from experimental studies and models interpreting the experimental studies. The general topic of how to calculate values of ΔG¯f,i0 as functions of pressure and temperature is discussed later in this paper.

In equation (3), the thermodynamic activity of each species in the reaction is a dimensionless quantity which is related to the actual concentration of the species in the reaction of interest by the choice of standard states.

It follows from equations (1)–(3) that 
ΔGr=iviΔG¯f,i0+2.303RTlogi(ai)vi.
(4)
The two terms on the right-hand side of equation (4) can be summarized as follows. The first term is the free energy of the reaction of interest in the standard state, and is a function of temperature and/or pressure only. The standard Gibbs free energies of the individual species in the reactions depend on the standard states chosen. The second term is the activities of the real species in the reaction of interest, and is a function of temperature, pressure and composition. The activities are directly relatable to actual physically measurable concentrations, depending on the specific standard states chosen.

Evaluation of equation (4) leads to a predictive capability: if ΔGr<0, the reaction should proceed as written. However, if ΔGr>0, the reaction should proceed in reverse.

For example, suppose we want to answer the question: will calcite dissolve in rainwater? It is important to realize that, in principle, many alternative reactions representing the dissolution reaction could be written; for example, with alternative species such as CO2(aq), CO2(g), HCO3 or CO32. The species to choose is the one that is most practical. Rainwater can be assumed to be in equilibrium with atmospheric CO2. Therefore, for convenience, we can write the reaction as 
CaCO3calcite+2H+Ca2++CO2(g)+H2O
(5)
for which we would need to evaluate the sign of ΔGr calculated from the equation 
ΔGr=ΔGr0+RTlnaCa2+aCO2(g)aH2OaCaCO3calciteaH+2.
(6)
We can evaluate the first term on the right-hand side of equation (6) if we have values of ΔG¯f,i0 for each of the species in equation (5) in their standard states. We can evaluate the second term if we know how standard states relate the activities to measurable concentrations in the actual samples of calcite and rainwater. For example, how is the aCaCO3calcite related to how much Ca versus Mg is in the calcite? And how is the aCa2+ related to the actual concentration of Ca2+ in the rainwater? Consequently, the definitions of the standard states to be used are of the utmost importance to addressing the question of whether the calcite should dissolve in the rainwater. However, before discussing standard states, it is important to recognize a special circumstance of equation (1) above. What happens when ΔGr=0?

Prediction of equilibrium reactions

When ΔGr=0 in equation (1), the reaction being described by the equation is said to be at equilibrium. The rate of the forward reaction is balanced by the rate of the back reaction so that there is no net change with time; that is, 
0=ΔGr0+2.303RTlogQ.
(7)
Standard states depend only on pressure and temperature. In these circumstances, at a given pressure and temperature, ΔGr0 in equation (7) is fixed. Consequently, the value of Q is fixed. This special value of Q is called the thermodynamic equilibrium constant and is denoted by K according to 
ΔGr0=2.303RTlogK.
(8)
It should be emphasized that K represents the special value of Q that depends only on the pressure and temperature under consideration. This result follows from the definitions of the standard states as depending only on the pressure and temperature of the reaction as will be defined below. Consequently, K is also a standard state quantity. To recognize this, K should have a superscript (e.g. it should be written K0), but by convention the superscript is omitted. Here, again, ΔGr0 represents the standard Gibbs free energy of the reaction given by 
ΔGr0=iviΔG¯f,i0
(9)
and 
K=i(ai)vi.
(10)
K is the thermodynamic equilibrium constant, which fixes a unique activity ratio for the species in the reaction at a given temperature and a specific pressure. It should be noted that the individual activities are not fixed; only the activity quotient is fixed.

Definitions of standard states

Minerals

The simplest possible case is a solid solution of the components MZ and NZ in which compositional variation occurs only on a single site in the structure; that is, (M,N)Z. For multi-site mixing, a summary has been given by Anderson (2005). For the example here, the activity of MZ, aMZ, is given by 
aMZ=λMZXMZ
(11)
where XMZ is a mole fraction of MZ in the solid solution and λMZ is a rational activity coefficient, a function of temperature, pressure and XMZ.

Here the standard state refers to unit activity of the pure component MZ at the temperature and pressure of interest and it is assumed that Raoult's Law is followed at high mole fractions of MZ; that is, it is assumed that λMZ1.0 as XMZ1.0, as depicted in Figure 3.

In other words, the mineral has a physically achievable standard state when it is pure. An example is quartz. In pure quartz, which for all practical purposes does exist, SiO2 is in its standard state; that is, 
aSiO2quartz=XSiO2quartz=1.
(12)
But even when the mineral is not pure, as in Figure 3, we can often assume that with high enough Xj → 1.0 the component MZ will behave as though it is experiencing ideal mixing when XMZ is in the region where λj1.0. For example, the component CaCO3 in impure calcite is not in its standard state.
But if calcite is only slightly impure, λCaCO31, so 
aCaCO3calciteXCaCO3calcite.
(13)

Water and mixed solvents

Water is such an important solvent that we will discuss it first, before discussing mixed solvent fluids. Under surficial conditions, water can contain a wide range of dissolved ionic solutes. The usual standard state chosen for water under such conditions is exactly analogous to that discussed above for minerals in which the activity is related to the mole fraction by the equation 
aH2O=λH2OXH2O
(14)
where XH2O is the mole fraction of H2O and λH2O is a rational activity coefficient (a function of temperature, pressure, and XH2O).

Here, the standard state refers to unit activity of the pure component H2O at the pressure and temperature of interest according to Raoult's Law; that is, λH2O1.0 as XH2O1.0.

In other words, pure liquid water is in its standard state and we can write 
aH2O=1.0.
(15)
Furthermore, we can often assume that even impure water is in the region where XH2O1.0 and λH2O1.0, and therefore behaves ideally; that is, that 
aH2OXH2O1.0.
(16)
Under surficial, crustal conditions, where the mole fraction of water in aqueous solutions is very high, this is a good approximation. Examples include rainwater, river water and shallow groundwaters. However, in general for seawater, and particularly for saline brines such as those associated with evaporites, or in oilfield brines, or saline magmatic fluids, equation (16) is insufficient.

More generally, for a component MZ in a multicomponent fluid, such as metamorphic fluids, the activity of H2O or other dissolved gases may be highly non-ideal. In these circumstances, the activity of MZ may be depicted as in Figure 4.

The standard state shown in Figure 4 is commonly used for fluid CO2 and other neutral molecules that can be thought of as mixing with H2O at elevated temperatures and pressures. It should be noted that the standard Gibbs free energy of formation for CO2 that should be used for this standard state must be that of pure fluid CO2, completely different from the standard free energy of formation of CO2 in aqueous solution that is required by the hypothetical 1.0 molal standard state commonly used for aqueous species (see below).

Gases

For gaseous species, the general equation for the relationship between the thermodynamic activity and the measurable partial pressure of a gas is given by 
aj=fjfj0=υjυj0PjPj0
(17)
where fj is the fugacity (with units of pressure) and fj0 is the fugacity of species j in the standard state. Each fugacity is related to a partial pressure of a gas by a fugacity coefficient, where υj is the fugacity coefficient for the real gas and υj0 is the fugacity coefficient for the gas when it is in its standard state.
The gas standard state in the present study refers to unit activity of the pure ideal gas at the temperature of interest and a pressure of 1.0 bar; that is, υj0 and pj0 are equal to 1.0 for the standard state. In other words, 
aj=fjfj0=υjPj.
(18)
A real gas cannot actually behave ideally at the temperature of interest and 1.0 bar. Therefore, the gas standard state is an example of a hypothetical standard state, one that cannot be physically achieved by a gas. It contrasts with the standard states chosen for minerals and water discussed above, which are physically achievable (at least in principle) when the minerals and the water are pure.
At ambient pressures, we can assume that the fugacity coefficient is unity and so 
aj=pj.
(19)
However, at elevated temperatures, the fugacity coefficients are necessary because they are a function of pressure, temperature, and the composition of a gas mixture.

Aqueous species

In the case of aqueous species of all types, the relationship between the thermodynamic activity and concentration is given by 
aj=γjmjγj0mj0
(20)
where γj and mj refer to the activity coefficient and molality of j, and γj0 and mj0 refer to the standard state values of these quantities (Anderson 2005). The most frequently used standard state refers to unit activity of a hypothetical 1.0 molal solution of j referenced to infinite dilution at the temperature and pressure of interest. Because of this definition, real aqueous species can never be in their standard state! In other words, the standard state for aqueous species is not physically achievable. The hypothetical 1.0 molal aqueous species standard state is depicted in Figure 5.
It can be seen in Figure 5 that at a molality of 1.0, the activity is 1.0, the slope of the dashed line is 1.0, and so the activity coefficient is 1.0 in the hypothetical 1.0 molal standard state. As a result, equation (20) becomes 
aj=γjmjγj0mj0=γjmj(1.0)(1.0)=γjmj
(21)
which is the commonly used expression for the relationship between the thermodynamic activity of an aqueous species and the molality, assuming use of the hypothetical 1.0 molal standard state.

Activity coefficients linking standard states to measurable concentrations in fluids

Much has been written about the calculation of aqueous activity coefficients (e.g. Anderson 2005). With aqueous fluids under surficial conditions or shallow crustal temperatures and pressures, the discussions in the literature come down to choices and recommendations of the approaches of Pitzer versus Helgeson (e.g. Anderson 2005). The focus here has always been primarily on how to treat the activity coefficients of aqueous electrolytes and their associated neutral species and ionic complexes. Little consideration has been given to the effects of large amounts of dissolved gases (e.g. CO2, CH4, N2). For metamorphic and deep Earth fluids, numerous variants of so-called COH-models of dissolved gases are used (e.g. Zhang & Duan 2009). A few studies have attempted to combine COH-models with the effects of dissolved electrolytes, and their associated neutral species and ionic complexes (Galvez et al. 2015, 2016). However, it is important to realize that, from a practical standpoint, the modelling of fluid–rock interactions in the deep Earth at elevated temperatures and pressures in fluids of complex composition requires some simplifying assumptions within the framework of a predictive model for computing activity coefficients. The only usable framework built to be as general, and predictive, as possible is the one described by Helgeson et al. (1981). Other treatments have too many adjustable parameters to be able to function in a predictive mode.

Evaluation of the activity coefficients of aqueous ionic species

For aqueous ionic species in solutions of a single 1:1 electrolyte (k) containing a monovalent cation (i) and anion (l), the individual ionic activity coefficient is related to the mean ionic activity coefficient of the electrolyte (γ¯±,k) by the equation 
γ¯±,k=(γ¯iγ¯l)0.5
(22)
which is based on equation (64) of Helgeson et al. (1981).
The individual ionic activity coefficient of any ionic species (j) with a charge Zj is represented by 
logγ¯j;P,T=AγZj2I¯0.51+akBγI¯0.5+bγ,kI¯+Γγ
(23)
which is based on equation (165) of Helgeson et al. (1981).
In equation (23), the first term on the right-hand side is referred to as the Debye–Hückel term. It represents long-range ionic interactions. In this term, Aγ and Bγ are functions of temperature, and the density and dielectric constant of pure water. The ionic strength of the solution (I¯) is given by 
I¯=12imiZi2
(24)
and ak represents an ion-size parameter for the kth salt in solution (see equations (124) and (125) of Helgeson et al. 1981).
The second term, referred to as the extended term, represents a combination of the ionic strength dependence of ion solvation and the short-range electrostatic interactions of ions. In this term, the parameter bγ,k is a function of temperature and pressure for the kth predominant salt in solution. Extensive description of the evaluation of this term has been given by Helgeson et al. (1981). The last term in equation (23) represents a conversion factor from the mole fraction to the molality concentration scale given by 
Γγ=log(1+0.0180153m)
(25)
where m* represents the sum of the molalities of all aqueous species in the solution.

Behaviour of electrolyte activity coefficients

To evaluate the extended term contribution in equation (23), a deviation function for the kth electrolyte (ργ,k;P,T) can be defined according to 
ργ,k;P,T=logγ¯k;P,TAγ|ZiZl|I¯0.51+akBγI¯0.5Γγ
(26)
that is, 
ργ,k;P,T=bγ,kI¯
(27)
where ργ,k;P,T is obtained from experimental data and plotted against the stoichiometric ionic strength of the solution. Linear behaviour according to equation (27) allows evaluation of bγ,k directly (Fig. 6). Extensive analysis of experimental data for dozens of salts detailed by Helgeson et al. (1981) has demonstrated that in solutions of a single electrolyte, the simplicity of equation (27) is valid to very high ionic strengths, something that seems to have been forgotten in the geochemical literature. Systematic departures from linearity as shown in Figure 6 for the NaOH electrolyte at the highest ionic strengths indicate the onset of aqueous complexing, in this instance the likely formation of the species NaOH0, which is part of the speciation of the overall NaOH electrolyte.

Further analysis for many electrolytes allowed the development of an equation of state characterization for bγ,k as a function of temperature and pressure. This equation of state was revised and parameterized for NaCl, building on the revised, HKF standard state equations (Oelkers & Helgeson 1990). However, none of the above treatments included an extended term contribution in equation (23) for the solvation dependence of the ions on changes in the solvent associated, for example, with large variations in the CO2 content of the fluid.

Evaluation of the activity coefficients of neutral aqueous species

The activity coefficient of the nth neutral aqueous species in a solution of the kth electrolyte can be expressed by 
logγ¯n;P,T=bγ,nI¯+Γγ
(28)
where the parameter bγ,n refers to the ionic strength dependence of the solvation of the neutral species and the short-range electrostatic interactions of ions with the neutral species. Although no specific parameterization of bγ,n was made by Helgeson et al. (1981), estimates of the temperature and pressure dependence of bγ,n were made subsequently for a number of 1:1 ion pairs (Oelkers & Helgeson 1991).

CO2 and other dissolved gases: solute or solvent?

In principle, equation (28) could be used not only for neutral ion pairs as described above, but also for dissolved gases such as CO2. The thermodynamic treatment of CO2 in water has typically taken one of two paths. For low temperature (i.e. surficial conditions) and geothermal systems, CO2 has been treated with the hypothetical 1.0 molal standard state, as for aqueous ions. In these circumstances, equation (28) could be used, provided that sufficient experimental data to calibrate the equation are available. At the higher temperatures and pressures of metamorphic fluids and deep Earth fluids in general, where the potential range of fluid compositions can be much wider, extending even to extremely CO2-rich fluids, a completely different approach has been taken. Here, in the COH-type fluid model, the fluid is treated as a number of dissolved gases (e.g. Zhang & Duan 2009) for which the standard states are the pure fluids at the temperature and pressure of interest; that is, analogous to that of H2O discussed above. Furthermore, the mole fraction concentration scale is used, as it is more practical when the range of fluid chemistry can approach pure CO2. However, recent ab initio molecular dynamics calculations (Pan & Galli 2016) have shown that at 1000 K and 100 kbar Na2CO3 solutions contain mainly ions, and that the predominant neutral species is not CO2 but instead is H2CO3. The discovery of the importance of H2CO3 will require a complete re-evaluation of the speciation and activity–concentration models for COH fluids over a wide range of temperatures and pressures from crustal to mantle conditions.

Nonequilibrium reactions at ambient conditions: prediction of calcite dissolution

Let us consider the dissolution of calcite in a natural water of some given composition at 25°C and 1 bar: 
CaCO3calcite+2H(aq)+Ca(aq)2++CO2(g)+H2O(l).
(29)
To predict whether or not calcite will dissolve in a water of given composition it is necessary to evaluate the sign of ΔGr using 
ΔGr=ΔGr0+2.303RTlogQ.
(30)
The sign of ΔGr tells us whether or not the reaction should proceed as written, or opposite to what is written, or whether the reaction is at equilibrium, according to 
ΔGr<0,reactionproceedsaswritten
(31)
 
ΔGr>0,reactionproceedsinreverse
(32)
 
ΔGr=0,reactionisatequilibrium.
(33)
The two terms on the right-hand side of equation (30) must be evaluated. The first of these is the standard Gibbs free energy of reaction, which is evaluated from the equation 
ΔGr=iviΔG¯f,i0
(34)
where viΔG¯f,i0 represents the number of moles of the ith species in the reaction multiplied by the standard Gibbs free energy of formation of the species from its elements in their defined reference states. It should be noted that equation (34) only involves standard state free energies, which depend only on temperature and pressure. Each standard state refers to a fixed chemical composition, regardless of the actual state of interest in the natural water referred to in equation (29).
Using values of ΔG¯f,i0 at 25°C and 1 bar (Berman 1988; Shock et al. 1997) we have 
ΔG¯f,calcite0=269880calmol1;ΔG¯f,CO2(g)0=94254calmol1;ΔG¯f,H2O0=56686calmol1;ΔG¯f,Ca2+0=132120calmol1;ΔG¯f,H+0=0.0calmol1.
Substituting these standard free energies into equation (34) permits evaluation of the standard Gibbs free energy of the reaction according to 
ΔGr0=ΔGf,Ca2+0+ΔGf,CO2(g)0+ΔGf,H2O0ΔGf,calcite02ΔGf,H+0=1321209425456686+2698800=13180calmol1.
(35)
It is important to remember that the value of ΔGr0 given by equation (35) depends only on temperature and pressure, here 25°C and 1.0 bar. It is independent of the chemistry of the water under consideration for reaction with the calcite and whether or not the calcite is a solid solution containing other chemical elements.

We now consider four different situations in which the chemistry of the water and the calcite are variables.

Pure calcite reacting with rainwater equilibrated with atmospheric carbon dioxide.

It is assumed that the concentration of Ca in the rainwater is 2.0 ppm, that the pH is 5.7, and that pCO2 is 10−3.5 bar.

The stipulation of pure calcite can be explicitly included in equation (29) in a way that indicates that the component CaCO3 is the only component in the calcite; that is, 
CaCO3purecalcite+2H+Ca2++CO2(g)+H2O.
(36)
Together with recognition of the dilute nature of rainwater, which allows us to assume that the thermodynamic activity of Ca2+ will be equal to its molality, equation (30) becomes 
ΔGr=ΔGr0+RTlnaCa2+aCO2(g)aH2OaCaCO3purecalciteaH+2=ΔGr0+2.303RTlogmCa2+pCO2(g)aH+2
(37)
where the aCaCO3purecalcite=1.0 because the component CaCO3 is in its standard state when the calcite is pure.
Now evaluating the activity quotient term in equation (37) results in 
logmCa2+pCO2(g)aH+2=log(mCa2+)+log(pCO2(g))2log(aH+)=log21000×40.083.5+2(5.7)=log(5.0×105)3.5+2(5.7)=4.33.5+11.4=3.6.
(38)
Overall, combining equations (35), (37), and (38) results in 
ΔGr=13180+2.303RT(3.6)=13180+2.303(1.987×298.15)(3.6)=13180+4912=8268calmol1.
(39)
Based on equation (31), the negative sign for the calculated value of ΔGr in equation (39) indicates that the reaction in equation (36) should proceed as written. Consequently, pure calcite should dissolve in this particular rainwater; that is, the rainwater is undersaturated with respect to pure calcite.

Calcite solid-solution (20% MgCO3) reacting with rainwater equilibrated with atmospheric carbon dioxide.

It is again assumed that the concentration of Ca in the rainwater is 2.0 ppm, that the pH is 5.7, and that pCO2 is 10−3.5 bar.

The stipulation of impure calcite can be explicitly included in equation (29) in a way that indicates that the component CaCO3 is not the only component in the calcite; that is, 
CaCO3impurecalcite+2H+Ca2++CO2(g)+H2O(l).
(40)
Proceeding as above we now obtain 
ΔGr=ΔGr0+RTlnaCa2+aCO2(g)aH2OaCaCO3calciteaH+2=ΔGr0+2.303RTlogmCa2+pCO2(g)0.8aH+2
(41)
where the aCaCO3calcite=0.8 by assuming that the CaCO3 component mixes ideally with the MgCO3 component.
Now evaluating the activity quotient term in equation (41) results in 
logmCa2+pCO2(g)0.8aH+2=log(mCa2+)+log(pCO2(g))2log(aH+)log(0.8)=log21000×40.083.5+2(5.7)+0.097=log(5.0×105)3.5+2(5.7)+0.097=4.33.5+11.4+0.097=3.7.
(42)
Overall, combining equations (35), (41) and (42) results in 
ΔGr=13180+2.303RT(3.7)=13180+2.303(1.987×298.15)(3.7)=13180+5408=8132calmol1.
(43)
It should be noted that in equation (43), the resulting free energy of the reaction is only 136 calories more positive than in equation (39). Consequently, the impure state of the calcite has only a very small effect on the free energy of the dissolution reaction. The negative sign for the calculated value of ΔGr in equation (43) indicates again that the reaction in equation (40) should proceed as written. Consequently, the impure calcite should also dissolve in this particular rainwater.

Pure calcite reacting with soil water equilibrated with pCO2=101.5 bar.

It is assumed that the concentration of Ca in the soil water is still 2.0 ppm, but that the pH is 7.7, owing to weathering reactions.

The stipulation of pure calcite is again included in the reaction: 
CaCO3purecalcite+2H+Ca2++CO2(g)+H2O(l).
(44)
Proceeding as in Case 1, we obtain 
ΔGr=ΔGr0+2.303RTlogmCa2+pCO2(g)aH+2
(45)
where now the soil water is so different in chemistry from the original rainwater that 
logmCa2+pCO2,gaH+2=log(mCa2+)+log(pCO2(g))2log(aH+)=log21000×40.081.5+2(7.7)=log(5.0×105)1.5+2(7.7)=4.31.5+15.4=9.6.
(46)
We note that the first term on the right-hand side of equation (45), the standard Gibbs free energy of reaction, will be identical to that in all the previous examples. The standard states for the species are the same in the rainwater and the soil water. Therefore 
ΔGr=13180+2.303RT(9.6)=13180+2.303(1.987×298.15)(9.6)=13180+13098=82calmol1.
(47)
The very small negative number for the calculated value of ΔGr in equation (47) is equal to zero within the uncertainties of all the standard state free energies used in the calculation. Consequently, the calcite should be in equilibrium with this soil water.

Pure calcite reacting with surface seawater equilibrated with the atmosphere at pCO2=103.5 bar.

It is assumed that the concentration of Ca is now 411 ppm, and that the pH is 8.1.

For pure calcite, we again have the reaction 
CaCO3purecalcite+2H+Ca2++CO2(g)+H2O(l)
(48)
and if we continue to assume that the aqueous activities can be approximated by molalities, which is not strictly valid for seawater, we obtain 
ΔGr=ΔGr0+2.303RTlogmCa2+pCO2(g)aH+2
(49)
where 
logmCa2+pCO2(g)aH+2=log(mCa2+)+log(pCO2(g))2log(aH+)=log4111000×40.083.5+2(8.1)=log(0.0103)3.5+2(8.1)=2.03.5+2(8.1)=10.7.
(50)
Therefore 
ΔGr=13180+2.303RT(10.7)=13180+2.303(1.987×298.15)(10.7)=13180+14599=+1419calmol1.
(51)
The positive sign for the calculated value of ΔGr in equation (51) now indicates that the reaction in equation (48) should proceed in the opposite direction to that written. Consequently, calcite should precipitate from surface seawater; that is, the seawater is supersaturated with respect to calcite.

Summary of how to predict a nonequilibrium chemical reaction

For a single chemical reaction, such as the dissolution of calcite according to 
CaCO3calcite+2H+Ca2++CO2(g)+H2O
(52)
the fundamental equation to be evaluated is 
ΔGr=ΔGr0+2.303RTlogQ
(53)
which requires evaluation of the standard Gibbs free energy of reaction (ΔGr0) and the activity quotient (Q). Several things should be clear from the above examples involving pure versus impure calcite and three different water chemistries representing rainwater, soil water and surface seawater, as follows.

  1. The numerical value of ΔGr0 was the same for all four examples. All the examples referred to the same temperature and pressure (i.e. 25°C and 1.0 bar). Therefore, the numerical value of ΔGr0 was the same for all four examples, even though the mineral chemistry and the water chemistry were varied. The reason is that the value of ΔGr0 depends on the standard state Gibbs free energies of formation of the species in the reaction, and the standard states depend only on temperature and pressure.

  2. The sign of ΔGr is what is important. The sign of ΔGr0 is not important.

  3. The only exception to the above rule occurs when all the species in a reaction are actually in their standard states; for example, a reaction containing two pure minerals, such as pure calcite reacting to pure aragonite, which can be represented by 
    CaCO3purecalciteCaCO3purearagonite.
    (54)

Another example would be any reaction consisting of pure minerals and water that is (traditionally) assumed to be pure or very close to pure; for example, the breakdown of talc according to 
Mg3Si4O10(OH)2puretalcMgSiO3pureenstatite+SiO2quartz+H2O(l).
(55)
For all the species in the reactions in equations (54) and (55), the activities are equal to unity because the minerals and the water are all assumed to be pure, end-member species. Therefore Q = 1.0 and hence 
ΔGr=ΔGr0.
(56)
In these circumstances, the sign of ΔGr0 becomes meaningful because it is the same as for ΔGr.

But this is a special circumstance.

  1. When aqueous ions or gases are in the reaction, Q is never equal to unity. Aqueous ion and gas species have hypothetical standard states; that is, they can never be in their standard state in the real reactions of interest.

  2. The activities of any aqueous species in the reaction must be evaluated. Calculation of Q for a reaction involving aqueous species requires that the activities of the specific species in the reaction be evaluated; for example, when we evaluated the equation 
    ΔGr=ΔGr0+RTlnaCa2+aCO2,gaH2OaCaCO3calciteaH+2=ΔGr0+2.303RTlogmCa2+pCO2(g)0.8aH+2
    (57)

we assumed that the only aqueous Ca-bearing species in the water was the Ca2+ ion. Therefore, the activity of the ion required for equation (57) could be set equal to the measured analytical value of Ca (mt,Ca); that is, 
aCa2+=mCa2+=mt,Ca.
(58)
Although the assumption represented by equation (58) is adequate for the examples above involving rainwater and soil water, it is not completely accurate for seawater and is totally inaccurate for brines or fluids at elevated temperatures and pressures. Under the latter conditions particularly, the molality of an ion is much lower than the total dissolved concentration of the element or elements making up that ion. The reason is that in brines and at elevated temperatures and pressures ions form complexes with other ions, and neutral species. More generally, we need to consider the aqueous speciation of each chemical element (Wolery 1983; Bethke 1996; Parkhurst & Appelo 1999; Anderson 2005). Even in our surface seawater example at ambient conditions, discussed above, the molality of the free (uncomplexed) Ca2+ ion is about 92% of the total dissolved Ca. In the cases of other important ions in seawater such as carbonate and sulphate, the percentages of the free ions relative to the total dissolved concentrations of the elements are much lower (Table 1). For example, the amount of free carbonate ion is only about 2.6% of the total dissolved carbon, the rest is complexed by H+, Na+, Mg2+ and Ca2+. For sulphate, the percentage of the free ion is about 55% of the total dissolved sulphur.

The need to take into account aqueous speciation means that, in general, we cannot just use the total dissolved concentrations of the elements as the molalities of individual ions when evaluating dissolution reactions such as in the examples discussed above. Equations such as equation (57) above should only be evaluated after carrying out a separate calculation in which a model of the full aqueous speciation is built. Based on the aqueous speciation model, the individual species’ molalities and thermodynamic activities can be used to assess the state of saturation of a given aqueous fluid or to assume equilibrium with one or more minerals to calculate solubilities.

Equilibrium aqueous speciation and solubility models at ambient conditions

Speciation of simple NaCl–HCl–H2O solutions

We start by assuming that the solution could contain the following eight aqueous species: Na+, Cl, H+, OH, NaCl0, HCl0, NaOH0 and H2O. Corresponding to these species there are eight thermodynamic activities and concentrations to know. For simplicity, if we assume that the activity of water is unity (i.e. aH2O=1) we are left with seven unknown activities and concentrations.

However, there are relationships between these species because of equilibrium in solution: 
Na++Cl=NaCl0
(59)
 
Na++OH=NaOH0
(60)
 
H++Cl=HCl0
(61)
 
H++OH=H2O.
(62)
These equilibria correspond to the following mass action equations: 
KNaCl0=aNaCl0aNa+aCl=γNaCl0γNa+γClmNaCl0mNa+mCl
(63)
 
KNaOH0=aNaOH0aNa+aOH=γNaOH0γNa+γOHmNaOH0mNa+mOH
(64)
 
KHCl0=aHCl0aH+aCl=γHCl0γH+γClmHCl0mH+mCl
(65)
 
KH2O=aH2OaH+aOH=1γH+γOH1mH+mOH
(66)
where KNaCl0, KNaOH0, KHCl0 and KH2O refer to the thermodynamic equilibrium constants for the reactions.

In addition to the mass action equations, we have mass balance and charge balance equations to use as constraints, as follows.

Mass balance equations: 
forNamt,Na+=mNa++mNaCl0+mNaOH0
(67)
 
forClmt,Cl=mCl+mNaCl0+mHCl0.
(68)
Charge balance equation: 
mCl+mOH=mNa++mH+.
(69)
Equations (63) –(69) constitute a total of seven equations in seven unknowns, which can be solved if we know the total sodium and the total chloride concentrations, and the K values and γ values. Solving the seven equations in seven unknowns results in the concentrations and activities of all seven species. Then the state of saturation of the fluid with respect to any minerals in the system can be evaluated using 
ΔGr=ΔGr0+2.303RTlogQ=ΔGr0+2.303RTlogΠi(ai)vi
(70)
for each possible mineral dissolution reaction.

Modelling at upper mantle conditions

I. Equilibrium constants over a wide range of temperatures and pressures

What is needed?

Equilibrium constants for reactions with all possible kinds of species are needed to compute aqueous speciation, solubility and chemical mass transfer models. Traditionally in aqueous geochemistry equilibrium constants involving aqueous species have been calculated using the computer code SUPCRT92. An example is shown below in Figure 7a for the dissociation constant of Mg(HCO3)+ where the theoretically calculated curve was fit to the experimental data (Siebert & Hostetler, 1977; Stefánsson et al. 2014). Extrapolation of such calculations using the equations of state for the standard Gibbs free energies of the individual aqueous species have been limited to 1000°C and 5.0 kbar.

In contrast, the DEW model can be used to calculate equilibrium constants involving aqueous species at temperatures greater than 1000°C and pressures up to 60 kbar. An example is shown in Figure 7b for the dissociation of aqueous silica Si(OH)40 to the monovalent silicate anion SiO(OH)3 (note that in the reaction shown in the figure water molecules have been subtracted). The calculated curves were fit to the experimental data points (Seward, 1974; Busey & Mesmer, 1977; Kessel et al. 2005b). The basis for calculations such as those shown in Figure 7a and b is presented below.

Any equilibrium constant can be calculated from the standard Gibbs free energy of the reaction according to the equation 
logK=ΔGr;P,To2.303RT.
(71)
Clearly, we need to be able to calculate standard Gibbs free energies of reactions (ΔG¯r;P,T0) at elevated temperatures and pressures. To do this we need the standard Gibbs free energies of each species in the reaction at elevated temperatures and pressures (ΔG¯f,j;P,T0); that is, 
ΔGr;P,To=jvjΔG¯f,j;P,T0.
(72)
As an example, the equilibrium between calcite and water may be written in terms of the dissolved CO20 species according to 
CaCO3calcite+2H+=Ca2++CO20+H2O
(73)
for which we can write 
ΔGr;P,To=ΔG¯Ca2+;P,T0+ΔG¯CO20;P,T0+ΔG¯H2O;P,T0ΔG¯CaCO3calcite;P,T02ΔG¯H+;P,T0.
(74)
It should be noted here that, in practice, the free energies on the right-hand side of equation (74) refer to apparent standard partial molal Gibbs free energies of formation from the elements for the jth species (Helgeson et al. 1978) where 
ΔG¯j;P,T0=ΔG¯f,j;Pr,Tr0+(G¯j;P,T0G¯j;Pr,Tr0).
(75)
It can be seen in equation (75) that the first term on the right-hand side refers to the true standard Gibbs free energy of formation of the jth species from the elements at the reference temperature and pressure. However, the second term on the right-hand side refers to the change in the standard partial molal free energy of the jth species alone with temperature and pressure. The corresponding changes for the elements that make up the jth species are not evaluated. By using equation (75) only the apparent standard Gibbs free energy of formation of the jth species is calculated (ΔG¯j;P,T0). Values of ΔG¯j;P,T0 from equation (75) are intended to be used in equations such as equation (74) because such equations refer to balanced chemical reactions (equation (73)) for which the changes in temperature and pressure for all the elements if used would simply cancel. Consequently, evaluation of equations (74) and (75) reduce to the issue of calculating the term G¯j;P,T0G¯j;Pr,Tr0, which can be carried out using the following equation: 
G¯j;P,T0G¯j;Pr,Tr0=S¯j;Pr,Tr0(TTr)+TrTC¯j;Pr0dTTTrTC¯j;Pr0dlnT+PrPV¯j;T0dP.
(76)
Here the standard partial molal entropy of the jth species is needed at the reference temperature and pressure (S¯j;Pr,Tr0), as well as the temperature dependence of the isobaric standard partial molal heat capacity (C¯j;Pr0), and the pressure dependence of the standard partial molal volume (V¯j;T0). Equations for evaluating these functions are described below for minerals and for aqueous species.

Mineral free energies at elevated pressures and temperatures

In the original SUPCRT92 model the thermodynamic properties of minerals were taken from Helgeson et al. (1978). However, the heat capacity function used for minerals does not extrapolate sensibly to very high temperature (Berman & Brown 1985) as indicated in Figure 8a. Furthermore, the volumes of the minerals (with the exception of quartz) were treated as constants referring to 25°C and 1.0 bar. The data for antigorite in Figure 8b clearly show that this is not the case. The SUPCRT92 model for the C¯j;Pr0(T) and V¯j;T0(P) of minerals was sufficiently accurate for shallow crustal temperatures and pressures, but conditions in the deep continental crust or upper mantle require more elaborate treatments of the heat capacities and volumes of minerals.

In the overall DEW model (Sverjensky et al. 2014a), the thermodynamic properties of minerals are calculated using the formulation for the C¯j;Pr0(T) and V¯j;T0(P) from Berman (1988) with a modification for placing the free energies and enthalpies of Na- and K-bearing silicates on an absolute basis (Sverjensky et al. 1991). These changes resulted in the code SUPCRT92b used in Figure 8a and b.

For the isobaric heat capacities of minerals, the equation 
Cp=k0+k1T0.5+k2T2+k3T3
(77)
is used, where k0, k1, k2 and k3 are empirical coefficients fit to experimental heat capacity data as a function of temperature. This heat capacity function for minerals extrapolates sensibly to very high temperatures of the order of 2500–3000°C. Furthermore, when experimental heat capacity data are not available, the coefficients can be estimated using a predictive scheme (Berman & Brown 1985).
For the volumes of minerals, the equation 
vP,TvPr,Tr=1+v1(PPr)+v2(PPr)2+v3(TTr)+v4(TTr)2
(78)
is used, where v1, v2, v3 and v4 are empirical coefficients fitted to experimental volume data.
Combining the integrals of equations (77) and (78) in equations (75) and (76) results in the overall equation for calculating the apparent standard Gibbs free energies of minerals at elevated pressure and temperature consistent with Berman (1988) according to 
ΔGP,T0=ΔGf;Pr,Tr0TSPr,Tr0+k0{(TTr)T(lnTlnTr)}+2k1{(T0.5Tr0.5)+T(T0.5Tr0.5)}k2(T1Tr1)+T2(T2Tr2)k3(T2Tr2)2T3(T3Tr3)+VPr,Tr0[v12v2(P2Pr2)+v23(P3Pr3)+{1v1+v2+v3(TTr)+v4(TTr)2}(PPr)].
(79)

Aqueous species at elevated temperatures and pressures

Here again, the fundamental information needed to evaluate equations (75) and (76) is the S¯Pr,Tr0, C¯Pr0(T) and V¯T0(P) of the species in the reactions. Aqueous species behave very differently from minerals as functions of temperature and pressure. Decades of experimental and theoretical work have gone into investigating the temperature and pressure dependence of the thermodynamic properties of aqueous species. In the present study, we focus on the Helgeson–Kirkham–Flowers (HKF) approach, which goes back to the foundational four-part series (Helgeson & Kirkham 1974a, b, 1976; Helgeson et al. 1981) and forms the basis for all subsequent revisions and extensions. The primary reason is that the HKF approach was built to have predictive capabilities, particularly for the standard state properties of aqueous species. Some background on the development of the HKF approach is summarized below.

For the standard state properties of species in water, the solvent water is the model substance. Therefore, we start with the properties of pure water as a function of temperature and pressure. At room temperature and pressure, the electronegativity of oxygen in the water molecule leads to a net negative charge on the oxygen and net positive charges on the hydrogens. Consequently, the water molecule has a large dipole moment. In turn, this leads to a structure for water through H-bonding (Fig. 9).

On a molecular scale, when ions are dissolved in water they exert an electric field on the water molecules. The behaviour of water in applied electric fields is therefore important for understanding the behaviour of salts in water. Water is a dielectric (i.e. an insulator) because application of an electric field only results in displacement of electrons within molecules. But this affects the orientation of the water molecules relative to the ions. The behaviour of water in electric fields is expressed by the very high dielectric constant of water (εH2O=78) at 25°C and 1 bar. The dielectric constant of pure water can be thought of as a measure of the degree of polarization of water molecules in an electric field. The very high dielectric constant at ambient conditions results in the strong solvation of ions by water molecules. However, the dielectric constant of water changes very strongly with temperature and pressure.

Dielectric constant of water at elevated temperatures and pressures

Experimental measurements of the dielectric constant of water are shown in Figure 10a and b (Heger et al. 1980), together with estimations based on a literature review (Fernández et al. 1997), and ab initio molecular dynamics calculations (Pan et al. 2013). The curves in Figure 10a represent a fit of the data with an empirical polynomial, whereas the curves in Figure 10b represent a semi-empirical theory calibrated with only two parameters (Franck et al. 1990; Sverjensky et al. 2014a). Clearly, the dielectric constant of water decreases strongly with temperature, and increases strongly with pressure when the pressure is in the tens of kilobars range.

The strong decrease of the dielectric constant of water with increasing temperature leads to strong complexing between oppositely charged ions. Even for salts that are strong electrolytes at ambient conditions, the changing dielectric properties of water lead to complexing (ion-pairing) at elevated temperatures and pressures (Figs 11a and b). Decreases in the dielectric constant of water lead to weaker solvation shells around cations and anions, which allow stronger electrostatic and chemical interactions between the oppositely charged ions, favouring the formation of complexes. The types of complexes that form in crustal hydrothermal fluids have been intensively studied in relation to the origins of hydrothermal ore deposits. However, in deep crustal and upper mantle fluids, much less information is available, particularly at high temperatures and pressures above 2.0 GPa. It is an area of intense current studies. We will return to this topic towards the end of the study.

Volume and heat capacity of water at elevated temperatures and pressures

The molar volume, or density, of pure water is a fundamental quantity that is a basis for the HKF approach. The DEW model uses an equation of state that combines a fit to experimental data with the results of molecular dynamics calculations to cover a very wide range of temperatures and pressures (Zhang & Duan 2005). It can be seen in Figure 12a that the Zhang & Duan equation is closely consistent with experimental data from 100 to 900°C and 1.0 to 9.0 kbar (Burnham et al. 1969). It is also in reasonable agreement with PVT data obtained from synthetic fluid inclusion studies (Brodholt & Wood 1994; Withers et al. 2000) up to about 1600°C and 40.0 kbar. In this regard, it is more accurate than the IAPWS (International Association for the Properties of Water and Steam) formulation (Wagner & Pruβ 2002), which did not consider the datasets shown in Figure 12a and b (Sverjensky et al. 2014a).

The strong increases in the molar volume of water with increasing temperature at 1.0–3.0 kbar shown in Figure 12a are magnified greatly at even lower pressures (Fig. 13a). It can be seen in Figure 13a that the increasingly steep changes of the volume with temperature as pressure is lowered to 0.5 kbar echoes the behaviour of the saturation curve as it heads up to the critical point of water. The extreme behaviour of the volume of water near its critical point is therefore manifested at higher temperatures and pressures in the supercritical fluid regime. However, the rapid changes with temperature and pressure die out by pressures of 10 kbar.

The striking behaviour of the molar volume of water near its critical point is also seen in the calculated values of the molar heat capacity of water in Figure 13b, with the additional complexity that at the lowest pressures in the supercritical region, the heat capacity not only increases steeply with temperature, but also maximizes (e.g. at 0.5 kbar) and then plunges steeply down.

The above properties of water in the vicinity of the critical point are reflected in the standard partial molal properties of dissolved electrolytes. It can be seen in Figure 14a and b that experimental measurements of the standard partial molal volume and heat capacity of strong electrolytes such as NaCl show extreme behaviour in the vicinity of the critical point of water. This behaviour for electrolytes contrasts dramatically with the behaviour of minerals illustrated above. It necessitated the development of equations of state that could describe large fluctuations in the standard partial molal properties of dissolved salts in water.

Equations of state for the standard partial molal volumes and heat capacities of aqueous species

An equation of state for describing the standard partial molal volumes of electrolytes was developed by Helgeson & Kirkham (1976). Three contributions to the standard partial molal volume of an ion in water were hypothesized: an intrinsic contribution (V¯i,j0) characteristic of the ion itself, a collapse contribution (ΔV¯c,j0) to describe the consequences of the alteration of solvent structure around the ion, and a solvation contribution (ΔV¯s,j0) to describe the consequences of orientation of water dipoles around the ion. The three contributions are summed as shown in the equation 
V¯j0=V¯i,j0intrinsic+ΔV¯c,j0collapse+ΔV¯s,j0solvation.
(80)
The intrinsic and the collapse terms can be summed to give the non-solvation volume (ΔV¯n,j0) according to 
ΔV¯n,j0=V¯i,j0+ΔV¯c,j0.
(81)
Consequently, 
V¯j0=ΔV¯n,j0nonsolvation+ΔV¯s,j0solvation.
(82)
A similar set of three contributions was subsequently proposed (Helgeson et al. 1981) for the standard partial molal heat capacity according to 
C¯P;j0=C¯P;j,i0intrinsic+ΔC¯P;j,c0collapse+ΔC¯P;j,s0solvation
(83)
 
=ΔC¯P;j,n0nonsolvation+ΔC¯P;j,s0solvation.
(84)
We now examine the specific functional forms proposed for these contributions to the standard partial molal volume and heat capacity of an ion. We start with the solvation contribution.

Born equation for solvation of an aqueous ion

The solvation contributions to the standard partial molal volumes and heat capacities in equations (71) and (73) are built on an application of the Born equation (Helgeson & Kirkham 1976). The electrical work associated with orientation of water dipoles around the jth ion can be expressed as the standard partial molal Gibbs free energy (ΔG¯s,j0) given by 
ΔG¯s,j0=ωj1εH2O1
(85)
where ωj represents the conventional Born solvation coefficient of the aqueous ion and εH2O represents the dielectric constant of pure water. Equation (73) represents a quantification of the importance of the dielectric constant of water discussed above. It involves only the dielectric constant of pure water because it is a standard state solvation contribution. In real fluids containing mixed solvents, a contribution from the dielectric constant of the mixed solvent should be used in the non-standard state, extended-term contributions to the aqueous activity coefficient; for example, an additional term in equation (23), as discussed above.
The conventional Born solvation coefficient for the jth ion is related to the absolute Born solvation coefficient (ωjabs.) by 
ωj=ωjabs.ZjωH+abs.
(86)
where ωH+abs.=0.5387×105 cal mol−1 refers to the absolute Born coefficient of the H+ ion (ωH+abs.) and Zj. refers to the charge on the ion.
The absolute Born solvation coefficient for the jth aqueous ion can in turn be related to a radius for the jth ion according to 
ωjabs.=ηZj2re,j1εH2O1
(87)
where η is a constant equal to 1.66027×105 Å cal mol−1 and re,j represents the effective electrostatic radius of the aqueous ion (Å). According to Helgeson & Kirkham (1976), re,j is an electrostatic parameter that represents the distance from the centre of the ion to the location where the water molecules around the ion have the properties of bulk water.
The beauty of this approach is that values of re,j for monatomic ions can be estimated from crystallographic radii. The effective electrostatic radius for the jth aqueous ion has been related to the crystallographic radius by the equation 
re,j=rx,j+|Zj|(kZ+g)
(88)
where rx,j represents the crystallographic radius of the jth ion (Shock & Helgeson 1988), Zj again represents the charge on the jth ion, kZ is equal to 0.0 for anions and 0.94 for cations based on the analysis of the temperature dependence of the standard state partial molal volumes of aqueous electrolytes (Helgeson & Kirkham 1976) and g(P, T) designates a solvent function of temperature and density (Shock et al. 1992).

Analysis of a wide variety of electrolyte standard partial molal properties, primarily as functions of temperature (Fig. 14a and b) together with analysis of the dissociation constant of NaCl0 (Fig. 15a) have led to the conclusion that g varies most strongly at temperatures and pressures in the vicinity of the critical point of water; for example, at temperatures of 250–500°C and pressures less than 2000 bar (Fig. 15b).

However, at T < 200°C or at P > 4.0 kbar, the properties of electrolytes and the dissociation constant of NaCl0 are consistent with the assumption that g = 0.0 (Shock et al. 1992). Under these conditions, equation (88) simplifies to 
re,j=rx,j+|Zj|kZ
(89)
and consequently ωj is independent of temperature and pressure for T < 200°C or P > 4.0 kbar. In particular, it follows that for T < 200°C or Psat conditions, where most of the experimental data for the standard partial molal properties of electrolytes have been obtained, we can differentiate the Born equation (equation (85)) assuming that ωj is independent of temperature and pressure to obtain the standard partial molal volumes, compressibilities and heat capacities of solvation from the equations 
ΔV¯s,j0=ωj;Pr,TrQP,T
(90)
 
Δκ¯s,j0=ωj;Pr,TrNP,T
(91)
 
ΔC¯P;s,j0=ωj;Pr,TrTXP,T
(92)
where QP,T, NP,T and XP,T are partial derivatives of εH2O given by 
QP,T=1εlnεPT
(93)
 
NP,T=1ε2lnεPTlnεPT2
(94)
 
XP,T=1ε2lnεTPlnεTP2.
(95)

Equations of state for the standard partial molal properties of an aqueous species at T < 200°C and Psat

The non-solvation contributions to equations (82) and (84) were first developed by analysis of experimental data for the temperature dependence of the standard partial molal properties of electrolytes using experimental data referring to T < 200°C and Psat and led to the first version of equations of state for the non-solvation contributions (Helgeson & Kirkham 1976; Helgeson et al. 1981). With the availability of higher temperature and higher pressure experimental data, these equations were subsequently revised (Tanger & Helgeson 1988) and can be written as follows for ions or neutral species.

For the non-solvation volume 
ΔV¯n;j0=σj+ξjTθ
(96)
where σ and ξ are temperature-independent coefficients characteristic of the jth ion and the pressure dependence is given by 
σj=a1,j+a2,jP+ψ
(97)
and 
ξj=a3,j+a4,jP+ψ
(98)
where θ and ψ are solvent constants with values of 228 K and 2600 bars, respectively. Taking the partial derivative with respect to pressure leads to the non-solvation compressibility given by 
Δκ¯n;j0=a2,j+a4,j(Tθ)1P+ψ2.
(99)
For the non-solvation heat capacity 
ΔC¯P,n;j0=c1,j+c2,j(Tθ)2.
(100)
Combining the non-solvation and the solvation contributions, the full equations of state for the standard partial molal volume, compressibility and heat capacity under conditions where ωPr,Tr;j is a constant can be written 
V¯j0=σj+ξjTθωPr,Tr;jQP,T
(101)
 
V¯j0=a1,j+a2,jψ+P+a3,j+a4,jψ+P1TθωPr,Tr;jQP,T
(102)
 
κ¯j0=a2,j+a4,j(Tθ)1P+ψ2+ωPr,Tr;jNP,T
(103)
and 
C¯P;j0=c1,j+c2,j(Tθ)2+ωPr,Tr;jTXP,T.
(104)

Relationships between the properties of ions and electrolytes

The properties of the individual ions cannot be measured directly. Instead, experimental measurements are made on dissolved electrolytes from which the individual ion properties can be obtained through the following approach. For the kth 1:1 electrolyte consisting of the ith cation and the lth anion, the standard partial molal volume of the electrolyte is, by definition, equal to the sum of the standard partial molal volumes of the ions; that is, 
V¯k0=V¯i0+V¯l0.
(105)
For the kth 2:1 electrolyte, 
V¯k0=V¯i0+2V¯l0
(106)
and so on. General expressions for all electrolytes have been given by Helgeson et al. (1981).
It follows from equation (105) and all similar equations that 
σk=σi+σl
(107)
and 
ξk=ξi+ξl
(108)
and that similar equations apply to all the equation of state coefficients.
Together with the H+ ion convention that the standard partial molal properties of formation from the elements at all temperatures and pressures are exactly equal to zero it follows that the equation of state coefficients of the H+ ion are also all exactly equal to zero. Because the standard partial molal properties of the aqueous HCl electrolyte can be written in terms of the sum of the properties of H+ and Cl, it follows, in the case of the volume of HCl, that 
V¯HCl0=V¯H+0+V¯Cl0=V¯Cl0(conventional).
(109)
In other words, the H+ ion convention results in the conventional properties of the individual Cl ion being exactly equal to those of the electrolyte HCl measured experimentally. More generally, the relationship between the conventional and absolute properties of the jth ion with a charge of Zj are expressed by the equation, in the case of volumes, as 
V¯j0=V¯j0abs.ZjV¯H+0abs.
(110)
where V¯H+0abs. represents the absolute standard partial molal volume of the H+ ion. Expressions analogous to those in equations (109) and (110) exist for all the individual standard partial molal properties of ions.

Regression equation for the volume of an aqueous electrolyte at T < 200°C and Psat

For the restricted temperatures and pressures being considered, experimental measurements of the kth electrolyte (V¯k0) as a function of temperature can be used to obtain ΔV¯n0(expt.) as a function of temperature using the equation 
ΔV¯n;k0=V¯k0(expt.)ΔV¯s;k0(calc.)
(111)
and predicted values of ΔV¯S0 from equation (90) and estimated values of ωPr,Tr;j. Values of ΔV¯n;k0 obtained in this way are not sensitive to the prediction of ΔV¯s;k0 because at such low temperatures ΔV¯s;k0 is small. Based on equation (96), we can then expect that ΔV¯n;k0(expt.) will show a linear dependence on (Tθ)1 according to the equation 
ΔV¯n;k0(expt.)=σk+ξkTθ.
(112)
Consistency with equation (112) allows regression to obtain values of σk and ξk for an electrolyte. Examples are shown in Figure 16 for NaCl and MgCl2 (Tanger & Helgeson 1988). Combination of the results of this analysis with a corresponding analysis of compressibility data, discussed next, allows the equation of state coefficients for the pressure dependence of σk and ξk to be obtained.

Regression equation for the compressibility of an aqueous electrolyte at T < 200°C and Psat

By analogy with the procedure for volumes described above, experimental measurements of κ¯0 as a function of temperature can be used to obtain Δκ¯n0(expt.) as a function of temperature using the equation 
Δκ¯n;k0(expt.)=κ¯k0(expt.)Δκ¯s;k0(calc.).
(113)
We can then expect that Δκ¯n0(expt.) will show a linear dependence on (Tθ)1 according to the equation 
Δκ¯n;k0=a2,k+a4,k(Tθ)1P+ψ2.
(114)
Consistency with equation (114) allows a linear regression to obtain values of a2,k and a4,k for an electrolyte. Examples are shown in Figure 17 for NaCl and MgCl2 (Tanger & Helgeson 1988).

The values of a2,k and a4,k obtained from the temperature dependence of the compressibility can then be used with the values of σk and ξk obtained previously from the temperature dependence of the volume, together with equations (97) and (98), to generate values of a1,k and a3,k. In this way, a complete set of equation of state coefficients is available for prediction of the standard partial molal volume at temperatures and pressures outside the experimentally studied range.

Calibration of the volume behaviour of an aqueous electrolyte at pressures greater than Psat

The pressure dependence of the equation of state coefficients σk and ξk was revised by Tanger & Helgeson (1988) to incorporate an asymptotic dependence on pressure with a low-pressure limiting value given by the parameter ψ=2600 bar. This value was obtained from the only available data at the time for the volumes of NaCl and K2SO4 at 25°C and pressures up to 10 kbar (Adams 1931, 1932). The regression curve for the NaCl data is shown in Figure 18. Experimental density data are now available from diamond anvil cell studies (Mantegazzi et al. 2012, 2013) at 100–400°C and pressures up to 4.0 GPa for several electrolytes. However, at the concentrations studied, 1.0 m and greater, analysis of the data requires a model for not only the standard state volumes of the electrolytes but also the non-standard state properties and those of the complexes.

Regression equation for the standard partial molal heat capacity of an aqueous electrolyte using data at T < 200°C and Psat

Again proceeding as above for the analysis of the volume and compressibility data at low temperatures, experimental measurements of C¯P;k0 as a function of temperature can be used to obtain ΔC¯P,n;k0(expt.) as a function of temperature using the equation 
ΔC¯P,n;k0=C¯P;k0(expt.)ΔC¯P,s;k0(calc.).
(115)
We can then expect that ΔC¯P,n;k0(expt.) will show a linear dependence on (Tθ)2 according to the equation 
ΔC¯P,n;k0(expt.)=c1,k+c2,k(TΘ)2.
(116)
Consistency with equation (116) allows regression to obtain values of c1,k and c2,k for the electrolytes NaCl and MgCl2 shown in Figure 19 (Tanger & Helgeson 1988).

The properties of many additional types of aqueous species including inorganic and organic salts and neutral species have been analysed in this way (Shock & Helgeson 1988, 1990; Shock et al. 1989, 1997; Amend & Helgeson 1997; Plyasunov & Shock 2001; Dick et al. 2006; LaRowe & Helgeson 2006). The results are incorporated into the datafiles for SUPCRT92 and the DEW model. It should be noted, however, that a full set of experimental data needed to derive all the equation of state coefficients is available for only a very restricted number of aqueous species. The rest all involve estimation techniques to obtain all the equation of state coefficients. Such techniques provide a unique predictive power to the overall HKF approach (see below).

Full equations of state for the standard partial molal volume and heat capacity of individual aqueous species over wide ranges of temperature and pressure

The regression procedures discussed above were applied at temperatures <200°C and pressures at or a few hundred bars above Psat to obtain equation of state coefficients for many aqueous electrolytes and neutral species. The equation of state coefficients so obtained can be used in the full equations for the standard partial molal properties of aqueous species over wide ranges of temperature and pressure. For example, the V¯j0 and C¯P;j0 of the jth aqueous species are given by 
V¯j0=a1,j+a2,jP+ψ+a3,j+a4,jP+ψ1TθωjQT,P+1ε1ωjPT
(117)
 
C¯P;j0=c1,j+c2,j(Tθ)22T(Tθ)3a3,j(PPr)+a4,jlnP+ψPr+ψ+ωjTX+2TYT,PωjTPT1ε12YT,P2TP
(118)
where the derivatives of ωj as a function of temperature and pressure depend on the g function discussed above, as well as the partial derivatives of the dielectric constant of water.

Equation for the apparent standard partial molal Gibbs free energy of formation of an aqueous ion as a function of temperature and pressure

The primary use of the equations of state for the standard partial molal heat capacity and volume is their integration to allow calculation and prediction of the standard partial molal free energies of individual aqueous species of all kinds, from simple, monatomic ions, to neutral species, complexes and aqueous organic molecules. The integrated free energy equation yielding the apparent standard partial molal Gibbs free energy of formation for an aqueous species is given below: 
ΔG¯f,j;P,T0=ΔG¯f,j;Pr,Tr0S¯j;Pr,Tr0(TTr)c1,jTlnTTrT+Trc2,j{1Tθ1TrθθTθTθ2lnTr(Tθ)T(Trθ)}+a1,j(PPr)+a2,jlnψ+Pψ+Pr+a3,j1TθlnPPrTθ+a4,j1Tθlnψ+Pψ+Pr+ωj1ε1+ωj;Pr,Tr1εPr,Tr1+ωj;Pr,Tr{YPr,Tr(TTr)
(119)
where the coefficients c1,j, c2,j, a1,j, a2,j, a3,j, a4,j and ωj are the seven equation of state coefficients for the aqueous species. Knowing the seven equation of state coefficients therefore allows prediction of apparent standard free energies over wide ranges of pressure and temperature. The beauty of this approach is that the equation of state coefficients can be obtained from regression of experimental low-temperature and -pressure volumes, compressibilities and heat capacities (i.e. at T < 200°C and Psat) and used to predict free energy changes at elevated temperatures and pressures. In this approach, uncertainties in the predicted free energies are minimized because the equation of state coefficients were calibrated using the derivative properties of the free energy; that is, volumes, compressibilities and heat capacities. Even when the uncertainties in the derivative properties are substantial, the integrated free energy changes are less sensitive to uncertainty because the character of the free energy function with temperature and pressure is always less pronounced than in the derivative properties, which show more extreme behaviour with temperature and pressure (Helgeson et al. 1981; Shock & Helgeson 1988). This point has not been sufficiently recognized in the subsequent literature.

For most electrolytes, however, experimental low-temperature and -pressure volumes, compressibilities and heat capacities as functions of temperature have not been measured. A similar situation applies to aqueous neutral inorganic and organic species, and all but a few aqueous metal complexes. Complete sets of equation of state coefficients from experimental V¯T0, κ¯T0 and C¯Pr,T0 are very few. Therefore, to be able to carry out geochemical modelling, we need ways to estimate the equation of state coefficients for aqueous species. A great strength of the HKF approach is that it was built with the need for estimation algorithms in mind.

Three situations can be imagined when there is a need for estimating equation of state coefficients for aqueous species, as follows: (1) partial sets of experimental data for V¯T0 and C¯Pr,T0 as functions of temperature are available, or (2) experimental data for V¯Pr0 and C¯Pr0 at 25°C only are available, or (3) no experimental data are available.

Estimation of equation of state coefficients of aqueous species

The basis of estimating the equation of state coefficients of aqueous species rests on empirical correlations with the standard partial molal for V¯Pr0 and C¯Pr0 at 25°C and 1 bar. The first step involves evaluating the solvation volume and heat capacity using equations (90) and (92) in which the Born solvation functions Q and X refer to 25°C and 1 bar: Q=0.05903×105bar1 and X=0.0309×105K2; that is 
ΔV¯s0=ωj;Pr,Tr(0.05903×105×41.84)cm3mol1
(120)
and 
ΔC¯P;s0=ωj;Pr,Tr(298.15)(0.0309×105)calmol1K1
(121)
where ωj;Pr,Tr is in calmol1.

Using estimated values of ωj;Pr,Tr allows estimation of the solvation volume and heat capacity.

Values of ωj;Pr,Tr can be obtained from predictive correlations with S¯j;Pr,Tr0 as can be seen in Figure 20.

The equations of the correlation lines in Figure 20a and b are as follows.

Aqueous inorganic or organic ions: for aqueous ions 
ωj;Pr,Tr=(1514.4)S¯j;Pr,Tr0+βZ
(122)
where for cations 
βZ=0.544×105Zjcalmol1
(123)
and for anions 
βZ=1.62×105Zjcalmol1.
(124)
In contrast, for neutral species that are typically crystalline at ambient conditions when not dissolved in water 
ωj;Pr,Tr=+0.30×105calmol1.
(125)
For neutral species that are typically gases or liquids at ambient conditions when not dissolved in water 
ωj;Pr,Tr=0.30×105calmol1.
(126)
It should be noted that the equations for neutral species are no more than a first approximation, to be used in the absence of experimental data for the temperature dependence of the standard partial molal volumes or heat capacities at temperatures greater than 200°C, or when no data are available for the temperature dependence of the equilibrium constant for a reaction involving the neutral species at pressures less than about 2.0 kbar and temperatures >300°C.

Estimation of the volume equation of state coefficients σ, a1, a2 and a4 for aqueous species

After analysing numerous sets of experimental data for a variety of aqueous electrolytes, neutral species and aqueous organic species, a comprehensive approach was taken to develop correlations for estimating volume equation of state coefficients (Shock & Helgeson 1988, 1990; Shock et al. 1989). After estimating the solvation volume at 25°C and 1.0 bar (ΔV¯s0) as described above, an experimental volume at 25°C and 1.0 bar is used to calculate the non-solvation volume (ΔV¯n0). The value of ΔV¯n0 is then used to estimate the parameter a1. This step is a key part of the approach for predicting aqueous species free energies at upper mantle pressures. It can seen in the full free energy equation (equation (119)) that the term a1(PPr) appears. At pressures of 20 000 bar and higher, this term dominates the change of ΔG¯f;P,T0 with pressure. It is consequently very important to have as accurate an estimate of a1 as possible. The correlation between a1 and ΔV¯n0 currently used in the DEW model differs from that originally proposed by Shock & Helgeson (1988), which was based only on data for Mg2+, Na+, K+, Cl and Br. Using constraints derived from an interpretation of Raman speciation data for water in equilibrium with aragonite at pressures of 20–60 kbar (Facq et al. 2014), separate correlations for aqueous ions with charges of 1 or 2 were developed, as well as a correlation for all neutral species and aqueous metal-complexes (Figs 21a and b).

The equations of the correlation lines in Figure 21a–d are as follows.

For aqueous inorganic or organic ions 
Z=+1or1,a1×10=(0.10871)ΔV¯n02.0754
(127)
 
Z=+2or2,a1×10=(0.23999)ΔV¯n03.5151.
(128)
For aqueous complexes or neutral inorganic or organic species 
a1×10=(0.19421)ΔV¯n01.5204.
(129)
For aqueous species of all kinds 
σ=1.11ΔV¯n0+1.8.
(130)
For aqueous species of all kinds 
a4=4.134a227790
(131)
where σ and ΔV¯n0 are in cm3 mol−1, a1 is in cal mol−1 bar−1, a2 is in cal mol−1 and a4 is in cal K mol−1.

The procedure is as follows. After estimation of a1 from Figure 21a or b, σ can be estimated using Figure 21c. The coefficient a3 is then calculated from the overall equation of state for the volume at 25°C and 1.0 bar, the value of V¯0, and the coefficients a1, a2, a4 and ω.

Estimation of the heat capacity equation of state coefficients c1 and c2 for aqueous species

Based on the analysis of a large amount of experimental data for the standard partial molal heat capacities of a wide range of aqueous species, a correlation between c2 and C¯Pr0 at 25°C and 1.0 bar was developed for aqueous ions, including AlO2, an abbreviation for Al(OH)4 (Shock & Helgeson 1988). The same correlation was found to apply to the inorganic neutral species NH03, H2S0, B(OH)30 and CH3OOH0 (Shock et al. 1989). However, other aqueous organic species follow very strong but different correlations. Examples included here are shown in Figure 22b: one line for alkanes and aromatic molecules, and another for straight chain alcohols (Shock & Helgeson 1990; Shock 1995; Plyasunov & Shock 2001).

Based on the above correlations we estimate values of c2 using the following equations.

For aqueous inorganic ions and neutral species, organic ions and metal complexes 
c2×104=(0.2037)C¯Pr03.0346.
(132)
For aqueous alkanes and aromatic hydrocarbons 
c2×104=(0.17859)C¯Pr04.7893.
(133)
For aqueous straight chain alcohols 
c2×104=(0.068979)C¯Pr02.2842
(134)
where c2 is in cal mol−1 K−1 and C¯Pr0 is in cal mol−1 K−1.

Using these correlations, values of c1 can be calculated from equation (104) and estimated values of the solvation heat capacity contribution can be calculated from equation (121) above.

Summary of predictive algorithms for aqueous ions; estimation of equation of state coefficients

The greatest strength of the HKF approach is its predictive power. This strength rests on the correlations for estimating equation of state coefficients for aqueous species of all kinds, which in turn rest on the analysis of a large amount of experimental data for the standard partial molal volumes, compressibilities and heat capacities as functions of temperature and pressure. When experimental values of volumes, compressibilities and heat capacities are known only at 25°C and 1 bar, the correlations summarized above allow estimation of all the equation of state coefficients of an aqueous species.

A summary diagram of the algorithm for estimating equation of state coefficients for aqueous species is given in Figure 23.

As described above, the estimation algorithm is useful for filling in gaps in experimental data or for estimating equation of state coefficients when no experimental data for volumes, compressibilities and heat capacities as functions of temperature are known. However, a very common situation that arises is the case where volumes, compressibilities and heat capacities even at 25°C and 1.0 bar are not known. This is generally the situation for aqueous metal complexes of all kinds, inorganic and organic. To handle this, estimation algorithms have been developed to generate estimates of the volumes, heat capacities and entropies at 25°C and 1.0 bar of many kinds of complex species (Shock et al. 1997; Sverjensky et al. 1997). Consequently, a completely predictive approach is possible.

At upper mantle pressures, however, new kinds of complexes can occur. To establish the appropriate volumes, heat capacities and entropies of these species requires analysis of the experimental data available, which are primarily solubility measurements, either of single minerals or of complex synthetic rocks. Experimental aqueous speciation data are an invaluable aid in developing realistic solubility and speciation models, but are few and far between at pressures greater than about 20 kbar. Theoretical molecular calculations identifying the nature of species at high pressures are similarly valuable, but also lacking for the most part. Nevertheless, advances can be made in our understanding of complexes in high-pressure fluids by taking advantage of the predictive algorithm summarized above. For example, if solubility data are known over a range of pressures and temperatures, the predictive algorithm in Figure 23 can be used to constrain values of the entropy, volume and heat capacity of a complex at 25°C and 1.0 bar by regression of the solubility data. In this way, only a small number of parameters need to be established by the regression. An example will be discussed below.

II. Analysis of experimental solubility data

Activity coefficient approximations

The full equation for the activity coefficient of the jth ion is 
logγ¯j;P,T=AγZj2I¯0.51+akBγI¯0.5+bγ,kI¯+Γγ.
(135)
Values of Aγ and Bγ in the first term on the right-hand side can be evaluated from the density and dielectric constant of water at elevated temperatures and pressures (Helgeson & Kirkham 1974a). However, the application of an extended term equation of state to calculating bγ,k as functions of temperature and pressure has not yet been developed for upper mantle pressures. Consequently, a first approximation to activity coefficients of ions can be made by setting bγ,k=0.0, resulting in the equation 
logγ¯j;P,T=AγZj2I¯0.51+akBγI¯0.5+Γγ.
(136)
Equation (136) has been used recently to analyse solubility measurements of synthetic rocks at high temperatures and pressures (Huang 2017).
For neutral species such as CO2, CO and CH4 in fluids at high temperatures and pressures, an activity coefficient is needed to account for the strong non-ideality apparent in the traditional CO2–H2O models (Duan & Zhang 2006). Here we present a first approximation obtained by transforming the Zhang & Duan CO2–H2O model to the hypothetical 1.0 molal standard state, resulting in a semi-empirical equation of form 
logγ¯n;P,T=bP,T+Γγ
(137)
where bP,T is a constant fitted to the transformed Zhang & Duan CO2–H2O model and Γγ again refers to the sum of the molalities of all aqueous species in the fluid (i.e. ions, neutral species and complexes).
Finally, the activity of water in fluids at high temperatures and pressures will be affected by all the species dissolved in the water, including ions and neutral species (e.g. CO2). The non-ideality of water in such fluids can be accounted for by the expression 
aH2O=λH2OXH2O.
(138)
For simplicity here, we use 
log(aH2O)=logXH2O=log55.5155.51+m
(139)
where m* again represents the sum of the molalities of all the aqueous species in the solution.

Hydration of aqueous complexes at high temperatures and pressures in fluids

By convention in the predictive algorithm for estimation of the equation of state parameters of aqueous species summarized above we use the stripped-down aqueous species (i.e. without water in them); for example, in place of the hydrated aqueous complex Mg(OH)20 we use MgO0, which differs only by 1.0 H2O and allows us to write the reaction 
MgO0+2H+=Mg2++H2O.
(140)
When the species MgO0 is used in an aqueous speciation model and when the activity of water is equal to 1.0, the above equation is exactly equivalent to 
Mg(OH)20+2H+=Mg2++2H2O.
(141)
However, when the activity of water is less than 1.0, the two reactions have a different dependence on the activity of water, which affects the results of aqueous speciation and solubility calculations. Similar considerations apply to all aqueous OH-bearing species.

Retrieval of metal-complex equilibrium constants and equation of state characterization

Metal complexes at upper mantle pressures

Ion-association is typically promoted by increasing temperature and opposed by increasing pressure, as discussed above in the section on the properties of water. The typical behaviour has been extensively investigated at crustal, hydrothermal temperatures and pressures (Sverjensky et al. 1997; Brugger et al. 2016). However, at upper mantle conditions, the pressure is so much higher than in the crust that it seems possible that perhaps dissociation would dominate. However, several lines of evidence suggest otherwise. Theoretical evidence from molecular models suggests that complex species are very important, partly because of dissociation of simpler species to provide ligands that can complex with metals in ways that are not important at crustal conditions. For example, the speciation of Mg2+ and SO42 in water has been studied over a wide range of conditions (Jahn & Schmidt 2010) and indicates a transition from Mg2+ and HSO4 ions at low densities to polymeric species involving completely dissociated sulphate complexed with Mg (Fig. 24).

Similarly, Raman spectroscopy of aqueous fluids in equilibrium with aragonite (Facq et al. 2014) reveals the expected strong dissociation of HCO3 to CO32 with increasing pressure (Fig. 25), and in Na2CO3 fluids, ab initio molecular dynamics modelling (Pan & Galli 2016) revealed the formation of NaHCO30, Na2HCO3+, NaCO3 and Na2CO30. The latter study also revealed that CO20 as a molecule is unimportant compared with H2CO30 at high temperatures and pressures.

Because the simple species at crustal conditions tend to dissociate at upper mantle conditions, other possible ligands at high pressures can be expected to be important in upper mantle fluids. For example, aqueous silica dissociates according to 
Si(OH)40=H++SiO(OH)3.
(142)
A reaction that has the same equilibrium constant is shown in Figure 26a. It can be seen that the logK increases dramatically with pressure, suggesting that the silicate anion ought to be a significant ligand with metals in upper mantle fluids.
Similarly, water itself dissociates strongly with increasing pressure (Fig. 26b) according to the reaction 
H2O=H++OH
(143)
indicating the availability of the OH ion as a potential complex-forming ligand.

And finally, aqueous organic acid anions such as acetate, formate and propionate are predicted to have true thermodynamic predominance fields at pressures greater than about 30 kbar (Fig. 26c).

All the above results strongly suggest that the formation of complexes involving ligands such as sulphate, carbonate, silicate, hydroxide and organic species might be possible in upper mantle fluids. Further support for this comes from an analysis of aqueous Al-speciation in upper mantle fluids, which is discussed below.

Example: solubility and complexing of Al in complex upper mantle fluids

An incredibly useful feature of the HKF approach used in the DEW model is the predictive algorithm discussed above relating the standard partial molal properties of an aqueous species to the equation of state coefficients. For most aqueous metal complexes, we do not, and probably will not, have measurements of V¯0, κ¯0 and C¯P0 as functions of temperature to calibrate the equation of state coefficients. However, the correlations in Figures 21 and 22 and accompanying equations link the partial molal properties to the equation of state coefficients c1, c2, a1, a2, a3, a4 and ω. It therefore becomes possible to regress equilibrium constants as functions of temperature and pressure to directly obtain values of V¯0, C¯P0 and S¯0 referring to 25°C and 1.0 bar. If the equilibrium constants are known over a wide enough temperature and pressure range (e.g. at least 300°C and 10 kbar) a unique, well-constrained set of V¯0, C¯P0 and S¯0 can be retrieved, plus all the equation of state coefficients linked by the correlation and equations above. An alternative approach would be to regress the equilibrium constants directly for all seven equation of state coefficients c1, c2, a1, a2, a3, a4 and ω. Typically, this will result in over-fitting the experimental data. Furthermore, the equation of state coefficients will not be consistent with the many dozens of sets that are based on experimental data. Consequently, this approach is not recommended.

An example of the recommended approach is shown in Figure 27a and b. The data shown for Al solubility in equilibrium with a synthetic K-free eclogite (Kessel et al. 2005b) were fitted with a speciation model that simultaneously fitted data for Na, Mg, Ca, Fe and Si (Huang 2017). Included in the model speciation for Al were Al3+, Al(OH)30 and Al(OH)4, the last two of which have been calibrated previously (Sverjensky et al. 2014a). Al(OH)30 was calibrated using solubility data for corundum at low temperatures and pressures (Pokrovskii & Helgeson 1995), as well as at elevated temperatures and pressures (Becker et al. 1983; Tropper & Manning 2007). However, none of these species proved to be significant contributors to the solubility of Al in the more complex fluid in equilibrium with the K-free eclogite. Instead, after some trial regressions, an Al–Si complex previously studied at Psat up to 300°C (Tagirov & Schott 2001) was found to account for the data well (Fig. 27a) according to the reaction 
Al(OH)3OSi(OH)3+4H+=Al3++Si(OH)40+3H2O
(144)
which is equivalent to the reaction shown in the figure as 
AlO2SiO2+4H+=Al3++SiO20+2H2O.
(145)

Furthermore, the equilibrium constants for this reaction, which were retrieved from solubilities at 40 and 50 kbar, appear to be closely consistent with the equilibrium constants reported by Tagirov & Schott (2001). The overall dataset was fitted with the predictive algorithm to retrieve values of V¯0, C¯P0 and S¯0 referring to 25°C and 1.0 bar, which is only a three-parameter fit over a temperature range of about 1000°C and a pressure range of 50 kbar. Most importantly, the same equilibrium constants provide excellent predictions of the solubility of Al in two different peridotitic fluids (Huang 2017). The fact that the peridotitic fluids have a very different chemistry from the eclogitic fluids provides strong support for the equilibrium constants of the Al–Si complex shown in Figure 27b.

Concluding remarks

The approach described above provides a quantitative, thermodynamic basis for modelling fluid–rock interactions under upper mantle conditions. It builds on an understanding of how to predict equilibrium and non-equilibrium reactions at ambient conditions, then it provides the means to extend predictive modelling from ambient conditions to upper mantle conditions. The principal data needed to do this are equilibrium constants for aqueous reactions and mineral–aqueous reactions at elevated temperatures and pressures, as well as estimations of the relationships between thermodynamic activities and measurable concentrations. Whereas approximations of the latter can be made, and will continue to be improved, the biggest source of uncertainty in the modelling of upper mantle fluids is the equilibrium constants for the principal species governing solubilities (Helgeson et al. 1981); that is, metal-complexes, which appear to involve inorganic ligands such as Cl, HCO3, CO32, SiO(OH)3, OH and SO42, as well as organic ligands such as HCOO, CH3OO and CH3CH2OO. Many of the complexes that are important in upper mantle fluids may not be important at crustal conditions.

The solubility of the major elements in upper mantle fluids is greatly enhanced because of complexes involving such ligands. The interpretation of mineral and rock solubilities is greatly facilitated by the predictive nature of the HKF approach embedded in the DEW model (http://www.dewcommunity.org/). By facilitating the development of equations of state for the complexes, the DEW model allows extrapolation of the properties of complexes derived from experimental temperatures, pressures and bulk compositions different from those used for calibration of the equations of state. A specific example of how this is done is discussed above. The practical value of the approach is that it is easy to implement, does not involve over-fitting the equation of state to experimental data, and that it affords extrapolations that can be experimentally tested or compared with geological observations.

However, it is still highly uncertain as to the precise nature of upper mantle ligands and complexes because of the dearth of fundamental experimental spectroscopic data on speciation in high-density aqueous fluids. Raman spectroscopy has been used to quantify aqueous carbonate and bicarbonate species (Facq et al. 2014, 2016; Schmidt 2014), metal–sulphate species (Frantz et al. 1994; Schmidt & Manning 2017), and silicate polymerization (Zotov & Keppler 2000, 2002; Mysen 2010; Hunt et al. 2011; Mysen et al. 2013). X-ray absorption near edge structure (XANES) and X-ray absorption fine structure (XAFS) studies have characterized aspects of Zr–silicate complexes (Louvel et al. 2013; Sanchez-Valle 2013; Sanchez-Valle et al. 2015). Nuclear magnetic resonance (NMR) studies are now being extended to high-density fluids with the potential to help identify the nature of metal complexes (Pautler et al. 2014; Ochoa et al. 2015; Augustine et al. 2017). Many more such studies are needed.

Theoretical ab initio molecular dynamics calculations are also starting to have a major impact on our ideas of speciation in upper mantle fluids. For example, it now appears that H2CO30 is much more important than CO20 in upper mantle fluids (Pan & Galli 2016). Polymeric Mg–SO4 species are more important than simple ions and ion-pairs in high-density fluids (Jahn & Schmidt 2010). These few examples complement the much larger set of studies carried out for crustal hydrothermal conditions (e.g. Brugger et al. 2016). But many more studies are needed to identify the nature of metal-complexing under upper mantle conditions.

Acknowledgements

The existence of this paper was catalysed by the invitation from M. L. Frezzotti (Università degli Studi di Milano–Bicocca, Milano, Italy) to present lectures as part of the School for PhD students on ‘Carbon forms, paths and processes in the Earth’, at Villa Grumello, Lake Como, Italy held in October 2017. The author is grateful for the help and support of the Johns Hopkins University and wishes to acknowledge helpful discussions with C. Manning (UCLA), M. Ghiorso, E. Shock (ASU), F. Huang, and J. Hao, and a detailed reading of the manuscript by J. Huang (Johns Hopkins University). S. Abelos (JHU) provided outstanding help with the figures and equations.

Scientific editing by Maria Luce Frezzotti

This is an Open Access article distributed under the terms of the Creative Commons Attribution 4.0 License (http://creativecommons.org/licenses/by/4.0/)