Thermal and Tidal Evolution of Uranus with a Growing Frozen Core

We used Car–Parrinello AIMD simulations to determine the physical properties governing thermal evolution and reexamine the viscosity of the superionic phase of H₂O. All of the AIMD simulations were performed using computer codes taken from the QUANTUM ESPRESSO package v. 6.1 (Giannozzi et al. 2009, 2017, 2020) and adopting the Perdew–Burke–Ernzerhof density functional and Hamann–Schluter–Chiang, as modified by Vanderbilt (Vanderbilt 1985) norm-conserving pseudopotentials (downloadable from http://fpmd.ucdavis.edu/potentials/index.htm), which are also accurate at planetary pressure–temperature conditions (Sun et al. 2015).

For the thermal conductivity, we use the Green–Kubo theory of linear response, leveraging recently discovered invariance principles in the theory and numerical computation of transport coefficients (Marcolongo et al. 2016; Grasselli & Baroni 2019) in multicomponent systems (Baroni et al. 2018; Bertossa et al. 2019) and the resulting cepstral analysis for the time series of the fluxes to be fed into the Green–Kubo formula (Ercole et al. 2017; Bertossa et al. 2019). The results are collected in Grasselli et al. (2020), where the reader can find a thorough description of the methods employed.

To obtain the isobaric specific heat capacity, C_P , and thermal expansion coefficient, α, as well as the isothermal compressibility, K_T , we ran AIMD simulations at a fixed number of particles N, pressure P, and temperature T (NPT). Simulation parameters were chosen as in Grasselli et al. (2020).

The isobaric specific heat and thermal expansion coefficient, as well as the isothermal bulk modulus, can be extracted in a single NPT simulation from fluctuations as

$$\begin{eqnarray}&&{C}_{P}=\displaystyle \frac{1}{{k}_{B}{T}^{2}M}{\left\langle {\rm{\Delta }}{H}^{2}\right\rangle }_{\mathrm{NPT}},\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\alpha =\displaystyle \frac{1}{{k}_{B}{T}^{2}}\displaystyle \frac{{\left\langle {\rm{\Delta }}H{\rm{\Delta }}V\right\rangle }_{\mathrm{NPT}}}{{\left\langle V\right\rangle }_{\mathrm{NPT}}},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{K}_{T}^{-1}=\displaystyle \frac{1}{{k}_{B}T}\displaystyle \frac{{\left\langle {\rm{\Delta }}{V}^{2}\right\rangle }_{\mathrm{NPT}}}{{\left\langle V\right\rangle }_{\mathrm{NPT}}},\end{eqnarray} \tag{ 3 }$$

where H is the enthalpy and M is the total mass of the system. We employed standard block analysis and error propagation to estimate the uncertainties on C_P , α, and K_T , which are acceptably small for the relatively long NPT simulations we ran (≥20 ps). By employing two NPT simulations at different temperatures T₁ and T₂ but the same P, a finite-difference method was also devised to obtain estimates of C_P and α with a lower statistical uncertainty:

$$\begin{eqnarray}&&{C}_{P}=\displaystyle \frac{1}{M}\displaystyle \frac{\left\langle {H}_{2}\right\rangle -\left\langle {H}_{1}\right\rangle }{{T}_{2}-{T}_{1}},\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&\alpha =\displaystyle \frac{2}{\left\langle {V}_{2}\right\rangle +\left\langle {V}_{1}\right\rangle }\displaystyle \frac{\left\langle {V}_{2}\right\rangle -\left\langle {V}_{1}\right\rangle }{{T}_{2}-{T}_{1}}.\end{eqnarray} \tag{ 5 }$$

We checked the convergence of the results on the parameters of the thermobarostat. The results are reported in Table 1, in which we also report the Grüneisen parameter $\gamma =V{\left(\partial P/\partial E\right)}_{V}=\alpha {K}_{T}/(\rho {C}_{P}-{\alpha }^{2}{K}_{T}T)$. Our values are in good agreement with previous AIMD simulations; from the internal energy (u) results reported in Table 5 of French et al. (2009), we find C_V = Δu/ΔT = 4.78 J (gK)^–1 at ρ = 2.5 g cm⁻³ and an average T = 2500 K, and in Supplementary Table 5 of Millot et al. (2018), we find C_V = 3.98 J (gK)^–1 at ρ = 3.753 g cm^{– 3} and an average T = 3500 K (using their values of energy at 3000 and 4000 K).

We examined the viscosity of H₂O by performing AIMD simulations at fixed N, volume V, and energy E (NVE) in liquid and superionic states in hydrostatic and sheared configurations.

The time integral of the autocorrelation function of the shear stress

$$\begin{eqnarray}&&{\eta }_{\mathrm{ACF}}=\displaystyle \frac{V}{{k}_{B}T}{\int }_{0}^{\infty }\left\langle {P}_{{xz}}(t){P}_{{xz}}(0)\right\rangle {dt}\end{eqnarray} \tag{ 6 }$$

does not represent the proportionality coefficient between the stress and the strain rate in viscoelastic materials, like crystals or superionic materials, nor does it dictate the convective motion of large masses of solid, crystalline material, which is instead governed by defect diffusion.

We ran two NVE simulations to extract η_ACF from Equation (6), one at 2500 K and 175 GPa, where the system is superionic, and one at approximately the same average temperature but a much lower pressure (50GPa), where the oxygen lattice is molten and the system is in the ionic liquid phase. We obtained η_ACF = 1.50 ± 0.08 mPa s for the superionic phase and η_ACF = 1.89 ± 0.11 mPa s for the liquid phase. Even if these values are consistent with those given in Redmer et al. (2011), we do not agree with the conclusion drawn there that "the superionic phase responds almost like a fluid." The similarity of η_ACF for the superionic and ionic liquid phases does not indicate their comparable behavior under external shear stress; if it were so, the stress induced by the instantaneous application of an external shear strain to the superionic ice system would decay to zero on a timescale related to the time τ ≈ 100 fs over which the autocorrelation function of the shear stress is extinguished. We verified that this is not the case by deforming the simulation cell so as to include nondiagonal components in the cell tensor (in particular, we introduced nonvanishing xz and zx components) and monitoring the behavior of the induced shear stress. On the timescale of τ, we observe a decay of the latter to a new equilibrium value (Figure 2). No approach to zero of P_xz was observed for the superionic phase, not even after several ps of simulation. On the contrary, in the case of ionic liquid water, P_xz vanishes at equilibrium. Therefore, in contrast to Redmer et al. (2011), we conclude that the superionic phase does not respond to an external shear stress like a fluid, and that it has a nonnegligible viscosity. We note that another recent paper also concluded that the viscosity of the superionic phase may be large compared with that of a fluid (Millot et al. 2019).

Zoom In Zoom Out Reset image size

In the presence of a growing, viscous, frozen core, thermal evolution is governed by the equations (Fortney & Nettelmann 2010; Cottaar & Buffett 2012)

$$\begin{eqnarray}&&{L}_{\mathrm{int}}={L}_{\mathrm{fluid}}+{L}_{\mathrm{core}}=4\pi {R}^{2}\sigma \left({T}_{\mathrm{eff}}^{4}-{T}_{\mathrm{eq}}^{4}\right),\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{L}_{\mathrm{fluid}}=-{\int }_{{M}_{c}}^{M}{C}_{P}\displaystyle \frac{\partial T}{\partial t}{dm},\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{L}_{\mathrm{core}}=-{\int }_{0}^{{M}_{c}}{C}_{P}\displaystyle \frac{\partial T}{\partial t}{dm}-{\int }_{S}{C}_{P}{\rm{\Delta }}T\displaystyle \frac{\partial c}{\partial t}{\rho }_{c}{dS},\end{eqnarray} \tag{ 9 }$$

where L_int is the total luminosity from the interior; L_fluid and ${L}_{\mathrm{core}}$ are the contributions from the fluid envelope and frozen core, respectively; R is the radius of the planet; σ is the Stefan–Boltzmann constant; T_eq is the luminosity due to thermalized and reradiated solar flux; T_eff is the temperature the planet would have in the absence of solar luminosity; M is the mass of the planet; M_c is the mass of the frozen core; C_P is the isobaric specific heat; c is the radius of the frozen core; S is the surface of the core; ρ_c is the density at the core radius; ΔT is the temperature contrast across the thermal boundary layer at the top of the core; ∂T/∂t is the cooling rate; and ∂c/∂t is the growth rate of the core.

The system of equations differs from the standard thermal evolution model for giant planets in the appearance of a thermal boundary layer at the top of the core (second term on the right-hand side of Equation (9)). As ∂c/∂t and ΔT are positive quantities, this equation shows that the presence of a thermal boundary layer reduces ${L}_{\mathrm{core}}$, as compared with a uniformly fluid planet, thereby storing heat in the deep interior and allowing the fluid envelope to cool more quickly. In the case of a purely fluid planet, ΔT → 0, the interior heat flux reduces to that of the adiabatic, inviscid case (Fortney & Nettelmann 2010): ${L}_{\mathrm{int}}=-{\int }_{0}^{M}{C}_{P}\tfrac{\partial T}{\partial t}{dm}$.

The radius of the core, c, is determined at any time by the intersection of the planetary temperature profile with the freezing curve of H₂O (Figure 1). The rate of growth is then (Gubbins et al. 2003)

$$\begin{eqnarray}&&\displaystyle \frac{\partial c}{\partial t}=-\displaystyle \frac{1}{\rho g\left({\rm{\Gamma }}-{{\rm{\Gamma }}}_{a}\right)}\displaystyle \frac{\partial {T}_{c}}{\partial t},\end{eqnarray} \tag{ 10 }$$

where Γ is the Clausius–Clapeyron slope of the phase boundary, and Γ_a = (T/P)∇ is the slope of the planetary temperature distribution evaluated at T_c and P_c , the temperature and pressure at the top of the core. We represent the phase boundary with the Simon equation

$$\begin{eqnarray}&&{T}_{m}={T}_{0}{\left(1+\displaystyle \frac{P-{P}_{0}}{a}\right)}^{b},\end{eqnarray} \tag{ 11 }$$

with the values of the parameters chosen to match the melting transition of Millot et al. (2018) and the temperature in the fluid envelope and solid core by

$$\begin{eqnarray}T(P,t)=\left\{\begin{array}{ll}{T}_{1}(t){\left(\displaystyle \frac{P}{{P}_{1}}\right)}^{{\rm{\nabla }}}, & P\lt {P}_{c}\\ {T}_{i}(t){\left(\displaystyle \frac{P}{{P}_{c}}\right)}^{{\rm{\nabla }}}, & P\geqslant {P}_{c}\end{array}\right.,\end{eqnarray} \tag{ 12 }$$

respectively, where T_i = T_c + ΔT, P_c is the pressure at the top of the core, P₁ = 1 bar, and the relationship between T₁ and T_eff is given by Guillot et al. (1995),

$$\begin{eqnarray}&&{T}_{1}={{BT}}_{\mathrm{eff}}^{1.244},\end{eqnarray} \tag{ 13 }$$

with B = 0.47529 reproducing the present-day observed value of T₁ = 76 K for T_eff = 59.1. We assume that the adiabatic gradient ∇ (away from thermal boundary layers) and the heat capacity are homogeneous. The assumed form of the temperature structure with homogeneous ∇ matches very well with multilayer fluid models in which the adiabatic gradient is given at every point by an assumed composition and equations of state (EOSs) of rock, ice, and gas components (Figure 3).

Zoom In Zoom Out Reset image size

The thermal evolution equations are closed by relating the core luminosity to the temperature contrast across the thermal boundary layer (Schubert et al. 1969),

$$\begin{eqnarray}&&{L}_{\mathrm{core}}=k{\left(\displaystyle \frac{\rho \alpha g}{{\mathrm{Ra}}_{{\rm{cr}}}\eta \kappa }\right)}^{1/3}{\left(\displaystyle \frac{h}{c}\right)}^{4/3}{\rm{\Delta }}{T}^{4/3},\end{eqnarray} \tag{ 14 }$$

which shows that the core luminosity increases with the thermal conductivity k and the temperature contrast. The other parameters are thermal expansivity α, thermal diffusivity κ, gravitational acceleration g, and viscosity η. The viscous scale height $h=-{\left(\partial \mathrm{ln}\eta /\partial r\right)}^{-1}$ and other quantities in Equation (14) are evaluated at the surface of the core with radius c. Here Ra_cr is the critical Rayleigh number for the onset of convection. This expression is appropriate for very deep convective systems, such as the core of Uranus, over which the viscosity varies by many orders of magnitude. We adopt the value Ra_cr = 30, corresponding to the limit h/c → 0 (Schubert et al. 1969).

For the viscosity, we adopt the homologous temperature relationship that has been widely used in studies of the thermal evolution of icy moons and is supported by experimental data (Moore 2006),

$$\begin{eqnarray}&&\eta ={\eta }_{0}\exp \left[A\left(\displaystyle \frac{{T}_{m}}{T}-1\right)\right],\end{eqnarray} \tag{ 15 }$$

which yields

$$\begin{eqnarray}&&h=\displaystyle \frac{{T}_{c}}{\rho {gA}\left({\rm{\Gamma }}-{{\rm{\Gamma }}}_{a}\right)}\end{eqnarray} \tag{ 16 }$$

for the viscous scale height. The value of the viscosity exponent A is set by experimental parameters summarized by Moore (2006). We note that our expression assumes that viscous flow is dominated by volume diffusion, and that the creep mechanism in the deep interior of ice giants could be different from that in moons. The inner envelope of Uranus is unlikely to be composed of pure H₂O (although H₂O may be the dominant component), and the presence of impurities may influence the crystal structure and melting temperature and therefore transport properties. A study performed on an impure H₂O system at the relevant conditions found that a 1:1 mix of water and ammonia had nearly the same melting point as pure H₂O (Bethkenhagen et al. 2015). While this result indicates that the melting temperature of pure H₂O may be a reasonable approximation to that of the ice component of giant planets, melting in H–C–N–O systems at multimegabar pressures is still little explored. We find little difference in transport properties between the two crystalline structures that have been found for superionic H₂O (Millot et al. 2019; Grasselli et al. 2020). Another recent study also finds little difference in the self-diffusion coefficient among different crystalline structures of superionic H₂O (Cheng et al. 2021).

We solve Equations (7)–(16) with a fourth-order Runge–Kutta scheme and values of physical parameters as specified in Table 2. As in the case of inviscid thermal evolution models (Guillot et al. 1995), our results are insensitive to initial conditions. We ignore the relatively minor contributions to luminosity from gravitational contraction (Hubbard et al. 1995). For planetary structure, we adopt the so-called "empirical" model of Helled et al. (2011). This model is specified by a sixth-order polynomial of density in radius and matches all relevant observational data. We compute other structural quantities from this polynomial, including pressure, gravitational acceleration, and mass as a function of radius. We have checked that our results are not sensitive to the details of planetary structure by using a very different empirical model that includes a dense core to 20% of the planet's radius and fits relevant observational data equally well (Kaspi et al. 2013). We find that computed cooling times differ by 10% between the two models.

We do not explicitly account for compositional layering in Uranus; we assume that over the range of depth in which the core forms in our models (out to two-thirds of the planet's radius), the composition is dominated by H₂O, consistent with detailed interior models (Nettelmann et al. 2013). We neglect the contrast in physical properties between superionic ice and a rocky innermost core because these have a small effect on thermal evolution and tidal response. Typical models of the interior of Uranus contain a rocky core that makes up less than 6% of the mass of the planet, corresponding to less than 20% of the radius (Scheibe et al. 2019; Nettelmann et al. 2013). Assigning to this innermost region a melting temperature as much as twice that of H₂O (and thus a higher viscosity according to Equation (15); Mazevet et al. 2019) and a shear modulus as much as 10 times that of H₂O changes the cooling time by less than 2% and Q by less than 3%.

We neglect the influence of the release of latent heat upon freezing of the core because this is poorly constrained and has a small influence on thermal evolution. Including the effects of latent heat, the energy balance (Equation (7)) becomes ${L}_{\mathrm{int}}={L}_{\mathrm{fluid}}+{L}_{\mathrm{core}}+{L}_{\mathrm{latent}}$, with L_latent = 4π c² ρ_c T_c Δs(∂c/∂t). In order to estimate the value of the specific entropy of melting Δs, we assume a typical high-pressure value for the entropy of melting (${k}_{{\rm{B}}}\mathrm{ln}2$ atom^–1; Stishov et al. 1973), where k_B is the Boltzmann constant, and recall that only the oxygen sublattice melts at the superionic-to-fluid transition, yielding ${\rm{\Delta }}s={k}_{{\rm{B}}}\mathrm{ln}2w$, where w is the molecular mass of H₂O. Adopting this value of Δs and including the latent heat term in our numerical calculations, we find that the cooling time of our fiducial case (Figure 2) is increased by 6%.

We characterize tidal dissipation by the quality factor

$$\begin{eqnarray}&&Q=\displaystyle \frac{\left|{k}_{2}\right|}{\mathrm{Im}({k}_{2})},\end{eqnarray} \tag{ 17 }$$

where k₂ is the (complex) tidal Love number. We compute the tidal Love number from the radial structure of the planet at each time step of our thermal evolution calculations by solving the set of six coupled differential equations that govern the viscoelastic response (Alterman et al. 1959; Sabadini & Vermeersen 2004),

$$\begin{eqnarray}&&{\dot{y}}_{1}=-\displaystyle \frac{2}{r}{y}_{1}+\displaystyle \frac{l(l+1)}{r}{y}_{2},\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&{\dot{y}}_{2}=-\displaystyle \frac{1}{r}{y}_{1}+\displaystyle \frac{1}{r}{y}_{2}+\displaystyle \frac{1}{\mu }{y}_{4},\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}\begin{array}{rcl}{\dot{y}}_{3} & = & \displaystyle \frac{4}{r}\left(\displaystyle \frac{3\mu }{r}-\rho g\right){y}_{1}-\displaystyle \frac{l(l+1)}{r}\left(\displaystyle \frac{6\mu }{r}-\rho g\right){y}_{2}\\ & & +\displaystyle \frac{l(l+1)}{r}{y}_{4}-\displaystyle \frac{\rho (l+1)}{r}{y}_{5}+\rho {y}_{6},\end{array}\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}\begin{array}{rcl}{\dot{y}}_{4} & = & -\displaystyle \frac{1}{r}\left(\displaystyle \frac{6\mu }{r}-\rho g\right){y}_{1}+\displaystyle \frac{2(2{l}^{2}+2l-1)\mu }{{r}^{2}}{y}_{2}\\ & & -\displaystyle \frac{1}{r}{y}_{3}-\displaystyle \frac{3}{r}{y}_{4}+\displaystyle \frac{\rho }{r}{y}_{5},\end{array}\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&{\dot{y}}_{5}=-4\pi G\rho {y}_{1}-\displaystyle \frac{l+1}{r}{y}_{5}+{y}_{6},\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&{\dot{y}}_{6}=-\displaystyle \frac{4\pi G\rho (l+1)}{r}{y}_{1}+\displaystyle \frac{4\pi G\rho l(l+1)}{r}{y}_{2}+\displaystyle \frac{l-1}{r}{y}_{6},\end{eqnarray} \tag{ 23 }$$

where l is the angular order of the deformation (l = 2 for tidal deformation), r is the radius, μ is the (complex) shear modulus, g is the gravitational acceleration, and G is the universal gravitational constant. The functions y_i , i = 1,6 are, respectively, the amplitude of the radial and tangential deformation, the radial and shear stresses, and the gravitational potential and its radial derivative. The Love number is

$$\begin{eqnarray}&&{k}_{l}=-{y}_{5}-{\left(r/R\right)}^{l}.\end{eqnarray} \tag{ 24 }$$

Fluid layers require special treatment because they do not transmit deformation or stress and couple to other layers only through the gravitational potential (Jara-Orue & Vermeersen 2011; Matsuyama et al. 2018). We have found that we can obtain accurate solutions through fluid layers, and indeed for entirely fluid planets, by the simple expedient of assigning a small but nonvanishing value of the shear modulus to fluid layers (μ = 10⁻⁴ ρ gr). We use the propagator matrix technique (Henning & Hurford 2014; Sabadini & Vermeersen 2004) to solve Equations (18)–(23).

We validate our code for planetary structures consisting of two homogeneous layers for which the analytical solution (Remus et al. 2015) and numerical solutions (Padovan et al. 2018) are available and for a suite of purely fluid Jupiter models (Gavrilov et al. 1976; Figures 4 and 5).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We relate the viscosity to tidal dissipation by approximating the core as a standard linear solid. The complex shear modulus (Nowick & Berry 1972)

$$\begin{eqnarray}&&\mu ={\mu }_{0}-\displaystyle \frac{\delta \mu }{1+i\omega {\tau }_{\epsilon }},\end{eqnarray} \tag{ 25 }$$

where

$$\begin{eqnarray}&&\delta \mu ={\mu }_{0}\left(1-\displaystyle \frac{{\mu }_{1}}{{\mu }_{0}+{\mu }_{1}}\right),\end{eqnarray} \tag{ 26 }$$

and the two characteristic timescales are

$$\begin{eqnarray}&&{\tau }_{\epsilon }=\displaystyle \frac{\eta }{{\mu }_{0}+{\mu }_{1}}\end{eqnarray} \tag{ 27 }$$

and

$$\begin{eqnarray}&&{\tau }_{\sigma }=\displaystyle \frac{\eta }{{\mu }_{1}}\end{eqnarray} \tag{ 28 }$$

and the quality factor

$$\begin{eqnarray}&&Q=\displaystyle \frac{1+{\omega }^{2}{\tau }_{\epsilon }{\tau }_{\sigma }}{\omega ({\tau }_{\sigma }-{\tau }_{\epsilon })},\end{eqnarray} \tag{ 29 }$$

where the viscosity is given by Equation (15). The unrelaxed shear modulus μ₀ is taken from AIMD simulations (Hernandez & Caracas 2016); we used the published values of the full elastic constant tensor to compute the Voigt–Reuss–Hill average shear modulus (Watt et al. 1976) and approximated the pressure and temperature dependence by

$$\begin{eqnarray}&&{\mu }_{0}(P,T)=101+P/2.41-(T-1650)/23.4\end{eqnarray} \tag{ 30 }$$

with pressure in GPa and temperature in K. We adopt a value of μ₁ = 60μ₀, similar to a recent study of Earth tides (Lau et al.2015).

We have chosen the standard linear solid as a balance between parametric simplicity and realistic description of the relevant physics. The standard linear solid is a rough approximation to the behavior of real materials, i.e., polycrystals with a variety of grain sizes, in part because of the predominance of a single relaxation time, which causes Q to be overly sensitive to frequency. The Maxwell model, to which the standard linear solid reduces in the limit μ₁/μ₀ → 0, is even simpler but shows unrealistically large attenuation at low frequency (Henning et al. 2009). Other more complex models have also been considered in the planetary literature, including the Burgers model (Henning et al. 2009) and the Andrade model (Castillo-Rogez et al. 2011), which agree better with experimental data where they are available but contain additional parameters that are unconstrained at present at the pressure–temperature conditions of the interior of Uranus.

We determine the evolution of the radii of satellite orbits from the values of k₂ and Q of Uranus determined at each time step of our thermal evolution model by numerically integrating (Tittemore & Wisdom 1989)

$$\begin{eqnarray}&&\displaystyle \frac{\dot{a}}{a}=3\displaystyle \frac{\left|{k}_{2}\right|}{Q}n\displaystyle \frac{m}{M}{\left(\displaystyle \frac{R}{a}\right)}^{5}\end{eqnarray} \tag{ 31 }$$

backward in time from present-day observed values (Jacobson 2014). Here a is the radius of the orbit, m is the mass of the moon, M is the mass of the planet, time-varying values of k₂ and Q are from our thermal evolution model, $n=\sqrt{{GM}/{a}^{3}}$ is the mean motion of the moon, and $\dot{a}={da}/{dt}$.