A New Desalination Pump Helps Define the pH of Ocean Worlds

Brine is much less dense than rock. Therefore, Earth's saline ocean stably floats on the oceanic lithosphere. For the case of a water-rich planet, this may not be so, due to the much lower density difference between that of brine and of the high-pressure water-ice polymorphs that compose the ocean's bottom surface. Therefore, gravity may impose some restrictions on the possible maximum salinity of an ocean world.

Journaux et al. (2013) examined the melting-point depression of ices VI and VII for four NaCl aqueous solution compositions of 0.01 molal, 1 molal, 2.5 molal, and 4 molal. They noticed that when the brine concentration was increased from 1 molal to 2.5 molal, there was an inversion in behavior, where ice VI became less dense than the brine and floated in the chamber. However, Ice VII was denser than the brine for all examined samples.

In Figure 6, we reproduce the phase diagram for the binary system H₂O–NaCl at room temperature from Adams (1931). If during planetary formation strong water–rock interactions at high temperature resulted in a primordial ocean saltier than 1.75 ± 0.75 molal, then ice VII, rather than ice VI, would form the gravitationally stable ocean floor. Fluctuations of density in the ocean would cause downward migration of some liquid parcels. In Figure 6, we give an example for the compositional evolution of such a liquid parcel (see red arrows). The sinking parcel will not change its composition until it reaches the liquidus, at which point grains of ice VI begin to form within the parcel. These grains have some NaCl dissolved within them according to the solvus conditions, not shown in the figure. However, most of the salt remains dissolved in the liquid parcel. Therefore, the now denser parcel continues to sink, its composition restrained to the liquidus. On the other hand, the ice VI grains are positively buoyant and thus migrate upward, crossing their melting condition and producing a less saline upper ocean layer. The liquid parcel sinks while following the liquidus until it reaches eutectic composition. Beyond that point, its salt content largely forms grains of halite, whose high density (see Figure 7) ensures they descend down to the ice VII layer. The ice mantle overturn then sediments these pure salt grains onto the ice–rock mantle interface.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

This process is likely very efficient in removing salt from the ocean, if indeed its initial salinity was high. When the ocean's average salinity drops below 1.75 ± 0.75 molal, the ice VI grains become negatively buoyant in the background ocean. Therefore, an ice VI layer begins to form above the ice VII layer. Consequently, the ocean desalination process as described here ceases. This gravitational upper limit on the oceanic salt concentration is however not very stringent, in the sense that it is higher than the concentration of NaCl in Earth's seawater, which is approximately 0.5 molal (in Earth's seawater, Na⁺ and Cl⁻ concentrations are 0.47 molal and 0.55 molal, respectively; see Zeebe & Wolf-Gladrow 2001). It is interesting to note that the salt concentration suggested for Titan (1.88 molal; Baland et al. 2014) falls within this range of allowed maximum concentration.

The criterion of gravitational stability of the ocean may prove to be even more restrictive. One ought to consider the proliferation of clathrate-forming chemical species, both dissolved in the ocean and embedded within its bottom surface ice layer (Levi et al. 2014, 2017). Thus, the ocean's gravitational stability requires the consideration of the densities of clathrates as well. We have compared between the densities of aqueous NaCl solutions (Lvov & Wood 1990) and of the CO₂ structure I clathrate hydrate (Levi et al. 2017). We find that at 273.15 K and 0.8 GPa, their densities become equal for 3.14 wt% (0.56 molal) of NaCl. This is about the value of Earth's seawater. Therefore, the richer the ocean floor is in CO₂ clathrate hydrates, the less likely it is that a water-planet ocean was ever saltier than Earth's oceans.

The constraints placed by gravity on the upper bound values of the ocean's salinity are not very limiting. In addition, whether these upper bound values can be reached depends on the outward transport of solutes along the evolutionary path of the ocean exoplanet. Such an exploration is beyond the scope of this work. However, it may be that the melting and solidification cycle discussed in the previous section can further constrain the ocean's salinity. This requires an examination of the high-pressure experimental data available for saltwater solutions.

Frank et al. (2006) studied the H₂O–NaCl system up to 30 GPa at room temperature using a diamond anvil cell. A pressure of 30 GPa would reflect the conditions prevailing at the bottom of an ice mantle of a 2 M_⊕ planet with an approximately 15% ice mass fraction (see the water-planet internal structure tables in Levi et al. 2014). These authors compressed two NaCl aqueous solutions: one of 1.6 mol% (0.9 molal) and the other of 3.2 mol% (1.9 molal) of NaCl. In the more dilute case, the diffraction lines could be indexed using ice VII alone. In the more solute-rich case, there were clear indications for the formation of halite. The authors concluded that, at room temperature, the many natural voids available within the ice VII structure could host Na⁺ and Cl⁻ ions with a solubility of 2.4 ± 0.8 mol% (1.4 ± 0.5 molal). Daniel & Manning (2008) give a brief report of their investigation of the H₂O–NaCl system, and their results do not contradict the high solubility found by Frank et al. (2006) and the analysis of this work.

Adams (1931) measured the volume of the liquid binary solution H₂O–NaCl for various compositions at room temperature. In his experiment, he investigated pressures up to the transition to ice VI. From the partial volumes he derived chemical potentials, and inferred the solubility of NaCl in liquid water in the pressure range of 1 bar to 1.6 GPa. At 0.3 GPa, he estimated the solubility to be 6.52 molal. Sawamura et al. (2007) yielded a solubility in liquid water of 6.61 molal of NaCl, at these same pressure and temperature conditions. Adams (1931) further reported that the solubility varies very little with pressure, and derived a value for the solubility of 6.3 molal at 1 GPa, i.e., at the bottom of a water-planet ocean.

When ionic solutes become embedded in the water-ice lattice, they ought to influence their surrounding water molecules. Based on their X-ray diffraction data and Raman spectroscopy of the O–H stretching frequency, Frank et al. (2006) suggested that water hydrogen atoms tend to align themselves closer toward the Cl⁻ ion and also experience partial ordering. The increased order, around the ionic solute, is manifested by the protons becoming equidistant from the oxygen atoms, thus resembling the structure of ice X. Therefore, for temperatures above room temperature, and in the stability regime of water ice VII, the protons might resist this ordering, with implications for the solubility of salt in high-pressure water-ice polymorphs. In order to test their findings, Frank et al. (2008) repeated their earlier experiments, though now for temperatures higher than room temperature. The authors prepared a 1.6 mol% (0.9 molal) aqueous solution of NaCl at room temperature, which again was compressed to different pressures up to about 30 GPa. As before, halite was not observed at room temperature. When the samples were heated to 500 K, halite diffraction lines began to appear and intensified with the increasing temperature. The appearance of halite beyond 500 K was reported to not be very sensitive to the pressure. The halite signature disappeared once more after the samples crossed the melting point for the given isobar, meaning the formed fluid was not saturated. The authors concluded that at 500 K the salt began to exsolve out of the ice. We refer the interested reader to Bove et al. (2015) for more details on the effects of the inclusion of salt on the structure of water ice.

From Figure 3, we see that for the case of the thin upper boundary layer, a temperature of 500 K is exceeded only within the lower boundary layer, separating the inner rock core from the ice mantle. For the case of the thick upper boundary layer, this temperature is exceeded at a pressure level of about 53 GPa along the ice mantle adiabat. These are likely upper bound values for the salt exsolution pressure level. As we have shown in Levi et al. (2014), an ice mantle enriched with filled ice would have a higher temperature than a pure-water ice mantle.

These experimental results indicate that during the solidification of ocean water, in the melting–solidification cycle explained above, salts in the form of hydrated ions from the ocean can become incorporated into the continuously forming ocean floor. In the process of the convective overturn of the ice mantle, thermal conditions in the downwelling arm of the convection cell are sufficient to promote the unmixing of the solid solution into pure ice VII and pure salt. In Figure 7, we plot the density of pure NaCl versus that for pure-water ice VII at room temperature and up to 30 GPa. It is clear that under these conditions, the NaCl salt grains are much denser than ice VII. Therefore, the salt grains formed following this exsolution process sediment out of the ice mantle and onto the rocky interior.

Water planets, as considered here, have a water mass fraction of a few tens of percent and pressures of tens of GPa at the bottom of their ice mantles. Temperatures at the rocky core and ice mantle boundary far exceed the conditions for the exsolution of salt (see Figure 3). Therefore, under such conditions, it is probable that salts from the rocky core cannot incorporate into the overlying ice mantle as dissolved ions, to be transported outward.

For the case of icy satellites, a mechanism that is often suggested for the outward transport of salt is pockets of brine (e.g., Choblet et al. 2017; Kalousová et al. 2018). For such small icy bodies, these pockets are positively buoyant, since their surrounding is a partially differentiated layer of ice and rock. Partial differentiation of ice from rock is not likely for a planet. Therefore, the outward migration of a pocket of brine from the rock–ice mantle boundary requires it to be less dense than the surrounding ice. From Goncharov et al. (2009), we know that the density of water ice at the range of 2–60 GPa is at most a factor of 1.15 higher than that of the pure fluid. The equation of state for high-pressure brine as a function of composition is not known. Here we approximate the contribution of dissolved NaCl to the solution density using low-pressure data from Lvov & Wood (1990). They find that adding more than about 20 wt% (i.e., an abundance of about 0.07) of NaCl to liquid water can increase the solution density by more than a factor of 1.15. Therefore, we tentatively assume that for a brine pocket to be positively buoyant, it should not contain a molar abundance of dissolved NaCl in excess of about 0.07.

A pocket of brine forms if local conditions coincide with melting conditions. The thermal profiles in the ice mantles of water planets, as described in Figure 3, show that the difference between the adiabatic temperature and the ice melting temperature increases substantially with increasing pressure. Therefore, in order to produce local melting conditions, added solutes ought to lower the melting temperature by hundreds of degrees. To test whether this is possible, we derive from Goncharov et al. (2009) the entropy of fusion for ice VII, for the pure-water system, from 2 GPa to 60 GPa. The use of the pure-water system data is explained in Appendix A.

For pressures higher than 45 GPa,

$$\begin{eqnarray}{S}_{f} & = & 3\times {10}^{-6}{P}^{3}-0.00046{P}^{2}+0.024P\\ & & -\,0.4[\mathrm{kJ}\,{{\rm{K}}}^{-1}\,{\mathrm{mol}}^{-1}],\end{eqnarray} \tag{ 13 }$$

and for pressures lower than 45 GPa

$$\begin{eqnarray}{S}_{f} & = & 2.83\times {10}^{-5}{P}^{2}+0.00082P+0.026\\ & & \times [\mathrm{kJ}\,{{\rm{K}}}^{-1}\,{\mathrm{mol}}^{-1}],\end{eqnarray} \tag{ 14 }$$

where P is in GPa. We use the following equation to estimate the melting-point depression due to the addition of salt (see Appendix A for the derivation):

$$\begin{eqnarray}-{S}_{f}{\rm{\Delta }}T & + & \tfrac{1}{2}\left(\displaystyle \frac{\partial {S}_{f}}{\partial T}\right){({\rm{\Delta }}T)}^{2}=k({T}_{0}-{\rm{\Delta }}T)\\ & & \times \mathrm{ln}({\gamma }_{{{\rm{H}}}_{2}{\rm{O}}}(1-{X}_{\mathrm{NaCl}}),)\end{eqnarray} \tag{ 15 }$$

where ΔT is the melting-point depression, k is Boltzmann's constant, T₀ is the undepressed melting temperature, ${\gamma }_{{{\rm{H}}}_{2}{\rm{O}}}$ is the activity coefficient of water, and X_NaCl is the abundance of salt dissolved into the fluid solution.

In Figure 8, we plot various iso-melting temperature depression curves. From the figure, we see that lowering the melting temperature by 400 K requires

$$\begin{eqnarray}&&{\gamma }_{{{\rm{H}}}_{2}{\rm{O}}}(1-{X}_{\mathrm{NaCl}})\lt 0.1,\end{eqnarray} \tag{ 16 }$$

which yields for our adopted maximum salt concentration of X_NaCl = 0.07 the following restriction on the activity coefficient,

$$\begin{eqnarray}&&{\gamma }_{{{\rm{H}}}_{2}{\rm{O}}}\lt 0.11.\end{eqnarray} \tag{ 17 }$$

Journaux et al. (2013) report on the melting-point depression of ice VII when formed from a solution with a 19wt% (or abundance of 0.067) of NaCl. We infer from their melting-point depression data that the activity coefficient of water varies from 0.76 to 0.85, when increasing the pressure from 3 GPa to 6 GPa. Frank et al. (2008) found for a 5 wt% NaCl–water solution a melting-point depression of 40 K at 25 GPa. We find that this corresponds to an activity coefficient for water of 0.74. We note that the latter activity coefficient will only allow for a maximum melting-point depression of the order of tens of Kelvin, for X_NaCl = 0.07. Therefore, both conditions, melting and the positive buoyancy of the brine pocket, are not met simultaneously.

Zoom In Zoom Out Reset image size

If salt is partitioned out of the ocean in the downwelling arm of the ice mantle overturn, while the upwelling arm is prevented from replenishing the ocean, a pump is established. This pump would tend to desalinate the ocean over time. The pump is illustrated in Figure 9.

Zoom In Zoom Out Reset image size

We turn to estimate the rate of oceanic desalination, as enforced by this pump. In steady state, the mass of the ocean, M_oc, is constant. Therefore, the rate of melting at spreading centers equals the rate at which ocean water solidifies (see ${\dot{M}}_{\mathrm{liq}}$ in Section 2). We separate the solidification into a fraction χ forming a new ocean floor and a fraction 1 − χ representing the macroscopic entrapment of ocean water in regions of plate bending. Hydrated ions of salt dissolved in the ocean are incorporated into the newly forming ocean floor with some fractionation factor, f₁. Direct entrapment where the plates bend would fractionate ions into the ocean floor with a fractionation factor of f₂.

The average salt concentration in the ocean is m_salt (in molal units). In Levi et al. (2017), we have shown that winds and tides are insufficient for efficiently mixing the deep oceans we study here. Therefore, the bottom of the ocean may be enriched in salt relative to the average ocean salinity. In Appendix B, we estimate this factor of enrichment to be at least 2.5. If an oceanic mass element near the ocean floor, dM_oc, solidifies, the number of moles of salt partitioned out of the ocean is

$$\begin{eqnarray}&&{{dn}}_{\mathrm{salt}}\,\approx \,2.5{{dM}}_{\mathrm{oc}}\chi {f}_{1}{m}_{\mathrm{salt}}+2.5{{dM}}_{\mathrm{oc}}(1-\chi ){f}_{2}{m}_{\mathrm{salt}}.\end{eqnarray} \tag{ 18 }$$

If n_salt is the total number of moles of salt dissolved in the ocean, we have

$$\begin{eqnarray}&&{n}_{\mathrm{salt}}={M}_{\mathrm{oc}}{m}_{\mathrm{salt}}.\end{eqnarray} \tag{ 19 }$$

Differentiating with respect to time and using the steady-state assumption for the ocean mass yield

$$\begin{eqnarray}&&\displaystyle \frac{d}{{dt}}{m}_{\mathrm{salt}}\approx -2.5\displaystyle \frac{{\dot{M}}_{\mathrm{liq}}}{{M}_{\mathrm{oc}}}{m}_{\mathrm{salt}}[\chi {f}_{1}+(1-\chi ){f}_{2}].\end{eqnarray} \tag{ 20 }$$

Taking t = 0 as the time when the ocean may have attained maximum salinity, as mandated by the criterion of gravitational stability, we can solve the last equation to obtain

$$\begin{eqnarray}&&\mathrm{ln}\left(\displaystyle \frac{{m}_{\mathrm{salt}}[\mathrm{mol}\,{\mathrm{kg}}^{-1}]}{0.56}\right)=-\displaystyle \frac{t}{{\tau }_{d}},\end{eqnarray} \tag{ 21 }$$

where the ocean desalination timescale is

$$\begin{eqnarray}&&{\tau }_{d}\approx \displaystyle \frac{1}{2.5[\chi {f}_{1}+(1-\chi ){f}_{2}]}\displaystyle \frac{{M}_{\mathrm{oc}}}{{\dot{M}}_{\mathrm{liq}}}.\end{eqnarray} \tag{ 22 }$$

The region where the upper boundary layer (see Figure 5) folds experiences both an enhanced compression and tension. The strain in the fold can be written as

$$\begin{eqnarray}&&\epsilon =\displaystyle \frac{y}{{R}_{c}},\end{eqnarray} \tag{ 23 }$$

where y is the normal distance from the neutral axis and R_c is the radius of curvature of the fold. Assuming that the neutral axis lies along the center of the boundary layer, the maximal value for y is

$$\begin{eqnarray}&&{y}_{\max }\,=\,\displaystyle \frac{{\delta }_{\mathrm{om}}+{\delta }_{\mathrm{int}}}{2}.\end{eqnarray} \tag{ 24 }$$

Thus, the differential stress in the fold is (Turcotte & Schubert 2002)

$$\begin{eqnarray}&&{\sigma }_{\mathrm{fold}}=\displaystyle \frac{Y}{1-{\nu }^{2}}\displaystyle \frac{y}{{R}_{c}}.\end{eqnarray} \tag{ 25 }$$

Here, Y is the Young's modulus of water ice VI. Single-crystal (Shimizu et al. 1996; Tulk et al. 1996) and polycrystalline (Shaw 1986) experiments assign the latter a value of about 18 GPa. For Poisson's ratio, ν, we adopt a value of 0.32 (Tulk et al. 1996).

The functional form for the radius of curvature, R_c, for various flow regimes is still being debated. For the purpose of studying Earth's internal dynamics, Crowley & O'Connell (2012) adopted a value of R_c/b = 0.1, where b is the depth of the convection cell. Rose & Korenaga (2011) and Korenaga (2010) reported that the value of R_c/b decreases when the Rayleigh number considered increases. The total mantle Rayleigh number for our case of study is larger than the value for Earth's mantle (Crowley & O'Connell 2012). Therefore, by adopting a value of 0.1, we may be overestimating the length scale of the radius of curvature. This results in a lower estimate for the stresses in the fold.

For the values given above, we find that the maximal stress in the fold is 10.3 GPa and 1.6 GPa for the thick and thin upper boundary solutions, respectively. On Earth, the intense differential stress at plate bending, during subduction, promotes both brittle failure and plastic creep (Billen & Gurnis 2005). This is often manifested as earthquakes, capable of cutting through a substantial part of Earth's ≈100 km lithosphere (Kikuchi & Kanamori 1995). If a similar activity occurs for the case we are studying, then ocean water my be forced into the downwelling oceanic floor. Consequently, the injected ocean water will follow the liquidus (see Figure 6) while increasing its density. It is probable that such a direct and macroscopic injection of ocean water would result in a fractionation factor of f₂ = 1, regardless of the identity of the hydrated ionic species.

However, away from regions of plate bending, the fractionation of ions from the ocean into the continuously forming ocean floor is microscopic in nature. Therefore, the fractionation factor f₁ ought to depend on the rate of ocean water solidification. If the solidification front advances fast enough, kinetics determines its value. As we have discussed above, the criterion for the gravitational stability of the ocean gives a maximum salinity of about 0.56 molal. This is less than the solubility of NaCl in high-pressure ice. Therefore, if solidification is fast, the ions dissolved in the liquid are easily incorporated as interstitials within the high-pressure ice without supersaturating it, and f₁ = 1. If, on the other hand, the solidification front advances slowly, thermodynamics controls the value of the fractionation factor. In this case, equilibrium is attained after a fraction 1 − f₁ of the dissolved ions diffuse out of the solidifying mass and remain in the ocean. For this case, we estimate the fractionation factor by taking the ratio between the solubility of Na⁺ and Cl⁻ in the phase composed of the ocean floor and the overlying pressurized liquid water.

The solidification front advances in the ocean with a velocity

$$\begin{eqnarray}&&{v}_{{sf}}\sim \displaystyle \frac{{\dot{M}}_{\mathrm{liq}}}{{\rho }_{\mathrm{oc}}4\pi {R}_{p}^{2}}\sim 1\,[\mathring{\rm A} \,{s}^{-1}],\end{eqnarray} \tag{ 26 }$$

where ρ_oc is the density of the ocean. This is a lower bound value for the solidification front velocity because we assume the solidification takes place on the entire ocean floor. The diffusion timescale of ions in water is

$$\begin{eqnarray}&&{t}_{\mathrm{dif}}\sim \displaystyle \frac{{l}^{2}}{{D}_{\mathrm{ion}}}.\end{eqnarray} \tag{ 27 }$$

The diffusion coefficient of NaCl ions in water is 10⁻⁵ cm² s⁻¹ (Yuan-Hui & Gregory 1974). Advancing a distance l = 100 Å takes 10⁻⁷ s, whereas it takes 100 s for the solidification front to cover the same distance. Therefore, the kinetics of solidification is slow enough to let f₁ attain its equilibrium value.

From the experimental data given above, we obtain for the fractionation factor a value of f₁ ≈ 0.22 ± 0.08. This value is derived by adopting the experimental solubility for NaCl in ice VII. The solubility is lower for the case of ice VI. The latter is the more appropriate water-ice structure for the case we are solving for, in which the ocean floor temperature is low. Journaux et al. (2017) report on the solubility of RbI in ice VI. They further report on the melting-point depression in the H₂O–RbI binary system. They find that the melting-point depression is similar to that in the H₂O–NaCl system. We derived the activity coefficients for NaCl and RbI in an aqueous solution, which can be translated to the osmotic coefficient for the solvent (water). The latter is easily related to the activity of water in the binary solution (Blandamer et al. 2005). We find that the activity of water differs by less than 1% between the two salts. Therefore, the similar melting-point depression can also be translated to a similar solubility. At about 10 K below the melting temperature, Journaux et al. (2017) find a fractionation factor f₁ ≈ 0.005. For 2–3 K below the melting temperature, they give an upper bound value for the fractionation factor of about 0.3. Since the latter is an upper bound value rather than a measurement of the solubility, we shall adopt here their reported value of 0.005 ± 0.002. We refer the interested reader to Journaux et al. (2017) for a discussion on the difficulty in measuring the fractionation factor close to melting conditions, as in the on-melt layer (see the discussion of the on-melt layer in Section 2).

In Levi et al. (2014), we have derived planetary radii for a 2 M_⊕ planet, composed of various ice mass fractions. For a 50% ice mass fraction, we calculated a planetary radius of R_p = 9480 km. For an ocean depth of 80 km, the ocean mass is M_oc ≈ 9 × 10²² kg. For the melting rates tabulated in the previous section (adopting the average value between the thin and thick boundary layer solutions), the desalination timescale is

$$\begin{eqnarray}&&{\tau }_{d}\sim \displaystyle \frac{4.3}{\chi ({f}_{1}-1)+1}\quad [\mathrm{Myr}].\end{eqnarray} \tag{ 28 }$$

If there is no macroscopic entrapment of oceanic water at regions of plate bending, then χ = 1. For a planet with an ocean floor temperature that can stabilize ice VII, we find τ_d ≈ 20 Myr. For an ocean floor composed of ice VI, we have τ_d ≈ 860 Myr. If macroscopic entrapment can account for 1% of the reprocessing of ocean water (i.e., χ = 0.99), then the latter reduces to τ_d ≈ 290 Myr. The median age of transiting planets' host stars is ∼5 Gyr (Bonfanti et al. 2016). Therefore, the desalination mechanism is efficient in removing salt out of the ocean, with fundamental implications for the habitability of our studied worlds.

In Figure 10, we plot the oceanic salinity versus time. After 5 Gyr, the maximum concentration is of the order of millimolal. Under these circumstances, we conclude that ocean worlds, which are likely rich in volatiles, have oceans that are probably poor in salts.

Zoom In Zoom Out Reset image size

The equilibrium relations between the carbonate species are (Zeebe & Wolf-Gladrow 2001)

$$\begin{eqnarray}&&{\mathrm{CO}}_{2}+{{\rm{H}}}_{2}{\rm{O}}\rightleftarrows {\mathrm{HCO}}_{3}^{-}+{{\rm{H}}}^{+},\end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray}&&{\mathrm{HCO}}_{3}^{-}\rightleftarrows {\mathrm{CO}}_{3}^{2-}+{{\rm{H}}}^{+}.\end{eqnarray} \tag{ 30 }$$

Here, CO₂ represents freely dissolved carbon dioxide plus a trace amount of true carbonic acid, H₂CO₃. For our pressure and temperature range of interest, H⁺ represents a hydrated proton, which is part of a larger complex (an unbonded proton below 20 GPa has a very short lifetime; Goncharov et al. 2009). Two such commonly discussed complexes are the Zundel (${{\rm{H}}}_{5}{{\rm{O}}}_{2}^{+}$) structure, where the proton is shared by a pair of neighboring water molecules, and the Eigen (${{\rm{H}}}_{9}{{\rm{O}}}_{4}^{+}$) structure, for which a hydronium (H₃O⁺) core is bonded to three water molecules (Markovitch & Agmon 2007). Marx et al. (1999) have shown that many complexes exist, between which the hydrated proton transfers. Therefore, the Zundel and Eigen structures should only be considered as limiting cases.

In equilibrium, the relations between the concentrations of the different carbonate species are

$$\begin{eqnarray}&&{K}_{1}(P,T)=\displaystyle \frac{{a}_{{\mathrm{HCO}}_{3}^{-}}{a}_{{{\rm{H}}}^{+}}}{{a}_{{\mathrm{CO}}_{2}}{a}_{{{\rm{H}}}_{2}{\rm{O}}}},\end{eqnarray} \tag{ 31 }$$

$$\begin{eqnarray}&&{K}_{2}(P,T)=\displaystyle \frac{{a}_{{\mathrm{CO}}_{3}^{2-}}{a}_{{{\rm{H}}}^{+}}}{{a}_{{\mathrm{HCO}}_{3}^{-}}},\end{eqnarray} \tag{ 32 }$$

where K₁ and K₂ are the first and second dissociation (i.e., ionization) constants of carbonic acid, respectively. a_i is the activity of the ith chemical species. The activity of a chemical species is its concentration (m_i, using here the molality scale) times its activity coefficient (γ_i). In the Debye–Hückel theory, which usually serves as the theoretical foundation for describing solutions containing ionic solutes, a natural parameter describing the solution is its ionic strength (Davidson 2003),

$$\begin{eqnarray}&&I\equiv \tfrac{1}{2}\displaystyle \sum _{i}{m}_{i}{z}_{i}^{2},\end{eqnarray} \tag{ 33 }$$

where z_i is the valence of the ith ionic species.

In order to estimate the activity of an ionic solute, one requires an estimation of its activity coefficient. The activity coefficient is the fraction of the concentration of a component that is readily active in the solution. As long as the concentration of dissolved ions is not too high, the main reason for an activity coefficient correction (i.e., correction for solution non-ideality) comes from electric shielding between unlike ions. The shielding effectively reduces the long-range interaction rendering some portion of the ions less "active". For electric shielding, the most important aspect is the concentration of ions of a given valence rather than their elemental composition. Therefore, electric shielding can approximately be described with the help of the ionic strength parameter alone. In the Debye–Hückel theory, only long-range interactions between the different ions in the solution are taken into account. This implies that the Debye–Hückel theory is a good approximation for rather dilute solutions. In this case, the activity coefficient is theoretically given by the following form (see p. 507 in Davidson 2003):

$$\begin{eqnarray}&&{\mathrm{log}}_{10}{\gamma }_{i}=-1.8246\times {10}^{6}{z}_{i}^{2}\displaystyle \frac{{(\epsilon T)}^{-3/2}{({\rho }_{s}I)}^{1/2}}{1+50.29d{(\epsilon T)}^{-1/2}{({\rho }_{s}I)}^{1/2}},\end{eqnarray} \tag{ 34 }$$

where and ρ_s are the dielectric constant and bulk density of the solvent, respectively, T is the temperature, and d is the ionic diameter in angstroms. Since only long-range interactions are considered in the Debye–Hückel theory, it is not applicable for solutions with high ionic strengths. It is found to be a good approximation as long as the ionic strength is less than 10^−2.3 molal (Zeebe & Wolf-Gladrow 2001).

Davies (1938) suggested that a very simple empirical correction should be added to the Debye–Hückel activity coefficient, yielding what is known as the Davies activity coefficient equation:

$$\begin{eqnarray}&&{\mathrm{log}}_{10}{\gamma }_{i}=-1.8246\times {10}^{6}{(\epsilon T)}^{-3/2}\,{z}_{i}^{2}\left(\displaystyle \frac{\sqrt{I}}{1+\sqrt{I}}-0.2I\right).\end{eqnarray} \tag{ 35 }$$

Comparing with experimentally deduced equilibrium constants of various salt solutions in water, Davies (1938) argued that Equation (35) has an accuracy of about 2% for ionic strengths up to 0.1 molal. In practice, the Davies activity coefficient is considered to be a good approximation for ionic strengths up to 0.5 molal (Zeebe & Wolf-Gladrow 2001). It is important to note that in the Davis equation the ionic diameter, d, is omitted and so ions of the same valence have equal activity coefficients.

Özsoy et al. (2008) measured the ionic composition of several rivers discharging into the Cilician basin. From their measurements at stations close to the rivers' sources, we can derive an average ionic strength of 0.005 molal. Measurements therefore show that for freshwater on Earth, the above equations for the activity coefficient are applicable.

As one considers solutions of higher ionic concentrations, two phenomena become increasingly important: ion pairing and ion complex formation. Ion pairing is when two ions of opposite charge are close enough to experience a strong coulombic interaction yet create no covalent bonds. An ion pair can be two ions that are in contact and share a hydration shell, or have a single solvent molecule in between them, or have separate hydration shells (Pluhařová et al. 2013). On the other hand, an ion complex forms when ions begin to share electrons, thus forming new ionic species. For example, Earth's seawater has on average an ionic strength of 0.7 molal (see the standard mean chemical composition of seawater in Appendix A in Zeebe & Wolf-Gladrow 2001). Indeed, ion complexes between OH⁻ or ${\mathrm{CO}}_{3}^{2-}$ and divalent or trivalent metals are dominant in seawater in comparison to free metal ions (Millero et al. 2009). However, this is somewhat dependent on the metal under consideration. Therefore, in addition to their dependence on the ionic strength, these higher ionic concentration phenomena depend on the type of element considered as well. Thus, activity coefficients that would adequately describe such phenomena could not depend on the ionic strength alone.

Guggenheim (1935) argued that the Debye–Hückel theory can be made applicable to electrolyte solutions of higher ionic strengths. He suggested adding a correction accounting for short-range interactions between neighboring pairs of unlike-charged ions, i.e., a second virial correction. This second virial correction would consist of adjustable parameters to be obtained by fitting to experimental data and would be dependent on the types of ions composing the existing pairs in the solution. The approach suggested by Guggenheim gave accurate activity coefficients, within experimental error, for ionic strengths up to 0.1 molal (Guggenheim & Turgeon 1955). A model for the activity coefficients that is commonly used for describing seawater is that of Pitzer (1973). In this model, the constant coefficients of the second virial correction of Guggenheim (1935) are taken to be functions of the ionic strength of the solution, and higher order correction terms are added to account for short-range interactions between triplets of ions as well. The Pitzer model is extensively used when modeling electrolyte solutions where ionic concentrations exceed 1 molal. The parameters for the Pitzer model are usually obtained by a fit to solubility experimental data (see Li & Duan 2011 and references therein). In this work, we solve for the speciation of the carbonate system in the oceans of water-rich planets and derive the resulting ionic strength. We then make sure our choice for the activity coefficient model is consistent with our results.

Before leaving the subject of the activity coefficient, one may question the applicability of the Debye–Hückel theory to highly pressurized solutions. Manning (2013) estimated theoretically the solubility of corundum in KOH and aqueous solutions at 1 GPa for various activity coefficient models, including the Davies equation. By comparing his predicted solubility to experimental data at this high pressure, the author was able to show that the extended Debye–Hückel activity theory, and the Davies equation specifically, was adequate for modeling high-pressure solubility. Manning et al. (2013) repeated the same procedure, and compared his predicted solubility for calcite in water with experimental data at 1 GPa, where, once again, the Davies equation proved adequate. We note that 1 GPa approximates the upper limit on the pressure expected at the bottom of our water-planet oceans.

In the ocean, DIC resides in the species CO₂, ${\mathrm{HCO}}_{3}^{-}$, and ${\mathrm{CO}}_{3}^{2-}$. The total concentration of the DIC is therefore (Zeebe & Wolf-Gladrow 2001)

$$\begin{eqnarray}&&\mathrm{DIC}\equiv [{\mathrm{CO}}_{2}]+[{\mathrm{HCO}}_{3}^{-}]+[{\mathrm{CO}}_{3}^{2-}],\end{eqnarray} \tag{ 36 }$$

where the square brackets denote concentration throughout this work.

Another governing equation that must be addressed is the charge neutrality requirement for the ocean. A variety of ions differing in their valences, concentrations, and elemental compositions are dissolved in Earth's oceans. This also ought to be the case for any ocean with a rocky floor. These introduce a higher level of complexity, because they can react with the carbonate system, for example, in the formation and solubility of calcite, or when ion-rich silicates form clays during silicate weathering.

When studying Earth's oceans, one can partially overcome this enormous complexity by introducing the concept of carbonate alkalinity, defined as (Zeebe & Wolf-Gladrow 2001)

$$\begin{eqnarray}&&\mathrm{CA}\equiv [{\mathrm{HCO}}_{3}^{-}]+2[{\mathrm{CO}}_{3}^{2-}].\end{eqnarray} \tag{ 37 }$$

Equations (31), (32), (36), and (37) represent four equations governing over six variables: $[{\mathrm{CO}}_{2}]$, $[{\mathrm{HCO}}_{3}^{-}]$, $[{\mathrm{CO}}_{3}^{2-}]$, $[{{\rm{H}}}^{+}]$, DIC, and CA. If one can measure any two of these variables, the system can be solved for (Millero et al. 2002). This is indeed the practice for studying Earth's oceans, where the DIC and the ocean's pH can be measured directly with relative ease (Zeebe & Wolf-Gladrow 2001). No such short cuts are available for exoplanets whose oceans float on a rocky floor.

For the case of Earth, most of the carbon is deposited in rocks (Williams & Follows 2011), therefore putting an upper limit on the allowed value for the DIC (as expected for a closed system). Considering the pH for Earth's ocean, one can solve for the carbonate system speciation and show that bicarbonate and carbonate are the dominant species in the DIC. In other words, the ocean makes use of the acidic nature of freely dissolved CO₂ to comply with the requirement of charge neutrality by transforming it to ionic species.

It is interesting to note that, because the DIC in Earth's oceans has an upper limit, our oceans can largely remain undersaturated with respect to freely dissolved CO₂. The consequence of this is the instability of SI CO₂ clathrate hydrate in the oceans on Earth. This may be a general conclusion for water worlds having a global ocean with a rocky bottom surface, if these oceans represent a closed system with respect to CO₂ (i.e., limited DIC) as well. However, in this study, we focus on ocean planets whose ocean floor is composed of high-pressure ice. Such planets can keep their oceans saturated in CO₂ and stabilize the appropriate clathrate hydrate (Levi et al. 2017).

As discussed in the previous section, the oceans we study here are likely poor in salts. Therefore, the charge neutrality of the ocean can be approximated as

$$\begin{eqnarray}&&0=[{\mathrm{HCO}}_{3}^{-}]+2[{\mathrm{CO}}_{3}^{2-}]+[{\mathrm{OH}}^{-}]-[{{\rm{H}}}^{+}].\end{eqnarray} \tag{ 38 }$$

Finally, when in equilibrium, the dissociation products of water obey

$$\begin{eqnarray}&&{K}_{w}(P,T)={a}_{{{\rm{H}}}^{+}}{a}_{{\mathrm{OH}}^{-}}.\end{eqnarray} \tag{ 39 }$$

Equations (31), (32), (36), (38), and (39) govern the speciation of the carbonate system in our studied water planets. The six variables of the system are $[{\mathrm{CO}}_{2}]$, $[{\mathrm{HCO}}_{3}^{-}]$, $[{\mathrm{CO}}_{3}^{2-}]$, $[{{\rm{H}}}^{+}]$, $[{\mathrm{OH}}^{-}]$, and the DIC. Therefore, if we know one of the variables, we can solve for the entire system. Because our studied oceans are saturated in CO₂, as mandated from being an open system, we can solve the system of governing equations as a function of $[{\mathrm{CO}}_{2}]$.

In this subsection, we discuss the first, K₁, and second, K₂, dissociation constants of carbonic acid and the dissociation constant for water, K_w. Experiments measuring K₁ and K₂ do not cover our full parameter space of interest. For example, Dedong & Zhenhao (2007) represented K₁ and K₂ dependence on the temperature and pressure using a polynomial equation. The polynomial coefficients, for the part representing the temperature dependence, were obtained by a fit to K₁ and K₂ experimental data. For the dissociation constants' dependence on pressure, the authors used experimental partial molar volumes and compressibility data. This procedure is reported to reproduce the available experimental data with good accuracy in the temperature range from 0 °C to 250 °C. This temperature range is adequate for the purpose of modeling water-rich planets. However, the applicability of the formulations of Dedong & Zhenhao (2007) have an upper limit on the pressure of 0.1 GPa. This upper limit is sufficient for modeling the carbonate system at the bottom of Earth's ocean. However, it is an order of magnitude lower than the pressure prevailing at the ocean floor for our studied water planets.

Experiments determining the equilibrium ionization constants at low temperatures (approximately 300 K) and for pressures exceeding 0.1 GPa are very few in number (see review of experiments in Dedong & Zhenhao 2007). The notable exceptions are the experiments of Ellis (1959), reaching up to 0.3 GPa, and Read (1975), reaching as high as 0.2 GPa, and both are only for K₁.

High-pressure (≥1 GPa) carbonic acid ionization constants are experimentally available. However, the parameter space probed by experiments is mostly guided by the conditions within planet Earth. At these relatively high pressures and within the Earth, one examines the deep crust or upper mantle. Because the thermal conductive profile rapidly increases the temperature with depth in Earth's crust, then at 1 GPa, temperatures fall in the range of 700–1300 K, depending on the type of tectonic environment (Jones et al. 1983). Therefore, efforts for evaluating K₁ and K₂ at high pressure have so far been focused on temperatures exceeding, by at least a factor of two, those that are of interest to us.

A better understanding of the carbon cycle on Earth requires quantifying the dissolution of calcite and aragonite in aqueous fluids that are released from the tectonic slab during its subduction. We note that Facq et al. (2014) hydrostatically compressed a sample of calcite immersed in double-distilled water using a heated diamond anvil cell at temperatures appropriate for describing cold subduction zones (573–673 K) to deep-crust pressures as high as 8 GPa.

In this subsection, we will make use of two very common methods for interpolating between and extrapolating over the experimentally determined ionization constants: the solvent density method and the revised Helgeson–Kirkham–Flowers (HKF) equation of state for electrolyte solutions (see review of methods in Dolejs 2013). Our approach will be to estimate the free parameters for each of the two methods by fitting to low-pressure (∼0.1 GPa) experimental data for the dissociation constants. Then, we use models to extrapolate to high pressure. We then test the extrapolations by their ability to predict the high-pressure experimental dissociation constants given in Facq et al. (2014). The perils of extrapolating beyond the experimental data without a sound theoretical model will also be shown and discussed in this subsection.

First, we examine the solvent density method. We note that it has several variants (see discussion in Dolejs 2013). The physical reasoning at the basis of this method is that the solvent (water) can be described as a continuous dielectric medium, that the free energy change in solvation of an ion in the solvent can be described to a good approximation with the aid of the solvent's dielectric constant, and that this dielectric constant is strongly correlated with the solvent's density. As an example, we adopt the formulation given in Dolejs & Manning (2010):

$$\begin{eqnarray}\mathrm{ln}K & = & \displaystyle \frac{{A}_{1}}{T}+{A}_{2}+{A}_{3}\mathrm{ln}T+{A}_{4}T+{A}_{5}\mathrm{ln}{\rho }_{w}\\ & \equiv & m(T)+{A}_{5}\mathrm{ln}{\rho }_{w},\end{eqnarray} \tag{ 40 }$$

where K is the ionization constant, ρ_w is the density of water, and T is the temperature. The form for m(T) comes from assuming that the neutral molecule breakdown followed by the hydration of the products can be modeled using the free energy derived from integration over heat capacity terms, which are assumed to be linear functions of the temperature. The term $\mathrm{ln}{\rho }_{w}$ is assumed to encapsulate the effects of electrostriction. Plotting the experimental isothermal data from Read (1975) and Ellis (1959) as $\mathrm{ln}{K}_{1}$ versus $\mathrm{ln}{\rho }_{w}$, we indeed find it to be linear to a good approximation. We also find that A₅ must be varied with the temperature, in order to fit the data. Adopting the parametrization for the dielectric constant of water from Sverjensky et al. (2014), we assign to A₅ the following form:

$$\begin{eqnarray}&&{A}_{5}(T)={b}_{1}T+{b}_{2}\sqrt{T}+{b}_{3}.\end{eqnarray} \tag{ 41 }$$

For each isothermal data set reported in Read (1975) and Ellis (1959), we obtain a value for m and A₅ using a linear fit to the data, thus obtaining m(T) and A₅(T). Then, using the least-squares method, we obtain the optimal values for the parameters A₁ through A₄ and b₁ through b₃. The values we obtain for the coefficients A₁ to A₄ match the discrete experimentally fitted values for m(T) with an absolute average deviation of 0.16%. The values we find for b₁ to b₃ produce the experimentally fitted values for A₅(T) with an absolute average deviation of 3.2%.

The formulation for K₁ given in Equations (40) and (41), with the fitted parameters, reproduces the data sets of Read (1975) and Ellis (1959) with an absolute average deviation of a few percent, reaching as high as 13% for the case of the 200.5 °C isotherm from Read 1975 (see the right panel in Figure 11). We have also tried to simultaneously reproduce the experimental values of K₁ along the water saturation vapor pressure curve from Stefánsson et al. (2013). Our solvent density model fit can reproduce the latter experimental data set with an absolute average deviation of 15.6% (see the left panel in Figure 11).

Zoom In Zoom Out Reset image size

Our results show that the solvent density model can provide a simple method for interpolating between data points to an accuracy of a few percent. Considering its simplicity, it is a valuable tool. In addition, we test this model as an extrapolation tool by comparing its predictions to available experimental data. In Table 2, we list the high-temperature and high-pressure experimental data from Facq et al. (2014) alongside extrapolated K₁ values using the solvent density model, for the same temperature and pressure conditions. At a pressure of 1 GPa, the extrapolated values, using the solvent density model, are good to a factor of a few. However, at higher pressures, there are deviations of orders of magnitude. Therefore, using this theoretical model to extrapolate on the experimental data would produce erroneous results.

A different theory is the HKF model. A neutral molecule dissolved in a medium can react with its surrounding molecules and form ions. For example, a CO₂ molecule in an aqueous fluid can react with water molecules and form bicarbonate. The result is an increased electrostatic field in the locality surrounding the ionic product. The increased electrostatic field deforms the surrounding dielectric solvent, a phenomenon termed "electrostriction" (Stratton 1941). Electrostriction is actually a sum of many interdependent phenomena that involve ion–solvent interaction and the structural deformation of the solvent shells surrounding the ion. Helgeson & Kirkham (1976) proposed a simplified way to describe the thermodynamic properties of electrolyte aqueous solutions. Their simplified model for describing electrostriction is the foundation of the widely used revised-HKF model.

Helgeson & Kirkham (1976) argued that every thermodynamic quantity of dissolved ions in a dielectric solvent can be broken down into a sum of three terms. Those are an ionic intrinsic term, a solvent structural term, and a direct ion–solvent interaction term. The last two terms combined represent the net effect of electrostriction. Taking the standard partial molal volume as an example of such a thermodynamic quantity, its value can be written as

$$\begin{eqnarray}&&{\bar{V}}^{\circ }={\bar{V}}_{i}^{\circ }+{\rm{\Delta }}{\bar{V}}_{e}^{\circ }={\bar{V}}_{i}^{\circ }+{\rm{\Delta }}{\bar{V}}_{c}^{\circ }+{\rm{\Delta }}{\bar{V}}_{s}^{\circ }.\end{eqnarray} \tag{ 42 }$$

Here, ${\bar{V}}_{i}^{\circ }$ represents an intrinsic ionic volume independent of the nature of the solvent. ${\rm{\Delta }}{\bar{V}}_{e}^{\circ }$ is the volume change associated with electrostriction, which is further broken down to a volume change from structural collapse in the solvent surrounding the dissolved ion, ${\rm{\Delta }}{\bar{V}}_{c}^{\circ }$, and a volume change associated with direct ion–solvent interactions, also called volume change of solvation, ${\rm{\Delta }}{\bar{V}}_{s}^{\circ }$.

Consider a neutral conducting sphere of radius $\hat{b}$ that is immersed in a medium with a dielectric constant . The work required to charge it to an ionic valence z is (Davidson 2003)

$$\begin{eqnarray}&&W=\tfrac{1}{2}\displaystyle \frac{{z}^{2}{q}^{2}}{\epsilon \hat{b}},\end{eqnarray} \tag{ 43 }$$

where q is the magnitude of fundamental charge. The free energy change associated with transferring such a charged sphere from a vacuum to a dielectric medium of dielectric constant is therefore

$$\begin{eqnarray}&&{\rm{\Delta }}{\bar{G}}_{s}^{\circ }=\tfrac{1}{2}\displaystyle \frac{{z}^{2}{q}^{2}}{\hat{b}}\left(\displaystyle \frac{1}{\epsilon }-1\right)\equiv \omega \left(\displaystyle \frac{1}{\epsilon }-1\right),\end{eqnarray} \tag{ 44 }$$

where ω is the Born coefficient. This last function is also known as the Born equation and was used in Helgeson & Kirkham (1976) to quantify the standard partial molal Gibbs free energy of solvation. Within the framework of the HKF model, derivatives of this formulation produce, for every thermodynamic quantity, the solvation contribution to the net electrostriction effect. This is contrary to other theoretical models where the last term is considered to represent the effect of electrostriction in its entirety (e.g., Benson & Copeland 1963). Thus, the standard partial molal volume change of solvation, the standard partial molal compressibility of solvation, and the standard partial molal isobaric heat capacity of solvation have the following forms (Helgeson & Kirkham 1976; Helgeson et al. 1981):

$$\begin{eqnarray}&&{\rm{\Delta }}{\bar{V}}_{s}^{\circ }=-\displaystyle \frac{\omega }{\epsilon }{\left(\displaystyle \frac{\partial \mathrm{ln}\epsilon }{\partial P}\right)}_{T},\end{eqnarray} \tag{ 45 }$$

$$\begin{eqnarray}&&{\bar{\kappa }}_{s}^{\circ }\equiv -{\left(\displaystyle \frac{\partial {\rm{\Delta }}{\bar{V}}_{s}^{\circ }}{\partial P}\right)}_{T}=\displaystyle \frac{\omega }{\epsilon }\left[{\left(\displaystyle \frac{{\partial }^{2}\mathrm{ln}\epsilon }{\partial {P}^{2}}\right)}_{T}-{\left(\displaystyle \frac{\partial \mathrm{ln}\epsilon }{\partial P}\right)}_{T}^{2}\right],\end{eqnarray} \tag{ 46 }$$

$$\begin{eqnarray}&&{\bar{C}}_{p,s}^{\circ }=\displaystyle \frac{\omega T}{{\epsilon }^{2}}\left[{\left(\displaystyle \frac{{\partial }^{2}\epsilon }{\partial {T}^{2}}\right)}_{P}-\displaystyle \frac{2}{\epsilon }{\left(\displaystyle \frac{\partial \epsilon }{\partial T}\right)}_{P}^{2}\right].\end{eqnarray} \tag{ 47 }$$

Helgeson & Kirkham (1976) examined the intrinsic and solvent structural collapse effects by subtracting from experimental values for the ionic standard partial molal volumes the theoretical solvation volume change effect. Repeating this algorithm with experimental standard partial molal compressibility data, they have formulated their results in the following form:

$$\begin{eqnarray}&&{\bar{V}}_{i}^{\circ }+{\rm{\Delta }}{\bar{V}}_{c}^{\circ }={a}_{1}+{a}_{2}P+{a}_{3}\displaystyle \frac{T}{T-{\rm{\Theta }}}+{a}_{4}P\displaystyle \frac{T}{T-{\rm{\Theta }}}.\end{eqnarray} \tag{ 48 }$$

Helgeson et al. (1981) further subtracted from experimental data for the standard partial molal isobaric heat capacities of electrolytes, at water saturation vapor conditions, the theoretical contribution from solvation. They then found for the intrinsic and solvent structural collapse contributions, at 1 bar, the following:

$$\begin{eqnarray}&&{\bar{C}}_{p,i}^{\circ }+{\bar{C}}_{p,c}^{\circ }={c}_{1}+{c}_{2}\displaystyle \frac{1}{T-{\rm{\Theta }}}.\end{eqnarray} \tag{ 49 }$$

The nature of the asymptotic behavior at the temperature Θ is still under debate. Although the asymptotic behavior in the HKF model was arrived at by using a simple fit to experimental data, there seems to be a real physical basis for this behavior. Several of the thermodynamic properties of water seem to diverge with decreasing temperature in the supercooled liquid regime, with a singularity in the temperature falling in the range from 223 K to 228 K (Torre et al. 2004). The nature of this behavior could be due to polymorphism of liquid water, where the singularity in the temperature could be indicative of a transition between a high-density and a low-density liquid (Mallamace et al. 2014). Another debated possibility is that the singularity in the temperature represents a state of dynamic self-arrest. Relaxation times have an Arrhenius-type behavior with some activation energy. For the case of water, the activation energy increases drastically with decreasing temperature in the supercooled regime and seems to turn infinite at some finite temperature, Θ, resulting in dynamic arrest (Hecksher et al. 2008). In the HKF model, a constant value of 228 K is adopted for Θ.

Tanger & Helgeson (1988) introduced two corrections to the formulations given above. They have formulated what is known as the revised-HKF equations of state. The asymptotic behavior can have a different exponent for different thermodynamic properties. On the grounds of a better fit to heat capacity experimental data, Tanger & Helgeson (1988) replaced Equation (49) with the following:

$$\begin{eqnarray}&&{\bar{C}}_{p,i}^{\circ }+{\bar{C}}_{p,c}^{\circ }={c}_{1}+{c}_{2}{\left(\displaystyle \frac{1}{T-{\rm{\Theta }}}\right)}^{2}.\end{eqnarray} \tag{ 50 }$$

Equation (48) was found to predict molal volumes that approximate experimental data well only up to 0.1 GPa. Therefore, using experimental data for the molal volumes of NaCl and K₂SO₄ to 1 GPa, Tanger & Helgeson (1988) revised Equation (48) to the following form:

$$\begin{eqnarray}{\bar{V}}_{i}^{\circ }+{\rm{\Delta }}{\bar{V}}_{c}^{\circ } & = & {a}_{1}+{a}_{2}\displaystyle \frac{1}{{\rm{\Psi }}+P}+{a}_{3}\displaystyle \frac{1}{T-{\rm{\Theta }}}\\ & & +{a}_{4}\displaystyle \frac{1}{T-{\rm{\Theta }}}\displaystyle \frac{1}{{\rm{\Psi }}+P},\end{eqnarray} \tag{ 51 }$$

where Ψ is a constant depending on the type of solvent. Hence, the standard partial molal Gibbs free energy of formation of an aqueous species, either ion or electrolyte, can be derived, yielding the following form:

$$\begin{eqnarray*}{\rm{\Delta }}{\bar{G}}_{f}^{\circ } & = & {\rm{\Delta }}{\bar{G}}_{f,r}^{\circ }-{\bar{S}}_{r}^{\circ }(T-{T}_{r})-{c}_{1}\left[T\mathrm{ln}\displaystyle \frac{T}{{T}_{r}}-T+{T}_{r}\right]+{a}_{1}(P-{P}_{r})\\ & & +{a}_{2}\mathrm{ln}\left(\displaystyle \frac{{\rm{\Psi }}+P}{{\rm{\Psi }}+{P}_{r}}\right)-{c}_{2}\left[\left\{\displaystyle \frac{1}{T-{\rm{\Theta }}}-\displaystyle \frac{1}{{T}_{r}-{\rm{\Theta }}}\right\}\displaystyle \frac{{\rm{\Theta }}-T}{{\rm{\Theta }}}\right.\\ & & -\left.\displaystyle \frac{T}{{{\rm{\Theta }}}^{2}}\mathrm{ln}\left(\displaystyle \frac{{T}_{r}(T-{\rm{\Theta }})}{T({T}_{r}-{\rm{\Theta }})}\right)\right]+\displaystyle \frac{1}{T-{\rm{\Theta }}}[{a}_{3}(P-{P}_{r})\\ & & \left.+\,{a}_{4}\mathrm{ln}\left(\displaystyle \frac{{\rm{\Psi }}+P}{{\rm{\Psi }}+{P}_{r}}\right)\right]+\omega \left(\displaystyle \frac{1}{\epsilon }-1\right)-{\omega }_{r}\left(\displaystyle \frac{1}{{\epsilon }_{r}}-1\right)\\ & & +\,{\omega }_{r}\displaystyle \frac{1}{{\epsilon }_{r}^{2}}{\left(\displaystyle \frac{\partial \epsilon }{\partial T}\right)}_{P={P}_{r}}(T-{T}_{r}),\end{eqnarray*} \tag{ 52 }$$

where the sub-index, r, refers the parameter to its value at the reference conditions ${P}_{r}=1\,\mathrm{bar}$ and T_r = 298.15 K. With the aid of the last equation, one can derive the standard partial molal Gibbs free energy of reaction for the first and second ionization constants of carbonic acid. This Gibbs free energy is related to the equilibrium constants K₁ and K₂ in the following way:

$$\begin{eqnarray}\mathrm{ln}{K}_{1}(P,T) & = & -\displaystyle \frac{1}{{RT}}({\rm{\Delta }}{\bar{G}}_{f,{\mathrm{HCO}}_{3}^{-}}^{\circ }+{\rm{\Delta }}{\bar{G}}_{f,{{\rm{H}}}^{+}}^{\circ }-{\rm{\Delta }}{\bar{G}}_{f,{\mathrm{CO}}_{2}}^{\circ }\\ & & -{\rm{\Delta }}{\bar{G}}_{f,{{\rm{H}}}_{2}{\rm{O}}}^{\circ }),\end{eqnarray} \tag{ 53 }$$

$$\begin{eqnarray}&&\mathrm{ln}{K}_{2}(P,T)=-\displaystyle \frac{1}{{RT}}({\rm{\Delta }}{\bar{G}}_{f,{\mathrm{CO}}_{3}^{-2}}^{\circ }+{\rm{\Delta }}{\bar{G}}_{f,{{\rm{H}}}^{+}}^{\circ }-{\rm{\Delta }}{\bar{G}}_{f,{\mathrm{HCO}}_{3}^{-}}^{\circ }),\end{eqnarray} \tag{ 54 }$$

where R is the ideal gas constant per mole. Usually the apparent Gibbs free energies of formation are defined with respect to the value for the hydrogen ion. Hence, the actual apparent standard partial Gibbs free energy of formation for ${{\rm{H}}}^{+}$ is by definition zero.

The revised-HKF model has become widely used in geochemical calculations. The parameters for this model, for a plethora of ions and electrolytes, are given in the literature (e.g., Shock & Helgeson 1988). Also, correlation techniques that help estimate the parameters needed for this model, even for materials for which relatively little experimental data exist (Shock & Helgeson 1988), were developed.

Investigating the possible ionic composition of aqueous fluids in Earth's deep crust and mantle, with the aid of the revised-HKF model, requires the dielectric constant for water at the appropriate conditions (see the discussion in Manning 2013). Therefore, most of the effort toward determining the value for the dielectric constant is concentrated on pressures and temperatures far exceeding those we are interested in (e.g Sverjensky et al. 2014). Here, we examine two published fitted models for the dielectric constant of water: that of Marshall (2008) and that of Fernández et al. (1997). Fernández et al. (1995) tabulated most of the then available experimental data for the dielectric constant of water, which he then used to build a parametrized model (Fernández et al. 1997). The latter was approved for release by the IAPWS (latest version 1997 September). In Figure 12, we compare between the models of Marshall (2008) and Fernández et al. (1997) by examining their ability to describe known experimental data. The model of Fernández et al. (1997) is a better fit and is adopted in this study. We note that the density for liquid water used in this work is from Wagner & Pruss (2002).

Zoom In Zoom Out Reset image size

If one adopts the revised-HKF model parameters given in the literature, for the different ions and electrolytes, one also has to adopt all other formulations that were used to produce these values. This could present a problem if some of the original formulations become outdated. For example, we have used a more recent formulation for the Gibbs free energy of water (Wagner & Pruss 2002) rather than the one used to constrain many of the tabulated revised-HKF parameters (see Oelkers et al. 1995 and references therein). Another common practice is to assign neutral molecules (e.g., CO₂) with a non-zero Born coefficient to aid in fitting the experimental data. This practice is in clear contradiction to the definition of the Born coefficient. We have decided not to adopt this approach. We believe that keeping the integrity of the physical foundation of the HKF model is an important prerequisite in assessing its ability to predict experimental data. Our approach required making small corrections to the published revised-HKF parameters so that the model could fit the experimental dissociation constants of Ellis (1959) and Read (1975). In Table 3, we list the revised-HKF parameters used in this study.

Table 3.  Revised-HKF Model Parameters

Revised-HKF Parameter	CO2	${\mathrm{HCO}}_{3}^{-}$	${\mathrm{CO}}_{3}^{-2}$
${\rm{\Delta }}{\bar{G}}_{f,r}^{\circ }$ [cal mol−1]	−91881	−142110	−128010
${\bar{S}}_{r}^{\circ }$ [cal mol−1 K−1]	27.2570	22.5600	−12.3025
a1 [cal mol−1 bar−1]	0.6563	0.7409	0.5742
a2 [cal mol−1]	747.00	111.55	500
a3 [cal K mol−1 bar−1]	2.81	1.23	−2.0
a4 [cal K mol−1]	−30900	−28300	−108000
c1 [cal mol−1 K−1]	33.600	18.116	27.1800
c2 [cal K mol−1]	237600	−261800	−600000
ωr [cal mol−1]	0	120960	460000

Note. List of parameter values adopted in this study to be used with the revised-HKF model.

Download table as:  ASCIITypeset image

The revised-HKF model fits the experimental data for K₁ along the water saturation vapor pressure curve (data from Stefánsson et al. 2013) with an absolute average deviation of 6.9% (see the left panel in Figure 11). The isothermal data sets from Ellis (1959) and Read (1975) for temperatures less than 200 °C are reproduced with an absolute average deviation of a few percent. For the higher temperatures, larger deviations occur, reaching as high as 30% for the 250 °C isotherm reported in Read (1975). This latter deviation can be further reduced by adopting various models for the Born coefficient's dependence on temperature and pressure, and indeed this is often the course of action adopted. However, we have adopted a constant Born coefficient. This is because models for the Born coefficient dependence on pressure and temperature are often constrained by a fit to experimental data. A Born coefficient constrained in this way is appropriate only for the fitted ion and is likely non-transferable.

For purposes of interpolating between experimental data points, both the solvent density model and the revised-HKF model present similar accuracies. Since solvent density models are much simpler, they could prove preferable. However, the test we are more interested in is that of extrapolation. Once again, we extrapolate, now using the revised-HKF model, whose parameters were fitted to the data of Ellis (1959) and Read (1975). As before we compare these extrapolations with the high-pressure data reported in Facq et al. (2014). We list the results in Table 2. One can see that for pressures below 2 GPa, the experimental data from Facq et al. (2014) differ from the extrapolated values by a multiplicative factor of no more than 1.45. Therefore, we suggest that an error of about 50% is adequate for our modeled values for K₁ at the pressures prevailing at the bottom of water-planet oceans.

Even for the highest experimental pressures reported in Facq et al. (2014), our extrapolated values using the revised-HKF model never differ from the experimental data by more than a factor of a few. Therefore, when it comes to extrapolation, the revised-HKF model is more reliable than the solvent density model, and is therefore used in this study to model both K₁ and K₂.

There is even less reported experimental data, at high pressure and low temperature, for the second ionization constant, K₂, than there is for K₁ (see Table 1 in Dedong & Zhenhao 2007). In order to extrapolate the experimental data for K₂ up to approximately 1 GPa, for our desired low temperatures, we again use the revised-HKF model. We constrain the revised-HKF model parameters by fitting to experimental data for K₂: along the saturation vapor pressure curve (Stefánsson et al. 2013) and room temperature data up to 0.1 GPa (Hershey et al. 1983). The adopted parameters are listed in Table 3. With the aid of these adjusted revised-HKF parameters, we extrapolate to high pressure and temperature, and compare with the measurements of Facq et al. (2014) for K₂. Here as well, we adopt the disagreement between the two as our probable error.

The revised-HKF model, with our adopted parameters, fits the experimental data for K₂ along the water saturation vapor pressure curve (Stefánsson et al. 2013) with an absolute average deviation of 3.18% (see Figure 13). The absolute average deviation between our fitted revised-HKF model and the data from Hershey et al. (1983) is 0.95%. In Table 2, we list the experimental values from Facq et al. (2014) for K₂, together with the values obtained by extrapolation using the revised-HKF model. When averaged over the entire experimental pressure range tested for in Facq et al. (2014), the disagreement between extrapolated and experimental values is approximately 25%. At 1 GPa, approximately the ocean bottom pressure, the revised-HKF model predicts values for K₂ that are as much as a factor of 2 smaller than the experimental data. This factor of 2 will be our assumed error, and we will test for its significance on the carbonate speciation and ocean pH below.

Zoom In Zoom Out Reset image size

The dissociation constant of water, K_w, is taken from Bandura & Lvov (2006). This model has been adopted by the IAPWS (publication 2007 August). These authors have modeled the Gibbs free energy of the water dissociation reaction by separating it into an ideal gas and a residual contribution. The residual contribution was estimated by a sum of three terms: the work required to form proper cavities in the water bulk, the ion–water potential of interaction, and the effect of the polarization of the water bulk. The last was modeled with the aid of the Born equation. The free parameters of their model were obtained by a fit to experimental data. The resulting model is able to reproduce the experimental data, with high accuracy, in the temperature range of 0 °C–800 °C, and pressures up to 1 GPa. This P−T range covers the conditions of interest to us.