## Abstract

Carbon is subducted to depths where metamorphism liberates water-bearing fluids. The C-bearing fluids facilitate partial melting of the upper mantle, generating magmas that may erupt as arc volcanics. Degassing of the magmas releases CO_{2} 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 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 2001*a*, 50*b*; Hacker *et al.* 2003*a*,32*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.* 2005*a*,52*b*, 2015; Dvir *et al.* 2011; Dvir & Kessel 2017; Tsay *et al.* 2017; Tumiati *et al.* 2017; Tiraboschi *et al.* 2018) and of aqueous speciation (Martinez *et al.* 2004; Mysen 2010; Louvel *et al.* 2013; Mysen *et al.* 2013; Facq *et al.* 2014, 2016) now clearly show that the fluids leaving subducting slabs must contain all the rock-forming elements and trace elements in addition to the dissolved gases mentioned above. All the above evidence emphasizes the need for comprehensive models of aqueous speciation, solubility and chemical mass transfer under subduction-zone conditions.

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

The dielectric constant of water is a vital parameter in the HKF approach because it appears in the Born equation for solvation that is part of the equation of state for computing the standard partial molal properties and the non-standard state contributions for aqueous species as functions of temperature and pressure (Helgeson *et al.* 1981). To remove the limitation of 5.0 kbar in the HKF approach, a characterization of the dielectric constant of water at greater pressures was developed by building on a semi-empirical equation for polar solvents (Franck *et al.* 1990) and *ab initio* molecular dynamics results to about 100 kbar (Pan *et al.* 2013). The new model for the dielectric constant of water (Fig. 2b) led to the development of the Deep Earth Water (DEW) model (Sverjensky *et al.* 2014*a*). 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.* 2014*a*) 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.* 2014*b*); 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 *f*O_{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, 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

*Q*is the thermodynamic activity quotient, which are defined by

*i*th species in the reaction and $\Delta G\xaff,i0$ represents the standard Gibbs free energy of formation of the

*i*th 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 $\Delta G\xaff,i0$ for a given species in a chosen standard state depend only on pressure and temperature. The values are derived from experimental studies and models interpreting the experimental studies. The general topic of how to calculate values of $\Delta G\xaff,i0$ as functions of pressure and temperature is discussed later in this paper.

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

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

## Prediction of equilibrium reactions

*Q*is fixed. This special value of

*Q*is called the thermodynamic equilibrium constant and is denoted by

*K*according to

*K*represents the special value of

*Q*that depends only on the pressure and temperature under consideration. This result follows from the definitions of the standard states as depending only on the pressure and temperature of the reaction as will be defined below. Consequently,

*K*is also a standard state quantity. To recognize this,

*K*should have a superscript (e.g. it should be written $K0$), but by convention the superscript is omitted. Here, again, $\Delta Gr0$ represents the standard Gibbs free energy of the reaction given by

*K*is the thermodynamic equilibrium constant, which fixes a unique activity ratio for the species in the reaction at a given temperature and a specific pressure. It should be noted that the individual activities are not fixed; only the activity quotient is fixed.

## Definitions of standard states

### Minerals

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 $\lambda MZ\u21921.0$ as $XMZ\u21921.0$, as depicted in Figure 3.

*X*→ 1.0 the component MZ will behave as though it is experiencing ideal mixing when X

_{j}_{MZ}is in the region where $\lambda j\u21921.0$. For example, the component $CaCO3$ in impure calcite is not in its standard state.

### Water and mixed solvents

_{2}O and $\lambda H2O$ is a rational activity coefficient (a function of temperature, pressure, and $XH2O$).

Here, the standard state refers to unit activity of the pure component H_{2}O at the pressure and temperature of interest according to Raoult's Law; that is, $\lambda H2O\u21921.0$ as $XH2O\u21921.0$.

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

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

### Gases

*j*in the standard state. Each fugacity is related to a partial pressure of a gas by a fugacity coefficient, where $\upsilon j$ is the fugacity coefficient for the real gas and $\upsilon j0$ is the fugacity coefficient for the gas when it is in its standard state.

### Aqueous species

*j*, and $\gamma j0$ and $mj0$ refer to the standard state values of these quantities (Anderson 2005). The most frequently used standard state refers to unit activity of a hypothetical 1.0 molal solution of

*j*referenced to infinite dilution at the temperature and pressure of interest. Because of this definition, real aqueous species can never be in their standard state! In other words, the standard state for aqueous species is not physically achievable. The hypothetical 1.0 molal aqueous species standard state is depicted in Figure 5.

## Activity coefficients linking standard states to measurable concentrations in fluids

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

### Evaluation of the activity coefficients of aqueous ionic species

*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 $(\gamma \xaf\xb1,k)$ by the equation

*et al.*(1981).

*j*) with a charge $Zj$ is represented by

*et al.*(1981).

*k*th salt in solution (see equations (124) and (125) of Helgeson

*et al.*1981).

*k*th 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

*m** represents the sum of the molalities of all aqueous species in the solution.

### Behaviour of electrolyte activity coefficients

*k*th electrolyte $(\rho \gamma ,k;P,T\u2217)$ can be defined according to

*et al.*(1981) has demonstrated that in solutions of a single electrolyte, the simplicity of equation (27) is valid to very high ionic strengths, something that seems to have been forgotten in the geochemical literature. Systematic departures from linearity as shown in Figure 6 for the $NaOH$ electrolyte at the highest ionic strengths indicate the onset of aqueous complexing, in this instance the likely formation of the species $NaOH0$, which is part of the speciation of the overall $NaOH$ electrolyte.

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

### Evaluation of the activity coefficients of neutral aqueous species

*n*th neutral aqueous species in a solution of the

*k*th electrolyte can be expressed by

*et al.*(1981), estimates of the temperature and pressure dependence of $b\gamma ,n$ were made subsequently for a number of 1:1 ion pairs (Oelkers & Helgeson 1991).

### $CO2$ and other dissolved gases: solute or solvent?

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

## Nonequilibrium reactions at ambient conditions: prediction of calcite dissolution

*i*th species in the reaction multiplied by the standard Gibbs free energy of formation of the species from its elements in their defined reference states. It should be noted that equation (34) only involves standard state free energies, which depend only on temperature and pressure. Each standard state refers to a fixed chemical composition, regardless of the actual state of interest in the natural water referred to in equation (29).

*et al.*1997) we have

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

### Pure calcite reacting with rainwater equilibrated with atmospheric carbon dioxide.

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

### Calcite solid-solution (20% MgCO_{3}) reacting with rainwater equilibrated with atmospheric carbon dioxide.

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

_{3}is not the only component in the calcite; that is,

### Pure calcite reacting with soil water equilibrated with $pCO2=10\u22121.5$ bar.

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

### Pure calcite reacting with surface seawater equilibrated with the atmosphere at $pCO2=10\u22123.5$ 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

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

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

The sign of $\Delta Gr$ is what is important. The sign of $\Delta Gr0$ 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$CaCO3purecalcite\u2192CaCO3purearagonite.$(54)

*Q*= 1.0 and hence

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$\Delta Gr=\Delta Gr0+RTlnaCa2+aCO2,gaH2OaCaCO3calciteaH+2=\Delta Gr0+2.303RTlogmCa2+pCO2(g)0.8aH+2$(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–H_{2}O solutions

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

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

*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

## Modelling at upper mantle conditions

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

### What is needed?

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

In contrast, the DEW model can be used to calculate equilibrium constants involving aqueous species at temperatures greater than 1000°C and pressures up to 60 kbar. An example is shown in Figure 7b for the dissociation of aqueous silica $Si(OH)40$ to the monovalent silicate anion $SiO(OH)3\u2212$ (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.* 2005*b*). The basis for calculations such as those shown in Figure 7a and b is presented below.

*j*th species (Helgeson

*et al.*1978) where

*j*th 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

*j*th species alone with temperature and pressure. The corresponding changes for the elements that make up the

*j*th species are not evaluated. By using equation (75) only the apparent standard Gibbs free energy of formation of the

*j*th species is calculated $(\Delta G\xafj;P,T0)$. Values of $\Delta G\xafj;P,T0$ from equation (75) are intended to be used in equations such as equation (74) because such equations refer to balanced chemical reactions (equation (73)) for which the changes in temperature and pressure for all the elements if used would simply cancel. Consequently, evaluation of equations (74) and (75) reduce to the issue of calculating the term $G\xafj;P,T0\u2212G\xafj;Pr,Tr0$, which can be carried out using the following equation:

*j*th species is needed at the reference temperature and pressure $(S\xafj;Pr,Tr0)$, as well as the temperature dependence of the isobaric standard partial molal heat capacity $(C\xafj;Pr0),$ and the pressure dependence of the standard partial molal volume $(V\xafj;T0)$. Equations for evaluating these functions are described below for minerals and for aqueous species.

### Mineral free energies at elevated pressures and temperatures

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

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

### Aqueous species at elevated temperatures and pressures

Here again, the fundamental information needed to evaluate equations (75) and (76) is the $S\xafPr,Tr0$, $C\xafPr0(T)$ and $V\xafT0(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 1974*a*, 38*b*, 1976; Helgeson *et al*. 1981) and forms the basis for all subsequent revisions and extensions. The primary reason is that the HKF approach was built to have predictive capabilities, particularly for the standard state properties of aqueous species. Some background on the development of the HKF approach is summarized below.

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

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

#### Dielectric constant of water at elevated temperatures and pressures

Experimental measurements of the dielectric constant of water are shown in Figure 10a and b (Heger *et al.* 1980), together with estimations based on a literature review (Fernández *et al*. 1997), and *ab initio* molecular dynamics calculations (Pan *et al.* 2013). The curves in Figure 10a represent a fit of the data with an empirical polynomial, whereas the curves in Figure 10b represent a semi-empirical theory calibrated with only two parameters (Franck *et al.* 1990; Sverjensky *et al.* 2014*a*). 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.* 2014*a*).

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

*et al.*1981) for the standard partial molal heat capacity according to

#### Born equation for solvation of an aqueous ion

*j*th ion can be expressed as the standard partial molal Gibbs free energy $(\Delta G\xafs,j0)$ given by

*j*th ion is related to the absolute Born solvation coefficient $(\omega jabs.)$ by

^{−1}refers to the absolute Born coefficient of the H

^{+}ion ($\omega H+abs.$) and $Zj.$ refers to the charge on the ion.

*j*th aqueous ion can in turn be related to a radius for the

*j*th ion according to

^{−1}and $re,j$ represents the effective electrostatic radius of the aqueous ion (Å). According to Helgeson & Kirkham (1976), $re,j$ is an electrostatic parameter that represents the distance from the centre of the ion to the location where the water molecules around the ion have the properties of bulk water.

*j*th aqueous ion has been related to the crystallographic radius by the equation

*j*th ion (Shock & Helgeson 1988),

*Z*again represents the charge on the

_{j}*j*th ion, $kZ$ is equal to 0.0 for anions and 0.94 for cations based on the analysis of the temperature dependence of the standard state partial molal volumes of aqueous electrolytes (Helgeson & Kirkham 1976) and

*g*(

*P*,

*T*) designates a solvent function of temperature and density (Shock

*et al.*1992).

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

*T*< 200°C or at

*P*> 4.0 kbar, the properties of electrolytes and the dissociation constant of $NaCl0$ are consistent with the assumption that

*g*= 0.0 (Shock

*et al.*1992). Under these conditions, equation (88) simplifies to

*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 $\omega j$ is independent of temperature and pressure to obtain the standard partial molal volumes, compressibilities and heat capacities of solvation from the equations

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

*j*th ion and the pressure dependence is given by

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

#### Relationships between the properties of ions and electrolytes

*k*th 1:1 electrolyte consisting of the

*i*th cation and the

*l*th anion, the standard partial molal volume of the electrolyte is, by definition, equal to the sum of the standard partial molal volumes of the ions; that is,

*k*th 2:1 electrolyte,

*et al.*(1981).

*j*th ion with a charge of $Zj$ are expressed by the equation, in the case of volumes, as

#### Regression equation for the volume of an aqueous electrolyte at *T* < 200°C and P_{sat}

*k*th electrolyte $(V\xafk0)$ as a function of temperature can be used to obtain $\Delta V\xafn0(expt.)$ as a function of temperature using the equation

#### Regression equation for the compressibility of an aqueous electrolyte at *T* < 200°C and P_{sat}

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

#### Calibration of the volume behaviour of an aqueous electrolyte at pressures greater than P_{sat}

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

#### Regression equation for the standard partial molal heat capacity of an aqueous electrolyte using data at *T* < 200°C and P_{sat}

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

*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\xafj0$ and $C\xafP;j0$ of the

*j*th aqueous species are given by

*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

*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\xafT0$, $\kappa \xafT0$ and $C\xafPr,T0$ are very few. Therefore, to be able to carry out geochemical modelling, we need ways to estimate the equation of state coefficients for aqueous species. A great strength of the HKF approach is that it was built with the need for estimation algorithms in mind.

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

### Estimation of equation of state coefficients of aqueous species

*Q*and

*X*refer to 25°C and 1 bar: $Q=0.05903\xd710\u22125bar\u22121$ and $X=0.0309\xd710\u22125K\u22122$; that is

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

Values of $\omega j;Pr,Tr$ can be obtained from predictive correlations with $S\xafj;Pr,Tr0$ 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 $\sigma $, $a1$, $a2$ and $a4$ for aqueous species

After analysing numerous sets of experimental data for a variety of aqueous electrolytes, neutral species and aqueous organic species, a comprehensive approach was taken to develop correlations for estimating volume equation of state coefficients (Shock & Helgeson 1988, 1990; Shock *et al.* 1989). After estimating the solvation volume at 25°C and 1.0 bar $(\Delta V\xafs0)$ as described above, an experimental volume at 25°C and 1.0 bar is used to calculate the non-solvation volume $(\Delta V\xafn0)$. The value of $\Delta V\xafn0$ is then used to estimate the parameter $a1$. This step is a key part of the approach for predicting aqueous species free energies at upper mantle pressures. It can seen in the full free energy equation (equation (119)) that the term $a1(P\u2212Pr)$ appears. At pressures of 20 000 bar and higher, this term dominates the change of $\Delta G\xaff;P,T0$ with pressure. It is consequently very important to have as accurate an estimate of $a1$ as possible. The correlation between $a1$ and $\Delta V\xafn0$ currently used in the DEW model differs from that originally proposed by Shock & Helgeson (1988), which was based only on data for $Mg2+$, $Na+$, $K+$, $Cl\u2212$ and $Br\u2212$. 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.

^{3}mol

^{−1}, $a1$ is in cal mol

^{−1}bar

^{−1}, $a2$ is in cal mol

^{−1}and $a4$ is in cal K mol

^{−1}.

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

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

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

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

^{−1}K

^{−1}and $C\xafPr0$ is in cal mol

^{−1}K

^{−1}.

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

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

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

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

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

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

## II. Analysis of experimental solubility data

### Activity coefficient approximations

*j*th ion is

*a*). However, the application of an extended term equation of state to calculating $b\gamma ,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\gamma ,k=0.0$, resulting in the equation

_{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}–H

_{2}O model to the hypothetical 1.0 molal standard state, resulting in a semi-empirical equation of form

_{2}–H

_{2}O model and $\Gamma \gamma $ again refers to the sum of the molalities of all aqueous species in the fluid (i.e. ions, neutral species and complexes).

_{2}). The non-ideality of water in such fluids can be accounted for by the expression

*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

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

#### Metal complexes at upper mantle pressures

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

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

*logK*increases dramatically with pressure, suggesting that the silicate anion ought to be a significant ligand with metals in upper mantle fluids.

^{−}ion as a potential complex-forming ligand.

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

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

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

An incredibly useful feature of the HKF approach used in the DEW model is the predictive algorithm discussed above relating the standard partial molal properties of an aqueous species to the equation of state coefficients. For most aqueous metal complexes, we do not, and probably will not, have measurements of $V\xaf0$, $\kappa \xaf0$ and $C\xafP0$ as functions of temperature to calibrate the equation of state coefficients. However, the correlations in Figures 21 and 22 and accompanying equations link the partial molal properties to the equation of state coefficients $c1$, $c2$, $a1$, $a2$, $a3$, $a4$ and $\omega $. It therefore becomes possible to regress equilibrium constants as functions of temperature and pressure to directly obtain values of $V\xaf0$, $C\xafP0$ and $S\xaf0$ 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\xaf0$, $C\xafP0$ and $S\xaf0$ can be retrieved, plus all the equation of state coefficients linked by the correlation and equations above. An alternative approach would be to regress the equilibrium constants directly for all seven equation of state coefficients $c1$, $c2$, $a1$, $a2$, $a3$, $a4$ and *ω*. Typically, this will result in over-fitting the experimental data. Furthermore, the equation of state coefficients will not be consistent with the many dozens of sets that are based on experimental data. Consequently, this approach is not recommended.

*et al.*2005

*b*) were fitted with a speciation model that simultaneously fitted data for Na, Mg, Ca, Fe and Si (Huang 2017). Included in the model speciation for Al were $Al3+$, $Al(OH)30$ and $Al(OH)4\u2212$, the last two of which have been calibrated previously (Sverjensky

*et al.*2014

*a*). $Al(OH)30$ was calibrated using solubility data for corundum at low temperatures and pressures (Pokrovskii & Helgeson 1995), as well as at elevated temperatures and pressures (Becker

*et al.*1983; Tropper & Manning 2007). However, none of these species proved to be significant contributors to the solubility of Al in the more complex fluid in equilibrium with the K-free eclogite. Instead, after some trial regressions, an Al–Si complex previously studied at

*P*

_{sat}up to 300°C (Tagirov & Schott 2001) was found to account for the data well (Fig. 27a) according to the reaction

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\xaf0$, $C\xafP0$ and $S\xaf0$ referring to 25°C and 1.0 bar, which is only a three-parameter fit over a temperature range of about 1000°C and a pressure range of 50 kbar. Most importantly, the same equilibrium constants provide excellent predictions of the solubility of Al in two different peridotitic fluids (Huang 2017). The fact that the peridotitic fluids have a very different chemistry from the eclogitic fluids provides strong support for the equilibrium constants of the Al–Si complex shown in Figure 27b.

## Concluding remarks

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

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

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

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

## Acknowledgements

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

*Scientific editing by Maria Luce Frezzotti*