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, 50b; Hacker et al. 2003a,32b; 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,52b, 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,38b, 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
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.
Evaluation of equation (4) leads to a predictive capability: if , the reaction should proceed as written. However, if , the reaction should proceed in reverse.
Prediction of equilibrium reactions
Definitions of standard states
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 as , as depicted in Figure 3.
Water and mixed solvents
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, as .
More generally, for a component MZ in a multicomponent fluid, such as metamorphic fluids, the activity of 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 and other neutral molecules that can be thought of as mixing with at elevated temperatures and pressures. It should be noted that the standard Gibbs free energy of formation for that should be used for this standard state must be that of pure fluid , completely different from the standard free energy of formation of in aqueous solution that is required by the hypothetical 1.0 molal standard state commonly used for aqueous species (see below).
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. , , ). 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
Behaviour of electrolyte activity coefficients
Further analysis for many electrolytes allowed the development of an equation of state characterization for as a function of temperature and pressure. This equation of state was revised and parameterized for , 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 content of the fluid.
Evaluation of the activity coefficients of neutral aqueous species
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 . The thermodynamic treatment of in water has typically taken one of two paths. For low temperature (i.e. surficial conditions) and geothermal systems, 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 -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 discussed above. Furthermore, the mole fraction concentration scale is used, as it is more practical when the range of fluid chemistry can approach pure . However, recent ab initio molecular dynamics calculations (Pan & Galli 2016) have shown that at 1000 K and 100 kbar solutions contain mainly ions, and that the predominant neutral species is not but instead is . The discovery of the importance of 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
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 is 10−3.5 bar.
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 is 10−3.5 bar.
Pure calcite reacting with soil water equilibrated with 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.
Pure calcite reacting with surface seawater equilibrated with the atmosphere at bar.
It is assumed that the concentration of Ca is now 411 ppm, and that the pH is 8.1.
Summary of how to predict a nonequilibrium chemical reaction
The numerical value of 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 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 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.
The sign of is what is important. The sign of is not important.
- 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(54)
But this is a special circumstance.
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.
- 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(57)
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: , , , , , , and . 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. ) we are left with seven unknown activities and concentrations.
In addition to the mass action equations, we have mass balance and charge balance equations to use as constraints, as follows.
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 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 to the monovalent silicate anion (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.
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 and 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 and 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.
Aqueous species at elevated temperatures and pressures
Here again, the fundamental information needed to evaluate equations (75) and (76) is the , and 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, 38b, 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 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
Born equation for solvation of an aqueous ion
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 (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).
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.
Relationships between the properties of ions and electrolytes
Regression equation for the volume of an aqueous electrolyte at T < 200°C and Psat
Regression equation for the compressibility of an aqueous electrolyte at T < 200°C and Psat
The values of and obtained from the temperature dependence of the compressibility can then be used with the values of and obtained previously from the temperature dependence of the volume, together with equations (97) and (98), to generate values of and . 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 and was revised by Tanger & Helgeson (1988) to incorporate an asymptotic dependence on pressure with a low-pressure limiting value given by the parameter bar. This value was obtained from the only available data at the time for the volumes of and 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
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
Equation for the apparent standard partial molal Gibbs free energy of formation of an aqueous ion as a function of temperature and pressure
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 , and 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 and as functions of temperature are available, or (2) experimental data for and at 25°C only are available, or (3) no experimental data are available.
Estimation of equation of state coefficients of aqueous species
Using estimated values of allows estimation of the solvation volume and heat capacity.
Values of can be obtained from predictive correlations with as can be seen in Figure 20.
The equations of the correlation lines in Figure 20a and b are as follows.
Estimation of the volume equation of state coefficients , , and 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 as described above, an experimental volume at 25°C and 1.0 bar is used to calculate the non-solvation volume . The value of is then used to estimate the parameter . 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 appears. At pressures of 20 000 bar and higher, this term dominates the change of with pressure. It is consequently very important to have as accurate an estimate of as possible. The correlation between and currently used in the DEW model differs from that originally proposed by Shock & Helgeson (1988), which was based only on data for , , , and . 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.
Estimation of the heat capacity equation of state coefficients and 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 and at 25°C and 1.0 bar was developed for aqueous ions, including , an abbreviation for (Shock & Helgeson 1988). The same correlation was found to apply to the inorganic neutral species NH03, , and (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 using the following equations.
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
Hydration of aqueous complexes at high temperatures and pressures in fluids
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 and in water has been studied over a wide range of conditions (Jahn & Schmidt 2010) and indicates a transition from and 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 to with increasing pressure (Fig. 25), and in fluids, ab initio molecular dynamics modelling (Pan & Galli 2016) revealed the formation of , , and . The latter study also revealed that as a molecule is unimportant compared with at high temperatures and pressures.
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 , and 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 , , , , , and . It therefore becomes possible to regress equilibrium constants as functions of temperature and pressure to directly obtain values of , and 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 , and 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 , , , , , 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.
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 , and 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.
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 , , , , and , as well as organic ligands such as , and . 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 is much more important than 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.
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