Thermodynamic modelling of fluids from surficial to mantle conditions

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 nearsurface 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 CO 2 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 CO 2 , CO, CH 4 , H 2 and H 2 O (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(Kessel et al. ,b, 2015Dvir 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. 2014Facq et al. , 2016 now clearly show that the fluids leaving subducting slabs must contain all the rockforming 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(Helgeson , 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, 1976Helgeson 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 1988Tanger & Helgeson 1988;Shock et al. 1989Shock et al. , 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 . 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 semiempirical 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(Facq et al. , 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 fO 2 , 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(Galvez et al. , 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   ) based on polynomial expressions of experimental data (Heger et al. 1980) and theoretical estimates (Pitzer 1983); (b) range of the dielectric constant in the DEW model (Pan et al. 2013;. MD, molecular dynamics.

349
Thermodynamic modelling of fluids in Earth 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 (DG r ) according to where DG 0 r is the standard Gibbs free energy of reaction and Q is the thermodynamic activity quotient, which are defined by and In equation (2), v i represents the stoichiometric reaction coefficient of the ith species in the reaction and D G 0 f ,i 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 D G 0 f ,i 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 D G 0 f ,i 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 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 DG r , 0, the reaction should proceed as written. However, if DG r . 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 CO 2(aq) , CO 2(g) , HCO À 3 or CO 2À 3 . The species to choose is the one that is most practical. Rainwater can be assumed to be in equilibrium with atmospheric CO 2 . Therefore, for convenience, we can write the reaction as for which we would need to evaluate the sign of DG r calculated from the equation We can evaluate the first term on the right-hand side of equation (6) if we have values of D G 0 f ,i 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 a calcite CaCO 3 related to how much Ca versus Mg is in the calcite? And how is the a Ca 2þ related to the actual concentration of Ca 2þ 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 DG r ¼ 0?

Prediction of equilibrium reactions
When DG r ¼ 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, Standard states depend only on pressure and temperature. In these circumstances, at a given pressure and temperature, DG 0 r 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 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 K 0 ), but by convention the superscript is omitted. Here, again, DG 0 r represents the standard Gibbs free energy of the reaction given by and 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.

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, a MZ , is given by where X MZ is a mole fraction of MZ in the solid solution and l MZ is a rational activity coefficient, a function of temperature, pressure and X MZ .

D. A. Sverjensky
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 l MZ ! 1:0 as X MZ ! 1: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, SiO 2 is in its standard state; that is, a quartz SiO 2 ¼ X quartz SiO 2 ¼ 1: But even when the mineral is not pure, as in Figure 3, we can often assume that with high enough X j → 1.0 the component MZ will behave as though it is experiencing ideal mixing when X MZ is in the region where l j ! 1:0. For example, the component CaCO 3 in impure calcite is not in its standard state. But if calcite is only slightly impure, l CaCO 3 % 1, so

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 where In other words, pure liquid water is in its standard state and we can write Furthermore, we can often assume that even impure water is in the region where X H 2 O ! 1:0 and l H 2 O ! 1:0, and therefore behaves ideally; that is, that 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 H 2 O 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 CO 2 and other neutral molecules that can be thought of as mixing with H 2 O at elevated temperatures and pressures. It should be noted that the standard Gibbs free energy of formation for CO 2 that should be used for this standard state must be that of pure fluid CO 2 , completely different from the standard free energy of formation of CO 2 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 where f j is the fugacity (with units of pressure) and f 0 j 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 y j is the fugacity coefficient for the real gas and y 0 j 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, y 0 j and p 0 j are equal to 1.0 for the standard state. In other words, A real gas cannot actually behave ideally at the temperature of interest and 1.0 bar. Therefore, the gas standard state is an example  The standard state is represented by the blue cross at X H2O ¼ 1:0; that is, the pure fluid is in its standard state. 351 Thermodynamic modelling of fluids in Earth 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 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 where g j and m j refer to the activity coefficient and molality of j, and g 0 j and m 0 j 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 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. CO 2 , CH 4 , N 2 ). 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(Galvez et al. , 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 ( g +,k ) by the equation which is based on equation (64) of Helgeson et al. (1981). The individual ionic activity coefficient of any ionic species ( j ) with a charge Z j is represented by log g j;P,T ¼ À 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 g and B g are functions of temperature, and the density and dielectric constant of pure water. The ionic strength of the solution ( I) is given by and a k 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 g,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 where m* represents the sum of the molalities of all aqueous species in the solution. The standard state is represented by the blue cross at m 0 j ¼ 1:0 on the tangent taken to the actual activity curve at very low molalities; that is, the standard state is not on the real activity curve. This standard state is hypothetical; it is not physically achievable by an aqueous species.

Behaviour of electrolyte activity coefficients
To evaluate the extended term contribution in equation (23), a deviation function for the kth electrolyte (r Ã g,k;P,T ) can be defined according to r Ã g,k;P,T ¼ log g k;P,T À A g jZ i Z l j I 0:5 1 þ a k B g I 0:5 À G g that is, where r Ã g,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 g,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 NaOH 0 , 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 g,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 CO 2 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 g n;P, where the parameter b g,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 g,n was made by Helgeson et al. (1981), estimates of the temperature and pressure dependence of b g,n were made subsequently for a number of 1:1 ion pairs (Oelkers & Helgeson 1991).
CO 2 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 CO 2 .
The thermodynamic treatment of CO 2 in water has typically taken one of two paths. For low temperature (i.e. surficial conditions) and geothermal systems, CO 2 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 CO 2 -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 H 2 O discussed above. Furthermore, the mole fraction concentration scale is used, as it is more practical when the range of fluid chemistry can approach pure CO 2 . However, recent ab initio molecular dynamics calculations (Pan & Galli 2016) have shown that at 1000 K and 100 kbar Na 2 CO 3 solutions contain mainly ions, and that the predominant neutral species is not CO 2 but instead is H 2 CO 3 . The discovery of the importance of H 2 CO 3 will require a complete re-evaluation of the speciation and activityconcentration 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: To predict whether or not calcite will dissolve in a water of given composition it is necessary to evaluate the sign of DG r using The sign of DG r 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 DG r , 0, reaction proceeds as written (31) DG r . 0, reaction proceeds in reverse (32) DG r ¼ 0, reaction is at equilibrium: 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 where v i D G 0 f ,i 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 Fig. 6. Examples of the linear behaviour of electrolytes over wide ranges of ionic strength, which can be described with a single parameter (from Helgeson et al. 1981, reprinted by permission of the American Journal of Science). The symbols represent values of b g,k for each electrolyte obtained from experimental measurements of single electrolyte activity coefficients log g k at 25°C and 1.0 bar used in equation (26). The line represents a one-parameter fit to the symbols according to equation (27).

353
Thermodynamic modelling of fluids in Earth 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 D G 0 f ,i at 25°C and 1 bar (Berman 1988;Shock et al. 1997 :0 cal mol À1 : Substituting these standard free energies into equation (34) permits evaluation of the standard Gibbs free energy of the reaction according to It is important to remember that the value of DG 0 r 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.
(1) 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 p CO 2 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 CaCO 3 is the only component in the calcite; that is, Together with recognition of the dilute nature of rainwater, which allows us to assume that the thermodynamic activity of Ca 2þ will be equal to its molality, equation (30) becomes where the a pure calcite ¼ log 2 1000 Â 40:08 À 3:5 þ 2(5:7) ¼ log(5:0 Â 10 À5 ) À 3:5 þ 2(5:7) ¼ À4:3 À 3:5 þ 11:4 ¼ 3:6: Overall, combining equations (35), (37), and (38) results in Based on equation (31), the negative sign for the calculated value of DG r 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.
It is again assumed that the concentration of Ca in the rainwater is 2.0 ppm, that the pH is 5.7, and that p CO 2 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 CaCO 3 is not the only component in the calcite; that is, Proceeding as above we now obtain where the a calcite CaCO 3 ¼ 0:8 by assuming that the CaCO 3 component mixes ideally with the MgCO 3 component. Now evaluating the activity quotient term in equation (41) results in ¼ log 2 1000 Â 40:08 À 3:5 þ 2(5:7) þ 0:097 ¼ log (5:0 Â 10 À5 ) À 3:5 þ 2(5:7) þ 0:097 ¼ À4:3 À 3:5 þ 11:4 þ 0:097 ¼ 3:7: Overall, combining equations (35), (41) and (42) 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 DG r 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.
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: Proceeding as in Case 1, we obtain where now the soil water is so different in chemistry from the original rainwater that ¼ log 2 1000 Â 40:08 À 1:5 þ 2(7:7) ¼ log (5:0 Â 10 À5 ) À 1:5 þ 2(7:7) ¼ À4:3 À 1:5 þ 15:4 ¼ 9:6: 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 DG r ¼ À13 180 þ 2:303RT (9:6) ¼ À13 180 þ 2:303(1:987 Â 298:15)(9:6) ¼ À13 180 þ 13 098 The very small negative number for the calculated value of DG r 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.

Summary of how to predict a nonequilibrium chemical reaction
For a single chemical reaction, such as the dissolution of calcite according to the fundamental equation to be evaluated is which requires evaluation of the standard Gibbs free energy of reaction (DG 0 r ) 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 DG 0 r 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 DG 0 r 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 DG 0 r 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 DG r is what is important. The sign of DG 0 r 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 Another example would be any reaction consisting of pure minerals and water that is (traditionally) assumed to be pure or very close to 355 Thermodynamic modelling of fluids in Earth pure; for example, the breakdown of talc according to 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 In these circumstances, the sign of DG 0 r becomes meaningful because it is the same as for DG r .
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 we assumed that the only aqueous Ca-bearing species in the water was the Ca 2þ ion. Therefore, the activity of the ion required for equation (57) could be set equal to the measured analytical value of Ca (m t,Ca ); that is, 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) Ca 2þ 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 þ , Mg 2þ and Ca 2þ . 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.

Speciation of simple NaCl-HCl-H 2 O solutions
We start by assuming that the solution could contain the following eight aqueous species: Na þ , Cl À , H þ , OH À , NaCl 0 , HCl 0 , NaOH 0 and H 2 O. 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. a H 2 O ¼ 1) we are left with seven unknown activities and concentrations.
However, there are relationships between these species because of equilibrium in solution: These equilibria correspond to the following mass action equations: where K NaCl 0 , K NaOH 0 , K HCl 0 and K H 2 O 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. for Na m t,Na þ ¼ m Na þ þ m NaCl 0 þ m NaOH 0 for Cl m t,Cl À ¼ m Cl À þ m NaCl 0 þ m HCl 0 : Charge balance equation: 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 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(HCO 3 ) þ 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) 0 4 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 Clearly, we need to be able to calculate standard Gibbs free energies of reactions (D G 0 r;P,T ) 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 (D G 0 f ,j;P,T ); that is, As an example, the equilibrium between calcite and water may be written in terms of the dissolved CO 0 2 species according to for which we can write 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 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 (D G 0 j;P,T ). Values of D G 0 j;P,T 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)

357
Thermodynamic modelling of fluids in Earth the issue of calculating the term G 0 j;P,T À G 0 j;P r ,T r , which can be carried out using the following equation: Here the standard partial molal entropy of the jth species is needed at the reference temperature and pressure ( S 0 j;P r ,T r ), as well as the temperature dependence of the isobaric standard partial molal heat capacity ( C 0 j;P r ), and the pressure dependence of the standard partial molal volume ( V 0 j;T ). 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 0 j;P r (T ) and V 0 j;T (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 0 j;P r (T ) and V 0 j;T (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 is used, where k 0 , k 1 , k 2 and k 3 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 is used, where v 1 , v 2 , v 3 and v 4 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

Aqueous species at elevated temperatures and pressures
Here again, the fundamental information needed to evaluate equations (75) and (76) is the S 0 P r ,T r , C 0 P r (T ) and V 0 T (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, 1976Helgeson 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 (1 H 2 O ¼ 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 semiempirical 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

359
Thermodynamic modelling of fluids in Earth 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. Fig. 11. Equilibrium dissociation constants decrease strongly with increasing temperature, and increase with increasing pressure: (a) experimental data points for NaCl based on an analysis of conductance data (Oelkers & Helgeson 1991), with theoretical curves extrapolated with the DEW model based on a fit to the data (Shock et al. 1992); (b) experimental data points for HCl based on solubility data (Ruaya & Seward, 1987;Sverjensky et al., 1991;Tagirov et al., 1997) fit with theoretical curves extrapolated with the DEW model.

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 The intrinsic and the collapse terms can be summed to give the non-solvation volume (D V 0 n,j ) according to Fig. 13. The molar volume and heat capacity of water as a function of temperature and pressure calculated by Helgeson & Kirkham (1974a, reprinted by permission of the American Journal of Science). The curves in (a) are based on fits to data from Burnham et al. (1969) consistent with the curves in Figure 12a and b. The curves are contoured in kbar.
Fig. 14. The standard partial molal volume and heat capacity of the NaCl electrolyte (and other similar salts) as a function of temperature at pressures just above the liquid-vapor saturation curve of water. The curves labelled 'present study' refer to Tanger & Helgeson (1988, reprinted by permission of the American Journal of Science), who fitted equations of state to the experimental data represented by the symbols. 361 Thermodynamic modelling of fluids in Earth A similar set of three contributions was subsequently proposed (Helgeson et al. 1981) for the standard partial molal heat capacity according to C 0 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 (D G 0 s,j ) given by where v j represents the conventional Born solvation coefficient of the aqueous ion and 1 H 2 O 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 (v abs: where v abs: H þ ¼ 0:5387 Â 10 5 cal mol −1 refers to the absolute Born coefficient of the H + ion (v abs: H þ ) and Z : j 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 v abs: where h is a constant equal to 1:66027 Â 10 5 Å cal mol −1 and r e,j represents the effective electrostatic radius of the aqueous ion (Å). According to Helgeson & Kirkham (1976), r e,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 r e,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 where r x,j represents the crystallographic radius of the jth ion (Shock & Helgeson 1988), Z j again represents the charge on the jth ion, k Z 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 NaCl 0 (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 NaCl 0 are consistent with the assumption that g = 0.0 (Shock et al. 1992). Under these conditions, equation (88) simplifies to r e,j ¼ r x,j þ jZ j jk Z and consequently v 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 P sat 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 v j is independent of temperature and pressure to obtain the standard partial molal volumes, compressibilities and heat capacities of where Q P,T , N P,T and X P,T are partial derivatives of 1 H 2 O given by Equations of state for the standard partial molal properties of an aqueous species at T < 200°C and P sat 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 P sat 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 where s and j are temperature-independent coefficients characteristic of the jth ion and the pressure dependence is given by and 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 For the non-solvation heat capacity 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 v P r ,T r;j is a constant can be written V and C 0 P;j ¼ c 1,j þ c 2,j (T À u) 2 þ v P r ,T r ;j TX P,T :

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 363 Thermodynamic modelling of fluids in Earth the standard partial molal volumes of the ions; that is, For the kth 2:1 electrolyte, 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 and 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 0 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 Z j are expressed by the equation, in the case of volumes, as V where V 0 abs: H þ represents the absolute standard partial molal volume of the H þ ion. Expressions analogous to those in equations (109) and (110) (96), we can then expect that D V 0 n;k (expt:) will show a linear dependence on (T À u) À1 according to the equation Consistency with equation (112) allows regression to obtain values of s k and j k for an electrolyte. Examples are shown in Figure 16 for NaCl and MgCl 2 (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 s k and j k to be obtained.

Regression equation for the compressibility of an aqueous electrolyte at T < 200°C and P sat
By analogy with the procedure for volumes described above, experimental measurements of k 0 as a function of temperature can be used to obtain D k 0 n (expt:) as a function of temperature using the equation We can then expect that D k 0 n (expt:) will show a linear dependence Fig. 17. The standard partial molal nonsolvation compressibilities of the NaCl and MgCl 2 electrolytes as functions of inverse temperature. The lines represent regression of the data to obtain the equation of state coefficients a 2,k and a 4,k (Tanger & Helgeson 1988, reprinted by permission of the American Journal of Science).
Consistency with equation (114) allows a linear regression to obtain values of a 2,k and a 4,k for an electrolyte. Examples are shown in Figure 17 for NaCl and MgCl 2 (Tanger & Helgeson 1988). The values of a 2,k and a 4,k obtained from the temperature dependence of the compressibility can then be used with the values of s k and j k obtained previously from the temperature dependence of the volume, together with equations (97) and (98), to generate values of a 1,k and a 3,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 P sat
The pressure dependence of the equation of state coefficients s k and j k was revised by Tanger & Helgeson (1988) to incorporate an asymptotic dependence on pressure with a low-pressure limiting value given by the parameter c ¼ 2600 bar. This value was obtained from the only available data at the time for the volumes of NaCl and K 2 SO 4 at 25°C and pressures up to 10 kbar (Adams 1931(Adams , 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(Mantegazzi et al. , 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 P sat Again proceeding as above for the analysis of the volume and compressibility data at low temperatures, experimental measurements of C 0 P;k as a function of temperature can be used to obtain D C 0 P,n;k (expt:) as a function of temperature using the equation We can then expect that D C 0 P,n;k (expt:) will show a linear Fig. 19. The standard partial molal nonsolvation heat capacities of the NaCl and MgCl 2 electrolytes as functions of the square of the inverse temperature. The lines represent regression of the data to obtain the equation of state coefficients c 1,k and c 2,k (Tanger & Helgeson 1988, reprinted by permission of the American Journal of Science).

365
Thermodynamic modelling of fluids in Earth dependence on (T À u) À2 according to the equation Consistency with equation (116) allows regression to obtain values of c 1,k and c 2,k for the electrolytes NaCl and MgCl 2 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 1988Shock et al. 1989Shock et al. , 1997Amend & Helgeson 1997;Plyasunov & Shock 2001;Dick et al. 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 P sat 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 0 j and C 0 P;j of the jth aqueous species are given by V where the derivatives of v 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: where the coefficients c 1,j , c 2,j , a 1,j , a 2,j , a 3,j , a 4,j and v 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 lowtemperature and -pressure volumes, compressibilities and heat capacities (i.e. at T < 200°C and P sat ) 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 0 T , k 0 T and C 0 P r ,T 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 0 T and C 0 P r ,T as functions of temperature are available, or (2) experimental data for V 0 P r and C 0 P r 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 0 P r and C 0 P r 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 Â 10 À5 bar À1 and X ¼ 0:0309 Â 10 À5 K À2 ; 366 D. A. Sverjensky that is D V 0 s ¼ Àv j;P r ,T r (0:05903 Â 10 À5 Â 41:84) cm 3 mol À1 (120) and D C 0 P;s ¼ v j;P r ,T r (298:15)(0:0309 Â 10 À5 ) cal mol À1 K À1 (121) where v j;P r ,T r is in cal mol À1 .
Using estimated values of v j;P r ,T r allows estimation of the solvation volume and heat capacity.
Values of v j;P r ,T r can be obtained from predictive correlations with S 0 j;P r ,T r 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 and for anions b Z ¼ 1:62 Â 10 5 Z j cal mol À1 : In contrast, for neutral species that are typically crystalline at ambient conditions when not dissolved in water v j;P r ,T r ¼ þ0:30 Â 10 5 cal mol À1 : For neutral species that are typically gases or liquids at ambient conditions when not dissolved in water v j;P r ,T r ¼ À0:30 Â 10 5 cal mol À1 : 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  (Shock & Helgeson 1988;Shock et al. 1989); (d) estimation of a 2 from a 4 for ions (Shock & Helgeson 1988).

367
Thermodynamic modelling of fluids in Earth 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 s, a 1 , a 2 and a 4 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 1988Shock et al. 1989). After estimating the solvation volume at 25°C and 1.0 bar (D V 0 s ) as described above, an experimental volume at 25°C and 1.0 bar is used to calculate the non-solvation volume (D V 0 n ). The value of D V 0 n is then used to estimate the parameter a 1 . 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 a 1 (P À P r ) appears. At pressures of 20 000 bar and higher, this term dominates the change of D G 0 f ;P,T with pressure. It is consequently very important to have as accurate an estimate of a 1 as possible. The correlation between a 1 and D V 0 n currently used in the DEW model differs from that originally proposed by Shock & Helgeson (1988), which was based only on data for Mg 2þ , 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 For aqueous species of all kinds For aqueous species of all kinds where s and D V 0 n are in cm 3 mol −1 , a 1 is in cal mol −1 bar −1 , a 2 is in cal mol −1 and a 4 is in cal K mol −1 .
The procedure is as follows. After estimation of a 1 from Figure 21a or b, s can be estimated using Figure 21c. The coefficient a 3 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 a 1 , a 2 , a 4 and v.
Estimation of the heat capacity equation of state coefficients c 1 and c 2 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 c 2 and C 0 P r at 25°C and 1.0 bar was developed for aqueous ions, including AlO À 2 , an abbreviation for Al(OH) À 4 (Shock & Helgeson 1988). The same correlation was found to apply to the inorganic neutral species NH 0 3 , H 2 S 0 , B(OH) 0 3 and CH 3 OOH 0 (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 c 2 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 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 highpressure 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.

Activity coefficient approximations
The full equation for the activity coefficient of the jth ion is log g j;P,T ¼ À Values of A g and B g 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 g,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 g,k ¼ 0:0, resulting in the equation log g j;P, T ¼ À A g Z 2 j I 0:5 1 þ a k B g I 0:5 þ G g : 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 CO 2 , CO and CH 4 in fluids at high temperatures and pressures, an activity coefficient is needed to account for the strong non-ideality apparent in the traditional CO 2 -H 2 O models (Duan & Zhang 2006). Here we present a first approximation obtained by transforming the Zhang & Duan CO 2 -  Figure 23 and accompanying equations allow estimation of all the equation of state coefficients c 1,j , c 2,j , a 1,j , a 2,j , a 3,j , a 4,j and v j .
where b P,T is a constant fitted to the transformed Zhang & Duan CO 2 -H 2 O model and G g 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. CO 2 ). The non-ideality of water in such fluids can be accounted for by the expression For simplicity here, we use log(a H 2 O ) ¼ logX H 2 O ¼ log 55:51 55: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) 0 2 we use MgO 0 , which differs only by 1.0 H 2 O and allows us to write the reaction When the species MgO 0 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) 0 2 þ 2H þ ¼ Mg 2þ þ 2H 2 O: 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 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 Mg 2þ and SO 2À 4 in water has been studied over a wide range of conditions (Jahn & Schmidt 2010) and indicates a transition from Mg 2þ and HSO À 4 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 HCO À 3 to CO 2À 3 with increasing pressure (Fig. 25), and in Na 2 CO 3 fluids, ab initio molecular dynamics modelling (Pan & Galli 2016) revealed the formation of NaHCO 0 3 , Na 2 HCO þ 3 , NaCO À 3 and Na 2 CO 0 3 . The latter study also revealed that CO 0 2 as a molecule is unimportant compared with H 2 CO 0 3 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) 0 4 ¼ H þ þ SiO(OH) À 3 : A reaction that has the same equilibrium constant is shown in Figure 26a. It can be seen that the logK increases dramatically with  Similarly, water itself dissociates strongly with increasing pressure (Fig. 26b) according to the reaction indicating the availability of the OH − ion as a potential complexforming 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 , k 0 and C 0 P 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 c 1 , c 2 , a 1 , a 2 , a 3 , a 4 and v. It therefore becomes possible to regress equilibrium constants as functions of temperature and pressure to directly obtain values of V 0 , C 0 P 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 0 P 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 c 1 , c 2 , a 1 , a 2 , a 3 , a 4 and v.
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 Al 3þ , Al(OH) 0 3 and Al(OH) À 4 , the last two of which have been calibrated previously (Sverjensky et al. 2014a). Al(OH) 0 3 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 P sat up to 300°C (Tagirov & Schott 2001) was found to account for the data well ( Fig. 27a) according to the reaction Al(OH) 3 OSi(OH) À 3 þ 4H þ ¼ Al 3þ þ Si(OH) 0 4 þ 3H 2 O (144) which is equivalent to the reaction shown in the figure as 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 0 P 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  (Kessel et al. 2005b) and model fit to retrieve log K of the Al-Si complex; (b) equilibrium constants retrieved from the solubility of eclogite at 4.0 and 5.0 GPa (Huang 2017) and published values at P sat (Tagirov & Schott 2001).