• 988 Accesses

• 1 Altmetric

• Metrics details

The SiO₂–H₂O binary system serves as a basis for understanding complex silicate-water systems. In this study, based on limited existing experimental data of solubility, we propose a new thermodynamic model for SiO₂–H₂O fluid by modifying the traditional non-random two-liquid model with a simplified polymerization reaction. This model is applicable from 773 K to the anhydrous quartz melting temperature and from 0.5 GPa to at least 2 GPa across the critical conditions. It can predict solid–liquid equilibrium and vapor–liquid equilibrium in good agreement with available experiments. The upper critical endpoint of the SiO₂–H₂O system is predicted to be at ~ 1.14 ± 0.18 GPa and 1344 ± 87 K. With the new model, we obtain a quantitative three-dimensional pressure–temperature–composition phase diagram of the SiO₂–H₂O fluid, which greatly facilitates the understanding of the complex phase behavior of this binary around the upper critical endpoint. In addition, since the model is based on the Gibbs free energy foundation, we further discuss the derived density variations of SiO₂–H₂O fluid along with its complex phase changes in typical geochemical processes.

The Earth's unique plate tectonic processes often involve a certain amount of water, making water–rock interactions in the Earth’s interiors important. Silicates are the dominant minerals in the crust and mantle, so thermodynamics of mixing silicate and water at high temperatures and pressures is critical to understanding deep fluid-rock interactions (Hack et al. 2007a, 2007b; Manning 2004; Tian et al. 2019). Unlike solution in water such as alkali salts, silicates generally have negligible solubilities in water at ambient conditions. With increase of temperature and pressure, however, their solubilities in water increase markedly, and so does the solubility of water in silicate melts. At sufficiently high temperature and pressure, aqueous fluids (labeled as V phase below) and hydrous silicate melts (labeled as L phase below) can be completely miscible and become supercritical fluids (SCF) (Boettcher and Wyllie 1969; Kennedy et al. 1962; Ni et al. 2017). Therefore, magma rising from the deep interiors would develop significantly altered geochemical and physical properties during their ascent, which may influence many geological processes such as volcanic venting and pegmatite formation (Mibe et al. 2002; Sowerby and Keppler 2002; Thomas 2000; Thomas and Davidson 2016).

As a prototype binary for the silicate-water systems, Fig. 1 outlines schematically the complex phase behavior in the SiO₂–H₂O system. From Fig. 1a, as the temperature rises, the concentration of SiO₂ in the aqueous fluid phase gradually increases. When the temperature reaches the fluid-saturated solidus (T₁), quartz begins to melt. With a sufficient amount of H₂O, hydrous melt and aqueous fluid coexist in vapor–liquid equilibrium (VLE). In cases of insufficient H₂O, quartz undergoes partial melting, resulting in solid–liquid equilibrium (SLE). The SLE ends at the anhydrous melting temperature (T₃). Between T₁ and T₂, aqueous fluid and hydrous melt coexist with increasing mutual solubility till their final miscibility.

Pressure has important effects on the phase boundaries. As the pressure increases, from P₁ to P₂ in Fig. 1a, the temperatures of the consolute point (T₂) and the hydrous melting point (T₁) decrease. The T–P variations of the consolute point and the hydrous melting point are schematically illustrated in Fig. 1b. These two curves eventually intersect at a point, which is called the upper critical endpoint (UCEP) (Hack et al. 2007b).

Data from existing experimental solubility measurements of SiO₂–H₂O can be roughly divided into two categories: SLE (the green ① and red ② curves in Fig. 1a) and VLE (the purple ④ and orange ⑤ curves in Fig. 1a). Available experimental data, with 875 data points, are summarized in Table 1. Most of these determinations were focused on the quartz solubility in aqueous solutions. Experimental data for hydrous silicate melts at pressures higher than 0.5 GPa are scarce. Only Holtz et al. (2000) measured several data points for water solubility in SiO₂-rich liquid, and there are no experimental data constraining the coexisting vapor phase. Experimental studies on the UCEP of the SiO₂–H₂O system have proven challenging. Kennedy et al. (1962) estimated it to be around 0.97 GPa and 1353 K, while some other studies concluded that the UCEP should be located at a pressure higher than 1.0 GPa (Mysen 1998; Stewart 1967). Nakamura (1974) discovered that SiO₂–H₂O fluids exhibit supercritical behavior at 1.5 GPa, thereby establishing the upper limit of the UCEP.

A few solubility models have been developed for SiO₂–H₂O fluids. Sverjensky et al. (2014) proposed the Deep Earth Water (DEW) model by refining the dielectric constant model at high pressures based on a well-behaved equation of state for water (Zhang and Duan 2005). The DEW model extends the temperature and pressure of the electrostatic model to 6 GPa and 1473 K. Other, empirical models based on linear relations with pure water density, have been proposed by a number of studies (Deng et al. 2021; Fournier and Potter 1982; Manning 1994; Shi et al. 2019; Wei et al. 2012), but these models are generally applied to aqueous fluids at low temperatures. Hunt and Manning (2012) developed a model (HM2012 model) based on the excess mixing properties in the SiO₂–H₂O fluids, which is applicable from 773 K to the anhydrous quartz melting temperature and from 0.5 to 2.0 GPa.

Acknowledging the potential disparity between local and overall mole fractions in the system, Wilson (1964) introduced the concept of local composition to the excess Gibbs free energy. This refinement facilitates characterizing equilibrium within polymer solutions.

The non-random two-liquid (NRTL) model introduces a correction factor based on the concept of local composition, which is advantageous in depicting the immiscibility in liquid systems (Renon and Prausnitz 1968). However, the traditional NRTL model is commonly used to study systems with organic compounds at low temperatures and pressures, and the interaction parameters are assumed to be pressure independent. To extend this model to the deep Earth environment, we modify its interaction parameters as a function of temperature and pressure. In addition, the NRTL model requires accounting for the interactions and local orderings between every two different species in the liquid. Therefore, the application of the traditional NRTL model to the silicate-water fluid systems seems to be unrealistic with too many parameters for fitting. Based on this, we propose a two-step model: in the first step, we quantify the mixing of SiO₂ and H₂O binary (i.e., no polymerization) with the traditional NRTL model; in the second step, using polymerization reactions to describe the generation of other new species. The two-step method has the advantage of greatly reducing the number of parameters that need to be fitted. It can accurately predict both SLE and VLE properties of SiO₂–H₂O and constrain the density that is still lacking from experimental measurements, thus paving the way to quantify various more complex silicate-water systems.

2.1 Fundamental two-step model for the fluid

We model the mixing of SiO₂ and H₂O in two steps. In the first step: we consider a binary system of N moles SiO₂–H₂O fluid with x_s mole fraction of SiO₂ and 1−x_s mole fraction of H₂O. In this step, we apply the traditional NRTL model (Additional file 1: Appendix A) and obtain the excess molar Gibbs free energy \({{g}_{NRTL}^{ex}}\) as follows: (Renon and Prausnitz 1968):

where R is the universal gas constant, τ₁₂ and τ₂₁ are dimensionless parameters to characterize the difference in interaction energy between different molecules pairs in the local composition, which are related to the interaction energy parameters G₁₂ and G₂₁ by: \({G}_{12}={\text{exp}}\left(-{\alpha }_{12}{\tau }_{12}\right)\) and \({G}_{21}={\text{exp}}\left(-{\alpha }_{21}{\tau }_{21}\right)\), where α₁₂ and α₂₁ (α₁₂ = α₂₁ = α) are the characteristic constants of the non-random mixture with an appropriate value between 0.2 and 0.47 (Prausnitz et al. 1998).

There is numerous experimental and computational evidence leading to the suggestion that SiO₂–H₂O system is dominated by monomers (Si (OH)₄) and dimers (Si₂O(OH)₆) at low temperatures and pressures (Mysen 2010; Zotov and Keppler 2000, 2002). However, with increasing temperature and pressure, other higher polymers may be present in significant concentrations (Doltsinis et al. 2007). It is possible to extend Eq. (1) and derive a multi-component NRTL model, but it would involve too many parameters. Alternatively, following the strategy by Hunt and Manning (2012), we further amend the free energy with a second step by considering the exchange of three types of oxygen (oxygen in free water molecules (H₂O), bridging oxygen (\({{\text{O}}}_{{\text{br}}}\)), and terminal hydroxyl oxygen of silanol groups (\({{\text{OH}}}_{{\text{tm}}}\))) in a simplified polymerization reaction (Gerya et al. 2005; Hunt and Manning 2012; Plyasunov and Zyubina 2021):

If we assume that the reaction (2) produces M moles\({{\text{OH}}}_{{\text{tm}}}\). Then the reaction involves 0.5M moles H₂O and 0.5M moles\({{\text{O}}}_{{\text{br}}}\). This system contains N(1−x_s)−0.5M moles H₂O and 2Nx_s−0.5M moles\({{\text{O}}}_{{\text{br}}}\). Therefore, for the effectively ternary system of H₂O, \({{\text{O}}}_{{\text{br}}}\) and\({{\text{OH}}}\), the system contains N(1 + x_s) moles oxygen-containing components and the concentrations of H₂O, \({{\text{O}}}_{{\text{br}}}\) and \({{\text{OH}}}_{{\text{tm}}}\) can be expressed as:\(\frac{N\left(1-{x}_{{\text{s}}}\right)-0.5M}{N\left(1+{x}_{{\text{s}}}\right)}\),\(\frac{2N{x}_{{\text{s}}}-0.5M}{N\left(1+{x}_{{\text{s}}}\right)}\),\(\frac{M}{N\left(1+{x}_{{\text{s}}}\right)}\), respectively. If we denote \({y}_{{{\text{OH}}}^{ }}=\frac{M}{N\left(1+{x}_{{\text{s}}}\right)}\), then H₂O and \({{\text{O}}}_{{\text{br}}}\) concentration can be simplified as \(\frac{1-{x}_{{\text{s}}}}{1+{x}_{{\text{s}}}}-0.5{y}_{{{\text{OH}}}}\) and\(\frac{2{x}_{{\text{s}}}}{1+{x}_{{\text{s}}}}-0.5{y}_{{{\text{OH}}}}\), respectively. Assuming that all the non-ideality has been accounted for with the NRTL model, the equilibrium constant K of Eq. (2) can be expressed as:

where \(\Delta {g}_{{\text{rec}}}^{\circ }\) is the standard molar Gibbs free energy change of Eq. (2) at any temperature and pressure.

Then, the molar excess Gibbs free energy of this SiO₂–H₂O binary can be expressed as:

Note that the \({y}_{{{\text{OH}}}}\) in Eq. (4) can be further expressed as a function of x_s and \(\Delta {g}_{{\text{rec}}}^{\circ }\) by solving Eq. (3). Now we have established a fundamental model for the SiO₂–H₂O fluid with only four hyper-parameters (α, \({\tau }_{12},{\tau }_{21} {\text{and}} \Delta {g}_{{\text{rec}}}^{\circ }\)).

2.2 Phase equilibrium calculations and the model regression

The key property for the phase equilibrium is the activity coefficient, which can be derived from excess Gibbs free energy by \(RT{\text{ln}}{\gamma }_{i}={g}^{{\text{ex}}}+\left(1-{x}_{i}\right)\frac{\partial {g}^{ex}}{\partial {x}_{i}}\) (Levine 2008; Ma 2001). From Eq. (4), we obtain the following expressions:

\({y}_{{{\text{OH}}}}\) can be solved by Eq. (3) and calculating its derivative to x_s, the \(\frac{\partial {y}_{{{\text{OH}}}^{ }}}{\partial {x}_{{\text{s}}}}\) in Eqs. (5–6) can be further expressed as a function of x_s and \(\Delta {g}_{{\text{rec}}}^{\circ }\).

For the solid–liquid equilibrium (SLE), the chemical potential of quartz and the liquid-phase SiO₂ are equal:

where \({\mu }_{{\text{quartz}}}^{\circ }\) is the standard Gibbs free energy of quartz (Berman 1988), \({\mu }_{{{\text{SiO}}}_{2},{\text{liq}}}^{\circ }\) is the standard Gibbs free energy of SiO₂ in the liquid phase (Mao et al. 2001), X_s is the determined solubility of SiO₂ in mole fraction.

For the vapor–liquid equilibrium (VLE) between aqueous fluid and hydrous silica melt, we have

where \({X}_{{\text{s}}}^{{\text{V}}}+{X}_{{{\text{H}}}_{2}{\text{O}}}^{{\text{V}}}=1\) and \({X}_{{\text{s}}}^{{\text{L}}}+{X}_{{{\text{H}}}_{2}{\text{O}}}^{{\text{L}}}=1\).

To describe phase equilibria over a broad range of temperature and pressure, while α is assumed to be constant, the remaining three hyper-parameters (\(\Delta {g}_{{\text{rec}}}^{\circ },{\tau }_{12},{\tau }_{21}\)) in the fundamental model are modelled as functions of temperature and pressure:

where T is temperature in Kelvin, P is pressure in bars, and \(\Delta {g}_{{\text{rec}}}^{\circ }\) is the standard molar Gibbs free energy change in J/mol.

From the available experimental solubility data, we obtained parameters of the equilibrium model with nonlinear least squares regressions and Bayesian optimization (Wang and Yang 2023). In this process, we tested the statistical significance of the fitting parameters, eliminated one parameter that was statistically insignificant, and finally obtained the fitting values of 13 independent parameters (Table 2).

3.1 Comparison with experimental data

As shown in Table 1, there are 875 experimental solubility data points available for the SiO₂–H₂O system, covering a wide range of temperature and pressure ranges (0.00003–2 GPa, 281.15–1674.15 K). Most of the data are in the low pressure-low temperature region (T  773 K or P  0.5 GPa). However, the purpose of our model is to characterize the fluid properties near the upper critical endpoint, so we focused on a more appropriate regime with pressure higher than 0.5 GPa. A total of 204 available experimental data points were collected (temperature from 773 K to quartz dry melting temperature and pressure from 0.5 to 2.0 GPa). It should be noted that certain experimental data points may not be reliable. For instance, Manning (1994) demonstrated that the data of Anderson and Burnham (1965) had quench-related issues leading to underestimation of quartz solubility, so our model does not consider these data. There were finally 133 experimental data used to build our model, some of the experimental data were utilized for parameter regression, while the remaining data were used to validate the fitting results (Table 3). Table 3 also summarizes the averaged absolute deviations (AAD) and averaged relative deviations (AARD) of the model predictions. Figure 2 illustrates the modeled solubilities and the experimental data. The overall AAD and AARD are 0.02 and 16%, respectively. Furthermore, except for several outlier data by Kennedy et al. (1962) and Holtz et al. (2000), the model gives good predictions of the measured solubilities accurate to 19%.

Figure 3 presents comparison of our and previous models for quartz solubility under various temperature or pressure conditions. The selected models include widely used density models or empirical models, namely F1982 (Fournier and Potter 1982), M1994 (Manning 1994), WDM2012 (Wei et al. 2012), SMH2019 (Shi et al. 2019), DWP2021 (Deng et al. 2021), G2005 (Gerya et al. 2005), DM2010 (Dolejs and Manning 2010); the electrostatic model DEW2014 (Sverjensky et al. 2014) and the solute–solvent interaction model HM2012 (Hunt and Manning 2012). Under constant temperature conditions, F1982 agrees well with the experimental data at low pressures. However, as the pressure increases, it deviates significantly from the experimental data and underestimates the solubility of quartz in water. At 873 K, G2005 agrees well with experimental data at low pressure, but it overestimates the solubility of quartz as the pressure increases (Fig. 3a). By contrast, both M1994 and DEW2014 models tend to overestimate the solubility of quartz. Additionally, models that demonstrate good performance in low-temperature regions (such as G2005, DM2010, WDM2012, DWP2021 and SMH2019) no longer maintain their accuracy when the temperature extends beyond 1173 K. Unlike density and electrostatic models, models based on interactions between molecules are easier to scale to higher temperatures. Therefore, over a wide range of temperatures and pressures, the HM2012 model and the two-step model are in good agreement with the experimental data (Fig. 3d). While these two models indeed show quite similar accuracies, as shown in Fig. 4, slight differences can be discriminated: from Fig. 4a, when the pressure is 1 GPa, HM2012 predicts that the system is already supercritical, while our model predicts that the miscible gap still exists (green line in Fig. 4a). When the pressure is 1.5 GPa and 2 GPa, the predicted results of the two models agree well with the available experimental data (Fig. 4b, c), the fluid is supercritical fluid, and the miscibility gap disappears.

As shown in Fig. 5a, we estimate the location of the upper critical endpoint to be around 1.14 ± 0.18 GPa and 1344 ± 87 K, while it is predicted by HM2012 to be about 0.93 GPa, 1340 K. The predictions of these two models for the upper critical endpoint as well as the hydrous melting curve temperatures agree with the measurements of Kennedy et al. (1962) (approximately 0.97 GPa and 1353 K). Besides these similarities, it is evident that the miscibility gap predicted by HM2012 is larger than that evaluated in the present study. A comparison of the X–T relationships (Fig. 5b) suggests that the two models show some differences in VLE: HM2012 tends to be more asymmetric with the critical curve leaning towards lower concentration of silicate components. As mentioned previously, experimental measurements of hydrous melt solubility for VLE are scarce, with only one Holtz et al. (2000) measurement collected at 0.6 GPa (Fig. 5b).

3.2 Phase diagram of the SiO₂–H₂O system

Our thermodynamic model leads to a self-consistent phase diagram (Fig. 6). Figure 5a, b corresponds to the T–P and X–T projections, respectively. With the P–X projection in Fig. 6, the VLE boundary changes at different temperatures and pressures also can be quantified.

Isopleths of SiO₂ solubility in the fluid across the critical region around the upper critical endpoint (UCEP) is shown in Fig. 7. As noticed by Hack et al (2007a, b), the aqueous fluid that resides on the most left side of the hydrous melting curve often involves small amount of silicate solute. The solubility of SiO₂ is generally less than 0.1. Nevertheless, in the vicinity of the hydrous melting curve, particularly at high pressures, the solubility of SiO₂ increases remarkably. On the other side of the hydrous melting curve, the situation is more complex with several phase regimes overlapping from this point of view. As indicated by the thin black solid lines, the concentration of SiO₂ in the fluid could be high when it is equilibrated with quartz.

The vapor–liquid equilibrium (VLE) regime between the hydrous melting curve and the critical curve greatly redefines the compositional range of the fluid. This VLE regime determines the lower limit of SiO₂ in the hydrous melt (gray isopleths), on the other hand, at the same pressure it accommodates much higher concentration of SiO₂ in the aqueous fluid (dashed isopleths) as compared with that at lower temperature (left side of the hydrous melting curve). As will be discussed in the next section, the variations of SiO₂ solubility would significantly affect the density of the fluid and therefore are important for the transportations of SiO₂–H₂O fluids.

4.1 Density of the SiO₂–H₂O fluid

Density is the key property that controls the physical state and transportation of fluid. Experimental evidence indicates that the density of SiO₂–H₂O fluids depends on the SiO₂ concentration and is correlated with species behavior (Mysen 2010). However, directly measuring fluid density in experiments is challenging. With the model proposed in this study, it is possible to constrain the density of the SiO₂–H₂O binary.

By calculating the partial derivative of the excess free energy \({g}_{{\text{tot}}}^{{\text{ex}}}\) in Eq. (4) to the pressure, we obtain the excess volume V^ex:

Then the volume of SiO₂–H₂O fluid can be expressed as:

where \({V}_{{{\text{H}}}_{2}{\text{O}}}^{\circ }\) and \({V}_{{{\text{SiO}}}_{2}}^{\circ }\) are the standard molar volumes of liquid H₂O and liquid SiO₂ as given by Zhang and Duan (2005) and Mao et al. (2001), respectively. The fluid density can be easily obtained from the volume:

where \({M}_{{{\text{H}}}_{2}{\text{O}}}\) and \({M}_{{{\text{SiO}}}_{2}}\) are the standard molar mass of liquid H₂O and liquid SiO₂.

However, experimental errors and model parameter errors may affect density predictions. We used the Monte Carlo method by deploying multiple sampling statistics within experimental and parameter uncertainty, yielding an uncertainty for the density of the fluid. The maximum uncertainty is 0.06 g/cm³. In Fig. 6, we show the predicted density of quartz-saturated and vapor-saturated SiO₂–H₂O fluid with colored contours. Obviously, the density of the fluid is mainly controlled by the concentration of dissolved SiO₂ (X_s). The effects of temperature and pressure seem to be minor, but they may indirectly regulate the density by changing the composition of the fluid, which may result in counterintuitive variation of density with temperature or pressure. This can be understood by relating it with the solubility profile as shown in Fig. 7: the aqueous fluid would dissolve much more solute as the temperature approaching the UCEP and this would greatly compensate the intrinsic decreasing density effect of temperature.

Figure 8 shows the comparison and analysis of the fluid density predicted by our two-step model (solid line) and the HM2012 model (dashed line). While the density of aqueous fluids is increased with pressure as expected, the abovementioned reverse trend of temperature can be observed at higher temperatures and pressures (Fig. 8a, 1223 K). For quartz-saturated hydrous melts, the two-step model is the opposite of what the HM2012 model predicts (dashed lines) showing that its density decreases with increasing pressure. All these features can also be found in Fig. 8b, which reveals the variation of density along various isobars. While those for quartz-saturated aqueous fluids reveal minor disparity, which explains the resulted close agreements between the two models as for the densities of aqueous fluids, the volume of saturated hydrous silicate melts show noticeable deviations especially under higher pressures. The HM2012 model predicts a more rapid decline in the volume of hydrous silicate melts with pressure, thus compensating for the decreased mass with less silica in the fluid at higher pressures (Fig. 8c, d). The two-step model proposed in this study, on the other hand, predicts a less perturbed volume and the density of the melts is therefore decreased at higher pressures. In summary, experimental density data for silicate solutions are urgently needed to clarify the deviations among the existing thermodynamic models.

4.2 Thermodynamics of water exsolution from the rising magma

The exsolution of hydrothermal fluids from magma controls the many aspects of eruption behavior of volcanoes. In the deep Earth, the immiscible transition state (transition process) fluid may cause deep source earthquakes (Wech et al. 2020). This fluid also is an important medium for the formation of many metallic mineral deposits (Hedenquist and Lowenstern 1994). However, the mechanism of magma differentiation of immiscible fluids (often named as “boiling”) is still unclear. This process occurred in the Earth’s interior is often classified into two types (Fig. 9a): (1) The pressure in the magma chamber is suddenly reduced due to tectonic collapse or the rapid rising of magma, so that the unsaturated fluid at depth quickly decreases the solubility of solute during the ascending path and causes the exsolution. This type of process is called “first boiling”, which involves essentially invariant temperature. (2) The unsaturated magma is temporarily stored at a certain depth and thus the pressure keeps unchanged. As the temperature decreases, the melt gradually crystallizes and the water concentration in the remaining melt gradually increases, eventually reaching the saturation of water and causing exsolution. This type of process is known as the “second boiling” (Candela 1997; Heinrich 2007). Considering that the SiO₂–H₂O system is too basic compared to the real complex silicate-water system, we only try to discuss the changes in the properties of the fluids due to the first and second boiling.

The schematic thermodynamic paths for the first boiling process is illustrated in Fig. 9b (dark red, light red and gray lines marked as ①, representing three different initial conditions as discussed below) and the second boiling process (purple line marked as ②).

For a fluid-unsaturated melt (dark red line in Fig. 9b), its density would gradually decrease as the pressure decreases before the boiling process begins (blue solid line in Fig. 9c), when the melt reaches the fluid saturation and begins to exsolve aqueous fluid (Fig. 9b as the path touches the VLE surface). With further pressure decrease, more dilute aqueous fluid is exsolved and the remaining melt becomes more concentrated in SiO_2. This process results in decreasing density for the exsolved aqueous fluid and increasing density for the residual melt respectively (green and orange curves in Fig. 9c). When the pressure is further reduced to the hydrous melting curve, the quartz begins to crystallize and the first boiling stops.

Magma is usually thought to be stored at its neutral buoyancy level, where the density of the magma is equivalent to that of the host rock (Ryan 1987). However, another hypothesis suggests that the magma storage depth is influenced by the water content (Annen et al. 2006; Zellmer et al. 2016). As shown in Fig. 9b, when the initial magma at the same temperature and pressure has a higher water content (the gray line in Fig. 9b), it will reach the VLE surface and start the boiling process at higher pressure (Fig. 9c). Therefore, wetter magma would boil more readily than dry magma at depth. Rasmussen et al, (2022) found magma storage depths consistent with predicted water degassing depths. Although the wet magma contains higher water content and has a lower density, the degassing process leads to an increase in viscosity, allowing it to stay deeper during ascent.

When the magma of the same composition and depth initially is slightly cooler (the light red line in Fig. 9b), the boiling process will be started at a higher pressure. In addition, as more water enters the aqueous fluid, resulting in a larger density difference between the exsolved aqueous fluid and the residual melt, this would drive the boiling process (Fig. 9d). With further decreasing of the pressure till the remaining magma reaches the hydrous melting curve, the boiling process would stop more quickly. Therefore, the rising magma with lower initial temperature would also finish the boiling process at a deeper location.

The purple line in Fig. 9b shows the second boiling process. The temporarily stored magma becomes denser as the temperature gradually decreases (blue curve in Fig. 9e). When the temperature reaches the liquidus curve (quartz-saturated), the melt begins to crystallize. As the temperature decreases further, more and more solids would crystallize and the water content in the residual melt increases gradually. This process results in decreasing of density (orange curve in Fig. 9e). When the temperature decreases to the hydrous melting curve, the second boiling starts until most of the SiO₂ in the fluid has crystallized (grey curve in Fig. 9e). Obviously, when the fluid becomes saturated during this process, the variation of its density should be similar to the trend as revealed in Fig. 9b.

The SiO₂–H₂O binary system is the basis for understanding complex silicate-water systems, but relevant high temperature and pressure and density data remains scarce. This study considers the vapor–liquid equilibrium region using the NRTL model as the basis and modifies it with aggregated reactions to establish a two-step model that can simultaneously predict SLE and VLE. The two-step model is fairly accurate and applicable from 773 K to the anhydrous quartz melting temperature and from 0.5 GPa to at least 2 GPa. Based on the model proposed in this study, we provide a quantitative 3D phase diagram of SiO₂–H₂O, which greatly helps our understanding of the complex phase relations of this binary across the critical conditions. The upper critical endpoint is predicted to be around 1.14 ± 0.18 GPa and 1344 ± 87 K.

Because the two-step model quantifies the Gibbs free energy of the SiO₂–H₂O fluid under different conditions, we observe nontrivial different trends of density along various P–T paths. The temperature and pressure are vital for the density through indirectly regulating the SiO₂ dissolved in the fluid. By further considering the complex VLE, such as the boiling processes when the volatile exsolution occurs in the rising magma, we find distinctive thermodynamics would be involved along different paths.

Finally, it should be noted that thermodynamic models inevitably rely on the experimental data. For the SiO₂–H₂O fluids concerned in this study, experimental determination of VLE and density of SiO₂–H₂O fluid would be needed in the future.

NRTL:

Non-random two-liquid

V:

Aqueous fluids

L:

Hydrous silicate melts

S:

Solid (quartz)

SLE:

Solid–liquid equilibrium

VLE:

Vapor–liquid equilibrium

SCF:

Supercritical fluids

UCEP:

Upper critical endpoint

AAD:

Averaged absolute deviations

AARD:

Averaged relative deviations

• Anderson GM, Burnham CW (1965) The solubility of quartz in supercritical water. Am J Sci 263:494–511

  Article  CAS  Google Scholar 

• Anderson GM, Burnham CW (1967) Reactions of quartz and corundum with aqueous chloride and hydroxide solutions at high temperatures and pressures. Am J Sci 265:12–27

  Article  CAS  Google Scholar 

• Annen C, Blundy JD, Sparks RSJ (2006) The genesis of intermediate and silicic magmas in deep crustal hot zones. J Petrol 47:505–539

  Article  CAS  Google Scholar 

• Berman RG (1988) Internally consistent thermodynamic data for minerals in the system Na₂O-K₂O-CaO-MgO-FeO-Fe₂O₃-Al₂O₃-SiO₂-TiO₂-H₂O-CO₂. J Petrol 29:445–522

  Article  CAS  Google Scholar 

• Boettcher AL, Wyllie PJ (1969) Phase relationships in system NaAlSiO₄-SiO₂-H₂O to 35 kilobars pressure. Am J Sci 267:875–909

  Article  CAS  Google Scholar 

• Candela PA (1997) A review of shallow, ore-related granites: textures, volatiles, and ore metals. J Petrol 38:1619–1633

  Article  CAS  Google Scholar 

• Crerar DA, Anderson GM (1971) Solubility and solvation reactions of quartz in dilute hydrothermal solutions. Chem Geol 8:107–122

  Article  CAS  Google Scholar 

• Deng WP, Wei Q, Liu X (2021) A new density model of quartz solubility in H₂O-CO₂-NaCl ternary systems up to high temperatures and high pressures. Geofluids 2021(4):1–23

  Google Scholar 

• Dolejs D, Manning CE (2010) Thermodynamic model for mineral solubility in aqueous fluids: theory, calibration and application to model fluid-flow systems. Geofluids 10:20–40

  Article  CAS  Google Scholar 

• Doltsinis NL, Burchard M, Maresch WV, Boese AD, Fockenberg T (2007) Ab initio molecular dynamics study of dissolved SiO₂ in supercritical water. J Theor Comput Chem 6:49–62

  Article  CAS  Google Scholar 

• Fournier RO, Potter RW (1982) An equation correlating the solubility of quartz in water from 25° to 900°C at pressures up to 10,000 bars. Geochim Cosmochim Acta 46:1969–1973

  Article  CAS  Google Scholar 

• Gerya TV, Maresch WV, Burchard M, Zakhartchouk V, Doltsinis NL, Fockenberg T (2005) Thermodynamic modeling of solubility and speciation of silica in H₂O-SiO₂ fluid up to 1300°C and 20 kbar based on the chain reaction formalism. Eur J Mineral 17:269–283

  Article  CAS  Google Scholar 

• Gunnarsson I, Arnorsson S (2000) Amorphous silica solubility and the thermodynamic properties of H₄SiO₄ degrees in the range of 0 °C to 350°C at P-sat. Geochim Cosmochim Acta 64:2295–2307

  Article  CAS  Google Scholar 

• Hack AC, Hermann J, Mavrogenes JA (2007a) Mineral solubility and hydrous melting relations in the deep earth: analysis of some binary A-H₂O system pressure-temperature-composition topologies. Am J Sci 307:833–855

  Article  CAS  Google Scholar 

• Hack AC, Thompson AB, Aerts M (2007b) Phase relations involving hydrous silicate melts, aqueous fluids, and minerals. Rev Mineralal Geochem 65:129–185

  Article  CAS  Google Scholar 

• Hedenquist JW, Lowenstern JB (1994) The role of magmas in the formation of hydrothermal ore-deposits. Nature 370:519–527

  Article  CAS  Google Scholar 

• Heinrich CA (2007) Fluid-fluid interactions in magmatic-hydrothermal ore formation, Fluid-Fluid Interactions, 363–387

• Hemley JJ, Montoya JW, Marinenko JW, Luce RW (1980) Equilibria in the system Al₂O₃-SiO₂-H₂O and some general implications for alteration/mineralisation processes. Econ Geol 75:210–228

  Article  CAS  Google Scholar 

• Holtz F, Roux J, Behrens H, Pichavant M (2000) Water solubility in silica and quartzofeldspathic melts. Am Mineral 85:682–686

  Article  CAS  Google Scholar 

• Hunt JD, Manning CE (2012) A thermodynamic model for the system SiO₂-H₂O near the upper critical end point based on quartz solubility experiments at 500–1100 degrees C and 5–20 kbar. Geochim Cosmochim Acta 86:196–213

  Article  CAS  Google Scholar 

• Kennedy GC (1950) A portion of the system silica-water. Econ Geol 45:629–653

  Article  CAS  Google Scholar 

• Kennedy GC, Wasserburg GJ, Heard HC, Newton RC (1962) The upper three-phase region in the system SiO₂–H₂O. Am J Sci 260:501–521

  Article  Google Scholar 

• Levine IN (2008) Physical chemistry, 6th edn. McGraw-Hill Education, New York

  Google Scholar 

• Ma ZT (2001) Thermodynamic description for concentrated metallic solutions using interaction parameters. Metall Mater Trans B Process Metall Mater Process Sci 32:87–103

  Article  Google Scholar 

• Mackenzi FT, Gees R (1971) Quartz: synthesis at earth-surface conditions. Science 173:533–535

  Article  Google Scholar 

• Manning CE (1994) The solubility of quartz in H₂O in the lower crust and upper mantle. Geochim Cosmochim Acta 58:4831–4839

  Article  CAS  Google Scholar 

• Manning CE (2004) The chemistry of subduction-zone fluids. Earth Planet Sci Lett 223:1–16

  Article  CAS  Google Scholar 

• Manning CE, Boettcher SL (1994) Rapid-quench hydrothermal experiments at mantle pressures and temperatures. Am Mineral 79:1153–1158

  CAS  Google Scholar 

• Mao HH, Sundman B, Wang ZW, Saxena SK (2001) Volumetric properties and phase relations of silica—thermodynamic assessment. J Alloy Compd 327:253–262

  Article  CAS  Google Scholar 

• Mibe K, Fujii T, Yasuda A (2002) Composition of aqueous fluid coexisting with mantle minerals at high pressure and its bearing on the differentiation of the Earth’s mantle. Geochim Cosmochim Acta 66:2273–2285

  Article  CAS  Google Scholar 

• Morey GW, Fournier RO, Rowe JJ (1962) The solubility of quartz in water in the temperature interval from 25° to 300° C. Geochim Cosmochim Acta 26:1029–1043

  Article  CAS  Google Scholar 

• Mysen BO (1998) Interaction between aqueous fluid and silicate melt in the pressure and temperatures regime of the Earth’s crust and upper mantle. Neues Jahrb Mineral Abh 172:227–244

  CAS  Google Scholar 

• Mysen BO (2010) Speciation and mixing behavior of silica-saturated aqueous fluid at high temperature and pressure. Am Mineral 95:1807–1816

  Article  CAS  Google Scholar 

• Nakamura Y (1974) The system SiO₂-H₂O-H₂ at 15 kbar. Carnegie Institution of Washington Year Book 73:259–263

• Newton RC, Manning CE (2000) Quartz solubility in H₂O-NaCl and H₂O-CO₂ solutions at deep crust-upper mantle pressures and temperatures: 2–15 kbar and 500–900 °C. Geochim Cosmochim Acta 64:2993–3005

  Article  CAS  Google Scholar 

• Newton RC, Manning CE (2008) Thermodynamics of SiO₂-H₂O fluid near the upper critical end point from quartz solubility measurements at 10 kbar. Earth Planet Lett 274:241–249

  Article  CAS  Google Scholar 

• Ni HW, Zhang L, Xiong XL, Mao Z, Wang JY (2017) Supercritical fluids at subduction zones: evidence, formation condition, and physicochemical properties. Earth Sci Rev 167:62–71

  Article  CAS  Google Scholar 

• Novgorodov PG (1977) Solubility of quartz in mixture in of H₂O-CO₂ and H₂O-NaCl and at 700°C and at pressure of 1.5 Kbar. Geokhimiy 8:1270–1273

  Google Scholar 

• Plyasunov AV, Zyubina TS (2021) Thermodynamic modeling of the solubility of quartz in water up to high temperatures and pressures. Earth Space Chem 5:941–959

  Article  CAS  Google Scholar 

• Prausnitz JM, Lichtenthaler RN, de Azevedo EG (1998) Molecular thermodynamics of fluid-phase equilibria. 2nd Edn. Prentice-Hall, Upper Saddle River, NJ, 1986

• Ragnarsdottir KV, Walther JV (1983) Pressure sensitive “silica geothermometer” determined from quartz solubility experiments at 250 °C. Geochim Cosmochim Acta 47:941–946

  Article  CAS  Google Scholar 

• Rasmussen DJ, Plank TA, Roman DC, Zimmer MM (2022) Magmatic water content controls the pre-eruptive depth of arc magmas. Science 375:1169–1172

  Article  CAS  Google Scholar 

• Renon H, Prausnitz JM (1968) Local compositions in thermodynamic excess functions for liquid mixtures. Aiche J 14:135–144

  Article  CAS  Google Scholar 

• Rimstidt JD (1997) Quartz solubility at low temperatures. Geochim Cosmochim Acta 61:2553–2558

  Article  CAS  Google Scholar 

• Ryan MP (1987) Neutral buoyancy and the mechanical evolution of magmatic systems. The Geochemical Society, Special Publication No. 1, 259–287

• Shi XL, Mao SD, Hu JW, Zhang J, Zheng JX (2019) An accurate model for the solubilities of quartz in aqueous NaCl and/or CO₂ solutions at temperatures up to 1273 K and pressures up to 20,000 bar. Chem Geol 513:73–87

  Article  CAS  Google Scholar 

• Shmulovich K, Graham C, Yardley B (2001) Quartz, albite and diopside solubilities in H₂O-NaCl and H₂O-CO₂ fluids at 0.5-0.9 GPa. Contrib Mineral Petrol 141:95–108

  Article  CAS  Google Scholar 

• Shmulovich KI et al (2006) Solubility of quartz in crustal fluids: experiments and general equations for salt solutions and H₂O-CO₂ mixtures at 400–800 °C and 0.1-0.9 GPa. Geofluids 6:154–167

  Article  CAS  Google Scholar 

• Siever R (1962) Silica solubility, 0°C-200°C, and the diagenesis of siliceous sediments. J Geol 70:127–150

  Article  CAS  Google Scholar 

• Sommerfeld RA (1967) Quartz solution reaction: 400°C–500°C, 1000 bars. J Geophys Res 72:4253

  Article  CAS  Google Scholar 

• Sowerby JR, Keppler H (2002) The effect of fluorine, boron and excess sodium on the critical curve in the albite-H₂O system. Contrib Mineral Petrol 143:32–37

  Article  CAS  Google Scholar 

• Stewart DB (1967) Four-phase curve in the system CaAl₂Si₂O₈-SiO₂-H₂O between 1 and 10 kilobars. Schweiz Mineral Petrogr Mitt 47:35–39

  CAS  Google Scholar 

• Sverjensky DA, Harriaon B, Azzolini D (2014) Water in the deep Earth: the dielectric constant and the solubilities of quartz and corundum to 60 kb and 1200°C. Geochim Cosmochim Acta 129:125–145

  Article  CAS  Google Scholar 

• Thomas R (2000) Determination of water contents of granite melt inclusions by confocal laser Raman microprobe spectroscopy. Am Mineral 85:868–872

  Article  CAS  Google Scholar 

• Thomas R, Davidson P (2016) Revisiting complete miscibility between silicate melts and hydrous fluids, and the extreme enrichment of some elements in the supercritical state—consequences for the formation of pegmatites and ore deposits. Ore Geol Rev 72:1088–1101

  Article  Google Scholar 

• Tian M, Katz RF, ReeJones DW (2019) Devolatilization of Subducting Slabs, Part I: Thermodynamic parameterization and open system effects. Geochem Geophys Geosyst 20:5667–5690

  Article  CAS  Google Scholar 

• Vanlier JA, Debruyn PL, Overbeek JTHG (1960) The solubility quartz. J Phys Chem 64:1675–1682

  Article  CAS  Google Scholar 

• Walther JV, Orville PM (1983) The extraction quench technique for determination of the thermodynamic properties of solute complexes: application to quartz solubility in fluid mixtures. Am Mineral 68:731–741

  CAS  Google Scholar 

• Wang H, Yang K (2023) Bayesian optimization. In: Brockhoff D, Emmerich M, Naujoks B, Purshouse R (eds) Many-criteria optimization and decision analysis: state-of-the-art, present challenges, and future perspectives. Springer International Publishing, Cham,  271–297

  Chapter  Google Scholar 

• Wang HM, Henderson G, Brenan J (2004) Measuring quartz solubility by in situ weight-loss determination using a hydrothermal diamond cell. Geochim Cosmochim Acta 68:5197–5204

  Article  CAS  Google Scholar 

• Watson EB, Wark DA (1997) Diffusion of dissolved SiO₂ in H₂O at 1 GPa, with implications for mass transport in the crust and upper mantle. Contrib Mineral Petrol 130:66–80

  Article  CAS  Google Scholar 

• Wech AG, Thelen WA, Thomas HM (2020) Deep long-period earthquakes generated by second boiling beneath Mauna Kea volcano. Science 368:775–779

  Article  CAS  Google Scholar 

• Wei Q, Duan ZH, Mao SD (2012) A thermodynamic model of quartz solubility in H₂O-CO₂-NaCl systems up to 1000 °C and 1.5GPa. Acta Petrol Sin 28:2656–2666

  CAS  Google Scholar 

• Weill DF, Fyfe WS (1964) The solubility of quartz in H₂O in the range 1000–4000 bars and 400–550°C. Geochim Cosmochim Acta 28:1243–1255

  Article  Google Scholar 

• Wilson GM (1964) Vapor-liquid equilibrium. X I. A new expression for the excess free energy of mixing. J Am Chem Soc 86:127–130

  Article  CAS  Google Scholar 

• Xie ZX, Walther JV (1993) Quartz solubilities in NaCl solutions with and without wollastonite at elevated temperatures and pressures. Geochim Cosmochim Acta 57:1947–1955

  Article  CAS  Google Scholar 

• Zellmer GF, Pistone M, Iizuka Y, Anderws BJ, Gómez-Tuena A, Straub SM, Cottrell E (2016) Petrogenesis of antecryst-bearing arc basalts from the Trans-Mexican Volcanic Belt: insights into along-arc variations in magma-mush ponding depths, H₂O contents, and surface heat flux. Am Mineral 101:2405–2422

  Article  Google Scholar 

• Zhang ZG, Duan ZH (2005) Prediction of the PVT properties of water over wide range of temperatures and pressures from molecular dynamics simulation. Phys Earth Planet Inter 149:335–354

  Article  CAS  Google Scholar 

• Zotov N, Keppler H (2000) In-situ Raman spectra of dissolved silica species in aqueous fluids to 900 °C and 14 kbar. Am Mineral 85:600–603

  Article  CAS  Google Scholar 

• Zotov N, Keppler H (2002) Silica speciation in aqueous fluids at high pressures and high temperatures. Chem Geol 184:71–82

  Article  CAS  Google Scholar 

Download references

Authors and Affiliations

1. Key Laboratory of Earth and Planetary Physics, Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing, 100029, China

   Hui Gao, Guoge Li & Zhigang Zhang

2. College of Earth and Planetary Sciences, University of Chinese Academy of Sciences, Beijing, 100049, China

   Hui Gao, Guoge Li & Zhigang Zhang

Contributions

ZGZ proposed the topic, conceived and designed the study. HG constructed the model and was a major contributor in writing the manuscript. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Zhigang Zhang.

Competing interests

We declare that we have no competing interests.

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Reprints and permissions