TY - JOUR
AU - Deschamps,, Frédéric
AB - SUMMARY Convection is an efficient process to release heat from planetary interiors, but its efficiency depends on the detailed properties of planetary mantles and materials. A property whose impact has not yet been studied extensively is the temperature dependence of thermal conductivity. Because thermal conductivity controls heat fluxes, its variations with temperature may alter heat transfer. Here, I assess qualitatively and quantitatively the influence of temperature-dependent thermal conductivity on stagnant lid convection. Assuming that thermal conductivity varies as the inverse of temperature |$(k \propto 1/T)$|⁠, which is the case for ice Ih, the main component of outer shells of solar System large icy bodies, I performed numerical simulations of convection in 3-D-Cartesian geometry with top-to-bottom viscosity and conductivity ratios in the ranges 105 ≤ Δη ≤ 108 and 1 ≤ Rk ≤ 10, respectively. These simulations indicate that with increasing Rk, and for given values of the Rayleigh number and Δη, heat flux is reduced by a factor Rk0.82, while the stagnant lid is thickening. These results have implications for the structures and thermal evolutions of large icy bodies, the impact of temperature-dependent conductivity being more important with decreasing surface temperature, Tsurf. The heat fluxes and thermal evolutions obtained with temperature-dependent conductivity are comparable to those obtained with constant conductivity, provided that the conductivity is fixed to its value at the bottom or in the interior of the ice shell, that is, around 2.0–3.0 W m−1 K−1, depending on the body. By contrast, temperature-dependent conductivity leads to thicker stagnant lids, by about a factor 1.6–1.8 at Pluto (Tsurf = 40 K) and a factor 1.2–1.4 at Europa (Tsurf = 100 K), and smaller interior temperatures. Overall, temperature-dependent thermal conductivity therefore provides more accurate descriptions of the thermal evolutions of icy bodies. Composition of the planets, Numerical modelling, Planetary interiors, Dynamics: convection currents, and mantle plumes 1 INTRODUCTION Like rocky planets, icy moons and dwarf planets are cooling down, and their radial structure is constrained by their thermal histories. These evolutions are themselves controlled by the mode and efficiency of heat transfer through the outer ice Ih shells of icy bodies. Several mechanisms may influence, oppose or delay this cooling, including a regular and sustainable production of heat, the reduction of the melting temperature at the bottom of the outer ice shell and the limitation of heat transfer through this shell, for instance related to its physical and material properties. Tidal dissipation may, in some cases, provide a regular source of heat within or at the bottom of the ice shell. The amount of heat released, and thus the evolution of the body, depends on the orbital properties of icy bodies and may vary with time (e.g. Tobie et al. 2003, 2005; Roberts & Nimmo 2008). The presence of antifreeze compounds, for example, ammonia or methanol, in the subsurface ocean has long been advocated to inhibit crystallization (e.g. Grasset & Sotin 1996; Spohn & Schubert 2003; Deschamps et al. 2010a). Other properties that may impact the evolution of icy bodies include the presence of a thin layer of clathrate hydrates acting as an insulator either at the top (Tobie et al. 2006) or at the bottom (Kamata et al. 2019) of the ice shell. This process has been advocated to explain the presence of a subsurface ocean within Titan and Pluto. The properties of outer ice shells of icy moons and dwarf planets allow thermal convection to operate or have operated within them. Compared to thermal conduction, convection is an efficient mechanism to release heat to the surface, but how efficient depends on the detailed mode of convection, which depends itself on the physical and material properties of the system. A side effect of both tidal heating and antifreeze compounds is precisely to alter heat transfer. If tidal dissipation occurs within the ice shell, the amount of heat that can be extracted from the ocean and transported to the surface is reduced, as shown by simulations of convection in mixed heated systems (e.g. Travis & Olson 1994; Deschamps et al. 2010b). Because it lowers the water liquidus, and thus the temperature at the bottom of the ice shell, the presence of antifreeze compounds in the ocean increases the bulk viscosity of the ice shell, reducing the vigour of convection (Deschamps & Sotin 2001). Another key factor impacting heat transfer is rheology. The viscosity of ice strongly depends on temperature, such that convection within ice shells may operate in the so-called stagnant lid regime. In this mode of convection, a rigid, thermally conductive lid forms at the top of the system, reducing the amount of heat that can be transported to the surface (e.g. Christensen 1984; Moresi & Solomatov 1995). The temperature dependence of thermal conductivity may further influence convection and its ability to transport heat towards the surface. The thermal conductivity of water ice depends on temperature following a Debye model (Slack 1980; Andersson & Suga 1994), implying that above the Debye temperature, which, for ice Ih, is equal to 226 K, it is well approximated by an inverse law of the temperature, T. Interestingly, this approximation provides a good description of ice Ih conductivity well below the Debye temperature, from 273 K down to about 35 K (Slack 1980; Andersson & Suga 1994). For 40 K ≤ T ≤ 260 K the measurements of Andersson & Suga (1994) fit very well along k = 566.8/T. In this range of temperature, which encompasses the expected temperature conditions within icy bodies, thermal conductivity increases by about a factor of 15. By contrast, up to about 0.21 GPa, the thermal conductivity of ice Ih depends only slightly on pressure (Andersson & Suga 1994). The amount of heat that can be extracted from the subsurface ocean is controlled by the thermal conductivity at the bottom of the ice Ih shell, which, if conductivity varies with temperature, should be lower than at the surface. This, in turn, may impact models of thermal evolution of icy moons. To the first order, heat flux limitation can be accounted for by prescribing a low conductivity throughout the system. However, fixing conductivity to its bottom or interior value does not capture all details of the flow, which may, in fine, alter the structure and thermal evolution models. In particular, Tobie et al. (2003) suggested that, in the case of stagnant lid convection, the conductive lid may be thicker if thermal conductivity depends on temperature. So far, the effects of temperature-dependent thermal conductivity have not been explored in details. Here, I perform numerical simulations of convection combining temperature-dependent viscosity (stagnant lid convection) and thermal conductivity. These simulations show that heat transfer in outer ice shells of icy bodies is affected by changes in conductivity with temperature, and that these changes should be accounted for to refine the thermal evolutions and structures of these bodies. 2 NUMERICAL MODEL To model the dynamics of outer ice shells, I performed numerical simulations of thermal convection in 3-D-Cartesian geometry for an incompressible, infinite Prandtl number fluid using StagYY (Tackley 2008). Both the fluid viscosity, η, and thermal conductivity, k, vary with temperature. The conservation equations of momentum, mass and energy are then $$\begin{eqnarray} {\boldsymbol \nabla} \overline{\overline \sigma } - {\boldsymbol \nabla} P = - \alpha \rho g\delta T{{\boldsymbol{e}}_{\boldsymbol{z}}} \end{eqnarray}$$(1) $$\begin{eqnarray} {\boldsymbol \nabla} \cdot {\boldsymbol{v}} = 0 \end{eqnarray}$$(2) and $$\begin{eqnarray} \rho {C_P}\frac{{\partial T}}{{\partial t}} = {\boldsymbol \nabla} \cdot \left( {k{\boldsymbol \nabla}T} \right) - \rho {C_P}{\boldsymbol{v}} \cdot {\boldsymbol \nabla} T + \rho H, \end{eqnarray}$$(3) where the elements of the deviatoric stress tensor, |$\overline{\overline \sigma } $|⁠, are |${\sigma _{ij}} = \eta ( {\partial {v_i}/{x_j} + \partial {v_j}/{x_i}} )$|⁠, P is the non-hydrostatic pressure, v the velocity, δT the temperature difference between the local temperature T and a reference temperature (here, the surface temperature Tsurf), ez the radial unit vector, α, ρ and CP the fluid thermal expansion, density and heat capacity (all assumed constant throughout the system), g the gravity acceleration, and H the internal heat production per unit of mass, which is here set to zero (i.e. the fluid is only heated from its bottom and cooled at its top). Numerical methods used to solve eqs (1)–(3) are detailed in Tackley (2008). Conservation equations are non-dimensionalized with the characteristic properties of the system, in particular the thickness of the fluid layer D and the super-adiabatic temperature jump across this layer, ΔT. The non-dimensional temperature and velocity are then given by |$\tilde{T} = ( {T - {T_{\mathrm{ surf}}}} )/{\rm{\Delta }}T$| and |${\boldsymbol{\tilde{v}}} = {\boldsymbol{v}}\kappa /D$|⁠, where Tsurf is the surface temperature and |$\kappa = k/\rho {C_P}$| is the fluid thermal diffusivity (for convenience, non-dimensional quantities are denoted by adding a tilde, ∼, to the corresponding symbols denoting dimensional quantities). In addition, the characteristic thermal conductivity is defined as its surface value, ksurf. In addition, the source term of the momentum equation, αρgδT, is replaced by the Rayleigh number, $$\begin{eqnarray} Ra = \frac{{\alpha \rho g{\rm{\Delta }}T{D^3}}}{{\eta \kappa }}, \end{eqnarray}$$(4) measuring the ratio between the buoyancy and viscous forces. Because viscosity and thermal conductivity are allowed to vary with temperature, there is no unique definition of Ra. For calculations, one may prescribe this number at a specific temperature. Here, I prescribed the surface Rayleigh number, Rasurf, calculated at the surface temperature, which remains constant throughout the duration of simulations. Viscosity varies with temperature following the Frank–Kamenetskii (FK) approximation, which is a linearized representation of Arrhenius-type of laws, and is given by $$\begin{eqnarray} \eta = {\eta _0}\exp \left[ { - {a_{\rm{\eta }}}\frac{{\left( {T - {T_0}} \right)}}{{{\rm{\Delta }}T}}} \right], \end{eqnarray}$$(5) where η0 and T0 are the reference viscosity and temperature (here, the surface values of these parameters, ηsurf and Tsurf), and aη a parameter that controls the amplitude of viscosity variations. The non-dimensional viscosity, |$\tilde{\eta } = \eta /{\eta _0}$|⁠, is then given as a function of |$\tilde{T}$| by $$\begin{eqnarray} \tilde{\eta } = {\rm{exp}}\left( { - {a_{\rm{\eta }}}\tilde{T}} \right). \end{eqnarray}$$(6) FK approximation is known to overestimate heat flux by up to 20–30 per cent (Reese et al. 1999; Harel et al. 2020). In addition, because different mechanisms may control the deformation of ice depending on the grain size, a composite viscosity law may be more appropriate to model ice shells dynamics (Harel et al. 2020). On another hand, FK approximation facilitates calculations and, most importantly, a large number of simulations is available in the literature and can be used, by comparison, to estimate the specific influence of temperature-dependent thermal conductivity. A key parameter controlling the flow regime is the top-to-bottom viscosity ratio, which, in the case of FK approximation and following eq. (5) is given by Δη = exp(aη). For viscosity ratios around 104 and larger convection operates in the so-called stagnant lid regime (e.g. Christensen 1984; Davaille & Jaupart 1993; Moresi & Solomatov 1995). In this mode of convection, a rigid (or stagnant) lid, in which heat is transported by conduction, forms at the top of the fluid. The stagnant lid may be viewed as the upward extension of the top thermal boundary layer (TBL), but because TBLs are, strictly speaking, regions where thermal instabilities are developing, I consider here that it does not belong to the top TBL. Following this point of view, the convective flow is confined beneath this lid and decomposes in a top (cold) TBL, a well-mixed adiabatic interior, and a bottom (hot) TBL. Most of the viscosity contrast is accommodated by the stagnant lid, and the convective sub-layer is nearly isoviscous, with viscosity ratio typically around 10. The viscosity of ice Ih strongly depends on temperature, and rheological laws imply viscosity ratios through outer ice shells of icy bodies much larger than 104 (Durham et al. 2010). If it develops within these shells, convection should thus operate in the stagnant lid regime. Thermal conductivity also varies with temperature. Based on laboratory experiments on ice Ih (Slack 1980; Andersson & Suga 1994), I assumed that this dependence follows an inverse law of the temperature, |$k \propto 1/T$|⁠. Knowing the value of thermal conductivity at a reference temperature, for instance the surface temperature Tsurf, the conductivity at temperature T is $$\begin{eqnarray} k\left( T \right) = {k_{\mathrm{ surf}}}\frac{{{T_{\mathrm{ surf}}}}}{T}, \end{eqnarray}$$(7) where ksurf is the conductivity at T = Tsurf. Taking ksurf as the characteristic value of the conductivity, the non-dimensional conductivity is |$\tilde{k} = k/{k_{\mathrm{ surf}}}$| and may be written as a function of the non-dimensional temperature |$\tilde{T}$| following $$\begin{eqnarray} \tilde{k} = \frac{{{T_{\mathrm{ surf}}}}}{{{T_{\mathrm{ surf}}} + \tilde{T}{\rm{\Delta }}T}} = \frac{1}{{1 + {R_T}\tilde{T}}} \end{eqnarray}$$(8) where $$\begin{eqnarray} {R_T} = \frac{{{\rm{\Delta }}T}}{{{T_{\mathrm{ surf}}}}} \end{eqnarray}$$(9) is the ratio between the dimensional super-adiabatic temperature jump, ΔT = (Tbot—Tsurf), and the dimensional surface temperature, Tsurf. With this definition, the ratio between the top |$(\tilde{T} = 0)$| and bottom |$(\tilde{T} = 1)$| thermal conductivities is simply given by $$\begin{eqnarray} {R_{\rm{k}}} = {T_{\mathrm{ bot}}}/{T_{\mathrm{ surf}}} = 1 + {R_T}. \end{eqnarray}$$(10) The input Rayleigh number, Rasurf, is defined with surface values of the viscosity and thermal conductivity. To describe convection in the convective sub-layer, it is convenient to define an effective Rayleigh number, Raeff, calculated with values of the viscosity and thermal conductivity at the temperature of the well-mixed interior (or interior temperature, for short), |${\tilde{T}_m}$|⁠. This temperature is an output of the simulations and is defined as the volume averaged temperature within the adiabatic region, that is, the region located between the top and bottom thermal boundary layers, where the horizontally averaged temperature profile is approximately constant with depth. Following eqs (6) and (8), Raeff is given as a function of the Rasurf by $$\begin{eqnarray} R{a_{\mathrm{ eff}}} = R{a_{\mathrm{ surf}}}\left( {1 + {R_T}{{\tilde{T}}_m}} \right){\rm{exp}}\left( {{a_{\rm{\eta }}}{{\tilde{T}}_m}} \right). \end{eqnarray}$$(11) A key observable, together with |${\tilde{T}_m},$| is the non-dimensional horizontally averaged heat flux, |${{\rm{\tilde{\Phi }}}_{\mathrm{ conv}}}$|⁠. In Cartesian geometry, and due to energy conservation, this heat flux is constant with depth and can be estimated either at the top or at the bottom of the domain, once quasi-stationary state is reached. For applications to natural systems, |${{\rm{\tilde{\Phi }}}_{\mathrm{ conv}}}$| is rescaled with the characteristic values of thermal conductivity, temperature and length, leading to a characteristic heat flux |${{\rm{\Phi }}_{\mathrm{ carac}}} = {k_{\mathrm{ ref}}}{\rm{\Delta }}T/D$|⁠, which is equivalent to the conductive heat flux for a Cartesian geometry system with homogeneous conductivity |${k_{\mathrm{ ref}}}$|⁠. The efficiency of convection is measured with the Nusselt number, Nu, defined as the ratio between the convective and conductive heat fluxes, Φconv and Φcond, or equivalently the ratio of their non-dimensional forms. A value of Nu lower than 1 indicates that the system is not animated by convection. For a system with constant conductivity and Cartesian geometry |${{\rm{\Phi }}_{\mathrm{ cond}}} = {{\rm{\Phi }}_{\mathrm{ carac}}}$| and Nu is equal to |${{\rm{\tilde{\Phi }}}_{\mathrm{ conv}}}$|⁠. If, however, conductivity varies as 1/T, the conductive profile of temperature is not linear (Section 3 and the Appendix), and the conductive heat flux is given by $$\begin{eqnarray} {{\rm{\Phi }}_{\mathrm{ cond}}} = \frac{{{k_{\mathrm{ surf}}}{T_{\mathrm{ surf}}}}}{D}{\rm{ln}}\left( {\frac{{{T_{\mathrm{ bot}}}}}{{{T_{\mathrm{ surf}}}}}} \right), \end{eqnarray}$$(12) which is smaller than the heat flux for a system with constant conductivity ksurf (Fig. A1). Its non-dimensional form (with respect to the characteristic heat flux, |${k_{\mathrm{ surf}}}{\rm{\Delta }}T/D$|⁠) writes $$\begin{eqnarray} {{\rm{\tilde{\Phi }}}_{\mathrm{ cond}}} = \frac{{{\rm{ln}}\left( {1 + {R_T}} \right)}}{{{R_T}}}, \end{eqnarray}$$(13) and is lower than 1 for RT > 1. For a given RT, convection operates if |${{\rm{\tilde{\Phi }}}_{\mathrm{ conv}}}$| is larger than the value of |${{\rm{\tilde{\Phi }}}_{\mathrm{ cond}}}$| (i.e.Nu > 1). Because eq. (13) monotonically decreases with RT (and thus Rk), the lower possible value of |${{\rm{\tilde{\Phi }}}_{\mathrm{ conv}}}$| decreases with increasing conductivity ratio. For instance, taking RT = 2 (Rk = 3) and RT = 9 (Rk = 10), the lowest possible |${{\rm{\tilde{\Phi }}}_{\mathrm{ conv}}}$| are 0.549 and 0.256, respectively. In all numerical simulations reported in this study, convection is well developed with Nu ranging from 3.0 to a bit more than 6.0 (Table 1). Figure 8. Open in new tabDownload slide Pluto Ice Ih shell properties at t = 4.55 Gyr as a function of the reference viscosity (ηref) and for pure water (left column) or an initial mix of water and 3.0 per cent NH3 (right column). The top, middle and bottom rows show the ice shell thickness, stagnant lid thickness and ice shell interior temperature, respectively. Plain curves indicate heat transfer by convection, and dashed curves heat transfer by conduction. Values of the physical properties used for calculations are listed in Table 2. For constant conductivity, two values of the characteristic conductivity are considered, ksurf and 2.6 W m−1 K−1. Figure 8. Open in new tabDownload slide Pluto Ice Ih shell properties at t = 4.55 Gyr as a function of the reference viscosity (ηref) and for pure water (left column) or an initial mix of water and 3.0 per cent NH3 (right column). The top, middle and bottom rows show the ice shell thickness, stagnant lid thickness and ice shell interior temperature, respectively. Plain curves indicate heat transfer by convection, and dashed curves heat transfer by conduction. Values of the physical properties used for calculations are listed in Table 2. For constant conductivity, two values of the characteristic conductivity are considered, ksurf and 2.6 W m−1 K−1. Figure 9. Open in new tabDownload slide Europa Ice Ih shell properties at t = 4.55 Gyr as a function of the reference viscosity (ηref) and for pure water (left column) or an initial mix of water and 3.0 per cent NH3 (right column). The top, middle and bottom rows show the ice shell thickness, stagnant lid thickness and ice shell interior temperature, respectively. Plain curves indicate heat transfer by convection, and dashed curves heat transfer by conduction. Values of the physical properties used for calculations are listed in Table 2. For constant conductivity, two values of the characteristic conductivity are considered, ksurf and 2.6 W m−1 K−1. Figure 9. Open in new tabDownload slide Europa Ice Ih shell properties at t = 4.55 Gyr as a function of the reference viscosity (ηref) and for pure water (left column) or an initial mix of water and 3.0 per cent NH3 (right column). The top, middle and bottom rows show the ice shell thickness, stagnant lid thickness and ice shell interior temperature, respectively. Plain curves indicate heat transfer by convection, and dashed curves heat transfer by conduction. Values of the physical properties used for calculations are listed in Table 2. For constant conductivity, two values of the characteristic conductivity are considered, ksurf and 2.6 W m−1 K−1. Figure 10. Open in new tabDownload slide Temperature at the bottom of the ice Ih shell for pure water (top row) and for an initial mix of water and 3.0 per cent NH3 (middle row), and ocean concentration in NH3 (bottom row) at t = 4.55 Gyr as a function of the reference viscosity (ηref), and for Pluto (left column) and Europa (right column). Plain curves indicate heat transfer by convection, and dashed curves heat transfer by conduction. Values of the physical properties used for calculations are listed in Table 2. For constant conductivity, two values of the characteristic conductivity are considered, ksurf and 2.6 W m−1 K−1. Figure 10. Open in new tabDownload slide Temperature at the bottom of the ice Ih shell for pure water (top row) and for an initial mix of water and 3.0 per cent NH3 (middle row), and ocean concentration in NH3 (bottom row) at t = 4.55 Gyr as a function of the reference viscosity (ηref), and for Pluto (left column) and Europa (right column). Plain curves indicate heat transfer by convection, and dashed curves heat transfer by conduction. Values of the physical properties used for calculations are listed in Table 2. For constant conductivity, two values of the characteristic conductivity are considered, ksurf and 2.6 W m−1 K−1. Figure 11. Open in new tabDownload slide Pluto (left column) and Europa (right column) ice Ih shells properties at t = 4.55 Gyr as a function of the initial fraction of ammonia |$( {x_{\mathrm{ NH3}}^{\mathrm{ init}}} )$|⁠. The top, middle and bottom rows show the ice shell thickness, stagnant lid thickness and ice shell interior temperature, respectively. Plain curves indicate heat transfer by convection, and dashed curves heat transfer by conduction. The reference viscosity (ηref) is set to 1014 Pa s, and the values of other physical properties are listed in Table 2. For constant conductivity, two values of the characteristic conductivity are considered, ksurf and 2.6 W m−1 K−1. Figure 11. Open in new tabDownload slide Pluto (left column) and Europa (right column) ice Ih shells properties at t = 4.55 Gyr as a function of the initial fraction of ammonia |$( {x_{\mathrm{ NH3}}^{\mathrm{ init}}} )$|⁠. The top, middle and bottom rows show the ice shell thickness, stagnant lid thickness and ice shell interior temperature, respectively. Plain curves indicate heat transfer by convection, and dashed curves heat transfer by conduction. The reference viscosity (ηref) is set to 1014 Pa s, and the values of other physical properties are listed in Table 2. For constant conductivity, two values of the characteristic conductivity are considered, ksurf and 2.6 W m−1 K−1. Table 1. Simulations of stagnant lid convection with temperature-dependent thermal conductivity. # . Rasurf . Δη . Rk . |${\tilde{T}_m}$| . |${\tilde{k}_m}$| . |${{\rm{\tilde{\Phi }}}_{\mathrm{ conv}}}$| . Nu . |$\mathrm{ rms}( {\tilde{v}} )$| . |${\tilde{v}_{\mathrm{ surf}}}$| . Raeff . |${\tilde{d}_{\mathrm{ lid}}}$| . |${\tilde{T}_{\mathrm{ lid}}}$| . 1 31.62 105 1.0 0.891 1.0 3.144 3.144 148.1 3.0 × 10−1 9.06 × 105 0.257 0.809 2 10.00 106 1.0 0.913 1.0 3.437 3.437 257.3 1.4 × 10−1 2.99 × 106 0.245 0.842 3 3.16 107 1.0 0.923 1.0 3.887 3.887 458.5 3.9 × 10−2 9.18 × 106 0.218 0.848 4 1.00 108 1.0 0.933 1.0 4.445 4.445 875.2 8.6 × 10−3 2.91 × 107 0.196 0.871 5 10.00 106 1.143 0.909 0.885 3.173 3.395 240.9 1.0 × 10−1 3.23 × 106 0.249 0.836 6 31.62 105 1.333 0.885 0.772 2.659 3.080 128.1 4.9 × 10−1 1.09 × 106 0.263 0.787 7 10.00 106 1.333 0.907 0.768 2.937 3.403 224.7 1.3 × 10−1 3.62 × 106 0.249 0.826 8 3.16 107 1.333 0.922 0.765 3.336 3.865 407.6 2.9 × 10−2 1.18 × 107 0.226 0.857 9 1.00 108 1.333 0.931 0.763 3.816 4.421 751.0 7.9 × 10−3 3.69 × 107 0.201 0.874 10 10.00 106 1.5 0.908 0.688 2.749 3.389 213.1 9.9 × 10−2 4.10 × 106 0.253 0.833 11 10.00 106 2.0 0.903 0.525 2.331 3.362 184.4 1.2 × 10−1 5.01 × 106 0.256 0.816 12 14.86 107 2.0 0.919 0.521 3.867 5.578 908.6 6.8 × 10−2 7.68 × 107 0.165 0.894 13 10.00 106 2.5 0.901 0.425 2.069 3.387 161.9 9.6 × 10−2 6.01 × 106 0.258 0.817 14 31.62 105 3.0 0.878 0.363 1.723 3.137 86.6 3.8 × 10−1 2.15 × 106 0.270 0.766 15 10.00 106 3.0 0.900 0.357 1.816 3.305 152.6 1.2 × 10−1 7.02 × 106 0.264 0.806 16 16.00 106 3.0 0.899 0.357 2.066 3.761 210.4 1.6 × 10−1 1.11 × 107 0.232 0.803 17 32.00 106 3.0 0.899 0.357 2.514 4.576 328.7 2.0 × 10−1 2.22 × 107 0.193 0.817 18 100.00 106 3.0 0.901 0.357 3.374 6.142 667.4 3.6 × 10−1 7.11 × 107 0.150 0.874 19 3.16 107 3.0 0.916 0.353 2.026 3.687 279.7 2.1 × 10−2 2.31 × 107 0.242 0.832 20 1.00 108 3.0 0.927 0.350 2.279 4.149 518.7 5.5 × 10−3 7.44 × 107 0.221 0.870 21 10.00 106 4.0 0.898 0.271 1.539 3.331 133.6 8.4 × 10−2 9.06 × 106 0.265 0.801 22 4.74 107 4.0 0.914 0.267 1.876 4.060 321.9 2.5 × 10−2 4.45 × 107 0.226 0.855 23 31.62 105 6.0 0.865 0.188 1.128 3.148 61.4 6.5 × 10−1 3.54 × 106 0.274 0.737 24 10.00 106 6.0 0.895 0.183 1.173 3.272 111.5 9.0 × 10−2 1.29 × 107 0.271 0.780 25 3.16 107 6.0 0.912 0.180 1.293 3.607 207.0 1.5 × 10−2 4.29 × 107 0.254 0.830 26 1.00 108 6.0 0.925 0.178 1.438 4.013 382.8 4.0 × 10−3 1.43 × 108 0.235 0.887 27 10.00 106 8.0 0.894 0.138 0.973 3.274 97.0 8.3 × 10−2 1.68 × 107 0.274 0.781 28 10.00 106 10.0 0.894 0.111 0.837 3.270 90.7 5.5 × 10−2 2.08 × 107 0.275 0.771 # . Rasurf . Δη . Rk . |${\tilde{T}_m}$| . |${\tilde{k}_m}$| . |${{\rm{\tilde{\Phi }}}_{\mathrm{ conv}}}$| . Nu . |$\mathrm{ rms}( {\tilde{v}} )$| . |${\tilde{v}_{\mathrm{ surf}}}$| . Raeff . |${\tilde{d}_{\mathrm{ lid}}}$| . |${\tilde{T}_{\mathrm{ lid}}}$| . 1 31.62 105 1.0 0.891 1.0 3.144 3.144 148.1 3.0 × 10−1 9.06 × 105 0.257 0.809 2 10.00 106 1.0 0.913 1.0 3.437 3.437 257.3 1.4 × 10−1 2.99 × 106 0.245 0.842 3 3.16 107 1.0 0.923 1.0 3.887 3.887 458.5 3.9 × 10−2 9.18 × 106 0.218 0.848 4 1.00 108 1.0 0.933 1.0 4.445 4.445 875.2 8.6 × 10−3 2.91 × 107 0.196 0.871 5 10.00 106 1.143 0.909 0.885 3.173 3.395 240.9 1.0 × 10−1 3.23 × 106 0.249 0.836 6 31.62 105 1.333 0.885 0.772 2.659 3.080 128.1 4.9 × 10−1 1.09 × 106 0.263 0.787 7 10.00 106 1.333 0.907 0.768 2.937 3.403 224.7 1.3 × 10−1 3.62 × 106 0.249 0.826 8 3.16 107 1.333 0.922 0.765 3.336 3.865 407.6 2.9 × 10−2 1.18 × 107 0.226 0.857 9 1.00 108 1.333 0.931 0.763 3.816 4.421 751.0 7.9 × 10−3 3.69 × 107 0.201 0.874 10 10.00 106 1.5 0.908 0.688 2.749 3.389 213.1 9.9 × 10−2 4.10 × 106 0.253 0.833 11 10.00 106 2.0 0.903 0.525 2.331 3.362 184.4 1.2 × 10−1 5.01 × 106 0.256 0.816 12 14.86 107 2.0 0.919 0.521 3.867 5.578 908.6 6.8 × 10−2 7.68 × 107 0.165 0.894 13 10.00 106 2.5 0.901 0.425 2.069 3.387 161.9 9.6 × 10−2 6.01 × 106 0.258 0.817 14 31.62 105 3.0 0.878 0.363 1.723 3.137 86.6 3.8 × 10−1 2.15 × 106 0.270 0.766 15 10.00 106 3.0 0.900 0.357 1.816 3.305 152.6 1.2 × 10−1 7.02 × 106 0.264 0.806 16 16.00 106 3.0 0.899 0.357 2.066 3.761 210.4 1.6 × 10−1 1.11 × 107 0.232 0.803 17 32.00 106 3.0 0.899 0.357 2.514 4.576 328.7 2.0 × 10−1 2.22 × 107 0.193 0.817 18 100.00 106 3.0 0.901 0.357 3.374 6.142 667.4 3.6 × 10−1 7.11 × 107 0.150 0.874 19 3.16 107 3.0 0.916 0.353 2.026 3.687 279.7 2.1 × 10−2 2.31 × 107 0.242 0.832 20 1.00 108 3.0 0.927 0.350 2.279 4.149 518.7 5.5 × 10−3 7.44 × 107 0.221 0.870 21 10.00 106 4.0 0.898 0.271 1.539 3.331 133.6 8.4 × 10−2 9.06 × 106 0.265 0.801 22 4.74 107 4.0 0.914 0.267 1.876 4.060 321.9 2.5 × 10−2 4.45 × 107 0.226 0.855 23 31.62 105 6.0 0.865 0.188 1.128 3.148 61.4 6.5 × 10−1 3.54 × 106 0.274 0.737 24 10.00 106 6.0 0.895 0.183 1.173 3.272 111.5 9.0 × 10−2 1.29 × 107 0.271 0.780 25 3.16 107 6.0 0.912 0.180 1.293 3.607 207.0 1.5 × 10−2 4.29 × 107 0.254 0.830 26 1.00 108 6.0 0.925 0.178 1.438 4.013 382.8 4.0 × 10−3 1.43 × 108 0.235 0.887 27 10.00 106 8.0 0.894 0.138 0.973 3.274 97.0 8.3 × 10−2 1.68 × 107 0.274 0.781 28 10.00 106 10.0 0.894 0.111 0.837 3.270 90.7 5.5 × 10−2 2.08 × 107 0.275 0.771 All simulations are performed on 3-D-Cartesian grids of 256 × 256 × 128 points, with x and y aspect ratios both equal to 4. Calculations with constant thermal conductivity (Rk = 1) are from Deschamps & Lin (2014). For consistency with other cases, the Rk = 1 simulations with Δη = 105 and Δη = 106, originally calculated on grids with 128 × 128 × 64 points, were smoothed on grids with 256 × 256 × 128 points. Listed parameters are the surface Rayleigh number, Rasurf, the top-to-bottom thermal viscosity ratio, Δη, the top-to-bottom thermal conductivity ratio, Rk, the average non-dimensional temperature and thermal conductivity of the well-mixed interior, |${\tilde{T}_m}$| and |${\tilde{k}_m}$|⁠, respectively, the horizontally averaged non-dimensional convective heat flux, |${{\rm{\tilde{\Phi }}}_{\mathrm{ conv}}}$|⁠, the Nusselt number, Nu, the root mean square velocity of the whole system, |$\mathrm{ rms}( {\tilde{v}} )$|⁠, the average surface velocity, |${\tilde{v}_{\mathrm{ surf}}}$|⁠, the effective Rayleigh number, Raeff, calculated with viscosity and thermal conductivity at temperature |${\tilde{T}_m}$|⁠, the non-dimensional thickness of the stagnant lid, |${\tilde{d}_{\mathrm{ lid}}}$|⁠, estimated following the method of Davaille & Jaupart (1993), and the temperature at the base of this lid, |${\tilde{T}_{\mathrm{ lid}}}$|⁠, deduced from eq. (21) with observed values of |${{\rm{\tilde{\Phi }}}_{\mathrm{ conv}}}$| and |${\tilde{d}_{\mathrm{ lid}}}$|⁠. Open in new tab Table 1. Simulations of stagnant lid convection with temperature-dependent thermal conductivity. # . Rasurf . Δη . Rk . |${\tilde{T}_m}$| . |${\tilde{k}_m}$| . |${{\rm{\tilde{\Phi }}}_{\mathrm{ conv}}}$| . Nu . |$\mathrm{ rms}( {\tilde{v}} )$| . |${\tilde{v}_{\mathrm{ surf}}}$| . Raeff . |${\tilde{d}_{\mathrm{ lid}}}$| . |${\tilde{T}_{\mathrm{ lid}}}$| . 1 31.62 105 1.0 0.891 1.0 3.144 3.144 148.1 3.0 × 10−1 9.06 × 105 0.257 0.809 2 10.00 106 1.0 0.913 1.0 3.437 3.437 257.3 1.4 × 10−1 2.99 × 106 0.245 0.842 3 3.16 107 1.0 0.923 1.0 3.887 3.887 458.5 3.9 × 10−2 9.18 × 106 0.218 0.848 4 1.00 108 1.0 0.933 1.0 4.445 4.445 875.2 8.6 × 10−3 2.91 × 107 0.196 0.871 5 10.00 106 1.143 0.909 0.885 3.173 3.395 240.9 1.0 × 10−1 3.23 × 106 0.249 0.836 6 31.62 105 1.333 0.885 0.772 2.659 3.080 128.1 4.9 × 10−1 1.09 × 106 0.263 0.787 7 10.00 106 1.333 0.907 0.768 2.937 3.403 224.7 1.3 × 10−1 3.62 × 106 0.249 0.826 8 3.16 107 1.333 0.922 0.765 3.336 3.865 407.6 2.9 × 10−2 1.18 × 107 0.226 0.857 9 1.00 108 1.333 0.931 0.763 3.816 4.421 751.0 7.9 × 10−3 3.69 × 107 0.201 0.874 10 10.00 106 1.5 0.908 0.688 2.749 3.389 213.1 9.9 × 10−2 4.10 × 106 0.253 0.833 11 10.00 106 2.0 0.903 0.525 2.331 3.362 184.4 1.2 × 10−1 5.01 × 106 0.256 0.816 12 14.86 107 2.0 0.919 0.521 3.867 5.578 908.6 6.8 × 10−2 7.68 × 107 0.165 0.894 13 10.00 106 2.5 0.901 0.425 2.069 3.387 161.9 9.6 × 10−2 6.01 × 106 0.258 0.817 14 31.62 105 3.0 0.878 0.363 1.723 3.137 86.6 3.8 × 10−1 2.15 × 106 0.270 0.766 15 10.00 106 3.0 0.900 0.357 1.816 3.305 152.6 1.2 × 10−1 7.02 × 106 0.264 0.806 16 16.00 106 3.0 0.899 0.357 2.066 3.761 210.4 1.6 × 10−1 1.11 × 107 0.232 0.803 17 32.00 106 3.0 0.899 0.357 2.514 4.576 328.7 2.0 × 10−1 2.22 × 107 0.193 0.817 18 100.00 106 3.0 0.901 0.357 3.374 6.142 667.4 3.6 × 10−1 7.11 × 107 0.150 0.874 19 3.16 107 3.0 0.916 0.353 2.026 3.687 279.7 2.1 × 10−2 2.31 × 107 0.242 0.832 20 1.00 108 3.0 0.927 0.350 2.279 4.149 518.7 5.5 × 10−3 7.44 × 107 0.221 0.870 21 10.00 106 4.0 0.898 0.271 1.539 3.331 133.6 8.4 × 10−2 9.06 × 106 0.265 0.801 22 4.74 107 4.0 0.914 0.267 1.876 4.060 321.9 2.5 × 10−2 4.45 × 107 0.226 0.855 23 31.62 105 6.0 0.865 0.188 1.128 3.148 61.4 6.5 × 10−1 3.54 × 106 0.274 0.737 24 10.00 106 6.0 0.895 0.183 1.173 3.272 111.5 9.0 × 10−2 1.29 × 107 0.271 0.780 25 3.16 107 6.0 0.912 0.180 1.293 3.607 207.0 1.5 × 10−2 4.29 × 107 0.254 0.830 26 1.00 108 6.0 0.925 0.178 1.438 4.013 382.8 4.0 × 10−3 1.43 × 108 0.235 0.887 27 10.00 106 8.0 0.894 0.138 0.973 3.274 97.0 8.3 × 10−2 1.68 × 107 0.274 0.781 28 10.00 106 10.0 0.894 0.111 0.837 3.270 90.7 5.5 × 10−2 2.08 × 107 0.275 0.771 # . Rasurf . Δη . Rk . |${\tilde{T}_m}$| . |${\tilde{k}_m}$| . |${{\rm{\tilde{\Phi }}}_{\mathrm{ conv}}}$| . Nu . |$\mathrm{ rms}( {\tilde{v}} )$| . |${\tilde{v}_{\mathrm{ surf}}}$| . Raeff . |${\tilde{d}_{\mathrm{ lid}}}$| . |${\tilde{T}_{\mathrm{ lid}}}$| . 1 31.62 105 1.0 0.891 1.0 3.144 3.144 148.1 3.0 × 10−1 9.06 × 105 0.257 0.809 2 10.00 106 1.0 0.913 1.0 3.437 3.437 257.3 1.4 × 10−1 2.99 × 106 0.245 0.842 3 3.16 107 1.0 0.923 1.0 3.887 3.887 458.5 3.9 × 10−2 9.18 × 106 0.218 0.848 4 1.00 108 1.0 0.933 1.0 4.445 4.445 875.2 8.6 × 10−3 2.91 × 107 0.196 0.871 5 10.00 106 1.143 0.909 0.885 3.173 3.395 240.9 1.0 × 10−1 3.23 × 106 0.249 0.836 6 31.62 105 1.333 0.885 0.772 2.659 3.080 128.1 4.9 × 10−1 1.09 × 106 0.263 0.787 7 10.00 106 1.333 0.907 0.768 2.937 3.403 224.7 1.3 × 10−1 3.62 × 106 0.249 0.826 8 3.16 107 1.333 0.922 0.765 3.336 3.865 407.6 2.9 × 10−2 1.18 × 107 0.226 0.857 9 1.00 108 1.333 0.931 0.763 3.816 4.421 751.0 7.9 × 10−3 3.69 × 107 0.201 0.874 10 10.00 106 1.5 0.908 0.688 2.749 3.389 213.1 9.9 × 10−2 4.10 × 106 0.253 0.833 11 10.00 106 2.0 0.903 0.525 2.331 3.362 184.4 1.2 × 10−1 5.01 × 106 0.256 0.816 12 14.86 107 2.0 0.919 0.521 3.867 5.578 908.6 6.8 × 10−2 7.68 × 107 0.165 0.894 13 10.00 106 2.5 0.901 0.425 2.069 3.387 161.9 9.6 × 10−2 6.01 × 106 0.258 0.817 14 31.62 105 3.0 0.878 0.363 1.723 3.137 86.6 3.8 × 10−1 2.15 × 106 0.270 0.766 15 10.00 106 3.0 0.900 0.357 1.816 3.305 152.6 1.2 × 10−1 7.02 × 106 0.264 0.806 16 16.00 106 3.0 0.899 0.357 2.066 3.761 210.4 1.6 × 10−1 1.11 × 107 0.232 0.803 17 32.00 106 3.0 0.899 0.357 2.514 4.576 328.7 2.0 × 10−1 2.22 × 107 0.193 0.817 18 100.00 106 3.0 0.901 0.357 3.374 6.142 667.4 3.6 × 10−1 7.11 × 107 0.150 0.874 19 3.16 107 3.0 0.916 0.353 2.026 3.687 279.7 2.1 × 10−2 2.31 × 107 0.242 0.832 20 1.00 108 3.0 0.927 0.350 2.279 4.149 518.7 5.5 × 10−3 7.44 × 107 0.221 0.870 21 10.00 106 4.0 0.898 0.271 1.539 3.331 133.6 8.4 × 10−2 9.06 × 106 0.265 0.801 22 4.74 107 4.0 0.914 0.267 1.876 4.060 321.9 2.5 × 10−2 4.45 × 107 0.226 0.855 23 31.62 105 6.0 0.865 0.188 1.128 3.148 61.4 6.5 × 10−1 3.54 × 106 0.274 0.737 24 10.00 106 6.0 0.895 0.183 1.173 3.272 111.5 9.0 × 10−2 1.29 × 107 0.271 0.780 25 3.16 107 6.0 0.912 0.180 1.293 3.607 207.0 1.5 × 10−2 4.29 × 107 0.254 0.830 26 1.00 108 6.0 0.925 0.178 1.438 4.013 382.8 4.0 × 10−3 1.43 × 108 0.235 0.887 27 10.00 106 8.0 0.894 0.138 0.973 3.274 97.0 8.3 × 10−2 1.68 × 107 0.274 0.781 28 10.00 106 10.0 0.894 0.111 0.837 3.270 90.7 5.5 × 10−2 2.08 × 107 0.275 0.771 All simulations are performed on 3-D-Cartesian grids of 256 × 256 × 128 points, with x and y aspect ratios both equal to 4. Calculations with constant thermal conductivity (Rk = 1) are from Deschamps & Lin (2014). For consistency with other cases, the Rk = 1 simulations with Δη = 105 and Δη = 106, originally calculated on grids with 128 × 128 × 64 points, were smoothed on grids with 256 × 256 × 128 points. Listed parameters are the surface Rayleigh number, Rasurf, the top-to-bottom thermal viscosity ratio, Δη, the top-to-bottom thermal conductivity ratio, Rk, the average non-dimensional temperature and thermal conductivity of the well-mixed interior, |${\tilde{T}_m}$| and |${\tilde{k}_m}$|⁠, respectively, the horizontally averaged non-dimensional convective heat flux, |${{\rm{\tilde{\Phi }}}_{\mathrm{ conv}}}$|⁠, the Nusselt number, Nu, the root mean square velocity of the whole system, |$\mathrm{ rms}( {\tilde{v}} )$|⁠, the average surface velocity, |${\tilde{v}_{\mathrm{ surf}}}$|⁠, the effective Rayleigh number, Raeff, calculated with viscosity and thermal conductivity at temperature |${\tilde{T}_m}$|⁠, the non-dimensional thickness of the stagnant lid, |${\tilde{d}_{\mathrm{ lid}}}$|⁠, estimated following the method of Davaille & Jaupart (1993), and the temperature at the base of this lid, |${\tilde{T}_{\mathrm{ lid}}}$|⁠, deduced from eq. (21) with observed values of |${{\rm{\tilde{\Phi }}}_{\mathrm{ conv}}}$| and |${\tilde{d}_{\mathrm{ lid}}}$|⁠. Open in new tab Simulations are performed in 3-D-Cartesian domains with a horizontal to vertical aspect ratio equal to 4 in both x and y directions, and a grid resolution of 256 × 256 × 128 points. The grid is vertically refined at the top and at the bottom of the box, which provides a better description of TBLs. The top and bottom boundaries are free slip and isothermal, and reflective boundary conditions are imposed on sidewalls. Initial conditions on temperature consist of small random perturbations superposed on a 1-D radial adiabatic profile with thin TBLs at the top and the bottom. Using this set-up, I performed 24 numerical simulations (Table 1). For comparison, I added 4 simulations from Deschamps & Lin (2014), which were conducted with a similar set-up, but with constant thermal conductivity. Explored ranges of top-to-bottom viscosity and thermal conductivity ratios are 105 ≤ Δη ≤ 108 and 1 ≤ Rk ≤ 10, respectively. Surface Rayleigh numbers are chosen in the range 1–102, leading to effective Rayleigh numbers (eq. 11) between 9.0 × 105 and 1.5 × 108. For these ranges of values the flow is time-dependent and reaches a quasi-stationary state, meaning that global output properties (e.g. the volume average and interior temperatures, and the horizontally averaged heat flux) oscillate around values that are constant in time. Output properties are estimated after the quasi-stationary phase has been reached, by time-averaging of each property over several oscillations. The averaging time-window is typically around 0.5–1.0 (in non-dimensional time), corresponding to 105–106 iterations (note that the number of iterations per time unit increases with the Rayleigh number). The non-dimensional heat fluxes and Nusselt number are calculated from the top and bottom temperature gradients at each time step and averaged over the time-window, while the interior temperature is obtained from the vertical average (in the well-mixed interior) of the profile of horizontally averaged temperature at each output frame (recorded with a non-dimensional time interval of 10−3, different from the iteration time step) and, again, averaged over the time-window. 3 FLOW AND THERMAL STRUCTURES All 28 experiments are well within the stagnant lid regime, as indicated by the low values of the non-dimensional surface velocity, |${\tilde{v}_{\mathrm{ surf}}}$|⁠, and mobility M, defined as the ratio between |${\tilde{v}_{\mathrm{ surf}}}$| and the root mean square (rms) velocity of the whole system. More precisely, in all cases |${\tilde{v}_{\mathrm{ surf}}}$| < 1 and M < 0.01, satisfying the criteria defined by Stein et al. (2013) for the presence of a stagnant lid at the top of the system. A more detailed analysis based on the horizontally averaged profile of vertically advected heat (Davaille & Jaupart 1993) shows that the thickness of the lid varies between 0.15 and 0.28 depending on the case (Section 4.3 and Table 1). Fig. 1 plots snapshots of temperature fields for different cases. The flow structure is very similar to that obtained for stagnant lid convection with constant thermal conductivity in 3-D-Cartesian and spherical geometries (e.g. Deschamps & Lin 2014; Yao et al. 2014), and organizes as network of plumes rising from the bottom thermal boundary layer up to the base of the stagnant lid. A closer examination indicates that, for given values of the surface Rayleigh number, Rasurf, and viscosity contrast, Δη, plumes get thinner as the thermal conductivity ratio, Rk, increases (Fig. 2). This suggests that convection in the well-mixed interior is more vigorous, in agreement with the observation that the effective Rayleigh number, Raeff, increases with Rk due to the decrease in the average conductivity |$\tilde{k}$| of this region (Table 1). Still for given values of Rasurf and Δη, the temperature of the well-mixed interior slightly decreases with increasing Rk (Table 1). More importantly, because thermal conductivity at the bottom of the system decreases, the basal heat flow entering the system is reduced. This strongly impacts heat transfer, outweighing the increase in the vigour of convection, and resulting in a sharp decrease in the convective heat flux, |${{\rm{\tilde{\Phi }}}_{\mathrm{ conv}}}$|⁠. For instance, taking Δη = 106 and Rasurf = 10, |${{\rm{\tilde{\Phi }}}_{\mathrm{ conv}}}$| drops by a factor 4 for Rk in the range 1–10. In addition, because less energy is available to entertain convection, the average flow velocity, measured with the rms non-dimensional velocity, |$\mathrm{ rms}( {\tilde{v}} )$|⁠, is reduced. Taking, again, Δη = 106 and Rasurf = 10, |$\mathrm{ rms}( {\tilde{v}} )$| drops by a factor 2.5 for Rk in the range 1–10. The horizontally averaged vertical component of the velocity, |${\tilde{v}_z}$|⁠, is also smaller, indicating that the reduction in |$\mathrm{ rms}( {\tilde{v}} )$| is not related to the thickening of the stagnant lid, and implying that the amount of vertically advected heat, |${\tilde{v}_z}\tilde{T}$|⁠, is also substantially reduced (Fig. 3). Figure 1. Open in new tabDownload slide Snapshots of the temperature field for cases with surface Rayleigh number Rasurf = 10, thermal viscosity ratio Δη = 106 and different values of the top-to-bottom thermal conductivity ratio, Rk. (a) Constant conductivity (case #2, Rk = 1). In other cases, thermal conductivity depends on temperature following eq. (8) with (b) RT = 1/3 (case #7, Rk = 4/3), (c) RT = 2 (case #15, Rk = 3) and (d) RT = 9 (case #28, Rk = 10). Figure 1. Open in new tabDownload slide Snapshots of the temperature field for cases with surface Rayleigh number Rasurf = 10, thermal viscosity ratio Δη = 106 and different values of the top-to-bottom thermal conductivity ratio, Rk. (a) Constant conductivity (case #2, Rk = 1). In other cases, thermal conductivity depends on temperature following eq. (8) with (b) RT = 1/3 (case #7, Rk = 4/3), (c) RT = 2 (case #15, Rk = 3) and (d) RT = 9 (case #28, Rk = 10). Figure 2. Open in new tabDownload slide Residual non-dimensional temperature, |$\delta {\tilde{T}_{\mathrm{ plume}}} = ( {\tilde{T} - {{\tilde{T}}_{\mathrm{ plume}}}} )\ $|⁠, relative to the plume threshold temperature defined as |${\tilde{T}_{\mathrm{ plume}}} = {\tilde{T}_z}\ + 0.3( {{{\tilde{T}}_{\mathrm{ max}}} - {{\tilde{T}}_z}} )$|⁠, where |${\tilde{T}_z}$| and |${\tilde{T}_{\max}}$| are the average and maximum temperatures at non-dimensional depth |$\tilde{z}$|⁠. With this definition, plumes conduits appear in orange and red. Cases represented are similar to those shown in plots (a) and (c) of Fig. 1, with surface Rayleigh number and viscosity ratio equal to Rasurf = 10 and Δη = 106 for the two cases, and thermal conductivity ratio, Rk, equal to 1 (case #2 in Table 1; left column) and 3 (case #15; right column). Figure 2. Open in new tabDownload slide Residual non-dimensional temperature, |$\delta {\tilde{T}_{\mathrm{ plume}}} = ( {\tilde{T} - {{\tilde{T}}_{\mathrm{ plume}}}} )\ $|⁠, relative to the plume threshold temperature defined as |${\tilde{T}_{\mathrm{ plume}}} = {\tilde{T}_z}\ + 0.3( {{{\tilde{T}}_{\mathrm{ max}}} - {{\tilde{T}}_z}} )$|⁠, where |${\tilde{T}_z}$| and |${\tilde{T}_{\max}}$| are the average and maximum temperatures at non-dimensional depth |$\tilde{z}$|⁠. With this definition, plumes conduits appear in orange and red. Cases represented are similar to those shown in plots (a) and (c) of Fig. 1, with surface Rayleigh number and viscosity ratio equal to Rasurf = 10 and Δη = 106 for the two cases, and thermal conductivity ratio, Rk, equal to 1 (case #2 in Table 1; left column) and 3 (case #15; right column). Figure 3. Open in new tabDownload slide Horizontally averaged profiles of the vertically advected heat flow (left plot in each panel) and temperature (right plot) for cases with surface Rayleigh number Rasurf = 10, viscosity ratio Δη = 106 and different values of the top-to-bottom thermal conductivity ratio, Rk. (a) Rk = 1 (constant thermal conductivity, case #2 in Table 1). (b) Rk = 4/3 (case #7). (c) Rk = 3 (case #15). (d) Rk = 10 (case #28). The grey areas denote the vertical extension of the stagnant lid. The dashed lines in plots of advected heat flow show the tangent to the point of inflexion, whose intersection with the origin axis defines the bottom of the lid. The dashed dark-red curves in temperature plots are determined assuming a conductive temperature profile in the stagnant lid and are calculated following either |$\tilde{T}\ ( {\tilde{z}} ) = {\tilde{T}_{\mathrm{ lid}}}\ z/{\tilde{d}_{\mathrm{ lid}}}$| (panel a) or eq. (19) (panels b–d) with values of |${{\rm{\tilde{\Phi }}}_{\mathrm{ conv}}}$| and |${\tilde{d}_{\mathrm{ lid}}}$| listed in Table 1, and values of |${\tilde{T}_{\mathrm{ lid}}}$|⁠. estimated from eq. (21). Figure 3. Open in new tabDownload slide Horizontally averaged profiles of the vertically advected heat flow (left plot in each panel) and temperature (right plot) for cases with surface Rayleigh number Rasurf = 10, viscosity ratio Δη = 106 and different values of the top-to-bottom thermal conductivity ratio, Rk. (a) Rk = 1 (constant thermal conductivity, case #2 in Table 1). (b) Rk = 4/3 (case #7). (c) Rk = 3 (case #15). (d) Rk = 10 (case #28). The grey areas denote the vertical extension of the stagnant lid. The dashed lines in plots of advected heat flow show the tangent to the point of inflexion, whose intersection with the origin axis defines the bottom of the lid. The dashed dark-red curves in temperature plots are determined assuming a conductive temperature profile in the stagnant lid and are calculated following either |$\tilde{T}\ ( {\tilde{z}} ) = {\tilde{T}_{\mathrm{ lid}}}\ z/{\tilde{d}_{\mathrm{ lid}}}$| (panel a) or eq. (19) (panels b–d) with values of |${{\rm{\tilde{\Phi }}}_{\mathrm{ conv}}}$| and |${\tilde{d}_{\mathrm{ lid}}}$| listed in Table 1, and values of |${\tilde{T}_{\mathrm{ lid}}}$|⁠. estimated from eq. (21). Another key difference, compared to simulations with constant thermal conductivity, is the thermal structure within the stagnant lid. In cases with constant thermal conductivity, the top part of the horizontally averaged profile of temperature is linear, as a result of the presence of a stagnant and thermally conductive lid at the top of the system. In cases with temperature-dependent conductivity, by contrast, the thermal gradient decreases at the top of the system, and the temperature profile is not linear (Fig. 3). A stagnant lid is still present and heat is still transported by conduction, but due to the top-to-bottom decrease in conductivity, satisfying the conservation energy (which, in Cartesian geometry, implies that the top and bottom heat fluxes are equal) requires an adjustment of the thermal gradient and temperature profile. Solving the heat conduction equation for a conductivity varying as 1/T (see the Appendix) indicates that in Cartesian geometry, and in the absence of internal heating, the conductive temperature profile through the whole domain is $$\begin{eqnarray} T\left( z \right) = {T_{\mathrm{ surf}}}{\rm{exp}}\left[ {\frac{z}{D}{\rm{ln}}\left( {\frac{{{T_{\mathrm{ bot}}}}}{{{T_{\mathrm{ surf}}}}}} \right)} \right], \end{eqnarray}$$(14) where Tbot is the bottom temperature, z the depth and D the thickness of the domain. The non-dimensional form of eq. (14) writes $$\begin{eqnarray} \tilde{T}\left( {\tilde{z}} \right) = \frac{1}{{{R_T}}}\left\{ {{\rm{exp}}\left[ {\tilde{z}{\rm{ln}}\left( {1 + {R_T}} \right)} \right] - 1} \right\}, \end{eqnarray}$$(15) where |${R_T} = {\rm{\Delta }}T/{T_{\mathrm{ surf}}}$|⁠. A closer examination of the temperature profiles shows that the top sections of these profiles fit well along a non-dimensional version of eq. (14) in which D and Tbot are replaced by the thickness of the stagnant lid and the temperature at its base, respectively (Section 4.3 and dashed dark red curves in Fig. 3). 4 HEAT FLUX, TEMPERATURE AND STAGNANT LID For planetary applications, in particular the modelling of thermal evolution of icy bodies, it is convenient to define scaling laws, which relate input parameters (e.g. Rayleigh number and viscosity ratio) to observables (e.g. average temperature and surface heat flux). This section derives scaling laws for fluids with temperature-dependent (1/T) thermal conductivity. 4.1 Temperature of the well-mixed interior Following previous studies (e.g. Davaille & Jaupart 1993; Moresi & Solomatov 1995; Deschamps & Sotin 2000), I assumed that the temperature jumps across the top and bottom TBLs are controlled by a viscous temperature scale, ΔTv, defined as $$\begin{eqnarray} \Delta {T_\mathrm{ v}} = {\left( { - \frac{1}{\eta }{{\left. {\frac{{\mathrm{ d}\eta }}{{\mathrm{ d}T}}} \right|}_{T = {T_m}}}} \right)^{ - 1}}. \end{eqnarray}$$(16) This approach implicitly assumes that the stagnant lid does not belong to the top TBL. Estimating the temperature jump across the top TBL therefore requires the knowledge of the temperature at the limit between the stagnant lid and the top TBL. In practice, this temperature is difficult to measure, and it is easier to relate the viscous temperature scale to the temperature jump across the bottom TBL given by |$( {1 - {{\tilde{T}}_m}} )$|⁠. Table 1 indicates that for given values of Rasurf and Δη the interior temperature decreases only slightly with increasing thermal conductivity ratio, Rk. These differences are however significant, that is, they have the same trend and are larger than the expected error bars (estimated from the time-variations of the average temperature), and they should thus be accounted for in a modified scaling law for |${\tilde{T}_m}$|⁠. Here, I implemented the scaling obtained by Deschamps & Lin (2014) for similar calculations with constant thermal conductivity, by adding the variable Rk to describe the influence of conductivity changes on temperature, leading to $$\begin{eqnarray} {\tilde{T}_m} = 1 - {a_\mathrm{ T}}\frac{{{R_k}^{{d_\mathrm{ T}}}}}{{\rm{\gamma }}}, \end{eqnarray}$$(17) where aT and dT are two constants that can be inferred by inverting the values of |${\tilde{T}_m}$| deduced from the simulations, and |${\rm{\gamma }} = \Delta T/\Delta {T_v}$| is the inverse of the non-dimensional viscous temperature scale. Note that in the case of FK approximation (eq. 3), γ is equal to the logarithmic viscosity ratio, aη = ln(Δη). Fixing the value of aT to that obtained by Deschamps & Lin (2014), aT = 1.23 ± 0.05, inversion of interior temperatures listed in Table 1 leads to dT = 0.1 ± 0.02, which, as illustrated in Fig. 4(a), fits the results of the numerical simulations very well. Figure 4. Open in new tabDownload slide Comparison between observed and modelled output properties. (a) Temperature of the well-mixed interior, |${\tilde{T}_m}$|⁠. Observed values are listed in Table 1, and modelled values are given by eq. (17) with aT = 1.23, and dT = 0.1. Error bars in observed and modelled |${\tilde{T}_m}$| are estimated from the time-variations of |${\tilde{T}_m}$| and from error bars in aT and dT, respectively. (b) Non-dimensional horizontally averaged heat flux, |${{\rm{\tilde{\Phi }}}_{\mathrm{ conv}}}$|⁠. Observed values are listed in Table 1, and modelled values are calculated by eq. (18) with a = 1.46, b = 0.27, c = 1.21 and d = 0.82. Error bars in observed and modelled |${{\rm{\tilde{\Phi }}}_{\mathrm{ conv}}}$| are estimated from time variations of |${{\rm{\tilde{\Phi }}}_{\mathrm{ conv}}}$| and from uncertainties on parameters of eq. (18), respectively. Note that because they are small (∼1 per cent of observed values), error bars on the observed |${{\rm{\tilde{\Phi }}}_{\mathrm{ conv}}}$| do not appear in the plot. (c) Stagnant lid thickness, |${\tilde{d}_{\mathrm{ lid}}}$|⁠. Observed values are measured from the profile of vertically advected heat and are listed in Table 1, and modelled values are calculated from eq. (23) with |${\tilde{T}_m}$| and |${{\rm{\tilde{\Phi }}}_{\mathrm{ conv}}}$| given by eqs (17) and (18), respectively. (d) Temperature at the base of the stagnant lid, |${\tilde{T}_{\mathrm{ lid}}}$|⁠. Observed values are deduced from eq. (19) together with values of |${{\rm{\tilde{\Phi }}}_{\mathrm{ conv}}}$| and |${\tilde{d}_{\mathrm{ lid}}}$| listed in Table 1, and modelled values are obtained from eq. (22) with |${\tilde{T}_m}$| given, again, by eq. (17). Figure 4. Open in new tabDownload slide Comparison between observed and modelled output properties. (a) Temperature of the well-mixed interior, |${\tilde{T}_m}$|⁠. Observed values are listed in Table 1, and modelled values are given by eq. (17) with aT = 1.23, and dT = 0.1. Error bars in observed and modelled |${\tilde{T}_m}$| are estimated from the time-variations of |${\tilde{T}_m}$| and from error bars in aT and dT, respectively. (b) Non-dimensional horizontally averaged heat flux, |${{\rm{\tilde{\Phi }}}_{\mathrm{ conv}}}$|⁠. Observed values are listed in Table 1, and modelled values are calculated by eq. (18) with a = 1.46, b = 0.27, c = 1.21 and d = 0.82. Error bars in observed and modelled |${{\rm{\tilde{\Phi }}}_{\mathrm{ conv}}}$| are estimated from time variations of |${{\rm{\tilde{\Phi }}}_{\mathrm{ conv}}}$| and from uncertainties on parameters of eq. (18), respectively. Note that because they are small (∼1 per cent of observed values), error bars on the observed |${{\rm{\tilde{\Phi }}}_{\mathrm{ conv}}}$| do not appear in the plot. (c) Stagnant lid thickness, |${\tilde{d}_{\mathrm{ lid}}}$|⁠. Observed values are measured from the profile of vertically advected heat and are listed in Table 1, and modelled values are calculated from eq. (23) with |${\tilde{T}_m}$| and |${{\rm{\tilde{\Phi }}}_{\mathrm{ conv}}}$| given by eqs (17) and (18), respectively. (d) Temperature at the base of the stagnant lid, |${\tilde{T}_{\mathrm{ lid}}}$|⁠. Observed values are deduced from eq. (19) together with values of |${{\rm{\tilde{\Phi }}}_{\mathrm{ conv}}}$| and |${\tilde{d}_{\mathrm{ lid}}}$| listed in Table 1, and modelled values are obtained from eq. (22) with |${\tilde{T}_m}$| given, again, by eq. (17). 4.2 Heat flux Thermal boundary layer analysis shows that the heat flux through a TBL scales as a power law of the Rayleigh number and of the temperature jump across the TBL (e.g. Moore & Weiss 1973). For stagnant lid convection, this implies that heat flux scales as the viscous temperature scale. The horizontally averaged non-dimensional heat flux, |${{\rm{\tilde{\Phi }}}_{\mathrm{ conv}}}$|⁠, may then be written as a function of the Rayleigh number and of the parameter γ, which controls the viscosity variations. Because heat flux is proportional to thermal conductivity, the amount of heat entering the system is partially controlled by the bottom thermal conductivity, and it is reasonable to assume that |${{\rm{\tilde{\Phi }}}_{\mathrm{ conv}}}$| also scales as this conductivity, that is, as the inverse of the top-to-bottom conductivity ratio, Rk. Following, again, Deschamps & Lin (2014), and adding the parameter Rk to account for top-to-bottom variations in thermal conductivity, |${{\rm{\tilde{\Phi }}}_{\mathrm{ conv}}}$| may be written $$\begin{eqnarray} {{\rm{\tilde{\Phi }}}_{\mathrm{ conv}}} = a\frac{{R{a_{\mathrm{ eff}}}^b}}{{{{\rm{\gamma }}^c}{R_k}^d}}, \end{eqnarray}$$(18) where Raeff is the effective Rayleigh number (eq. 11), and a, b, c and d are constants that can be deduced by inversion of observed |${{\rm{\tilde{\Phi }}}_{\mathrm{ conv}}}$|⁠. Fixing the values of a, b and c to those obtained by Deschamps & Lin (2014), that is,a = 1.46 ± 0.06, b = 0.270 ± 0.004 and c = 1.21 ± 0.03, inversion of the values of |${{\rm{\tilde{\Phi }}}_{\mathrm{ conv}}}$| listed in Table 1 using a nonlinear generalized inverse method (Tarantola & Valette 1982) leads to d = 0.82 ± 0.01. Fig. 4(b) shows that within error bars, eq. (18) together with these parameters values fit the values of |${{\rm{\tilde{\Phi }}}_{\mathrm{ conv}}}$| obtained by simulations very well. 4.3 Thickness of the stagnant lid Following a method developed by Davaille & Jaupart (1993), I inferred the thickness of the stagnant lid from the intersection between the tangent at the point of inflexion of the horizontally averaged profile of vertically advected heat, |${\tilde{v}_z}\tilde{T}$|⁠, and the origin axis |$({\tilde{v}_z}\tilde{T} = 0)$|⁠. Left plots in Fig. 3 illustrate this method, and the obtained thicknesses of the stagnant lid, |${\tilde{d}_{\mathrm{ lid}}}$|⁠, are reported in Table 1. Because in the stagnant lid heat is transported by conduction, the horizontally averaged profile of temperature in this region should be described by eq. (14), provided that D and Tbot are replaced by the thickness of the conductive lid, dlid, and the temperature at its base, Tlid. With these changes, the non-dimensional temperature profile writes $$\begin{eqnarray} \tilde{T}\left( {\tilde{z}} \right) = \frac{1}{{{R_\mathrm{ T}}}}\left\{ {{\rm{exp}}\left[ {\frac{{\tilde{z}}}{{{{\tilde{d}}_{\mathrm{ lid}}}}}{\rm{ln}}\left( {1 + {R_\mathrm{ T}}{{\tilde{T}}_{\mathrm{ lid}}}} \right)} \right] - 1} \right\}. \end{eqnarray}$$(19) Eq. (14) is of course not valid for cases with constant conductivity and should be replaced by a linear increase with depth, |$T( z ) = {T_{\mathrm{ surf}}} + \Delta Tz/D$|⁠, whose non-dimensional form applied to the stagnant lid is |$\tilde{T}( {\tilde{z}} ) = {\tilde{T}_{\mathrm{ lid}}}\tilde{z}/{\tilde{d}_{\mathrm{ lid}}}$|⁠. The horizontally averaged non-dimensional heat flux is obtained by deriving eq. (19) with respect to depth and multiplying this derivative by the non-dimensional conductivity. In Cartesian geometry, changes in the thermal gradient are compensated by changes in conductivity, such that the non-dimensional heat flux is constant with depth and equal to |${{\rm{\tilde{\Phi }}}_{\mathrm{ conv}}}$|⁠. An expression of this heat flux may be obtained from the thermal gradient at the surface |$(\tilde{z} = 0)$|⁠, where the non-dimensional conductivity is equal to 1. One then gets $$\begin{eqnarray} {{\rm{\tilde{\Phi }}}_{\mathrm{ conv}}} = \frac{{{\rm{ln}}\left( {1 + {R_T}{{\tilde{T}}_{\mathrm{ lid}}}} \right)}}{{{R_T}{{\tilde{d}}_{\mathrm{ lid}}}}}. \end{eqnarray}$$(20) This further provides a definition of |${\tilde{T}_{\mathrm{ lid}}}$| based on the measured values of |${{\rm{\tilde{\Phi }}}_{\mathrm{ conv}}}$| and |${\tilde{d}_{\mathrm{ lid}}}$|⁠, $$\begin{eqnarray} {\tilde{T}_{\mathrm{ lid}}} = \frac{{{\rm{exp}}\left( {{R_T}{{{\rm{\tilde{\Phi }}}}_{\mathrm{ conv}}}{{\tilde{d}}_{\mathrm{ lid}}}} \right) - 1}}{{{R_T}}}. \end{eqnarray}$$(21) Dashed red curves in Fig. 3 plot the conductive temperature profiles calculated from eq. (19) together with the values of |${\tilde{T}_{\mathrm{ lid}}}$| obtained from eq. (21) and the values of |${{\rm{\tilde{\Phi }}}_{\mathrm{ conv}}}$| and |${\tilde{d}_{\mathrm{ lid}}}$| listed in Table 1. These profiles are in very good agreement with the observed temperature profiles throughout the stagnant lid, except at its base, where the transition to the top thermal boundary layer occurs. This agreement validates a posteriori the method used to estimate the stagnant lid thickness. Eq. (20) can further be used to build a scaling law for the thickness of the stagnant lid. Assuming that the fluid layer beneath the lid is nearly isoviscous, and considering that the top TBL is only the bottom part of the top conductive layer (i.e. excluding the stagnant lid), the temperature jump across the top TBL should be equal to that in the bottom TBL, given by |$( {1 - {{\tilde{T}}_m}} )$|⁠. Following this approach, the value of |${\tilde{T}_{\mathrm{ lid}}}$| is $$\begin{eqnarray} {\tilde{T}_{\mathrm{ lid}}} = 2{\tilde{T}_m} - 1. \end{eqnarray}$$(22) Combining eqs (20) and (22), the thickness of the stagnant lid writes $$\begin{eqnarray} {\skew{4}\tilde{d}_{\mathrm{ lid}}} = \frac{{{\rm{ln}}\left[ {1 + {R_\mathrm{ T}}\left( {2{{\tilde{T}}_m} - 1} \right)} \right]}}{{{R_\mathrm{ T}}{{{\rm{\tilde{\Phi }}}}_{\mathrm{ conv}}}}}, \end{eqnarray}$$(23) where |${\tilde{T}_m}$| and |${{\rm{\tilde{\Phi }}}_{\mathrm{ conv}}}$| can be obtained using eqs (17) and (18). Again, eq. (23) is not valid for a system with constant thermal conductivity. In this later case, |${\tilde{d}_{\mathrm{ lid}}}$| is given by |$( {2{{\tilde{T}}_m} - 1} )/{{\rm{\tilde{\Phi }}}_{\mathrm{ conv}}}$|⁠. Despite large error bars, eq. (23) provides a very good description of the stagnant lid thickness obtained from profiles of advected heat (Fig. 4c). The agreement between the temperature at the base of the stagnant lid modelled with eq. (22) and those estimated from eq. (21) with observed |${{\rm{\tilde{\Phi }}}_{\mathrm{ conv}}}$| and |${\tilde{d}_{\mathrm{ lid}}}$| is less good but still reasonable (Fig. 4d), indicating that the system beneath the stagnant lid is close to, but not perfectly isoviscous. 5 IMPLICATIONS FOR THERMAL EVOLUTIONS OF ICY BODIES The results discussed in Sections 3 and 4 have important implications on the thermal evolutions of large solar System icy bodies, which I underline in this section. The purpose here is not to provide detailed descriptions of these evolutions, as several important complexities are not included in my modelling, but instead to assess quantitatively the role played by temperature-dependent thermal conductivity. A property shared by many, if not all, large icy bodies is the presence of a subsurface liquid ocean beneath an outer shell of low pressure water ice (phase Ih; e.g. Hussman et al. 2007). This property has been hypothesized in the early 1970s (Lewis 1971) as a consequence of the phase diagram of water ice. Since then, spacecraft missions brought observational evidences for the existence of such oceans in Europa (Khurana et al. 1998), Ganymede (Kivelson et al. 2002) and Titan (Lorenz et al. 2008). The presence of subsurface oceans may extend to smaller moons, in particularly Enceladus (Thomas et al. 2016), and to dwarf planets, including Ceres (McCord & Sotin 2005) and Pluto (Nimmo et al. 2016). In Titan and Ganymede, high pressure ice (phases III, V and VI) layers may surround the rocky core in which case the subsurface ocean would be inserted in between ice Ih and high pressure ice shells. The survival of subsurface oceans and their thickness is controlled by the thermal history of each body. Complete crystallization may be prevented by several factors, most importantly tidal dissipation within the ice layers and/or the subsurface ocean and the silicate core (e.g. Tobie et al. 2003, 2005; Choblet et al. 2017), the presence of antifreeze compounds in the initial ocean (e.g. Grasset & Sotin 1996; Spohn & Schubert 2003), which lowers the temperature at the bottom of the outer ice shells and reduce indirectly the vigour of convection within this shell (Deschamps & Sotin 2001), and the presence of insulating layers either at the top of the ice Ih shell (Tobie et al. 2006) or at the interface between the outer ice shell and the subsurface ocean (Kamata et al. 2019). Because it modifies the heat transfer through outer ice shells and the thickness of the stagnant lid, temperature-dependent thermal conductivity may further affect the evolution of icy bodies. 5.1 Modelling For applications to icy bodies, I assumed that the viscosity of ice Ih is described by $$\begin{eqnarray} \eta \left( T \right) = {\eta _{\mathrm{ ref}}}\exp\left[ {\frac{E}{{R{T_{\mathrm{ ref}}}}}\left( {\frac{{{T_{\mathrm{ ref}}}}}{T} - 1} \right)} \right], \end{eqnarray}$$(24) where E is the activation energy, R the ideal gas constant and ηref the reference viscosity at temperature Tref. The reference viscosity is not well constrained. Close to the melting point, that is, for Tref equal to the bottom temperature, Tbot, a range of values based on polar ice sheet creep is 1013–1015 Pa s (Montagnat & Duval 2000). Here, I fixed ηref to 1014 Pa s to compare shells properties (Section 5.2), and I considered this parameter as a free parameter to compare thermal evolutions (Section 5.3). Activation energy is better constrained, with values in the range 49–60 kJ mol-1 depending on the creep regime (Durham et al. 2010), and around 60 kJ mol-1 for atomic diffusion (Weertman 1983). Here, I used E = 60 kJ mol-1 in all calculations. Again, ice Ih rheology under icy moons conditions is likely more complex than the diffusion creep mechanism assumed in eq. (24) (see discussion in Section 6), but it is reasonable to think that the impact of temperature-dependent conductivity on ice shell dynamics follows a similar trend for different rheologies. Following eqs (17) and (24), the viscous temperature scale is $$\begin{eqnarray} \Delta {T_\mathrm{ v}} = \frac{{R{T_m}^2}}{E}, \end{eqnarray}$$(25) such that rescalling eq. (17) leads to a degree 2 polynomial of the temperature Tm, whose positive solution is $$\begin{eqnarray} {T_m} = \frac{E}{{2{a_T}{R_k}^{{d_T}}R}}\left( {\sqrt {1 + 4{T_{\mathrm{ bot}}}{a_\mathrm{ T}}{R_k}^{{d_\mathrm{ T}}}R/E} - 1} \right), \end{eqnarray}$$(26) where Rk = Tbot/Tsurf (eq. 10), aT = 1.23 and dT = 0.1 (Section 4.1). In the case of constant thermal conductivity, Rk is fixed to 1 independently of the value of the bottom temperature, Tbot. This bottom temperature is given either by the liquidus of water or, if impurities are present, by the liquidus of the water-impurities mixture. Impurities act as an antifreeze and may include ammonia (NH3), methanol (CH3OH), and salts (e.g. magnesium sulphate, MgSO4). Here, I more specifically considered ammonia, which has been observed in the water plumes of Enceladus (Waite et al. 2009). In the case of Europa, magnesium sulphate may further be an important compound of the ocean (Vance et al. 2018). Qualitatively, however, the impact on the evolution of the icy bodies is similar whatever the nature of the impurities. For instance, Vilella et al. (2020) pointed out that the impact of 30 per cent MgSO4 on the liquidus is equivalent to that of 3.5 per cent NH3. On another hand, it should be noted that different compositions may impact physical properties of the ocean, in particular its density. Adding 30 per cent MgSO4 would increase density by about 150 kg m−3, while 3.5 per cent NH3 would reduce it. Details on the calculation of the water–ammonia system liquidus can be found in Deschamps & Sotin (2001). Importantly, up to the eutectic composition (equal to 32.2 vol per cent in the case of NH3), only water crystalizes and impurities (here ammonia) are left in the subsurface ocean, such that their concentrations increase and their impacts are more pronounced as the ice shell thickens. Heat fluxes are obtained by rescaling the heat flux scaling law (eq. 18) with the characteristic heat flux, |${{\rm{\Phi }}_{\mathbf{ carac}}} = {k_{\mathrm{ ref}}}{\rm{\Delta }}T/D$|⁠, where kref is the characteristic thermal conductivity. Note that because the effective Rayleigh number, Raeff (eq. 11), depends on thermal conductivity (since thermal diffusivity is given by κ = k/ρCP), heat flux and the quantities derived from it further depend on the choice of kref. For temperature-dependent conductivity, numerical simulations (Section 3) implicitly assume that the kref used in Φcarac is equal to the surface conductivity, ksurf, while Raeff is calculated with the conductivity at temperature Tm, as implied by eq. (11). By contrast, for constant conductivity kref is not a priori imposed by numerical simulations, and is not only used to define Φcarac, but also Rasurf and thus Raeff, since RT = 1. To model the thermal evolutions of icy moons, most studies (e.g. Grasset & Sotin 1996; Hussman et al. 2002; Tobie et al. 2003; Mitri & Showman 2005; Běhounková et al. 2010) used values in the range 2.0–3.0 W m−1 K−1, corresponding to the conductivity at the temperature of the well mixed interior or at the bottom of the shell. These values are smaller than the surface conductivity by a factor 2–5, depending on the surface temperature. For comparison with these studies, I fixed kref to 2.6 W m−1 K−1 (Grasset & Sotin 1996; Deschamps & Lin 2014) in constant conductivity calculations, corresponding to the conductivity at a temperature of ∼ 220 K (Andersson & Suga 1994). To illustrate the potential effect of temperature-dependent conductivity, I also did additional calculations with kref = ksurf (Table 2). Accounting for the shells curvature, measured with the ratio between the inner and outer radii, f, the basal and surface heat fluxes write $$\begin{eqnarray} {{\rm{\Phi }}_{\mathrm{ surf}}} = f{{\rm{\Phi }}_{\mathrm{ carac}}}{{\rm{\tilde{\Phi }}}_{\mathrm{ conv}}} \end{eqnarray}$$(27) and $$\begin{eqnarray} {{\rm{\Phi }}_{\mathrm{ bot}}} = {{\rm{\Phi }}_{\mathrm{ carac}}}{{\rm{\tilde{\Phi }}}_{\mathrm{ conv}}}/f. \end{eqnarray}$$(28) Table 2. Europa, Pluto and materials properties. Parameter . Symbol . Unit . Value/expression . Europa . Pluto . Ice Ih properties Density ρI kg m−3 920 Thermal expansion αI K−1 1.56 × 10−4 Thermal conductivity k W m−1 K−1 566.8/T Heat capacity Cp J kg−1 K−1 7.037T + 185 Thermal diffusivity κI m2 s−1 k/ρICp Latent heat of fusion LI kJ kg−1 284 Reference bulk viscosity ηref Pa s 1012–1015 Activation energy E kJ mol−1 60 Liquid water properties Density ρI kg m−3 1000 Thermal expansion αw K−1 3.0 × 10−4 Heat capacity Cw J kg−1 K−1 4180 Silicate core properties Density ρc kg m−3 3300 Thermal diffusivity κc m2 s−1 10−6 Europa/Pluto properties Total radius R Km 1561 1188 Core radius rc Km 1400 870 Gravity acceleration g m s−2 1.31 0.62 Surface temperature Tsurf K 100 40 Surface thermal conductivity ksurf W m−1 K−1 5.7 14.2 Parameter . Symbol . Unit . Value/expression . Europa . Pluto . Ice Ih properties Density ρI kg m−3 920 Thermal expansion αI K−1 1.56 × 10−4 Thermal conductivity k W m−1 K−1 566.8/T Heat capacity Cp J kg−1 K−1 7.037T + 185 Thermal diffusivity κI m2 s−1 k/ρICp Latent heat of fusion LI kJ kg−1 284 Reference bulk viscosity ηref Pa s 1012–1015 Activation energy E kJ mol−1 60 Liquid water properties Density ρI kg m−3 1000 Thermal expansion αw K−1 3.0 × 10−4 Heat capacity Cw J kg−1 K−1 4180 Silicate core properties Density ρc kg m−3 3300 Thermal diffusivity κc m2 s−1 10−6 Europa/Pluto properties Total radius R Km 1561 1188 Core radius rc Km 1400 870 Gravity acceleration g m s−2 1.31 0.62 Surface temperature Tsurf K 100 40 Surface thermal conductivity ksurf W m−1 K−1 5.7 14.2 All data for ice Ih and liquid water properties are similar to that used by Kirk and Stevenson (1987, see references therein), except for thermal conductivity, which is from Andersson and Suga (1994), bulk viscosity, which is a free parameter with possible range of values extended from Montagnat & Duval (2000) estimates, and the activation energy, which is taken from the intermediate regime of Durham et al. (2010). Open in new tab Table 2. Europa, Pluto and materials properties. Parameter . Symbol . Unit . Value/expression . Europa . Pluto . Ice Ih properties Density ρI kg m−3 920 Thermal expansion αI K−1 1.56 × 10−4 Thermal conductivity k W m−1 K−1 566.8/T Heat capacity Cp J kg−1 K−1 7.037T + 185 Thermal diffusivity κI m2 s−1 k/ρICp Latent heat of fusion LI kJ kg−1 284 Reference bulk viscosity ηref Pa s 1012–1015 Activation energy E kJ mol−1 60 Liquid water properties Density ρI kg m−3 1000 Thermal expansion αw K−1 3.0 × 10−4 Heat capacity Cw J kg−1 K−1 4180 Silicate core properties Density ρc kg m−3 3300 Thermal diffusivity κc m2 s−1 10−6 Europa/Pluto properties Total radius R Km 1561 1188 Core radius rc Km 1400 870 Gravity acceleration g m s−2 1.31 0.62 Surface temperature Tsurf K 100 40 Surface thermal conductivity ksurf W m−1 K−1 5.7 14.2 Parameter . Symbol . Unit . Value/expression . Europa . Pluto . Ice Ih properties Density ρI kg m−3 920 Thermal expansion αI K−1 1.56 × 10−4 Thermal conductivity k W m−1 K−1 566.8/T Heat capacity Cp J kg−1 K−1 7.037T + 185 Thermal diffusivity κI m2 s−1 k/ρICp Latent heat of fusion LI kJ kg−1 284 Reference bulk viscosity ηref Pa s 1012–1015 Activation energy E kJ mol−1 60 Liquid water properties Density ρI kg m−3 1000 Thermal expansion αw K−1 3.0 × 10−4 Heat capacity Cw J kg−1 K−1 4180 Silicate core properties Density ρc kg m−3 3300 Thermal diffusivity κc m2 s−1 10−6 Europa/Pluto properties Total radius R Km 1561 1188 Core radius rc Km 1400 870 Gravity acceleration g m s−2 1.31 0.62 Surface temperature Tsurf K 100 40 Surface thermal conductivity ksurf W m−1 K−1 5.7 14.2 All data for ice Ih and liquid water properties are similar to that used by Kirk and Stevenson (1987, see references therein), except for thermal conductivity, which is from Andersson and Suga (1994), bulk viscosity, which is a free parameter with possible range of values extended from Montagnat & Duval (2000) estimates, and the activation energy, which is taken from the intermediate regime of Durham et al. (2010). Open in new tab Strictly speaking, eqs (17) and (18), which are used to estimate Tm (eq. 26) and Φsurf (eq. 27), are valid in 3-D-Cartesian geometry and should be modified for spherical systems, as ice shell curvature may slightly alter the interior temperature and heat fluxes (Yao et al. 2014). In the case of bodies like Europa, the ratio between the inner and outer radii of the shell, f, is large, around 0.90, and the impact of curvature should be limited. For thicker shells, as expected in the case of Pluto, Titan, or Ganymede, with f in the range 0.7–0.8, curvature may have a more significant impact. In the case of stagnant lid convection with constant conductivity, Yao et al. (2014) showed that dividing the second term of the right-hand side of eq. (17) by f1.5 gives good fits to the non-dimensional interior temperature provided that f is large enough, typically ≥ 0.7. Yao et al. (2014) further suggested that the top heat flux is well described by multiplying eq. (18) (with Rk = 1) by f0.22. A reappraisal of these simulations however indicates that for f ≥ 0.7 it is not needed to correct heat fluxes for the effect of curvature, that is, the top heat flux is well described by eq. (18). Additional calculations in full spherical geometry are however needed to verify that these geometrical effects are valid for stagnant lid convection with temperature-dependent thermal conductivity. Finally, the thickness of the conductive lid is either given by $$\begin{eqnarray} {d_{\mathrm{ lid}}} = \frac{{f{k_{\mathrm{ ref}}}\left( {2{T_m} - {T_{\mathrm{ bot}}} - {T_{\mathrm{ surf}}}} \right)}}{{{{\rm{\Phi }}_{\mathrm{ surf}}}}}. \end{eqnarray}$$(29) for constant conductivity, with Φsurf being given by eq. (27), or rescaled from eq. (23) for temperature-dependent conductivity, leading to $$\begin{eqnarray} {d_{\mathrm{ lid}}} = \frac{{f{k_{\mathrm{ surf}}}{T_{\mathrm{ surf}}}}}{{{{\rm{\Phi }}_{\mathrm{ surf}}}}}{\rm{ln}}\left( {\frac{{2{T_m} - {T_{\mathrm{ bot}}}}}{{{T_{\mathrm{ surf}}}}}} \right), \end{eqnarray}$$(30) where Φsurf is, again, given by eq. (27). In the following subsections, I apply this model to assess the influence of temperature-dependent thermal conductivity on the properties and evolutions of Europa and Pluto ice shells. The physical properties of these two bodies are listed in Table 2 together with the properties of ice Ih. 5.2 Ice shell properties For an ice shell of given thickness, the bottom temperature, Tbot, depends only the initial composition of the ocean and on the physical properties of the body. For reference, Fig. 5 shows Tbot as a function of the ice shell thickness for several initial compositions (pure water or initial fraction of NH3, |$x_{\mathrm{ NH3}}^{\mathrm{ init}}$|⁠) of the ocean. By definition, the top-to-bottom thermal conductivity ratio, Rk, is proportional to Tbot (eq. 10) and can thus be deduced from this temperature. Because the liquidus temperature of water decreases with pressure, Rk diminishes as the ice layer thickens. The top-to-bottom drop in thermal conductivity is further reduced by the presence of impurities, which also lowers Tbot, and this effect is strongly amplified as the ice shell thickens (i.e. as the concentration of impurities in the subsurface ocean increases). For instance, at the conditions of Pluto, conductivity within a 200 km thick ice shell drops by a factor 6.6 for a pure water ocean, and 6.3 for |$x_{\mathrm{ NH3}}^{\mathrm{ init}} = 3.0\ \mathrm{ per\, cent} $|⁠. For a 300 km thick ice shell, and assuming a pure water ocean, conductivity still decreases by a factor 6.4, while if ammonia is present the eutectic composition is reached and conductivity drops by only a factor 4.5. Finally, it is worth noting that because the surface temperature is smaller at Pluto than at Europa, conductivity changes and their potential effects are stronger at Pluto. Figure 5. Open in new tabDownload slide Bottom temperature (top row) and abundance of ammonia (NH3, bottom row) as a function of the thickness of ice Ih layer and for different input compositions (colour code). Two cases are considered, corresponding to conditions at Pluto (left column) and Europa (right column). Material and physical properties of ice Ih, Pluto and Europa are listed in Table 2. The fraction of NH3 and the bottom temperature reaches maximum and minimum values, respectively, for the eutectic composition, equal to a fraction of NH3 equal to 32.2 per cent. The top-to-bottom thermal conductivity ratio, Rk, is proportional to the bottom temperature (eq. 10), and its value is indicated on the right scale of top row plots. Figure 5. Open in new tabDownload slide Bottom temperature (top row) and abundance of ammonia (NH3, bottom row) as a function of the thickness of ice Ih layer and for different input compositions (colour code). Two cases are considered, corresponding to conditions at Pluto (left column) and Europa (right column). Material and physical properties of ice Ih, Pluto and Europa are listed in Table 2. The fraction of NH3 and the bottom temperature reaches maximum and minimum values, respectively, for the eutectic composition, equal to a fraction of NH3 equal to 32.2 per cent. The top-to-bottom thermal conductivity ratio, Rk, is proportional to the bottom temperature (eq. 10), and its value is indicated on the right scale of top row plots. Fig. 6 plots the bottom heat flux, Φbot, and stagnant lid thickness, dlid, as a function of the ice shell thickness, D, at the conditions of Pluto and Europa, and assuming a pure water composition and a reference viscosity of 1014 Pa s, corresponding to the median value estimated from Montagnat & Duval (2000). If too thin, the ice layer is not animated by convection, and heat is transported by conduction. Temperature-dependent conductivity slightly delays the onset of convection within the ice shell, that is, convection starts operating for slightly thicker shells. The bottom heat flux, Φbot, obtained for temperature-dependent conductivity is close to that predicted by constant conductivity with a characteristic conductivity, kref, fixed to 2.6 W m−1 K−1. By contrast, compared to constant conductivity calculations with kref = ksurf, temperature-dependent thermal conductivity results in a sharp reduction of Φbot. This drop increases with decreasing surface temperature, and is therefore more pronounced at Pluto, where Φbot is reduced by about a factor 5 compared to a shell with constant thermal conductivity, than at Europa, where it is reduced by only a factor 2.5. The thickness of the stagnant lid is strongly impacted by temperature-dependent conductivity, the exact effect depending on surface temperature. For Europa (Tsurf = 100 K), the lid obtained with temperature-dependent conductivity is thicker than that calculated for constant conductivity with kref = 2.6 W m−1 K−1 by a factor 1.4, in good agreement with the estimates of Tobie et al. (2003), and for Pluto (Tsurf = 40 K), it is thicker by a factor 1.8. Additional calculations show that for Tsurf = 60 K and the physical parameters of Europa, the ratio between the values of dlid calculated with temperature-dependent and constant conductivities is about 1.6, somewhat smaller than the ratio estimated by Tobie et al. (2003). The stagnant lid predicted by constant conductivity calculations thickens as kref increases, but remains thinner than that obtained with temperature-dependent conductivity even for kref = ksurf. Finally, and as expected from numerical experiments (Table 1), the temperature of the well-mixed interior is not significantly affected by temperature-dependent conductivity, and is only reduced by 1–3 K (plots c and f in Fig. 6). Figure 6. Open in new tabDownload slide Bottom convective heat flux (top row, plain curves), stagnant lid thickness (middle row, plain curves) and interior temperature (bottom row) as a function of the ice shell thickness, D, and for conditions at Pluto (left column) and Europa (right column). Dashed lines in plots (a) and (d) correspond to the conductive heat flux. Convective heat flux and stagnant lid thickness are drawn only for values of D at which convection operates. If convection does not operate, the heat flux is equal to the conductive heat flux, the entire layer is conductive, and the interior temperature is much smaller, around 156 K (for constant conductivity) and 104 K (for temperature-dependent conductivity) in the case of Pluto, and 186 and 165 K in the case of Europa. For constant conductivity, two values of the characteristic conductivity are considered, ksurf (Table 2) and 2.6 W m−1 K−1. This does not affect the interior temperature. Ocean composition is pure water. Calculations are made with the values of physical parameters listed in Table 2 and a reference viscosity ηref equal to 1014 Pa s. Note the different scales for Pluto and Europa. Figure 6. Open in new tabDownload slide Bottom convective heat flux (top row, plain curves), stagnant lid thickness (middle row, plain curves) and interior temperature (bottom row) as a function of the ice shell thickness, D, and for conditions at Pluto (left column) and Europa (right column). Dashed lines in plots (a) and (d) correspond to the conductive heat flux. Convective heat flux and stagnant lid thickness are drawn only for values of D at which convection operates. If convection does not operate, the heat flux is equal to the conductive heat flux, the entire layer is conductive, and the interior temperature is much smaller, around 156 K (for constant conductivity) and 104 K (for temperature-dependent conductivity) in the case of Pluto, and 186 and 165 K in the case of Europa. For constant conductivity, two values of the characteristic conductivity are considered, ksurf (Table 2) and 2.6 W m−1 K−1. This does not affect the interior temperature. Ocean composition is pure water. Calculations are made with the values of physical parameters listed in Table 2 and a reference viscosity ηref equal to 1014 Pa s. Note the different scales for Pluto and Europa. Adding impurities in the ocean (here, ammonia) reduces the melting temperature at the bottom of the ice shell, which has two main consequences. First, it delays the crystallization of the outer ice layer; and second, it increases its bulk viscosity, reducing the vigour of convection and the efficiency of heat transfer. This, in turn, delays the crystallization of the ice shell. Importantly, as the ice layer thickens the fraction of ammonia in the subsurface ocean increases, and its impact is enhanced. As a result, variations of Φbot and dlid as a function of D are more complex than in the case of a pure water ocean (Fig. 7). In particular, Φbot and dlid vary more sharply as the ice layer thickens. For ice layers thicker than about 220 km in the case of Pluto and 120 km in the case of Europa, a small increment in the ice layer thickness results in a large increase in the stagnant lid thickness. Finally, if the ice shell is too thick, convection stops, and heat is transported, again, by conduction. The impact of temperature-dependent conductivity is, however, globally unchanged, that is, heat flux is either similar or strongly reduced, depending on whether kref is equal to 2.6 W m−1 K−1 or ksurf, and the stagnant lid obtained with temperature-dependent conductivity is, again, thicker than that predicted by constant conductivity, by a about a factor 1.6 and 1.2 in the cases of Pluto and Europa, respectively, and for kref = 2.6 W m−1 K−1. Figure 7. Open in new tabDownload slide Bottom convective heat flux (top row, plain curve) and stagnant lid thickness (bottom row, plain curve) as a function of the ice shell thickness, D, and for conditions at Pluto (left column) and Europa (right column). Dashed lines in plots (a) and (c) correspond to the conductive heat flux. Convective heat flux and stagnant lid thickness are drawn only for values ofD at which convection operates. If convection does not operate, the heat flux is equal to the conductive heat flux, and the entire layer is conductive. For constant conductivity, two values of the characteristic conductivity are considered, ksurf (Table 2) and 2.6 W m−1 K−1. Ocean composition is a mixture of water and ammonia (NH3), with initial volume fraction of NH3 equal to 3.0 per cent. Calculations are made with the values of physical parameters listed in Table 2 and a reference viscosity ηref equal to 1014 Pa s. Note the different scales for Pluto and Europa. Figure 7. Open in new tabDownload slide Bottom convective heat flux (top row, plain curve) and stagnant lid thickness (bottom row, plain curve) as a function of the ice shell thickness, D, and for conditions at Pluto (left column) and Europa (right column). Dashed lines in plots (a) and (c) correspond to the conductive heat flux. Convective heat flux and stagnant lid thickness are drawn only for values ofD at which convection operates. If convection does not operate, the heat flux is equal to the conductive heat flux, and the entire layer is conductive. For constant conductivity, two values of the characteristic conductivity are considered, ksurf (Table 2) and 2.6 W m−1 K−1. Ocean composition is a mixture of water and ammonia (NH3), with initial volume fraction of NH3 equal to 3.0 per cent. Calculations are made with the values of physical parameters listed in Table 2 and a reference viscosity ηref equal to 1014 Pa s. Note the different scales for Pluto and Europa. 5.3 Thermal evolution The present day radial structure of icy bodies may be deduced from appropriate thermal evolution modelling. Here, I followed the approach of Grasset & Sotin (1996), which calculates the evolution of ice layers thicknesses based on an energy balance accounting for the production of heat in the silicate core, the cooling of the ocean and the crystallization of ice shells. Tidal heating, by contrast, is not included, but is certainly an important ingredient, at least in the case of Europa. Pluto and Europa are not large enough to host high pressure ices, such that the inner radius of the outer ice Ih shell, rbot, can be calculated by solving the energy conservation equation at the boundary between this shell and the subsurface ocean. Energy conservation at this boundary then writes $$\begin{eqnarray} && -\hspace{-10pt} \frac{{\mathrm{ d}{r_{\mathrm{ bot}}}}}{{\mathrm{ d}t}}\left[ {{\rho _\mathrm{ w}}{C_\mathrm{ w}}\left( { - \frac{{\partial {T_{\mathrm{ ad}}}}}{{\partial r}} + \frac{{\partial {T_{\mathrm{ bot}}}}}{{\partial r}}} \right)\frac{{\left( {{r_{\mathrm{ bot}}}^3 - {r_\mathrm{ c}}^3} \right)}}{3} - {\rho _\mathrm{ I}}{L_\mathrm{ I}}{r_{\mathrm{ bot}}}^2} \right] \\ && \quad = {r_{\mathrm{ bot}}}^2{{\rm{\Phi }}_{\mathrm{ bot}}} - {r_\mathrm{ c}}^2{{\rm{\Phi }}_\mathrm{ c}} \end{eqnarray}$$(31) where t is time, Tbot and Φbot are the temperature and heat flux at the bottom of the ice layer, given by the liquidus of the ocean and by eq. (28), respectively, rc is the core radius, Φc the heat flux at the top of the core, ρw and Cw the liquid water density and heat capacity, ρI and LI the ice Ih density and latent heat of fusion, and Tad, the adiabatic temperature in the ocean, given by $$\begin{eqnarray} {T_{\mathrm{ ad}}}\left( r \right) = {T_{\mathrm{ bot}}}\left( {{r_{\mathrm{ bot}}}} \right)\left[ {1 - \frac{{{\alpha _\mathrm{ w}}}}{{{\rho _\mathrm{ w}}{C_\mathrm{ w}}}}{\rho _\mathrm{ I}}g\left( {r - {r_{\mathrm{ bot}}}} \right)} \right], \end{eqnarray}$$(32) with αw being the thermal expansion of liquid water. Within the silicate core, heat is assumed to be produced by the decay of four radiogenic elements, 40 K, 232Th, 235U and 238U. The heat flux at the top of the core is then calculated following Kirk and Stevenson (1987) by $$\begin{eqnarray} {{\rm{\Phi }}_\mathrm{ c}} = 2\sqrt {\frac{{{\kappa _\mathrm{ c}}t}}{\pi }} {\rho _\mathrm{ c}}\mathop \sum \limits_{i = 1}^4 {C_{0,i}}{H_i}\frac{{\left[ {1 - {\rm{exp}}\left( { - {\lambda _i}t} \right)} \right]}}{{{\lambda _i}t}}, \end{eqnarray}$$(33) where κc and ρc are the thermal diffusivity and density of the silicate core, and the subscript i refers to the four radiogenic elements, whose properties are listed in Table 3. I solved eq. (31) up to t = 4.55 Gyr using an adaptative stepsize control Runge–Kutta method (Press et al. 1992), and assuming an initial ice Ih thickness equal to 10 km together with the material and physical properties listed in Table 2. Because the reference viscosity ηref is a sensitive parameter but is poorly constrained, I performed calculations for values of ηref in the range 1012–1015 Pa s, corresponding to an extended range of the values estimated by Montagnat & Duval (2000). Figs 8 and 9 show the final (t = 4.55 Gyr) bottom radius, interior temperature and stagnant lid thickness as a function of the reference viscosity, ηref, and for two possible initial compositions, while Fig. 10 plots the final bottom temperature and concentration in NH3. Fig. 11 illustrates the influence of the initial composition of the ocean on the ice shell thickness, interior temperature and stagnant lid thickness for ηref = 1014 Pa s. Table 3. Properties of radioelements. Element . Decay constant, λ . Heat release, H . Initial abundance, C0 . . (yr−1) . (W kg−1) . (ppb) . 40K 5.4279 × 10−10 2.917 × 10−5 738.0 232Th 4.9405 × 10−11 2.638 × 10−5 38.7 235U 9.8485 × 10−10 5.687 × 10−4 5.4 238U 1.5514 × 10−10 9.465 × 10−5 19.9 Element . Decay constant, λ . Heat release, H . Initial abundance, C0 . . (yr−1) . (W kg−1) . (ppb) . 40K 5.4279 × 10−10 2.917 × 10−5 738.0 232Th 4.9405 × 10−11 2.638 × 10−5 38.7 235U 9.8485 × 10−10 5.687 × 10−4 5.4 238U 1.5514 × 10−10 9.465 × 10−5 19.9 Open in new tab Table 3. Properties of radioelements. Element . Decay constant, λ . Heat release, H . Initial abundance, C0 . . (yr−1) . (W kg−1) . (ppb) . 40K 5.4279 × 10−10 2.917 × 10−5 738.0 232Th 4.9405 × 10−11 2.638 × 10−5 38.7 235U 9.8485 × 10−10 5.687 × 10−4 5.4 238U 1.5514 × 10−10 9.465 × 10−5 19.9 Element . Decay constant, λ . Heat release, H . Initial abundance, C0 . . (yr−1) . (W kg−1) . (ppb) . 40K 5.4279 × 10−10 2.917 × 10−5 738.0 232Th 4.9405 × 10−11 2.638 × 10−5 38.7 235U 9.8485 × 10−10 5.687 × 10−4 5.4 238U 1.5514 × 10−10 9.465 × 10−5 19.9 Open in new tab For a pure water ocean, and assuming a temperature-dependent conductivity, a subsurface ocean is still present after 4.55 Gyr of evolution if ηref exceeds 1014 Pa s (left columns in Figs 8 and 9). For constant conductivity, the evolution of the ice layer depends on the assumed value of the characteristic conductivity, kref. If this parameter is fixed to the surface conductivity, ksurf, crystallization is completed whatever ηref, and for both Pluto and Europa. By contrast, for kref = 2.6 W m−1 K−1 an ocean is still present for ηref ≥ 6.0 × 1013 Pa s in the case of Pluto, and ηref ≥ 1014 Pa s in the case of Europa. Given the value of ηref, the ice layer is thicker than that obtained with a temperature-dependent conductivity in case of Pluto, and slightly thinner in the case of Europa. Still for kref = 2.6 W m−1 K−1, and as expected, the stagnant lid predicted by constant conductivity is thinner than that obtained with temperature-dependent conductivity. If kref is fixed to ksurf, the stagnant lid is, again, slightly thinner. For larger values of ηref, the ice shell is thicker than that obtained with temperature-dependent conductivity, leading to thicker lids. The interior temperature calculated with constant conductivity does not depend on the value of the conductivity. As long as the ice shell thicknesses are similar (e.g. in the case of Pluto, for ηref ≤ 5.0 × 1014 Pa s), interior temperature is therefore unchanged whatever kref, and its value is slightly larger than that obtained with temperature-dependent conductivity. The presence of ammonia in the subsurface ocean prevents full crystallization for all the tested values of ηref (right columns in Figs 8 and 9; middle and bottom rows in Fig. 10). The ice shell thicknesses, bottom and interior temperatures, and abundances of NH3 at t = 4.55 Gyr obtained for temperature-dependent conductivity are comparable to those predicted by constant conductivity with kref = 2.6 W m−1 K−1. Remarkably, in the case of Europa, the ice shell thicknesses predicted by these two sets of calculations are very close to each other whatever the reference viscosity (Figs 8 and 9) and initial ocean composition (Fig. 11). Setting kref to values smaller (larger) than 2.6 W m−1 K−1 underestimates (overestimates) the ice shell thickness. For Pluto, a good agreement with the ice shell thickness predicted by temperature-dependent conductivity is obtained for kref = 3.0 W m−1 K−1. As expected, the thickness of the stagnant lid obtained with temperature-dependent conductivity, is thicker than that predicted by constant conductivity with kref = 2.6 W m−1 K−1. Note that, in the case of Pluto, part of the difference is simply related to the fact that the ice shell thickness obtained with temperature-dependent conductivity is slightly larger, implying a thicker stagnant lid (Figs 7b and d). Constant conductivity with kref = ksurf leads to systems in which convection is weak and close to shut-off. In the case of Pluto, convection even stops for ηref ≥ 2.0 × 1014 Pa s and |$x_{\mathrm{ NH3}}^{\mathrm{ init}} = 3.0\,\mathrm{ per\, cent} $| (Fig. 10), or for ηref = 1014 Pa s and |$x_{\mathrm{ NH3}}^{\mathrm{ init}}$| ≥ 3.5 per cent (Fig. 11). The ice shell is systematically thicker (i.e. the subsurface ocean is thinner) than if temperature-dependent conductivity is accounted for, around 250 km for Pluto and 90–140 km for Europa, depending on ηref. This implies high concentration in NH3 and low bottom and interior temperatures (Fig. 10), increasing in turn the shell viscosity and the top-to-bottom viscosity ratio, and leading, in fine, to very thick stagnant lids. Finally, it is interesting to note that, for cases with oceans composed of a water-ammonia mixture, the interior temperature depends very few on the initial concentration in ammonia (Fig. 11). The reduction in the liquidus temperature triggered by the increase in the fraction of ammonia is compensated by the thinning of the ice shell, which is also related to the enrichment in ammonia. For similar reasons, the thickness of the stagnant lid is only slightly impacted by the initial fraction of ammonia. 6 CONCLUDING DISCUSSION The numerical simulations performed in this study indicate that variations of thermal conductivity with temperature have a strong impact on heat transfer by thermal convection. In systems composed of material with conductivity varying as the inverse of temperature |$(k \propto 1/T)$|⁠, like ice Ih, and animated by stagnant lid convection, heat transfer is reduced by a factor Rk0.82, where Rk is the top-to-bottom conductivity ratio (eq. 8). In addition, the stagnant lid that develops at the top of the system thickens. This, in turn, alters heat transfer through the outer ice layers of solar System icy bodies, and the thermal evolution of these objects. For a given value of the bottom temperature, and because Rk is proportional to the ratio between the surface and bottom temperatures, the impact of temperature-dependent conductivity intensifies with decreasing surface temperature. Compared to constant thermal conductivity calculations with conductivity fixed to its surface value, heat transfer is reduced by about a factor 5 in the case of Pluto (Tsurf = 40 K), and by a factor 2.5 for Europa, where the surface temperature is larger (Tsurf = 100 K). Interestingly, calculations with constant conductivity predict heat fluxes comparable to those obtained with temperature-dependent conductivity, provided that the thermal conductivity is fixed to its value in the well-mixed interior or at the bottom of the shell, that is, around 2.0–3.0 W m−1 K−1. Thermal evolutions reconstructed in studies using such values (e.g. Grasset & Sotin 1996; Hussman et al. 2002; Tobie et al. 2003; Mitri & Showman 2005; Běhounková et al. 2010) should thus be mostly unaffected by the fact they do not account for the temperature dependence of thermal conductivity. However, strong differences in stagnant lid thickness remain, the stagnant lid predicted by temperature-dependent conductivity being thicker than that obtained with constant conductivity, even if conductivity is fixed to its surface value. Therefore, a detailed evolution of the icy shells structure still requires to account for temperature-dependent conductivity. The viscosity laws used in numerical simulations and their applications to icy bodies (eqs 5 and 24) implicitly assume that the deformation of ice Ih follows a diffusion creep mechanism. Under icy moons conditions, ice Ih rheology is likely more complex and would be better described by a composite viscosity taking into account, in particular, the influence of grain size (Harelet al. 2020). Note, however, that polar ice sheet creep suggests a Newtonian behaviour for low |$(\dot{\varepsilon } \le {10^{ - 11}}{\mathrm{ s}^{ - 1}})$| strain rates (Montagnat & Duval 2000). Compared to the FK approximation, the scaling laws built for a composite rheology predict smaller heat fluxes and thicker stagnant lids (Harel et al. 2020). Because temperature-dependent thermal conductivity acts on the transported heat flux but not on rheology, it is reasonable to assume that the effects observed with FK rheology are also valid with more complex rheologies, at least in their trend, that is, a reduction of the heat flux and a thickening of the stagnant lid. If Harel et al. (2020) conclusions still hold for temperature-dependent conductivity, the calculations of Section 5, quantifying the impact of temperature-dependent conductivity on heat flux and stagnant lid thickness, may be considered as upper and lower bounds, respectively, that is, efficiency of heat transfer would be further reduced and the crystallization of the subsurface ocean further delayed with a composite rheology. Detailed calculations should however be conducted to quantify more precisely the combined effects of complex rheology and temperature-dependent thermal conductivity. As demonstrated in many studies (e.g. Grasset & Sotin 1996; Deschamps & Sotin 2001; Spohn & Schubert 2003), the presence of antifreeze compounds delays or stops the crystallization of the ice layer. As long as the eutectic composition is not reached, only water ice crystallizes. Once this composition is reached impurities crystallize together with water. This would result in a compositional change at the bottom of the shell, which may, in turn, alter the thermal conductivity at the bottom of the ice layer, further impacting the ocean crystallization. For instance, experimental data on methanol show that the presence of 10 wt per cent methanol reduces the thermal conductivity of ice by a factor 2 (Hsieh & Deschamps 2015). Similarly, the conductivity of methane clathrates is lower than that of water ice by about a factor 5–10, and the presence of a thin clathrates layer at the top or at the bottom of the ice shell may act as an efficient insulator (Tobie et al. 2006; Kamata et al. 2019). Another important factor limiting the crystallization of subsurface oceans is a regular production of heat, such as that provided by tidal dissipation within icy bodies (e.g. Tobie et al. 2003, 2005; Roberts & Nimmo 2008). Tidal dissipation within the ocean or within the core increases the heat flux available at the bottom of the ice shell. Tidal dissipation within the ice shell itself reduces the plumes strength and the convective heat flux at the bottom of the shell (e.g. Travis & Olson 1994; Deschamps et al. 2010b). In both cases removing heat from the subsurface ocean is more challenging, in particular if tidal heating is combined with stagnant lid convection and temperature-dependent thermal conductivity, which both reduces the conductive heat flux. Temperature-dependent conductivity may, in particular, lower the threshold value of dissipated heat for which the bottom heat flux at the bottom of the ice shell turns negative (i.e. the ice layer heats up the ocean and thins). In addition, because it increases the ratio between internal and basal sources of heat and leads to thicker stagnant lids, temperature-dependent conductivity may favour the existence of pockets of partially melted material close to the surface (Vilella et al. 2020), an ingredient that is essential to explain past and recent cryo-volcanism at the surfaces of icy bodies. Beside the limitations discussed in this section, temperature-dependent thermal conductivity certainly impacts the evolution of large icy bodies. Because the thermal conductivity of silicated rocks also varies as 1/Tn, with n depending on the mineralogical composition, it may further play a, somewhat smaller, role in heat transfer within rocky mantles. For iron-bearing minerals, n is likely around 1/2 (Klemens et al. 1962), thus moderating the top-to-bottom thermal conductivity ratio. Temperature-dependent thermal conductivity may however impact heat transfer in the mantles of the Moon and Mercury. Note that in the case of Mercury, the high surface temperature could limit the top-to-bottom reduction in thermal conductivity. For larger bodies, such as Mars and Venus, pressure effects should be accounted for. Laboratory experiments then indicate that thermal conductivities of the main lower-mantle minerals, bridgmanite and ferropericlase, increase with depth (e.g. Ohta et al. 2012; Goncharov et al. 2015; Hsieh et al. 2017, 2018), and that the combined effects of temperature and pressure along an adiabat lead to an increase of thermal conductivity by about a factor 2 at the bottom of the Earth's mantle (Hsieh et al. 2018). If convection operates, the top-to-bottom temperature jump is super-adiabatic, and the bottom conductivity should be reduced compared to the adiabatic case. For instance, assuming a surface temperature of 700 K (typical of Venus), a super-adiabatic jump of 2000 K, and n = 1/2, the thermal conductivity should be divided by 2, thus cancelling the pressure effect (which, in the case of Venus, should be similar to that of the Earth's mantle). Additional calculations, including simulations of convection with temperature and pressure dependent thermal conductivity, are however needed to assess more precisely the impact of heterogeneous thermal conductivity on heat transfer within rocky mantles. ACKNOWLEDGEMENTS The research presented in this paper was supported by the National Science Council of Taiwan (MOST) under grants 107–2116-M-001–010 and 108–2116-M-001–017. The author is grateful to Shijie Zhong, to an anonymous colleague and to the editor (Gaël Choblet) for their detailed comments and suggestions that helped improving the first version of this paper. The code used in this work to perform numerical simulations is not publicly available but was thoroughly described in Tackley (2008). The modifications made to the code are detailed in this article. All the data presented in this article was generated with this code and are available on demand to the author (frederic@earth.sinica.edu.tw), as well as the data used to generate the figures. REFERENCES Andersson O. , Suga H., 1994 . Thermal conductivity of the Ih and XI phases of ice , Phys. Rev. B , 50 , 6583 – 6588 . 10.1103/PhysRevB.50.6583 Google Scholar Crossref Search ADS WorldCat Běhounkova M. , Tobie G., Choblet G., Čadek O., 2010 . Coupling mantle convection and tidal dissipation: applications to Enceladus and Earth‐like planets , J. geophys. Res. , 115 , doi:10.1029/2009JE003564 . 10.1029/2009JE003564 Google Scholar OpenURL Placeholder Text WorldCat Crossref Choblet G. , Tobie G., Sotin C., Běhounková M., Čadek O., Postberg, Souček O., 2017 . Powering prolonged hydrothermal activity inside Enceladus , Nat. Astron. , 12 , 841 – 847 . 10.1038/s41550-017-0289-8 Google Scholar Crossref Search ADS WorldCat Christensen U.R. , 1984 . Heat transport by variable viscosity convection and implications for the Earth's thermal evolution , Phys. Earth planet. Inter. , 35 , 264 – 282 . 10.1016/0031-9201(84)90021-9 Google Scholar Crossref Search ADS WorldCat Davaille A. , Jaupart C., 1993 . Transient high-Rayleigh-number thermal convection with large viscosity variations , J. Fluid Mech. , 253 , 141 – 166 . 10.1017/S0022112093001740 Google Scholar Crossref Search ADS WorldCat Deschamps F. , Lin J.-R., 2014 . Stagnant lid convection in 3D-Cartesian geometry: scaling laws and applications to icy moons and dwarf planets , Phys. Earth planet. Inter. , 229 , 40 – 54 . 10.1016/j.pepi.2014.01.002 Google Scholar Crossref Search ADS WorldCat Deschamps F. , Mousis O., Sanchez-Valle C., Lunine J.I., 2010a . The role of methanol on the crystallization of Titan's primordial ocean , Astrophys. J. , 724 , 887 – 894 . 10.1088/0004-637X/724/2/887 Google Scholar Crossref Search ADS WorldCat Deschamps F. , Sotin C., 2000 . Inversion of two-dimensional numerical convection experiments for a fluid with a strongly temperature-dependent viscosity , Geophys. J. Int. , 143 , 204 – 218 . 10.1046/j.1365-246x.2000.00228.x Google Scholar Crossref Search ADS WorldCat Deschamps F. , Sotin C., 2001 . Thermal convection in the outer shell of large icy satellites , J. geophys. Res. , 106 , 5107 – 5121 . 10.1029/2000JE001253 Google Scholar Crossref Search ADS WorldCat Deschamps F. , Tackley P.J., Nakagawa , 2010b . Temperature and heat flux scalings for isoviscous thermal convection in spherical geometry , Geophys. J. Int. , 182 , 137 – 154 . Google Scholar OpenURL Placeholder Text WorldCat Durham W.B. , Prieto-Ballesteros O., Goldsby D.L., Kargel J.S., 2010 . Rheological and thermal properties of icy minerals , Space Sci. Rev. , 153 , 273 – 298 . 10.1007/s11214-009-9619-1 Google Scholar Crossref Search ADS WorldCat Goncharov A.F. et al. , 2015 . Experimental study of thermal conductivity at high pressures: Implications for the deep Earth's interior , Phys. Earth planet. Inter. , 247 , 11 – 16 . 10.1016/j.pepi.2015.02.004 Google Scholar Crossref Search ADS WorldCat Grasset O. , Sotin C., 1996 . The cooling rate of a liquid shell in Titan's interior , Icarus , 123 , 101 – 112 . 10.1006/icar.1996.0144 Google Scholar Crossref Search ADS WorldCat Harel L. , Dumoulin C., Choblet G., Tobie G., Besserer J., 2020 . Scaling of heat transfer in stagnant lid convection for the outer ice shells of icy moons: influence of rheology , Icarus , 338 , 113448 . 10.1016/j.icarus.2019.113448 Google Scholar Crossref Search ADS WorldCat Hsieh W.-P. , Deschamps F., 2015 . Thermal conductivity of H2O-CH3OH mixtures at high pressures: implications for the dynamics of icy super-Earths outer shells , J. geophys. Res. , 120 , 1697 – 1707 . 10.1002/2015JE004883 Google Scholar Crossref Search ADS WorldCat Hsieh W.-P. , Deschamps F., Okuchi T., Lin J.F., 2017 . Reduced lattice thermal conductivity of Fe-bearing bridgmanite in Earth's deep mantle , J. geophys. Res. , 122 , 4900 – 4917 . 10.1002/2017JB014339 Google Scholar Crossref Search ADS WorldCat Hsieh W.-P. , Deschamps F., Okuchi T., Lin J.F., 2018 . Effects of iron on the lattice thermal conductivity and dynamics of Earth's deep mantle , Proc. Natl. Acad. Sci. USA , 115 , 4099 – 4104 . 10.1073/pnas.1718557115 Google Scholar Crossref Search ADS WorldCat Hussmann H. , Sotin C., Lunine J.I., 2007 . Interiors and evolution of icy satellites , in Planets and Moons , Treatise on Geophysics, Vol. 10, pp. 509 – 539 ., ed. Spohn T., Elsevier . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Hussmann H. , Spohn T., Wieczerkowshi K., 2002 . Thermal equilibrium states of Europa's ice shell: implications for internal ocean thickness and surface heat flow , Icarus , 156 , 143 – 151 . 10.1006/icar.2001.6776 Google Scholar Crossref Search ADS WorldCat Kamata S. , Nimmo F., Sekine Y., Kuramoto K., Nogushi N., Kimura J., Tani A., 2019 . Pluto's ocean is caped and insulated by gas hydrates , Nat. Geosci. , 12 , 407 – 410 . 10.1038/s41561-019-0369-8 Google Scholar Crossref Search ADS WorldCat Khurana K.K. , Kivelson M.G., Stevenson D.J., Schubert G., Russell C.T., Walker R.J., Polanskey C., 1998 . Induced magnetic field as evidence for subsurface ocean in Europa and Callisto , Nature , 395 , 777 – 780 . 10.1038/27394 Google Scholar Crossref Search ADS PubMed WorldCat Kirk R.L. , Stevenson D.J., 1987 . Thermal evolution of a differentiated Ganymede and implications for surface features , Icarus , 69 , 91 – 134 . 10.1016/0019-1035(87)90009-1 Google Scholar Crossref Search ADS WorldCat Kivelson M.G. , Khurana K.K., Volwerk M., 2002 . The permanent and inductive magnetic moments of Ganymede , Icarus , 157 , 507 – 522 . 10.1006/icar.2002.6834 Google Scholar Crossref Search ADS WorldCat Klemens P.G. , White G.K., Tainsh R.J., 1962 . Scattering of lattice waves by point defects , Phil. Mag. , 7 , 1323 – 1335 . 10.1080/14786436208213166 Google Scholar Crossref Search ADS WorldCat Lewis J.S. , 1971 . Satellites of the outer planets: their physical and chemical nature , Icarus , 15 , 174 – 185 . 10.1016/0019-1035(71)90072-8 Google Scholar Crossref Search ADS WorldCat Lorenz R.D. et al. , 2008 . Titan's rotation reveals an internal ocean and changing zonal winds , Science , 319 , 1649 – 1651 . 10.1126/science.1151639 Google Scholar Crossref Search ADS PubMed WorldCat McCord T.B. , Sotin C., 2005 . Ceres: evolution and current state , J. geophys. Res. , 110 , E05009 . Google Scholar OpenURL Placeholder Text WorldCat Mitri G. , Showman A.P., 2005 . Convective-conductive transitions and sensitivity of a convecting ice shell to perturbations in heat flux and tidal-heating rate: implications for Europa , Icarus , 177 , 447 – 460 . 10.1016/j.icarus.2005.03.019 Google Scholar Crossref Search ADS WorldCat Montagnat M. , Duval P., 2000 . Rate controlling processes in the creep of polar ice, influence of grain boundary migration associated with recrystallization , Earth planet. Sci. Lett. , 183 , 179 – 186 . 10.1016/S0012-821X(00)00262-4 Google Scholar Crossref Search ADS WorldCat Moore D.R. , Weiss N.O., 1973 . Two-dimensional Rayleigh–Bénard convection , J. Fluid Mech. , 58 , 289 – 312 . 10.1017/S0022112073002600 Google Scholar Crossref Search ADS WorldCat Moresi L.-N. , Solomatov V.S., 1995 . Numerical investigation of 2D convection with extremely large viscosity variations , Phys. Fluids , 7 , 2154 – 2162 . 10.1063/1.868465 Google Scholar Crossref Search ADS WorldCat Nimmo F. et al. , 2016 . Reorientation of Sputnik Planitia implies a subsurface ocean on Pluto , Nature , 540 , 94 – 96 . 10.1038/nature20148 Google Scholar Crossref Search ADS PubMed WorldCat Ohta K. , Yagi T., Taketoshi N., Hirose K., Komabayashi T., Baba T., Ohishi Y., Hernlund J., 2012 . Lattice thermal conductivity of MgSiO3 perovskite and post-perovskite at the core–mantle boundary , Earth planet. Sci. Lett. , 349–350 , 109 – 115 . 10.1016/j.epsl.2012.06.043 Google Scholar Crossref Search ADS WorldCat Press W.H. , Flannery B.P., Teukolsky S.A., Vetterling W.T., 1992 . Numerical Recipes , 2nd edn, Cambridge Univ. Press , pp. 701 – 725 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Reese C.C. , Solomatov V.S., Moresi L.N., 1999 . Non-Newtonian stagnant lid convection and magmatic resurfacing on Venus , Icarus , 139 , 67 – 80 . 10.1006/icar.1999.6088 Google Scholar Crossref Search ADS WorldCat Roberts J. H. , Nimmo F., 2008 . Tidal heating and the long-term stability of a subsurface ocean on Enceladus , Icarus , 194 , 675 – 689 . 10.1016/j.icarus.2007.11.010 Google Scholar Crossref Search ADS WorldCat Slack G.A. , 1980 . Thermal conductivity of ice , Phys. Rev. B , 22 , 3065 – 3071 . 10.1103/PhysRevB.22.3065 Google Scholar Crossref Search ADS WorldCat Spohn T. , Schubert G., 2003 . Oceans in the icy Galilean satellites of Jupiter? . Icarus , 161 , 456 – 467 . 10.1016/S0019-1035(02)00048-9 Google Scholar Crossref Search ADS WorldCat Stein C. , Lowman J.P., Hansen U., 2013 . The influence of mantle internal heating on lithospheric mobility: implications for super-Earths , Earth planet. Sci. Lett. , 361 , 448 – 459 . 10.1016/j.epsl.2012.11.011 Google Scholar Crossref Search ADS WorldCat Tackley P.J. , 2008 . Modelling compressible mantle convection with large viscosity contrasts in a three-dimensional spherical shell using the yin-yang grid , Phys. Earth planet. Inter. , 171 , 7 – 18 . 10.1016/j.pepi.2008.08.005 Google Scholar Crossref Search ADS WorldCat Tarantola A. , Valette B., 1982 . Generalized nonlinear inverse problems solved using the least square criterion , Rev. Geophys. Space Phys. , 20 , 219 – 232 . 10.1029/RG020i002p00219 Google Scholar Crossref Search ADS WorldCat Thomas P.C. , Tajeddine R., Tiscareno M.S., Burns J.A., Joseph J., Loredo T.J., Helfeinstein P., Porco C., 2016 . Enceladus's measured physical libration requires a global subsurface ocean , Icarus , 264 , 37 – 47 . 10.1016/j.icarus.2015.08.037 Google Scholar Crossref Search ADS WorldCat Tobie G. , Choblet G., Sotin C., 2003 . Tidally heated convection: constraints on Europa's ice shell thickness , J. geophys. Res. , 108 , doi:10.1029/2003JE002099 . 10.1029/2003JE002099 Google Scholar OpenURL Placeholder Text WorldCat Crossref Tobie G. , Lunine J.I., Sotin C., 2006 . Episodic outgassing as the origin of atmospheric methane on Titan , Nature , 440 , 61 – 64 . 10.1038/nature04497 Google Scholar Crossref Search ADS PubMed WorldCat Tobie G. , Mocquet A., Sotin C., 2005 . Tidal dissipation within large icy satellites: applications to Europa and titan , Icarus , 177 , 534 – 549 . 10.1016/j.icarus.2005.04.006 Google Scholar Crossref Search ADS WorldCat Travis B. , Olson P., 1994 . Convection with internal sources and turbulence in the Earth's mantle , Geophys. J. Int. , 118 , 1 – 19 . 10.1111/j.1365-246X.1994.tb04671.x Google Scholar Crossref Search ADS WorldCat Vance S.D. , Panning M.P., Stähler S., 2018 . Geophysical investigations of habitability in ice-covered ocean worlds , J. geophys. Res. , 123 , 180 – 205 . 10.1002/2017JE005341 Google Scholar Crossref Search ADS WorldCat Vilella K. , Choblet G., Tsao W.E., Deschamps F., 2020 . Tidally heated convection and the occurrence of melting in icy satellites: application to Europa , J. geophys. Res. , 125 , e2019JE006248 , doi:10.1029/2019JE006248 . 10.1029/2019JE006248 Google Scholar Crossref Search ADS WorldCat Waite J.-H. , Lewis W.S., Maggee B.A., 2009 . Liquid water on Enceladus from observations of ammonia and 40Ar in the plumes , Nature , 460 , 487 – 490 . 10.1038/nature08153 Google Scholar Crossref Search ADS WorldCat Weertman J. , 1983 . Creep deformation of ice , Annu. Rev. Earth Planet. Sci. , 11 , 215 – 240 . 10.1146/annurev.ea.11.050183.001243 Google Scholar Crossref Search ADS WorldCat Yao C. , Deschamps F., Lowman J.P., Sanchez-Valle C., Tackley P.J., 2014 . Stagnant-lid convection in bottom-heated thin 3-D spherical shells: influence of curvature and implications for dwarf planets and icy moons , J. geophys. Res. , 119 , 1895 – 1913 . 10.1002/2014JE004653 Google Scholar Crossref Search ADS WorldCat APPENDIX A: RADIAL PROFILES OF TEMPERATURE AND HEAT FLUX FOR TEMPERATURE-DEPENDENT THERMAL CONDUCTIVITY Experimental data show that the thermal conductivity of ice Ih depends on temperature following an inverse law with a very good approximation (e.g. Slack 1980; Andersson & Suga (1994). Taking the conductivity ksurf at surface temperature Tsurf as reference, the conductivity at temperature T writes $$\begin{eqnarray} k \left( T \right) = {k_{\mathrm{ surf}}} \frac{{{T_{\mathrm{ surf}}}}}{T}. \end{eqnarray}$$(A1) This dependence modifies the radial profiles of temperature and heat flux for a system in which heat is transported by conduction. A1 Cartesian geometry In Cartesian geometry, the steady state 1-D conduction equation is $$\begin{eqnarray} \frac{\partial }{{\partial z}} \left( {k\frac{{\partial T}}{{\partial z}}} \right) = 0, \end{eqnarray}$$(A2) where z is depth. Replacing k with its expression in eq. (A1), the conduction equation is $$\begin{eqnarray} \frac{\partial }{{\partial z}} \left( {{k_{\mathrm{ surf}}}\frac{{{T_{\mathrm{ surf}}}}}{T}\frac{{\partial T}}{{\partial z}}} \right) = 0. \end{eqnarray}$$(A3) Double integration of eq. (A3) with respect to z leads to $$\begin{eqnarray} \ln \left( T \right) = \frac{{{c_1}}}{{{k_{\mathrm{ \mathrm{ surf}}}}{T_{\mathrm{ surf}}}}} z + {c_2} \end{eqnarray}$$ and $$\begin{eqnarray} T\left( z \right) = {T_0}\ \exp\left( {\frac{{{c_1}}}{{{k_{\mathrm{ surf}}}{T_{\mathrm{ surf}}}}}z} \right), \end{eqnarray}$$(A4) where c1 and T0 = ln(c2) are integration constants. As usual, these constants may be deduced from the values of the temperature and its derivative at the interfaces of the system. At the surface, z = 0, T = Tsurf, and therefore T0 = Tsurf. Constant c1 may be defined either from the bottom temperature Tbot, or from the surface heat flux, Φsurf. In the case of ice Ih shells, the bottom temperature is known a priori because it is, by definition, equal to the liquidus of water (or of the water-impurities mixture). For a shell of thickness zbot = D, and reminding that T0 = Tsurf, eq. (A4) gives $$\begin{eqnarray} {c_1} = \frac{{{k_{\mathrm{ surf}}}{T_{\mathrm{ surf}}}}}{D}{\rm{\ ln}}\left( {\frac{{{T_{\mathrm{ bot}}}}}{{{T_{\mathrm{ surf}}}}}} \right), \end{eqnarray}$$(A5) and thus $$\begin{eqnarray} T\left( z \right) = {T_{\mathrm{ surf}}} \exp\left[ {\frac{z}{D}{\rm{ln}}\left( {\frac{{{T_{\mathrm{ bot}}}}}{{{T_{\mathrm{ surf}}}}}} \right)} \right]. \end{eqnarray}$$(A6) Strictly speaking the liquidus temperature is given as a function of pressure, which at depth D is given by P = ρicegD, where g is gravity acceleration and ρice ice density. Heat flux simply derives from eq. (A6) following $$\begin{eqnarray} {\rm{\Phi }}\left( z \right) = k \frac{{\partial T}}{{\partial z}} = {k_{\mathrm{ surf}}}\ \frac{{{T_{\mathrm{ surf}}}}}{T}\frac{{\partial T}}{{\partial z}}, \end{eqnarray}$$ which leads to $$\begin{eqnarray} {\rm{\Phi }}\left( z \right) = \frac{{{k_{\mathrm{ surf}}}{T_{\mathrm{ surf}}}}}{D}{\rm{\ ln}}\left( {\frac{{{T_{\mathrm{ bot}}}}}{{{T_{\mathrm{ surf}}}}}} \right). \end{eqnarray}$$(A7) As expected, in Cartesian geometry Φ is constant with depth. A comparison with constant c1 (eq. A5) then indicates that c1 = Φ. If surface heat flux, rather than bottom temperature, is known, the temperature profile may then be written $$\begin{eqnarray} T\left( z \right) = {T_{\mathrm{ surf}}}\exp\left( {\frac{{{{\rm{\Phi }}_{\mathrm{ surf}}}}}{{{k_{\mathrm{ surf}}}{T_{\mathrm{ surf}}}}}z} \right). \end{eqnarray}$$(A8) Eq. (A8) is useful to determine the radial profile of temperature in the top part of systems animated by stagnant lid convection, where a rigid, thermally conductive lid is developing at the top of the system. It allows, in particular, to derive the thickness of the conductive lid from the surface heat flux (Section 4.3). A2 Spherical geometry A similar reasoning can be followed in spherical geometry, except that the steady state conduction equation is given by $$\begin{eqnarray} \frac{1}{{{r^2}}}\frac{\partial }{{\partial r}}\ \left( {{r^2}k\frac{{\partial T}}{{\partial r}}} \right) = \ 0, \end{eqnarray}$$(A9) where r is radius. First and second integrations with respect to r, after replacing k with its expression in eq. (A1) gives $$\begin{eqnarray} \frac{1}{T} \frac{{\partial T}}{{\partial r}} = \frac{{{c_1}}}{{{k_{\mathrm{ surf}}}{T_{\mathrm{ surf}}}}}\ \frac{1}{{{r^2}}}, \end{eqnarray}$$ and $$\begin{eqnarray} \ln \left( T \right) = \ - \frac{{{c_1}}}{{{k_{\mathrm{ surf}}}{T_{\mathrm{ surf}}}}}\frac{1}{r} + {c_2}. \end{eqnarray}$$(A10) At surface, r = R, where R is the total radius, T = Tsurf, and therefore |${c_2} = \ln ( {{T_{\mathrm{ surf}}}} )\ + \frac{{{c_1}}}{{{k_{\mathrm{ surf}}}{T_{\mathrm{ surf}}}R}}$|⁠. Temperature is then given by $$\begin{eqnarray} T\left( r \right) = {T_{\mathrm{ surf}}}\exp\left[ {\frac{{{c_1}}}{{{k_{\mathrm{ surf}}}{T_{\mathrm{ surf}}}R}}\left( {1 - \frac{R}{r}} \right)} \right]. \end{eqnarray}$$(A11) Again, the constant c1 can be either defined from temperature at the bottom of the shell (i.e. at radius rbot), or from the surface heat flux (note that in spherical geometry, heat flux depends on radius). Assuming that the bottom temperature is known, one gets $$\begin{eqnarray} {c_1} = - \frac{{{k_{\mathrm{ surf}}}{T_{\mathrm{ surf}}}}}{D}f{R^2}{\rm{ln}}\left( {\frac{{{T_{\mathrm{ bot}}}}}{{{T_{\mathrm{ surf}}}}}} \right), \end{eqnarray}$$(A12) where |$f\ = {r_{\mathrm{ bot}}}/R$| is the ratio between the inner and outer radii of the shell, and D = (R—rbot) = R(1–f) is the ice shell thickness. The radial temperature profile is then $$\begin{eqnarray} T\left( r \right) = {T_{\mathrm{ surf}}}\exp\left[ { - \frac{{Rf}}{D}\left( {1 - \frac{R}{r}} \right){\rm{ln}}\left( {\frac{{{T_{\mathrm{ bot}}}}}{{{T_{\mathrm{ surf}}}}}} \right)} \right]. \end{eqnarray}$$(A13) The heat flux at radius r derives form eq. (A13) following |${\rm{\Phi }}( r ) = - k\partial T/\partial r$|⁠, leading to $$\begin{eqnarray} {\rm{\Phi }}\left( r \right) = \frac{{{k_{\mathrm{ surf}}}{T_{\mathrm{ surf}}}}}{D}\ f{\left( {\frac{R}{r}} \right)^2}{\rm{ln}}\left( {\frac{{{T_{\mathrm{ bot}}}}}{{{T_{\mathrm{ surf}}}}}} \right). \end{eqnarray}$$(A14) Its surface value is given by $$\begin{eqnarray} {{\rm{\Phi }}_{\mathrm{ surf}}} = \frac{{{k_{\mathrm{ surf}}}{T_{\mathrm{ surf}}}}}{D}\ f{\rm{ln}}\left( {\frac{{{T_{\mathrm{ bot}}}}}{{{T_{\mathrm{ surf}}}}}} \right), \end{eqnarray}$$(A15) and can be used to define the constant c1 (eq. A12), which is then given by c1 = −ΦsurfR2. Note also that the heat flux at any radius can be deduced from its surface value by |${\rm{\Phi }}( r ) = {{\rm{\Phi }}_{\mathrm{ surf}}}{R^2}/{r^2}$|⁠. Finally, as in Cartesian geometry, the radial temperature profile may be defined as a function of the surface heat flux following $$\begin{eqnarray} T\left( r \right) = {T_{\mathrm{ surf}}}\exp\left[ { - \frac{{{{\rm{\Phi }}_{\mathrm{ surf}}}R}}{{{k_{\mathrm{ surf}}}{T_{\mathrm{ surf}}}}}\left( {1 - \frac{R}{r}} \right)} \right]. \end{eqnarray}$$(A16) A3 Conductive profiles Table A1 summarizes the relationships for radial profiles of temperature and heat flux in Cartesian and spherical geometries and for a thermal-conductivity given by eq. (A1), together with the corresponding relationships for constant conductivity. Fig. A1 compares the radial conductive profiles of temperature and heat flux obtained with these expressions for two cases with different surface temperature and ice shell thickness, Tsurf = 40 K and D = 300 km, representative of Pluto, and Tsurf = 100 K and D = 150 km, representative of Europa. In both cases, there is no internal heating, and the bottom temperature is set to Tbot = 250 K, a value close to the liquidus of water at a pressure of 0.21 GPa. The heat flux profiles for constant conductivity are calculated either with the surface value of the thermal conductivity (plain curves) or with k = 2.6 W m−1 K−1 (dashed curves), a value typical of the expected conductivity within the interior of icy bodies outer ice shells. Clearly, temperature-dependent conductivity has a strong impact on the conductive thermal structure and heat flux. Compared to the constant conductivity case, temperature is smaller by up to ∼40 K and ∼15 K in the cases with Tsurf = 40 K and Tsurf = 100 K, respectively. If thermal conductivity is fixed to its surface value, the heat flux is reduced by a factor 3 in the case with Tsurf = 40 K, and by a factor 2 in the case with Tsurf = 100 K. By contrast, the heat flux obtained for temperature-dependent conductivity is slightly larger, but comparable to that obtained with constant conductivity and k = 2.6 W m−1 K−1. Figure A1. Open in new tabDownload slide Radial profiles of temperature (left column) and heat flux (right column) for a purely conductive heat transfer with either constant (k-cte) or temperature-dependent (k-var) thermal conductivity. For heat flux with constant conductivity, two values of the thermal conductivity are used: the surface conductivity ksurf listed in Table 2 (plain curves), and 2.6 W m−1 K−1 (dashed curves). In temperature-dependent conductivity cases, conductivity varies as 1/T (eq. A1). Calculations are made for surface temperature and ice shell thickness Tsurf = 40 K and D = 300 km (top row), and Tsurf = 100 K and D = 150 km (bottom row), and two geometries are considered, Cartesian (C) or spherical (sph). In all calculations, the rate of internal heating is set to 0, and the bottom temperature, Tbot, is fixed to 250 K. Relationships used to calculate these profiles are listed in Table A1. Figure A1. Open in new tabDownload slide Radial profiles of temperature (left column) and heat flux (right column) for a purely conductive heat transfer with either constant (k-cte) or temperature-dependent (k-var) thermal conductivity. For heat flux with constant conductivity, two values of the thermal conductivity are used: the surface conductivity ksurf listed in Table 2 (plain curves), and 2.6 W m−1 K−1 (dashed curves). In temperature-dependent conductivity cases, conductivity varies as 1/T (eq. A1). Calculations are made for surface temperature and ice shell thickness Tsurf = 40 K and D = 300 km (top row), and Tsurf = 100 K and D = 150 km (bottom row), and two geometries are considered, Cartesian (C) or spherical (sph). In all calculations, the rate of internal heating is set to 0, and the bottom temperature, Tbot, is fixed to 250 K. Relationships used to calculate these profiles are listed in Table A1. Table A1. Relationships for radial profiles of temperature and heat flux for a conductive system. Quantity . Geometry . Expression . k = constant Temperature Cartesian |${T_{\mathrm{ surf}}} + \Delta T\frac{z}{D}$| - Spherical |${T_{\mathrm{ surf}}} - \Delta T\frac{R}{D}f( {1 - \frac{R}{r}} )$| Heat flux Cartesian |$k\frac{{\Delta T}}{D}$| - Spherical |$k\frac{{\Delta T}}{D}f{( {\frac{R}{r}} )^2}$| |$k\ ( T ) = {k_{\mathrm{ surf}}}{T_{\mathrm{ surf}}}/T$| Temperature Cartesian |${T_{\mathrm{ surf}}}\exp[ {\frac{z}{D}{\rm{ln}}( {\frac{{{T_{\mathrm{ bot}}}}}{{{T_{\mathrm{ surf}}}}}} )} ]$| - Spherical |${T_{\mathrm{ surf}}}\exp[ { - \frac{R}{D}f( {1 - \frac{R}{r}} ){\rm{ln}}( {\frac{{{T_{\mathrm{ bot}}}}}{{{T_{\mathrm{ surf}}}}}} )} ]$| Heat flux Cartesian |$\frac{{{k_{\mathrm{ surf}}}{T_{surf}}}}{D}{\rm{ln}}( {\frac{{{T_{\mathrm{ bot}}}}}{{{T_{\mathrm{ surf}}}}}} )$| - Spherical |$\frac{{{k_{\mathrm{ surf}}}{T_{\mathrm{ surf}}}}}{D}f{( {\frac{R}{r}} )^2}{\rm{ln}}( {\frac{{{T_{\mathrm{ bot}}}}}{{{T_{\mathrm{ surf}}}}}} )$| Quantity . Geometry . Expression . k = constant Temperature Cartesian |${T_{\mathrm{ surf}}} + \Delta T\frac{z}{D}$| - Spherical |${T_{\mathrm{ surf}}} - \Delta T\frac{R}{D}f( {1 - \frac{R}{r}} )$| Heat flux Cartesian |$k\frac{{\Delta T}}{D}$| - Spherical |$k\frac{{\Delta T}}{D}f{( {\frac{R}{r}} )^2}$| |$k\ ( T ) = {k_{\mathrm{ surf}}}{T_{\mathrm{ surf}}}/T$| Temperature Cartesian |${T_{\mathrm{ surf}}}\exp[ {\frac{z}{D}{\rm{ln}}( {\frac{{{T_{\mathrm{ bot}}}}}{{{T_{\mathrm{ surf}}}}}} )} ]$| - Spherical |${T_{\mathrm{ surf}}}\exp[ { - \frac{R}{D}f( {1 - \frac{R}{r}} ){\rm{ln}}( {\frac{{{T_{\mathrm{ bot}}}}}{{{T_{\mathrm{ surf}}}}}} )} ]$| Heat flux Cartesian |$\frac{{{k_{\mathrm{ surf}}}{T_{surf}}}}{D}{\rm{ln}}( {\frac{{{T_{\mathrm{ bot}}}}}{{{T_{\mathrm{ surf}}}}}} )$| - Spherical |$\frac{{{k_{\mathrm{ surf}}}{T_{\mathrm{ surf}}}}}{D}f{( {\frac{R}{r}} )^2}{\rm{ln}}( {\frac{{{T_{\mathrm{ bot}}}}}{{{T_{\mathrm{ surf}}}}}} )$| Here, ΔT = (Tbot − Tsurf) is the bottom-to-top temperature jump, where Tsurf and Tbot are the surface and bottom temperature, and D is the thickness of the shell. In Cartesian geometry, z is depth, and in spherical geometry, r is radius, R the total radius, and |$f\ = {r_{\mathrm{ bot}}}/R\ $| the ratio between the inner and outer radii of the shell. k is the thermal conductivity, and its surface value is ksurf. Derivation of relationships for temperature-dependent conductivity are detailed in the Appendix. Open in new tab Table A1. Relationships for radial profiles of temperature and heat flux for a conductive system. Quantity . Geometry . Expression . k = constant Temperature Cartesian |${T_{\mathrm{ surf}}} + \Delta T\frac{z}{D}$| - Spherical |${T_{\mathrm{ surf}}} - \Delta T\frac{R}{D}f( {1 - \frac{R}{r}} )$| Heat flux Cartesian |$k\frac{{\Delta T}}{D}$| - Spherical |$k\frac{{\Delta T}}{D}f{( {\frac{R}{r}} )^2}$| |$k\ ( T ) = {k_{\mathrm{ surf}}}{T_{\mathrm{ surf}}}/T$| Temperature Cartesian |${T_{\mathrm{ surf}}}\exp[ {\frac{z}{D}{\rm{ln}}( {\frac{{{T_{\mathrm{ bot}}}}}{{{T_{\mathrm{ surf}}}}}} )} ]$| - Spherical |${T_{\mathrm{ surf}}}\exp[ { - \frac{R}{D}f( {1 - \frac{R}{r}} ){\rm{ln}}( {\frac{{{T_{\mathrm{ bot}}}}}{{{T_{\mathrm{ surf}}}}}} )} ]$| Heat flux Cartesian |$\frac{{{k_{\mathrm{ surf}}}{T_{surf}}}}{D}{\rm{ln}}( {\frac{{{T_{\mathrm{ bot}}}}}{{{T_{\mathrm{ surf}}}}}} )$| - Spherical |$\frac{{{k_{\mathrm{ surf}}}{T_{\mathrm{ surf}}}}}{D}f{( {\frac{R}{r}} )^2}{\rm{ln}}( {\frac{{{T_{\mathrm{ bot}}}}}{{{T_{\mathrm{ surf}}}}}} )$| Quantity . Geometry . Expression . k = constant Temperature Cartesian |${T_{\mathrm{ surf}}} + \Delta T\frac{z}{D}$| - Spherical |${T_{\mathrm{ surf}}} - \Delta T\frac{R}{D}f( {1 - \frac{R}{r}} )$| Heat flux Cartesian |$k\frac{{\Delta T}}{D}$| - Spherical |$k\frac{{\Delta T}}{D}f{( {\frac{R}{r}} )^2}$| |$k\ ( T ) = {k_{\mathrm{ surf}}}{T_{\mathrm{ surf}}}/T$| Temperature Cartesian |${T_{\mathrm{ surf}}}\exp[ {\frac{z}{D}{\rm{ln}}( {\frac{{{T_{\mathrm{ bot}}}}}{{{T_{\mathrm{ surf}}}}}} )} ]$| - Spherical |${T_{\mathrm{ surf}}}\exp[ { - \frac{R}{D}f( {1 - \frac{R}{r}} ){\rm{ln}}( {\frac{{{T_{\mathrm{ bot}}}}}{{{T_{\mathrm{ surf}}}}}} )} ]$| Heat flux Cartesian |$\frac{{{k_{\mathrm{ surf}}}{T_{surf}}}}{D}{\rm{ln}}( {\frac{{{T_{\mathrm{ bot}}}}}{{{T_{\mathrm{ surf}}}}}} )$| - Spherical |$\frac{{{k_{\mathrm{ surf}}}{T_{\mathrm{ surf}}}}}{D}f{( {\frac{R}{r}} )^2}{\rm{ln}}( {\frac{{{T_{\mathrm{ bot}}}}}{{{T_{\mathrm{ surf}}}}}} )$| Here, ΔT = (Tbot − Tsurf) is the bottom-to-top temperature jump, where Tsurf and Tbot are the surface and bottom temperature, and D is the thickness of the shell. In Cartesian geometry, z is depth, and in spherical geometry, r is radius, R the total radius, and |$f\ = {r_{\mathrm{ bot}}}/R\ $| the ratio between the inner and outer radii of the shell. k is the thermal conductivity, and its surface value is ksurf. Derivation of relationships for temperature-dependent conductivity are detailed in the Appendix. Open in new tab © The Author(s) 2020. Published by Oxford University Press on behalf of The Royal Astronomical Society. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model)
TI - Stagnant lid convection with temperature-dependent thermal conductivity and the thermal evolution of icy worlds
JO - Geophysical Journal International
DO - 10.1093/gji/ggaa540
DA - 2020-12-09
UR - https://www.deepdyve.com/lp/oxford-university-press/stagnant-lid-convection-with-temperature-dependent-thermal-vM0ROmQaFd
SP - 1870
EP - 1889
VL - 224
IS - 3
DP - DeepDyve
ER -