We present equations of state relevant to conditions encountered in ramp and multiple-shock compression experiments of water. These experiments compress water from ambient conditions to pressures as high as about 14 GPa and temperatures of up to several hundreds of Kelvin. Water may freeze into ice VII during this process. Although there are several studies on the thermodynamic properties of ice VII, an accurate and analytic free energy model from which all other properties may be derived in a thermodynamically consistent manner has not been previously determined. We have developed such a free energy model for ice VII that is calibrated with pressure-volume-temperature measurements and melt curve data. Furthermore, we show that liquid water in the pressure and temperature range of interest is well-represented by a simple Mie-Grüneisen equation of state. Our liquid water and ice VII equations of state are validated by comparing to sound speed and Hugoniot data. Although they are targeted towards ramp and multiple-shock compression experiments, we demonstrate that our equations of state also behave reasonably well at pressures and temperatures that lie somewhat beyond those found in the experiments.

## I. INTRODUCTION

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.^{1} One of the high-pressure polymorphs that is of great interest to planetary science is ice VII.^{2} 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.^{3} 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 $\mu $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^{15}) 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.^{3} 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,^{21} which has been calibrated to closely match properties up to 1 GPa. There are also extensions of IAPWS-95 to higher pressures.^{22} 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^{23} 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.^{24} 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.

## II. LIQUID WATER

### A. The Mie-Grüneisen equation of state

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^{30} and ALE3D,^{31} 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:^{32}

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\u2009=\u2009(\u2202E/\u2202T)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 $\Gamma $ is a dimensionless quantity that is defined^{24} 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 $\Gamma \u2009=\u2009\Gamma (V)$. Carrying out the integration,

We may therefore make the identification $\u2009f(V)\u2009=\u2009\Gamma CV/V$. The pressure *P*(*V*, *T*) is now completely determined if we assume that $\Gamma (V)$ is known and that *P*_{isotherm}(*V*) has been fitted to *P*(*V*) data at *T* = *T*_{ref}. Specific forms of $\Gamma $ 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=\u2212(\u2202F/\u2202V)T$ and $S=\u2212(\u2202F/\u2202T)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

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

in which

### B. Constant heat capacity approximation

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^{−4} 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 7*R*, in which *R* is the gas constant. For instance, according to the IAPWS-95 EOS, *C*_{V} varies between 5*R* and 9*R* in the temperature range 300–973 K at 0.1 GPa, while it varies more narrowly between 5.7*R* and 7.5*R* in the same temperature range at 1 GPa. We assume *C*_{V} to be a constant equal to 7*R*, as indicated by the solid horizontal line in the figure. This is very close to *C*_{V} of 7.1*R* chosen by Gurtman *et al.*^{25} 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^{38} up to 24 GPa. We verify this claim in Sec. II D. Lyzenga *et al.*^{27} 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.7*R* is necessary. Since we are interested in pressures far below 50 GPa, we make the choice of *C*_{V} = 7*R*.

### C. Mie-Grüneisen representation and pressure-volume-temperature behavior

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 $\Gamma $. Once *C*_{V} and $\Gamma $ 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 $\Gamma $ must necessarily be independent of *T* if *C*_{V} is a constant, determining the functional form of $\Gamma $ boils down to finding $\Gamma (V)$.

Figure 2 illustrates our $\Gamma =\Gamma (V)$ relation, along with those from previously published Mie-Grüneisen equations of state^{25–29} for liquid water. The $\Gamma (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 $\Gamma (V)$ to Hugoniot data at higher pressures. Figure 2 shows that there is significant variation in the form of $\Gamma $, despite the fact that all of the studies demonstrate good agreement with *PVT* data. This suggests that $\Gamma $ is somewhat insensitive to the *PVT* data so that more constraints are necessary to uniquely determine $\Gamma (V)$. One such constraint on $\Gamma $ comes from properties at ambient conditions, where the value of $\Gamma $ 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 $\alpha $, *B*, and *V* from the work of Kell^{39} and *C*_{V} from IAPWS-95 to find that $\Gamma =0.1489$ at ambient conditions. This is comparable to the value of $\Gamma \u2009=\u20090.11$ reported by Boehler and Kennedy,^{40} who obtained $\Gamma \u2009=\u2009(BS/T)(\u2202T/\u2202P)S$ directly by measuring the isentropic derivative $(\u2202T/\u2202P)S$ in a cylindrical piston apparatus and multiplying their result with the isentropic bulk modulus *B*_{S}. We represent $\Gamma $ as a fourth-order polynomial in *V*,

where the coefficients and the volume *V*_{ambient} are presented in Table I. Our form for $\Gamma (V)$ meets three criteria: (1) produces the expected value of $\Gamma =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.*^{25} which we have presented in Fig. 2.

$\Gamma 0$ . | $\Gamma 1$ (mol/cm^{3})
. | $\Gamma 2$ (mol^{2}/cm^{6})
. | $\Gamma 3$ (mol^{3}/cm^{9})
. | $\Gamma 4$ (mol^{4}/cm^{12})
. | V_{ambient} (cm^{3}/mol)
. |
---|---|---|---|---|---|

0.1489 | −0.2615 | $\u22126.5703\xd710\u22122$ | $\u22121.5903\xd710\u22122$ | $\u22121.1839\xd710\u22123$ | 18.073 02 |

$\Gamma 0$ . | $\Gamma 1$ (mol/cm^{3})
. | $\Gamma 2$ (mol^{2}/cm^{6})
. | $\Gamma 3$ (mol^{3}/cm^{9})
. | $\Gamma 4$ (mol^{4}/cm^{12})
. | V_{ambient} (cm^{3}/mol)
. |
---|---|---|---|---|---|

0.1489 | −0.2615 | $\u22126.5703\xd710\u22122$ | $\u22121.5903\xd710\u22122$ | $\u22121.1839\xd710\u22123$ | 18.073 02 |

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\u2032$ 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*_{3} 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^{−4} GPa = 1 bar. Since $Pisotherm(Vref)\u2009=\u20090\u2009\u2248\u200910\u22124$ GPa, *B*_{ref} is expected to be very similar to the value of *B* at ambient pressure, which Kell^{39} reports to be 2.22 GPa. Thus, the only true adjustable parameter in *P*_{isotherm} is $Bref\u2032$. Table II lists the value of $Bref\u2032$, along with those for the other parameters of our EOS.

T_{ref} (K)
. | F_{ref} (kJ/mol)
. | C_{V}/R
. | V_{ref} (cm^{3}/mol)
. | B_{ref} (GPa)
. | $Bref\u2032$ . |
---|---|---|---|---|---|

300 | 0 | 7 | 18.07384 | 2.22 | 7.6 |

T_{ref} (K)
. | F_{ref} (kJ/mol)
. | C_{V}/R
. | V_{ref} (cm^{3}/mol)
. | B_{ref} (GPa)
. | $Bref\u2032$ . |
---|---|---|---|---|---|

300 | 0 | 7 | 18.07384 | 2.22 | 7.6 |

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 $\Gamma (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 $\Gamma $ becomes zero around *V* = 18.57 cm/mol^{3} (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\u22652.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,

### D. Sound speed and Hugoniot

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+\alpha 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.*^{43} and Sanchez-Valle *et al.*,^{36} 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*_{0} and *P*_{0} represent the internal energy and pressure of the undisturbed material ahead of the shock, whose state is characterized by its volume *V*_{0} and temperature *T*_{0}. 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\u22487R$ can closely reproduce the results from that study. In addition, our calculations match well with the low-pressure measurements of Lysne.^{44} Although not shown in the figure, we have also verified that our EOS agrees with the particle speeds measured by Lysne.

## III. ICE VII

### A. Equation of state from a general pressure-volume-temperature function

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 $\Gamma (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.

### B. Pressure-volume-temperature behavior

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 $\eta $, *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.*^{19} 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.*,^{46} 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.*,^{50} Fei *et al.*,^{51} and Frank *et al.*^{52} The volumes reported by Bezacier *et al.*^{45} lie between the extremes represented by the two aforementioned sets of studies.

T_{ref} (K)
. | V_{ref} (cm^{3}/mol)
. | B_{ref} (GPa)
. | $Bref\u2032$ . | $\eta $ (GPa) . | V_{thermal} (cm^{3}/mol)
. | T_{thermal} (K)
. |
---|---|---|---|---|---|---|

300 | 13.7 | 5.8 | 9.9 | −2.2 | 14.9 | 130 |

T_{ref} (K)
. | V_{ref} (cm^{3}/mol)
. | B_{ref} (GPa)
. | $Bref\u2032$ . | $\eta $ (GPa) . | V_{thermal} (cm^{3}/mol)
. | T_{thermal} (K)
. |
---|---|---|---|---|---|---|

300 | 13.7 | 5.8 | 9.9 | −2.2 | 14.9 | 130 |

As we have mentioned in the Introduction, French and Redmer^{3} 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.*,^{51} 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 $\alpha $ is represented as

in which $\alpha 0$, $\alpha 1$, and $\xi $ are constants. A wide range of values have been published for $\alpha $ at ambient conditions. It has been reported to be $\alpha \u2009=\u20090.48\u2009\xd7\u200910\u22124$ K^{−1}, $0.6\u2009\xd7\u200910\u22124$ K^{−1}, $1.16\u2009\xd7\u200910\u22124$ K^{−1}, and $15.0\u2009\xd7\u200910\u22124$ K^{−1} by Frank *et al.*,^{52} Fei *et al.*,^{51} Bezacier *et al.*,^{45} and Sugimura *et al.*,^{19} respectively. According to our EOS, $\alpha $ at ambient conditions is about $2.45\xd710\u22124$ K^{−1}, which lies between that of Bezacier *et al.* and Sugimura *et al.* We have calculated this result by rearranging (8) to obtain $\alpha =(\u2202P/\u2202T)V/B$. The large variation in the reported ambient value of $\alpha $ 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.

### C. Free energy, melt curve, and heat capacity

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 $\Gamma $ 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)=\u222bVVrefPisothermdV\u2032$. 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.*^{20} 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.*^{37} 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\u2009=\u2009\u2211k=02FkTrefk\u2009=\u20094.1109\u2009kJ/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

F_{0} (kJ/mol)
. | F_{1} [(kJ/mol)/K]
. | F_{2} [(kJ/mol)/K^{2}]
. |
---|---|---|

−5.1461 | $4.1527\xd710\u22122$ | $\u22123.5569\xd710\u22125$ |

F_{0} (kJ/mol)
. | F_{1} [(kJ/mol)/K]
. | F_{2} [(kJ/mol)/K^{2}]
. |
---|---|---|

−5.1461 | $4.1527\xd710\u22122$ | $\u22123.5569\xd710\u22125$ |

Some properties along the melt curve are shown in Figs. 7(b)–7(d). Goncharov *et al.*^{59} have applied X-ray diffraction to estimate the density of both phases. Asahara *et al.*^{58} 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 $\Delta V$ is an important consideration for simulations of phase transitions under ramp and multiple-shock compression. The entropy change $\Delta 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.*^{58} and Ahart *et al.*^{60} 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^{3} 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.*^{51} 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 $\theta D$ through *C*_{V}, which is given by^{63}

They report the Debye temperature $\theta D$ to be 1470 K. The integral in (44) has an exact analytical expression,^{64}

where $Lin(x)\u2009=\u2009\u2211k=1\u221exk/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,^{61} 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.*^{62} 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.*^{47} 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+\alpha 2TBV$ and $(\u2202P/\u2202T)V=\alpha B$. Asahara *et al.* find *V* by inverting the Vinet fit of Sugimura *et al.*, compute $\alpha $ from (38) using the parameters in the work of Fei *et al.*,^{51} and take *B*_{S} from their own work. Values for *B* are chosen from four different studies.^{47,51,65,66}*C*_{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.^{51} 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 54*R* [25 (J/g)/K], whereas the heat capacity is less than 10*R* for all of the studies shown in Fig. 8.

### D. Sound speed and Hugoniot

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.*^{62} 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.*^{47} 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.*^{60} 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.^{3} 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^{61} 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*_{0} = 100 K and *T*_{0} = 263 K. Their study does not indicate the initial volumes *V*_{0}, but from a data table we have obtained through private communication,^{67} we have found *V*_{0} to be 19.318 cm^{3}/mol (0.932 g/cm^{3}) and 19.575 cm^{3}/mol (0.920 g/cm^{3}) for the *T*_{0} = 100 K and *T*_{0} = 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*_{0}. Instead, we obtain *E*_{0} from the five-phase tabular EOS presented in the work of Senft and Stewart.^{68} 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^{−4} GPa (1 bar) and *T* = 273 K ($0\xb0$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*_{0} 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*_{0} and *T*_{0}. Note that because French and Redmer^{3} do not present a liquid water EOS, we cannot apply the above procedure to determine their energy offset with respect to the tabular EOS.

. | V_{0} (cm^{3}/mol)
. | P_{0} (GPa)
. | E_{0} (kJ/mol)
. |
---|---|---|---|

T_{0} = 100 K | 19.318 | 10^{−4} | 6.531 |

T_{0} = 263 K | 19.575 | 10^{−4} | 10.717 |

. | V_{0} (cm^{3}/mol)
. | P_{0} (GPa)
. | E_{0} (kJ/mol)
. |
---|---|---|---|

T_{0} = 100 K | 19.318 | 10^{−4} | 6.531 |

T_{0} = 263 K | 19.575 | 10^{−4} | 10.717 |

Figure 10 compares our results for the ice VII Hugoniot with experimental data for this phase from the work of Stewart and Ahrens.^{61} 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*_{0} = 100 K data are clearly not consistent with the static measurements, but on first glance, it appears as though their *T*_{0} = 263 K data might be consistent.

If we were to try adjusting our *P*(*V*, *T*) model to fit the Stewart and Ahrens *T*_{0} = 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 [$(\u2202P/\u2202V)T$ would need to become very small] for $V\u2009>\u200911$ cm^{3}/mol. In other words, ice VII would need to undergo an unrealistic, abrupt softening. Furthermore, the very small $(\u2202P/\u2202V)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^{3}/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.*^{45}). 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,^{61} and it does not present diffraction measurements. We encourage more studies on the shock properties of ice VII to help resolve the discrepancies.

### E. Liquid Hugoniot, two-phase isentropes, and the melt curve

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 $\Delta 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.^{8} 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.^{72}

## IV. CONCLUSIONS

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.^{61} 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.

## ACKNOWLEDGMENTS

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.

## REFERENCES

_{2}O ice VII

_{2}O: Thermodynamic functions of the phase transitions of high pressure ices

_{2}O), and hydrogen (H

_{2})

_{2}O ice to 126 GPa and implications for hydrogen-bond symmetrization: Synchrotron x-ray diffraction measurements and density-functional calculations

_{2}O-VII and D

_{2}O-VII

_{2}O-ice to 128 GPa (1.28 Mbar)

_{2}O-ice VII to 20 GPa

_{2}O to 80 GPa using the melting curve, bulk modulus, and thermal expansivity of ice VII

_{2}O melting to 22 GPa and 900 K

_{2}O to 90 GPa measured in a laser-heated diamond cell

_{2}O: Another ice phase and its melting curve

_{2}O melting to 26 GPa

_{2}O ice

_{0.6}Fe

_{0.4})O

_{2}O ice VII