Thermal Conductivity of Water at Extreme Conditions

Measuring the thermal conductivity (κ) of water at extreme conditions is a challenging task, and few experimental data are available. We predict κ for temperatures and pressures relevant to the conditions of the Earth mantle, between 1,000 and 2,000 K and up to 22 GPa. We employ close to equilibrium molecular dynamics simulations and a deep neural network potential fitted to density functional theory data. We then interpret our results by computing the equation of state of water on a fine grid of points and using a simple model for κ. We find that the thermal conductivity is weakly dependent on temperature and monotonically increases with pressure with an approximate square-root behavior. In addition, we show how the increase of κ at high pressure, relative to ambient conditions, is related to the corresponding increase in the sound velocity. Although the relationships between the thermal conductivity, pressure and sound velocity established here are not rigorous, they are sufficiently accurate to allow for a robust estimate of the thermal conductivity of water in a broad range of temperatures and pressures, where experiments are still difficult to perform.


I. INTRODUCTION
Water under pressure (P) at high temperature (T) is an important constituent of the continental crust of the Earth, 1 and of the interiors of ice giants, e.g., Uranus and Neptune, 2 as well as many exoplanets. 3,4 Hence the characterization of the heat transport properties of water at extreme conditions is central to Earth and planetary sciences. 5−7 For example, understanding water heat transport may help explain the remarkably low luminosity of Uranus 8 as well as derive models for the core erosion processes in Jupiter. 9 However, the thermal conductivity (κ) of water at high P and T (HPT) is poorly known and difficult to measure.
Direct measurements at extreme conditions are challenging not only because of the reactivity of water but also for the errors that may be introduced in the experiments by convection and radiation processes. 10 Measurements for liquid water are available for P < ∼3.5 GPa and T < ∼1,000 K 11,12 and for ice up to 22 GPa, 13 below 1,000 K.
Similar to experiments, simulations of heat transport in water at HPT are challenging. Reliable empirical force fields are not available and so far first-principles molecular dynamics (FPMD) simulations based on density functional theory (DFT) have been mostly limited to structural, vibrational and electronic properties, 14−19 due to the long simulation times and large unit cells usually required to investigate transport properties. However, recently, French 20 conducted FPMD simulations of water at extreme conditions to obtain its thermal conductivity, although the heat current was approximated by fitting pair potentials. In addition, thanks to important theoretical advances, 21,22 ab initio calculations of the thermal conductivity of water using linear response and the Green−Kubo (GK) formalism 23−26 have been conducted at both ambient 21 and extreme conditions, 27 but only for pressures higher than 33 GPa. However, analytical expressions for the energy density and flux required in GK calculations are not be easily available for sophisticated DFT functionals; 28 furthermore, despite novel noise-reduction methods, 29,30 long simulation times of several hundreds of picoseconds for systems with several hundred atoms are required to obtain converged results for the thermal conductivity, making first-principles simulations a rather demanding task. Hence a computational framework avoiding the explicit calculation of the heat flux and allowing for long simulation times is desirable, to explore the thermal conductivity of water in a wide range of conditions.
Here we use the sinusoidal approach to equilibrium molecular dynamics (SAEMD) method recently proposed for fluids 31 with a deep neural network potential (DP), 32−34 allowing for long simulation times with relatively large cells. The DP interatomic potential, trained on first-principles data, can accurately describe interatomic interactions at a cost slightly higher than that of classical force fields, but much lower than FPMD. We compute the thermal conductivity of water for 1,000 < T < 2,000 K and 1.0 < ρ < 1.86 g/cm 3 , namely, at conditions relevant to the Earth mantle. At these conditions, we confirm that water is a liquid by monitoring the mean-squared displacement of atoms, and we found results consistent with the water phase diagram reported in ref 35. We then interpret our results by computing the equation of state (EOS) of water on a fine grid of points and using a simple model derived from our EOS results and the computed values of κ. We find that at the conditions studied here κ increases relative to ambient conditions, is weakly dependent on temperature and monotonically increases with pressure with an approximate square-root behavior.
The rest of the paper is organized as follows. The methods used here to compute the thermal conductivity and equation of state are described in the next section, followed by a presentation of our results and finally by our conclusions.

IIA. Thermal Conductivity Calculations.
We investigated the thermal conductivity of water at high pressure and temperature by carrying out molecular dynamics simulations with a deep neural network potential 32−34 and the LAMMPS code. 36, 37 The potential was trained with the DeepMD-kit package 32 using ice and water structures from low temperature and pressure to about 2,400 K and 50 GPa. The training data were obtained from density functional theory calculations using the strongly constrained and appropriately normed (SCAN) meta-GGA exchange−correlation functional. 38 More details can be found in ref 35.
Specifically, we used the SAEMD method, 31 which allowed us to avoid the calculation of the heat flux, and we computed the thermal conductivity of the liquid from its response to a perturbation. This perturbation is a nonhomogeneous constant temperature profile T(x, y, z), maintained by a thermostat: where L is the length of the simulation cell chosen to represent the system and ΔT is the difference between the maximum and the minimum temperature. During MD simulations we monitored how much energy the thermostat is providing to the system and computed the thermal conductivity from the solution of the heat equation: where q is the heat generation rate per unit volume from the thermostat.
We carried out eight simulations: (i) one at ambient conditions, at T = 300 K and ρ = 1 g/cm 3 ; (ii) three calculations at T = 1,000 K and ρ ∈ [1.2, 1.57, 1.86] g/cm 3 ; (iii) four calculations at T = 2,000 K and ρ ∈ [1.0, 1.2, 1.57, 1.86] g/cm 3 . We do not report calculations for ρ = 1 g/cm 3 and T = 1,000 K as it was difficult to properly converge our simulations due to the presence of large fluctuations in the heat generation rate q. We used ΔT = 10, 30, and 100 K for calculations at T = 300, 1,000 and 2,000 K, respectively. For each combination of density and temperature, we performed 20 independent runs, over which we averaged the amount of energy transferred to the system to compute the thermal conductivity. We used a cubic cell containing 512 water molecules, which was large enough to obtain approximately converged results, as previously verified. 31 For example, at T = 1,000 K and ρ = 1.57 g/cm 3 SAEMD simulations with 512 molecules yield a slight underestimation of the thermal conductivity of ∼5%, compared to the extrapolated value to infinite size.
For each independent run, we equilibrated the system for 3 × 10 5 steps, followed by production runs of 10 × 10 5 steps. We used a time step of 0.2 fs and collected data for a total of 4 ns for ρ = 1.0 and 1.2 g/cm 3 , and we used a time step of 0.25 fs and collected data for a total of 5 ns for ρ = 1.57 and 1.86 g/cm 3 .

IIB. Equation of State Calculations.
We also carried out equation of state calculations by considering 90 T−ρ conditions on an evenly spaced 9 × 10 mesh, for 1,000 < T < 2,000 K (9 grid points) and 1.0 < ρ < 1.9 g/cm 3 (10 grid points). At each T−ρ condition, we performed MD simulations using the DP potential in the NVT ensemble with a time step of 0.2 fs and a cubic cell containing 128 water molecules. For each MD simulation, we equilibrated the system for 20 ps, followed by a production run of 54 ps. In order to test finite size effects, we compared total energies and pressures obtained when using cubic cells of 128 and 512 water molecules at ρ = 1.2 g/cm 3 and T from 1,000 to 2,000 K; the relative differences in total energy are <0.1% and those in pressure are <1%, which are attributed to statistical errors. However, at T = 1,000 K and ρ = 1.8 and 1.9 g/cm 3 we found that the system did not exhibit a diffusive behavior, when using 128 water molecules in our cell. Hence we discarded the results of these simulations and we used a larger cell (512 water molecules) at T = 1,000 K and ρ = 1.86 g/cm 3 , where the system did behave as a liquid; for this simulation we used a time step of 0.2 fs and equilibrated the system for 30 ps, followed by a production run of 120 ps.
At each T−ρ condition we computed the total energy (E), the pressure (P) and the water dissociation ratio, obtained by using a cutoff distance for O−H bonds of 1.25 Å. We then interpolated E(T, ρ), P(T, ρ) and the water dissociation ratio over the whole parameter range considered here, by using the Gaussian process regression method as implemented in the sklearn package. 39 We used the radial basis function kernel with independent length scales for T and ρ. The hyper-parameters of the model were obtained by maximizing the log-marginal-likelihood.
Based on the interpolated functions, which are differentiable, we computed additional properties of the system. In particular, we computed the constant volume heat capacity per atom (C V ) as where E(T, ρ) is the total energy per atom. Further, we computed the constant pressure heat capacity per atom (C P ) 40 as ( , ) The Journal of Physical Chemistry B pubs.acs.org/JPCB Article where m is the average mass per atom. We also obtained the adiabatic index as γ(T, ρ) = C P (T, ρ)/C V (T, ρ), and computed the sound velocity C S 41 as We calculated all the properties described above on a dense 100 × 40 mesh, for 1,000 < T < 2,000 K (100 point grid) and 1.0 < ρ < 1.9 g/cm 3 (40 point grid).

IIIA. Computed Thermal Conductivity.
Our computed values of the thermal conductivity κ at extreme conditions are summarized in Table 1. We also present results at ambient conditions for comparison. In Figure 1, we show κ of water at extreme conditions as a function of the density ( Figure 1A) and pressure ( Figure 1B).
We start by comparing our results at ambient conditions to those of previous studies and experiments. The calculated value of 0.81 W/mK at 300 K and 1 g/cm 3 , agrees relatively well with that obtained via spectral analysis of the energy flux in NVE simulations with 128 water molecule cells, 28 as expected since both studies used the DP potential trained on a SCANgenerated data set. Based on the finite-size scaling study reported in ref 31 using empirical potentials, we expect our results to represent an underestimate of the data one would obtain for infinite sizes (possibly up to 15%). The overall overestimate from simulations compared to the experimental value (0.609 W/mK 42,43 ) may be due to the neglect of nuclear quantum effects and to errors introduced by the SCAN functional. We note that when using the DP model at the SCAN level of theory, the freezing temperature of water is ∼310 K. 44 At 300 K, water described by the SCAN functional is sluggish and solid-like; hence it is not surprising that the thermal conductivity at 300 K is overestimated with this functional. Nuclear quantum effects have been shown to affect several properties of water at ambient conditions. 45 Specifically, the heat capacity C V can be overestimated in classical MD simulations due to the activation of, e.g., high-frequency intramolecular motions. As suggested in recent studies based on the Green− Kubo method, 46,47 this overestimation of C V may also lead to an overestimation of the thermal conductivity, as found in our work. However, the effect of the quantum nuclear motion on the thermal conductivity computed by the SAEMD method remains to be established and will be the topic of future investigations.
We now turn to analyzing the dependence of κ on the temperature, density ( Figure 1A) and pressure ( Figure 1B). At the densities studied here, we find that the thermal conductivity increases slightly with T from 1,000 to 2,000 K. Incidentally, the κ for water at 1 g/cm 3 and 2,000 K is almost the same as that computed at ambient conditions. Consistent with experimental data at lower temperature and pressure, 12 and with high pressure studies of ice, 13 our simulations show an increase of the thermal conductivity as the density and pressure are increased. In addition, our results are consistent with the simulation reported in ref 20, where the authors found that the thermal conductivity in the T−ρ range investigated in our work is approximately independent of temperature. Remarkably, we find that a squareroot function = A P captures rather well the dependence of κ on P, at both 1,000 and 2,000 K (A is a parameter almost constant as a function of T, between 1,000 and 2,000 K). We quantify the relative error (RE) as where the fitted A (≈0.56) was used, and κ is the thermal conductivity from SAEMD simulations. We find the average RE over the 7 data points is 10.3%; the maximum RE is 21.9%. In order to interpret the temperature, pressure and density dependence of κ, particularly the relation = A P found in our simulation, we employed a simple model, described below.
IIIB. Model to Interpret Simulations. Numerous models have been proposed in the literature to describe the thermal conductivity of liquids. 48−51 Here we use a simple expression of κ encompassing several of these models: where C V is a model-dependent prefactor, k B is the Boltzmann constant, C S is the sound velocity, and δ = n −1/3 is the We also report the standard deviation error Δκ. and that the speed of energy transfer is approximately equal to the sound velocity C S ; the energy is transferred step by step, between neighboring molecules separated by a distance δ.
We extend the use of eq 7 to interpret the results of the thermal conductivity of water computed at extreme conditions. It should be noted that at HPT water may dissociate. 15,16,18,19,35,52 Hence, we first verified whether the use of eq 7, derived for simple liquids with no dissociating units, is at least approximately justified. We computed the ratio of dissociated water molecules in our samples, as shown in Figure 2A. We found that even at the highest T and ρ studied here, less than ∼15% of molecules were dissociated and therefore we expect that the dissociation of water molecules is not a major factor affecting thermal transport at the conditions considered in our work. Hence the use of the model of eq 7 appears to be reasonable to interpret our simulation results for HPT water.
Similar to previous studies, we use = m ( / )  Figure 2B,C. However, we do not have well-defined methods to compute C T ( , ) V , especially at HPT. In previous studies, C V was approximated by the specific heat per molecule (C V or C P ), 49−51 e.g., This may be a good approximation for liquids, including water, at near ambient conditions, where the major contribution to the heat capacity comes from intermolecular interactions. However, at extreme conditions and high temperature the contributions of intramolecular vibrations cannot be ignored. Therefore, we would expect a serious error in our estimate of κ if we used C V as C V /k B in eq 7. Hence here we treat C V as a parameter that we fit using the computed κ at high P and T (the 7 data points in Table 1).
As shown in Figure 3, we obtain a reasonable linear fit of κ versus k B C S /δ 2 , from which we determine C V 1.8. We quantify the RE as where the fitted C V (≈1.8) was used, and κ is the thermal conductivity from SAEMD simulations. We find that the average RE over the 7 data points is 9.5%; the maximum RE is 17.7%. The reasonable error found here indicates that water dissociation is unlikely to affect heat conduction in HPT water, in the T−ρ range considered in our work. However, while dissociation remains limited, proton conduction via Grotthus like mechanisms might play a role in determining heat transport. This aspect has not been studied in detail here and also for this reason we chose to fit the C V parameter to simulation data. To show qualitatively the difference between C V at ambient and extreme conditions, in Figure 3 we plot two lines corresponding to a value of C V equal to 3 and ≈1.8. When using 49,53 the measured value of κ at ambient conditions can be correctly predicted. The smaller C V (≈1. 8) found at HPT appears to be consistent with the presence of a disrupted hydrogen bonded network and a small fraction of dissociated water molecules at extreme conditions, leading to a decrease in the energy transfer between adjacent molecules, relative to ambient conditions. In addition, we computed the heat capacity of water at HPT conditions (see Methods), as shown in Figure 2D. We found that under the conditions considered here the heat capacity is smaller than that at ambient conditions, in line with the measurements reported in ref 54; further, at HPT, intra-  The Journal of Physical Chemistry B pubs.acs.org/JPCB Article molecular modes are activated and more energy is stored in those modes than at ambient conditions. Such energy is not expected to contribute to intermolecular energy transfer and then to heat conduction. All of these factors are expected to contribute to the decrease of the parameter C V at extreme conditions, as found in our simulations. We expect that treating C V as a function of T and ρ, instead of a fitting parameter (Figure 3), would increase the accuracy of the model (eq 7) in describing the thermal conductivity at HPT.
Using the model (eq 7) with the determined C V (≈1.8), we predicted the thermal conductivity in the whole T−ρ range. Our results are shown in Figure 4. Based on the fitting error ( Figure 3 and eq 8), the average RE of our predicted κ should be approximately 10% and the maximum RE ∼ 20%. We note that due to finite-size effects, our prediction here may also be slightly underestimated.
Finally, by substituting = m ( / ) 0 1/3 into eq 7, we obtain the relation κ ∝ ρ 2/3 C S (T, ρ). As shown in Figure 2C, C S increases slightly with T but significantly with ρ, which, according to the model (eq 7), leads to the same dependence found in our simulations for κ.
We note that C S is related to the derivative of P (see Methods); hence an analytical formula for P(T, ρ) is desirable to derive the relation between C S and P. To this end, we fit our interpolated function P(T, ρ) using the Benedict−Webb−Rubin (BWR) equation: 55,56 where m 0 is the mass of a water molecule and B 0 , B 1 , C 0 , C 1 , and D are fitting parameters. We have ignored the exponential term and terms higher than (1/T) 0 for simplicity. Our fitting data, denoted as P-DP, are evenly spaced over a 100 × 40 mesh; at each grid point T−ρ, a value of P is obtained from the interpolated function P(T, ρ). We optimized the parameters entering the BWR equation, and we show the computed pressure (P-BWR) in Figure 5, as well as the respective contributions. Interestingly, the BWR equation accurately describes the interpolated function P (T, ρ), with a small rootmean-square-error of ∼0.07 GPa. We find that the third term is dominant; the contributions of first, second and sixth terms are smaller than that of the cubic one. We note an approximate cancellation between the sum of the positive first and sixth terms and the negative second term; the third term alone is of similar magnitude to the total pressure (P-BWR). Hence, for simplicity, we assume P(T, ρ) ∼ C(T)ρ 3 , where C(T) refers to C 0 T + C 1 in eq 9.
Knowing the dependence of the sound velocity on pressure and a form of the pressure as a function of temperature and density, we can now obtain an approximate dependence of κ on the pressure:~× where γ is the adiabatic index (see Methods). In the T−ρ range studied here, we find that and ρ 1/6 are in the range of ∼(1.0, 1.1), i.e., nearly constant; as a result, we obtain P . Although the square-root relation P found here is not rigorously proven, it is a simple and useful functional relationship to approximately predict κ when P is measured in the range investigated in our work.

IV. CONCLUSIONS
By carrying out SAEMD simulations with the DP potential, we computed the thermal conductivity of water at high temperatures, 1,000 < T < 2,000 K and 1.0 < ρ < 1.86 g/cm 3 , at conditions relevant to the Earth mantle. We found that the thermal conductivity depends weakly on the temperature and increases monotonically with the density and pressure, reaching values approximately 4 times larger than that at ambient conditions at the highest density point, indicating a more efficient heat energy transport under pressure than at ambient conditions. We showed that a simple model (eq 7) can satisfactorily describe the thermal conductivity of water at extreme conditions, and using such a model, we provided predictions of the thermal conductivity in a broad range of density and temperature. Our results indicate that the heat is transferred roughly at the speed of sound over nearest-neighbor intermolecular distances and that the heat conduction mechanism is not significantly affected by water dissociation, when the proportion of dissociated molecules remains smaller than 15%. Our simulations and the model used here to interpret them indicate that an increased sound velocity and density at extreme conditions are responsible for a larger thermal conductivity in HPT water than at ambient conditions. Numerically, we identified a square-root relationship between the thermal conductivity and the pressure of the system.  We show values for several T and ρ conditions. The contributions from virial terms (first, second, third, and sixth) to P-BWR are also shown. We plot a dashed line, y = x as a guide to the eye, showing that P-BWR and P-DP values are close to each other. The P-BWR and P-DP are obtained on a dense grid, and the parameters in the BWR model are optimized (see text).
The Journal of Physical Chemistry B pubs.acs.org/JPCB Article Although this relationship is not rigorous, it can be useful to estimate the thermal conductivity at T−ρ conditions similar to those studied here, since its direct measurement may be more difficult than that of the pressure. Our study provides both insights and useful data on transport properties and equation of states of water at high temperature and pressure, which may be useful in planetary and geosciences.