Free energy models for ice VII and liquid water derived from pressure, entropy, and heat capacity relations

Water is one of the most remarkable substances known. It displays a rich polymorphism in the solid phase, where as many as 17 different ice phases have been discovered so far.¹ One of the high-pressure polymorphs that is of great interest to planetary science is ice VII.² The unit cell of this phase consists of a body-centered cubic (bcc) oxygen structure with hydrogens located asymmetrically along the diagonals of the bcc cell.³ There is satellite evidence which suggests that ice VII exists in the subsurface of the large moons of Jupiter and Saturn.^4,5 Outside of the Solar System, both ice VII and liquid water may be present in Gliese 1214 b, which is a super-Earth that is believed to contain a vast ocean.^6,7 Evolutionary models of these planetary bodies rely on an accurate description of the thermodynamics of ice VII and liquid water, including the phase transition from one form to the other.

Over the past 15 years, several authors^8–14 have investigated the kinetics of freezing from liquid water to ice VII through ramp compression and multiple-shock compression experiments. In their setups, a thin sample of water (typically between 10 and 200 $μ$m in width) is sandwiched between two thicker windows made of materials like silica or sapphire. The water is compressed from ambient conditions, where it is a liquid, to pressures as high as about 14 GPa and temperatures of up to several hundreds of Kelvin. The final state corresponds to conditions that lie well within the ice VII region of the phase diagram. The compression is achieved by means of a current-generated magnetic field (the Sandia Z machine¹⁵) for ramp compression^11,12 and a gas gun-accelerated impactor for multiple-shock compression.^8–10,13 In either case, the phase transition from liquid to ice VII is detected by observing the decrease in pressure along the portion of the rear window that is adjacent to the water sample. The pressure decreases during the freezing process because ice VII is denser than liquid, and the reduction in volume relieves some of the pressure along the window. This illustrates why in order to develop simulation codes that can successfully reproduce experimental diagnostics it is important to have an equation of state (EOS) that can accurately predict thermodynamic properties such as the change in volume upon freezing.

There are many studies on the thermodynamic properties of ice VII and liquid water at conditions relevant to the ramp and multiple-shock compression experiments. Fortunately, the relatively low pressures attained in the experiments allow us to avoid some of the still-unresolved complications that have been observed at higher pressures (e.g., formation of a dynamically disordered ice VII phase above 40 GPa).^16–20 Nevertheless, the development of an equation of state for ice VII that covers the conditions of interest in our study is still challenging. To date, we have found only one analytic free energy model for this phase.³ In addition to ice VII, this model by French and Redmer is also applicable to two other solid phases: ice X and dynamically disordered ice VII. However, while their model may give a good overall representation for the wide range of pressures where the three ice phases may exist, it is not designed to be highly accurate in a localized range of pressures, such as those encountered in the compression experiments. In addition to the lack of a reliable free energy model in the regime of interest, there are also no direct measurements of properties like the heat capacity of ice VII. Instead, past studies have inferred the heat capacity from other data using correlations that are not self-consistent. In the case of liquid water, the gold-standard equation of state is IAPWS-95,²¹ which has been calibrated to closely match properties up to 1 GPa. There are also extensions of IAPWS-95 to higher pressures.²² The various flavors of IAPWS-95 are, however, all based on free energy models with many parameters (up to 58), which may make them computationally prohibitive and less readily implementable.

The primary goal of this work is to develop an accurate and analytic equation of state (EOS) for ice VII. Since one of our main interests is to use this EOS to simulate phase transitions under ramp and multiple-shock compression, we have also developed an EOS for liquid water. Ideally, the equations of state would give us the free energy, and all other properties would be derived from it by applying thermodynamic identities. The formulation of such a free energy model typically involves either the use of sophisticated theories from statistical thermodynamics²³ or fitting to results from ab initio quantum simulations. Clearly, both approaches require specialized knowledge and/or computational resources that are not easily attainable. But even if they were, the resulting free energy model may still yield less than satisfactory agreement with experimental data.

We follow an alternative approach that is simpler yet leads to good agreement with experiments. We formulate expressions for three derivatives of the free energy—the pressure, entropy, and heat capacity—and integrate these quantities in a thermodynamically consistent way to obtain the free energy. The EOS for both phases involves a decomposition of the pressure into an isothermal part and a thermal part. For the liquid phase, we treat the heat capacity as being a constant so that the resulting equation of state is the Mie-Grüneisen EOS.²⁴ There are already some published Mie-Grüneisen equations of state for liquid water,^25–29 and our work complements the existing studies by comparing to heat capacity data that justifies the assumption of constant heat capacity for the regime of interest here. In addition, we calculate the sound speed along various isotherms to further confirm the applicability of the Mie-Grüneisen EOS to liquid water. Unlike in the past studies,^25–29 the Grüneisen parameter in our work is constrained by data at ambient conditions. Our free energy model for ice VII is also derived from an expression for the pressure, but it is more complicated because we do not treat the heat capacity of this phase as being a constant. The heat capacity is instead fit to reproduce melt temperature data. We show that our ice VII EOS provides good agreement with other properties along the melt curve. It is able to calculate the sound speed more accurately than the model by French and Redmer at the conditions of interest. We also compare our predictions for the Hugoniot of both ice VII and liquid water with experimental results.

As stated in the Introduction, we wish to decompose the pressure P = P(V,T) as a function of the volume V and temperature T into two terms,

where P_isotherm(V) is the pressure along an isotherm at an arbitrary reference temperature T_ref and P_thermal(V, T) is the thermal contribution to the pressure. We favor this type of formulation for P because hydrodynamic simulation codes for high-pressure applications, such as Ares³⁰ and ALE3D,³¹ tend to assume this functional form for the pressure. Our aim is to treat (1) as an equation of state from which we derive all other thermodynamic properties. If we make the additional assumption that the isochoric heat capacity C_V is a constant, the resulting equation of state is the Mie-Grüneisen EOS. We show in Sec. II B that it is justifiable to treat C_V of liquid water as a constant in the pressure and temperature ranges relevant to the ramp and multiple-shock compression experiments.^8–14 Although the Mie-Grüneisen EOS is rather simple, we present a brief derivation of it in this section for the sake of completeness. In addition, many of the steps in the derivation are also applicable to—and thus help to elucidate—the derivation of the more general EOS described in Sec. III A, where C_V may depend on both volume and temperature.

We begin by briefly reviewing a few basic thermodynamic relations used throughout this study. The first and second laws applied to a closed and simple (i.e., free from the influence of external fields) one-component system yield the following fundamental equation:³² 

where E is the internal energy and S is the entropy. From (2), we obtain

which leads to

where the isochoric heat capacity C_V is defined as $CV=(∂E/∂T)V$⁠. From (3), we have

Equation (4) is a relation between P and C_V that must be obeyed by any thermodynamically consistent equation of state. Furthermore, (2) and (3) together require that

The Grüneisen parameter $Γ$ is a dimensionless quantity that is defined²⁴ as

in which

is the coefficient of thermal expansion and

is the isothermal bulk modulus. The Grüneisen parameter gives an indication of how pressure varies with temperature,

There are many well-established forms for P_isotherm, such as the Vinet or Birch-Murnaghan equation,^24,33 so the main task is in deriving the functional form of P_thermal. Equation (4) in the case of constant C_V implies that the thermal pressure must have the functional form,

where f(V) depends only on volume. To determine what f(V) is, we integrate (8) with respect to temperature. From (8), it is clear that in order for P_thermal to have the form in (9), the Grüneisen parameter shall be independent of temperature so that $Γ=Γ(V)$⁠. Carrying out the integration,

We may therefore make the identification $f(V)=ΓCV/V$⁠. The pressure P(V, T) is now completely determined if we assume that $Γ(V)$ is known and that P_isotherm(V) has been fitted to P(V) data at T = T_ref. Specific forms of $Γ$ and P_isotherm applicable to liquid water are presented in Sec. II C. In principle, knowledge of P(V, T) provides a sufficient representation of the equation of state, since all other quantities may be obtained from P by applying thermodynamic relations. Towards this goal, we derive expressions for the Helmholtz energy F(V, T) and the entropy S(V, T) that are consistent with P(V, T). Once expressions for F and S have been established, other thermodynamic properties such as the internal energy E(V, T) = F + TS, enthalpy H(V, T) = F + TS + PV, and Gibbs energy G(V, T) = F + PV are straightforward to evaluate.

Since $P=−(∂F/∂V)T$ and $S=−(∂F/∂T)V$⁠, F and S are given by

and

where V_ref is a reference volume. The only unknown quantity in (11) and (12) is F(V_ref, T). In order to ascertain what this function is, another expression for F or S is needed. This expression is provided by (5), which may be integrated to get

By equating (12) with (13), we make the following identifications:

The former is obtained by grouping together the terms in (12) and (13) that depend only on temperature, while the latter is obtained by grouping together the terms that depend only on volume. Combining (13) and (15) leads to the complete expression for S,

Integrating (14), we have

where F_ref = F(V_ref, T_ref). Substituting (17) into (11) yields the full expression for F,

in which

A defining feature of the Mie-Grüneisen equation of state is that the isochoric heat capacity C_V is a constant. The behavior of C_V of liquid water is illustrated in Fig. 1(a). The figure shows that C_V depends on both pressure (volume) and temperature. Strictly speaking, it should not be modeled as being a constant. However, it can be approximated as a constant for the purposes of the ramp and multiple-shock compression experiments.^8–14 These experiments take water from ambient conditions (10⁻⁴ GPa = 1 bar and 300 K) through a loading path that is isentropic in the case of ramp compression and quasi-isentropic in the case of multiple-shock compression. Figure 1(b) shows a representative example of the liquid water isentrope. If the phase transition from liquid to ice VII is initiated by homogeneous nucleation, the experiments have found that freezing occurs when the pressure reaches around 7 GPa.^11–13 Due to the (quasi-)isentropic nature of the loading path, the temperature below this pressure remains relatively low. After the phase transition goes to completion, the pressure continues to rise to a final pressure of 8–14 GPa, depending on the conditions of the experiment.

Figure 1(a) suggests that over the pressure and temperature range stated above, C_V of liquid water may be approximated as a constant of about 7R, in which R is the gas constant. For instance, according to the IAPWS-95 EOS, C_V varies between 5R and 9R in the temperature range 300–973 K at 0.1 GPa, while it varies more narrowly between 5.7R and 7.5R in the same temperature range at 1 GPa. We assume C_V to be a constant equal to 7R, as indicated by the solid horizontal line in the figure. This is very close to C_V of 7.1R chosen by Gurtman et al.²⁵ for their Mie-Grüneisen EOS. They have found that this value of C_V, when combined with a properly chosen functional form for the Grüneisen parameter (which we discuss in Sec. II C), closely reproduces the liquid water Hugoniot curve of Walsh and Rice³⁸ up to 24 GPa. We verify this claim in Sec. II D. Lyzenga et al.²⁷ have found, however, that in order to match their temperature measurements along the Hugoniot between 50 and 80 GPa, a Mie-Grüneisen EOS with a larger heat capacity of C_V = 8.7R is necessary. Since we are interested in pressures far below 50 GPa, we make the choice of C_V = 7R.

In addition to choosing an appropriate value for C_V, the other main task in developing a Mie-Grüneisen EOS is determining the functional form of the Grüneisen parameter $Γ$⁠. Once C_V and $Γ$ are both established, the pressure P = P(V, T)—as well as other properties like the Helmholtz energy F(V, T) and entropy S(V, T)—may be derived in a systematic manner, as we have demonstrated in Sec. II A. Since $Γ$ must necessarily be independent of T if C_V is a constant, determining the functional form of $Γ$ boils down to finding $Γ(V)$⁠.

Figure 2 illustrates our $Γ=Γ(V)$ relation, along with those from previously published Mie-Grüneisen equations of state^25–29 for liquid water. The $Γ(V)$ functional forms presented by the earlier studies are obtained from two different sources: (1) by fitting a finite difference approximation of (8) to discrete P(V) points along different isotherms and (2) by fitting $Γ(V)$ to Hugoniot data at higher pressures. Figure 2 shows that there is significant variation in the form of $Γ$⁠, despite the fact that all of the studies demonstrate good agreement with PVT data. This suggests that $Γ$ is somewhat insensitive to the PVT data so that more constraints are necessary to uniquely determine $Γ(V)$⁠. One such constraint on $Γ$ comes from properties at ambient conditions, where the value of $Γ$ is well-established because all of the quantities on the right-hand side of (6) have been measured to a high degree of precision. We take values for $α$⁠, B, and V from the work of Kell³⁹ and C_V from IAPWS-95 to find that $Γ=0.1489$ at ambient conditions. This is comparable to the value of $Γ=0.11$ reported by Boehler and Kennedy,⁴⁰ who obtained $Γ=(BS/T)(∂T/∂P)S$ directly by measuring the isentropic derivative $(∂T/∂P)S$ in a cylindrical piston apparatus and multiplying their result with the isentropic bulk modulus B_S. We represent $Γ$ as a fourth-order polynomial in V,

where the coefficients and the volume V_ambient are presented in Table I. Our form for $Γ(V)$ meets three criteria: (1) produces the expected value of $Γ=0.1489$ at V = V_ambient, (2) yields good agreement with PVT data in the pressure and temperature range of interest (to be demonstrated shortly), and (3) satisfies the Hugoniot-fitted results obtained by Gurtman et al.²⁵ which we have presented in Fig. 2.

Equation (1) requires that we fit P(V) data along an isotherm at an arbitrary reference temperature T_ref, which we take to be 300 K. We choose to model the pressure P_isotherm(V) = P(V, T_ref) along the isotherm with the third-order Birch-Murnaghan equation,^24,33 which is traditionally expressed as

where B_ref is the isothermal bulk modulus at a reference volume V_ref and $Bref′$ is the pressure derivative of B_ref at V_ref. This equation may also be expressed as a pseudo-virial expansion,

in which the coefficients a_5/3, a_7/3, and a₃ are constants given by

The pseudo-virial form in (22) is more useful than the traditional form for practical purposes because it can more readily be differentiated and integrated with respect to volume. We adjust V_ref so that P_isotherm(V_ambient) = 10⁻⁴ GPa = 1 bar. Since $Pisotherm(Vref)=0≈10−4$ GPa, B_ref is expected to be very similar to the value of B at ambient pressure, which Kell³⁹ reports to be 2.22 GPa. Thus, the only true adjustable parameter in P_isotherm is $Bref′$⁠. Table II lists the value of $Bref′$⁠, along with those for the other parameters of our EOS.

In summary, the pressure P(V, T) is given by (1), where P_isotherm(V) and P_thermal(V, T) are represented by (22) and (10), respectively, and $Γ(V)$ is computed from (21). Figure 3 compares our EOS results with published data along eight isotherms at 300 K, 373 K, 473 K, …, and 973 K. Our EOS does not match the data at volumes much greater than V_ambient because $Γ$ becomes zero around V = 18.57 cm/mol³ (all isotherms converge at this volume), and it becomes negative for volumes greater than this value. However, from the discussion in Sec. II B [see, e.g., the isentrope in Fig. 1(b)], it is clear that the low pressure and high temperature conditions that correspond to these large volumes are not relevant to the ramp and multiple-shock compression experiments. We achieve good agreement for the PVT conditions that are traversed by these experiments. Along the 300 K and 373 K isotherms, the average deviation in the volume is 0.35%, and the maximum deviation is 1.5%. For the other six (higher temperature) isotherms, these numbers are 0.29% and 1.6%, respectively, if we examine only the data points for which $P≥2.5$ GPa.

The Helmholtz energy F(V, T) is given by (18)–(20) and the entropy S(V, T) is obtained from (16). In these equations,

One property we analyze to validate our EOS is the sound speed c_s, which is given by

where M is the molecular weight and the isentropic bulk modulus B_S is calculated from

Equation (24) may be obtained by combining (8) with the relations B_S = (C_P/C_V)B and $CP=CV+α2TBV$⁠, where C_P is the isobaric heat capacity. Figure 4 shows our results for c_s along three isotherms at 300 K, 473 K, and 673 K. Since we are not interested in the low pressure and high temperature regime, we do not examine the behavior of c_s below 0.5 GPa at 473 K and below 1 GPa at 673 K. However, we do consider the entire pressure range shown in the figure for the 300 K isotherm. The average deviation over all the data points is 2.8%, and the maximum deviation is 12.7%. Most of the departure from the published data occurs at 673 K. In fact, at this temperature, there is noticeable disagreement above 6 GPa between the studies of Decremps et al.⁴³ and Sanchez-Valle et al.,³⁶ even though both studies involve the same technique (Brillouin scattering). Our results at these higher pressures lie between those from the two studies.

Another way to validate our equation of state is to examine the Hugoniot, which is the locus of states that satisfy the following Rankine-Hugoniot equations:

The three equations express the conservation of mass, momentum, and energy, respectively, across a shock front. Here, u_p is the particle speed and u_s is the shock speed, while E₀ and P₀ represent the internal energy and pressure of the undisturbed material ahead of the shock, whose state is characterized by its volume V₀ and temperature T₀. Figure 5 compares our prediction for the Hugoniot with experimental measurements. We find good agreement with data from the work of Walsh and Rice, which is not surprising because as we have mentioned at the end of Sec. II B, Gurtman et al. have reported that a Mie-Grüneisen EOS with $CV≈7R$ can closely reproduce the results from that study. In addition, our calculations match well with the low-pressure measurements of Lysne.⁴⁴ Although not shown in the figure, we have also verified that our EOS agrees with the particle speeds measured by Lysne.

The constant C_V approximation in the case of the Mie-Grüneisen equation of state greatly simplifies the derivation of F(V, T) and S(V, T) from P(V, T). Nevertheless, we show in this section that it is straightforward to derive F and S from P even in the general case where C_V = C_V(V, T). Let us assume that we have fit pressure-volume-temperature data to a function P(V, T) of the form in (1). The Helmholtz energy and entropy are then given by

We again use (5), and this requires an expression for C_V. Since the heat capacity and the pressure are related according to (4), we have

Evaluating the temperature part of the double integral, (31) simplifies to

Just like in Sec. II A, we equate the two expressions for entropy [which in this case are (29) and (32)] and group together the terms that depend on temperature only to obtain

We integrate (33) to get

Combining (1) with (28) and (35), we may again represent F in form (18), where the free energy F_isotherm(V) is given by (19) and the thermal contribution F_thermal(V,T) is

One may verify that if C_V is a constant so that P_thermal must be of the form in (10), (34) and (36) simplify to the Mie-Grüneisen expressions (16) and (20), respectively. We note that C_V(V_ref,T) is analogous to the Grüneisen parameter $Γ(V)$ in that both functions can be determined only by fitting to experimental data; they cannot be established a priori. Section III C describes how C_V(V_ref, T) of ice VII may be determined by fitting it to temperature measurements along the melt curve.

The analysis in Sec. III A presumes that a functional form of P(V, T) has been provided. We represent the pressure of ice VII with (1), where the isothermal part P_isotherm(V) at T_ref is again described by the third-order Birch-Murnaghan equation (22). The thermal contribution P_thermal(V, T) is modeled as

where $η$⁠, V_thermal, and T_thermal are constants. The parameters of our ice VII equation of state are listed in Table III. These values yield generally good agreement with PVT data, though as indicated in Fig. 6 there are appreciable differences among some of the studies, so it is not possible to closely match all the published results. Sugimura et al.¹⁹ discuss some possible explanations (e.g., presence of non-hydrostatic environments and differences in pressure scale) for these discrepancies. The studies may be divided into two groups: for a given pressure, the volumes measured by Loubeyre et al.,⁴⁶ Sugimura et al.,^19,47 and two studies which are not shown in the figure^48,49 are significantly smaller than the volumes obtained by Hemley et al.,⁵⁰ Fei et al.,⁵¹ and Frank et al.⁵² The volumes reported by Bezacier et al.⁴⁵ lie between the extremes represented by the two aforementioned sets of studies.

As we have mentioned in the Introduction, French and Redmer³ have developed a free energy model for ice VII. To the best of our knowledge, theirs is the only such analytic model published to date. The free energy in their study is composed of three parts: (1) the electronic ground-state energy, (2) the nuclear ground-state energy, and (3) a thermal contribution from the motion of the nuclei. Each of these three parts is fit to results from their density functional theory molecular dynamics simulations. They have tried several different exchange-correlation functionals to compute the electronic ground-state energy, and they report that the Heyd-Scuseria-Ernzerhof functional gives the best agreement with PVT data. Figure 6 shows that their model with this functional tends to overpredict the volume (underpredict the density) below 15 GPa. For these pressures, their 300 K isotherm lies to the right of nearly all the 450 K data.

With the exception of Fei et al.,⁵¹ which we discuss more in Sec. III C, past experimental studies do not decompose the pressure into an isothermal and a thermal part as in (1). Rather, they fit volume measurements at 300 K to an isothermal expression (such as the Vinet or Birch-Murnaghan equation) and find the volume at other temperatures by integrating (7) with respect to temperature, where the thermal expansion coefficient $α$ is represented as

in which $α0$⁠, $α1$⁠, and $ξ$ are constants. A wide range of values have been published for $α$ at ambient conditions. It has been reported to be $α=0.48×10−4$ K⁻¹, $0.6×10−4$ K⁻¹, $1.16×10−4$ K⁻¹, and $15.0×10−4$ K⁻¹ by Frank et al.,⁵² Fei et al.,⁵¹ Bezacier et al.,⁴⁵ and Sugimura et al.,¹⁹ respectively. According to our EOS, $α$ at ambient conditions is about $2.45×10−4$ K⁻¹, which lies between that of Bezacier et al. and Sugimura et al. We have calculated this result by rearranging (8) to obtain $α=(∂P/∂T)V/B$⁠. The large variation in the reported ambient value of $α$ is a reflection of a few different factors: (1) experimental uncertainties in the volume measurements, (2) the choice of the 300 K isothermal EOS and the subsequent extrapolation to ambient pressure, (3) the fact that the various studies have examined different pressure ranges in their volume measurements.

Equation (36) shows that C_V(V_ref, T), which is the portion of the heat capacity that depends on temperature only, affects the Helmholtz energy through the double integral

We have stated in Sec. III A that C_V(V_ref, T) is an adjustable quantity, similar to how we have freedom in modeling the volumetric dependence of $Γ$ in the case of the Mie-Grüneisen equation of state. Thus, (39) is also an adjustable quantity. We now describe how (39) and the constant F_ref may be calibrated with data along the liquid water-ice VII melt curve. First, we temporarily set F_ref = 0 so that $Fisotherm(V)=∫VVrefPisothermdV′$⁠. Since the Gibbs energy G = F + PV of the liquid and ice VII phases is equal to each other along the melt curve, this equilibrium condition places a constraint on the Helmholtz energy F of the ice phase. Consequently, this condition also fully specifies how (39) must vary with temperature because we have already established P(V, T) of the ice phase in Sec. III B and all contributions to F besides (39) are determined by the pressure. Figure 7(a) shows a number of melt curves reported in the literature. Most of the melt curves are consistent with each other at lower pressures. The studies presented in the figure^37,51–57 and the review by Dunaeva et al.²⁰ have proposed various explanations for the observed discrepancies at higher pressures.

We have found that our EOS can exactly reproduce any of the melt curves in Fig. 7(a) if we represent (39) as a fourth-order polynomial in temperature, but this causes C_V(V_ref, T) to become unrealistically large with increasing temperature. A good compromise may be achieved if we instead represent (39) with a second-order polynomial; this produces a much more realistic heat capacity and still allows us to closely match any desired melt curve. We fit the second-order polynomial to the melt curve of Datchi et al.³⁷ because we have found that choosing this melt curve leads to the best agreement with sound speed data. We discuss more about our results for the sound speed in Sec. III D. Thus, we have

where the coefficients F_k are listed in Table IV. Since we have zeroed out F_ref in F_isotherm, this reference free energy is given by $Fref=∑k=02FkTrefk=4.1109kJ/mol$⁠. In summary, we model the Helmholtz energy of ice VII as

where the pressure is given by the expressions in (22) and (37). All other properties may be derived from (41) in a thermodynamically consistent manner. For example, the entropy is

Some properties along the melt curve are shown in Figs. 7(b)–7(d). Goncharov et al.⁵⁹ have applied X-ray diffraction to estimate the density of both phases. Asahara et al.⁵⁸ have performed a combination of X-ray diffraction and Brillouin scattering to obtain their liquid density measurements. Our predictions for the liquid density lie within the uncertainties of the measurements in both studies. Our ice VII densities are smaller than those of Goncharov et al., which implies that if their results were included in Fig. 6, they would be located quite a bit to the left of the results from Fei et al. Nevertheless, if the uncertainties in their ice measurements are similar to those in their liquid measurements (they do not present error bars for the ice phase), our results would lie within the uncertainties. Since the French and Redmer model tends to underpredict the density for pressures below 15 GPa, their ice VII density curve would lie below ours, even further away from the results of Goncharov et al. We have explained in the Introduction that an accurate prediction of the volume change $ΔV$ is an important consideration for simulations of phase transitions under ramp and multiple-shock compression. The entropy change $ΔS$ of melting is calculated from (16) and (42). The high-pressure results in Fig. 7(c) may be thought of as an extrapolation of our low-pressure results, where we are more confident because of greater availability of data. At lower pressures, there is close agreement for the liquid sound speed among all of the studies in Fig. 7(d), including ours. At higher pressures, our results lie between those of Asahara et al.⁵⁸ and Ahart et al.⁶⁰ The average deviation and maximum deviation for the liquid sound speed data in the figure are 5.2% and 14.5%, respectively.

We see from (30) that because we have modeled the thermal pressure as being linear in temperature, the heat capacity C_V must be independent of volume so that C_V = C_V(V_ref, T),

The solid line in Fig. 8 illustrates the temperature dependence of our C_V. The heat capacity from the French and Redmer EOS³ also increases linearly with temperature over the temperature range shown in the figure, as depicted by the two representative isochores. We have explored introducing a volume dependence to our C_V by modeling the thermal pressure with the following expression:

Because this expression for pressure has a nonlinear dependence on the temperature, the resulting C_V is a function of both volume and temperature. However, we have found that C_V(V, T) is negative for certain values of V and T, and so it is not physically valid at these conditions. As a result, we elect to keep (37) so that C_V is given by (43). These expressions are sufficient for our purposes because they are relatively simple, yet they lead to good agreement with both melt curve and sound speed data for the conditions of interest.

Because ice VII exists as a high-pressure polymorph, its heat capacity has not been measured. Two studies^51,62 have inferred the heat capacity from other measurements, but the resulting expressions are not self-consistent. The main focus of Fei et al.⁵¹ is on fitting their PVT data with (38), but they briefly mention that their data may also be modeled with a thermal pressure of the form

where the internal energy E is calculated from the Debye temperature $θD$ through C_V, which is given by⁶³ 

They report the Debye temperature $θD$ to be 1470 K. The integral in (44) has an exact analytical expression,⁶⁴ 

where $Lin(x)=∑k=1∞xk/kn$ is the polylogarithm function of order n. In addition to (44), Fei et al. have deduced the following expression for the isobaric heat capacity C_P from their melt temperature measurements:

in which a, b, and c are constants. Stewart and Ahrens,⁶¹ which we discuss more in Sec. III D, have developed a Debye model to closely match C_P in (45). Figure 8 shows that C_P of Fei et al. is more than three times greater than their C_V. Since C_P is expected to be only slightly larger than C_V, we conclude that (44) is not consistent with (45).

Asahara et al.⁶² have measured the sound speed along the 300 K isotherm, and they obtained a correlation for the isentropic bulk modulus B_S vs. pressure by fitting their sound speed results with the Vinet equation of Sugimura et al.⁴⁷ The following equation is used to determine C_P vs. pressure along the isotherm:

This equation may be derived by rearranging (24) and substituting the relations $CP=CV+α2TBV$ and $(∂P/∂T)V=αB$⁠. Asahara et al. find V by inverting the Vinet fit of Sugimura et al., compute $α$ from (38) using the parameters in the work of Fei et al.,⁵¹ and take B_S from their own work. Values for B are chosen from four different studies.^47,51,65,66C_P obtained using B from three of the studies^47,65,66 seem to agree fairly well with each other, though they are still not fully consistent. At lower pressures, they deviate significantly from the C_P curve obtained using B from the fourth study.⁵¹ Asahara et al. offer no opinion on which source is most credible, leaving the reader with ambiguity regarding the correct value for C_P. Furthermore, their heat capacities seem to be quite large. For example, their expression produces a $CP$ at ambient pressure and 300 K that is greater than 54R [25 (J/g)/K], whereas the heat capacity is less than 10R for all of the studies shown in Fig. 8.

We calculate the sound speed in ice VII from (23) and (24). Results are shown in Fig. 9. In order to reconstruct the sound speed measurements of Asahara et al.⁶² along the 300 K isotherm [the dotted curve in Fig. 9(a)], we reverse their procedure described above. That is, we combine their B_S vs. pressure correlation with the Vinet fit of Sugimura et al.⁴⁷ to compute the sound speed c_s from B_S and V by applying (23). The ice VII sound speeds calculated from our EOS tend to lie below theirs (especially at higher pressures), which is also the trend that we have observed for the melt-curve liquid sound speeds in Fig. 7(d). Nevertheless, there is fairly close agreement throughout the entire pressure range. Figure 9(b) depicts the temperature dependence of the sound speed at three different pressures. At 14 GPa, the sound speeds measured by Ahart et al.⁶⁰ are noticeably slower than those measured by Asahara et al. Again, this is consistent with what one might expect from the liquid sound speeds in Fig. 7(d). The conditions of interest in this study cover only a small portion of the much wider range of conditions considered by French and Redmer.³ As expected, our equation of state more closely matches experimental data than their model does over the relatively narrow range of interest.

Experimentally measuring the Hugoniot is a far more complicated matter for ice VII than it is for liquid water because ice VII is a non-quenchable phase. First of all, one may need to cool the system in order to experimentally obtain a Hugoniot that starts from ambient pressure and probes ice VII states. Stewart and Ahrens⁶¹ report to have measured shock Hugoniot curves that pass through the ice VII region of the phase diagram by using liquid nitrogen to achieve this cooling. They have shocked ice at two different initial temperatures: T₀ = 100 K and T₀ = 263 K. Their study does not indicate the initial volumes V₀, but from a data table we have obtained through private communication,⁶⁷ we have found V₀ to be 19.318 cm³/mol (0.932 g/cm³) and 19.575 cm³/mol (0.920 g/cm³) for the T₀ = 100 K and T₀ = 263 K experiments, respectively. One complication is that the ice phase present in the initial states of the Stewart and Ahrens experiments is not ice VII, but it is presumed to be ice I_h. Thus, we cannot naïvely apply our ice VII EOS to compute the internal energy E₀. Instead, we obtain E₀ from the five-phase tabular EOS presented in the work of Senft and Stewart.⁶⁸ This tabular EOS is designed to reproduce the shock measurements of Stewart and Ahrens. The table is offset in energy from our EOS by some arbitrary constant, and we set the value of this quantity to reproduce the expected freezing point at ambient pressure. That is, we define the energy offset so that the Gibbs energy of ice I_h returned by the Senft and Stewart tabular EOS at P = 10⁻⁴ GPa (1 bar) and T = 273 K (⁠$0°$C) is equal to the Gibbs energy computed from our liquid water EOS at these conditions. See Table V for a summary of the Hugoniot initial states. The E₀ values listed in the table are found by adding the energy offset to the internal energy returned by the tabular EOS for the corresponding V₀ and T₀. Note that because French and Redmer³ do not present a liquid water EOS, we cannot apply the above procedure to determine their energy offset with respect to the tabular EOS.

Figure 10 compares our results for the ice VII Hugoniot with experimental data for this phase from the work of Stewart and Ahrens.⁶¹ They state that a few different phases of water, including ice VII, are sampled by the Hugoniot curves in their study. Our results do not agree well with the experimental data. For instance, Stewart and Ahrens have found that at a pressure of 2.2 GPa, u_p is 0.6 km/s and u_s lies between 1.8 and 2.3 km/s, while our EOS predicts that u_p is about 0.9 km/s at this pressure and u_s is about 2.5 km/s. The reasons for this discrepancy may be understood from Fig. 11, which compares the shock measurements of Stewart and Ahrens to static measurements of the ice VII PVT behavior. The figure shows that their T₀ = 100 K data are clearly not consistent with the static measurements, but on first glance, it appears as though their T₀ = 263 K data might be consistent.

If we were to try adjusting our P(V, T) model to fit the Stewart and Ahrens T₀ = 263 K data in an attempt to ultimately reproduce their reported ice VII Hugoniot, we would encounter additional problems that would make the fulfillment of this goal impossible. First of all, in order to have P(V, T) match both static and shock data, the P(V) isotherms must taper off significantly [$(∂P/∂V)T$ would need to become very small] for $V>11$ cm³/mol. In other words, ice VII would need to undergo an unrealistic, abrupt softening. Furthermore, the very small $(∂P/∂V)T$ at larger volumes means that our reference volume V_ref (the volume at zero pressure and T = T_ref) would have to be very large, say 20 cm³/mol. None of the T_ref = 300 K isotherms published for ice VII have V_ref of anywhere near that value (see Table II from the work of Bezacier et al.⁴⁵). Moreover, we have pointed out that the Fei et al. heat capacity C_P in (45) on which the heat capacity in the work of Stewart and Ahrens is based is very different from C_V in (44) that Fei et al. infer from their PVT measurements. All of the aforementioned factors together suggest that there is indeed some level of inconsistency between the static and the shock data, assuming that the heat capacity model (i.e., Hugoniot temperatures) in the work of Stewart and Ahrens is accurate. It would be difficult, if not impossible, to develop a pressure model that satisfies both sets. Now, even if we could somehow constrain our P(V, T) to agree with both static and shock data, we would also have to adjust our C_V(V_ref, T) to obtain the same internal energy as Stewart and Ahrens do in order to match their Hugoniot. Considering how different our heat capacity is from theirs, it is very unlikely that after doing so, we would obtain a reasonable melt curve. All of this discussion presumes that the many difficulties which inevitably arise in cryogenic shock experiments (e.g., creation of voids and other defects) are well controlled in the work of Stewart and Ahrens. Given the inconsistencies between the static and the shock data, we choose to calibrate our ice VII EOS with the former because the static measurements represent a compilation of results from many studies, several of which have performed X-ray diffraction to confirm that they obtained ice VII. In contrast, the shock data come from only one study,⁶¹ and it does not present diffraction measurements. We encourage more studies on the shock properties of ice VII to help resolve the discrepancies.

Some interesting issues may be examined by plotting the liquid Hugoniot, plus a few isentropes that branch off from the Hugoniot, together with the melt curve (Fig. 12). The principal isentrope is a two-phase isentrope in the following sense: to the left (right) of the melt curve, the isentrope represents the liquid (ice VII) phase. At the melt curve, the isentrope traverses upward until it reaches conditions along the melt curve where the ice VII phase has the same entropy as the liquid. Isentropes 2 and 3 are also two-phase isentropes; together they give an idea of how one might experimentally probe high-pressure, high-temperature states of ice VII. According to our EOS, there is a transition point in which isentropes that branch off from the Hugoniot above a pressure of 6.9 GPa (corresponding to a temperature of 631 K) will remain in the liquid phase. Thus, isentrope 4, which branches off at 8 GPa, is a liquid-only isentrope.

Figure 12 shows that for pressures between 2.5 and 6 GPa, the liquid Hugoniot from our EOS resides below our melt curve. As stated in the caption, we refer to the Hugoniot in the figure as the principal Hugoniot because its initial state is at ambient conditions (P = 1 bar, T = 300 K). The fact that the principal Hugoniot seems to cross the melt curve means that it might be possible to form ice VII with a single shock initiated at ambient conditions, which is a long-standing controversy that has still not been resolved.^8,69–71

While ice VII formation on the principal Hugoniot may be possible, we cannot definitively conclude that it will occur. First of all, we note that the density change from liquid to ice VII is rather significant [roughly 10%, see Fig. 7(b)]. It seems unlikely that such a large volume change could occur for freezing under shock compression without it being obvious in Hugoniot measurements. Furthermore, there is experimental uncertainty in the precise determination of both the melt curve and the Hugoniot, an important consideration given that the field of ice stability along the Hugoniot predicted by our EOS is quite small. It is possible that the Hugoniot in reality does not actually cross the melt curve. Even if the two curves cross each other, ice VII formation may still not be observable with the principal Hugoniot because of issues relating to metastability and kinetics. It is likely that the Hugoniot only barely enters the ice VII region, meaning that the driving force $ΔG$ for the phase transition to ice VII will be low. Consequently, the time duration needed for the phase transition to reach completion may be longer than the time duration of a single-shock experiment, which typically lasts only a few microseconds.⁸ Thus, even if the principal Hugoniot crosses the melt curve, the water may remain in a metastable liquid state. These considerations illustrate why it is important to examine not only the thermodynamics of phase transitions but also the kinetics—a topic that is the focus of our current efforts.⁷² 

We have presented equations of state for ice VII and liquid water at conditions relevant to ramp and multiple-shock compression experiments.^8–14 In our approach, we formulate an expression for the pressure, and we relate it to two other derivatives of the free energy—the entropy and the heat capacity—to obtain the free energy. This is essentially the opposite of the conventional approach, where a free energy model is developed first and all other properties are derived from it. Experimental data on the free energy tend to be far less commonly available than data on its derivatives, particularly the pressure and the heat capacity. Thus, our approach to build up an equation of state by starting with the free energy derivatives is more natural and practical. The procedure we have described in Sec. III A for doing this is very general, and it may therefore be applied to develop equations of state for virtually every other condensed-phase or gaseous single-component substance.

Our ice VII and liquid water equations of state involve only polynomials and logarithms, allowing all thermodynamic properties to be computed analytically in an efficient manner. Despite their simplicity, the equations of state are also accurate. We have examined the pressure-volume-temperature behavior, the sound speed, and the Hugoniot, plus various properties along the melt curve, and the average deviation from experimental measurements for nearly all of these properties is around 5% or less. The only exception is the ice VII Hugoniot, where we have not found good agreement with the work of Stewart and Ahrens.⁶¹ However, we have pointed out that the pressure-density results in their study, which is the only experimental work that reports data on the ice VII Hugoniot, do not seem to be consistent with static measurements obtained from many other studies performed in diamond anvil cells. We have demonstrated how our equations of state may be utilized to help design experiments that generate isentropes which probe high-pressure, high-temperature states of ice VII. We have also briefly examined the possibility of ice VII formation along the path of the principal liquid Hugoniot. Although our equations of state are targeted towards ramp and multiple-shock compression experiments, we have demonstrated that they also behave reasonably well at pressures and temperatures that lie somewhat beyond those found in the experiments.

This work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract No. DE-AC52-07NA27344. We are grateful to S. T. Stewart for sending us experimental data, and we thank M. French, S. Hamel, R. G. Kraus, M. A. Millot, and R. Redmer for helpful discussions.