# Thermal evolution of lava planets

Mahesh Herath, Charles-Édouard Boukaré, Nicolas B. Cowan

► To cite this version:

Mahesh Herath, Charles-Édouard Boukaré, Nicolas B. Cowan. Thermal evolution of lava planets. Monthly Notices of the Royal Astronomical Society, 2024, 535, pp.2404-2414. 10.1093/mn-ras/stae2431. insu-05081544

HAL Id: insu-05081544

https://insu.hal.science/insu-05081544v1

Submitted on 23 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Thermal evolution of lava planets

Mahesh Herath,<fnref n="1" />,<fnref n="2" />* Charles-Édouard Boukaré<fnref n="3" />,<fnref n="4" /> and Nicolas B. Cowan † <fnref n="1" />,<fnref n="2" />,<fnref n="5" />*

Accepted 2024 October 15. Received 2024 October 15; in original form 2024 May 14

## ABSTRACT

Rocky planets typically form with a transient magma ocean. Lava planets, however, maintain a permanent day-side magma ocean. The extent of this magma ocean depends on the planet's thermal history. We present numerical simulations of the thermal history of tidally locked lava planets over 10 billion years, starting from a completely molten mantle. The day-side surface temperature is fixed at 3000 K, while the night-side surface temperature cools by thermal radiation. We consider planets with radii of 1.0 and 1.5R<sub>⊕</sub>; super-Earths have shallower steady-state magma oceans due to their greater gravity. The night-side begins crystallizing within a few thousand years, fully solidifying in 800 Myr, in the absence of tidal heating or day-night heat transport. We find that a mushy night-side can persist if at least 20 per cent of absorbed stellar power is transferred from the day to night hemisphere through magma currents, which would be feasible at a viscosity of 10<sup>-3</sup> Pa s. Maintaining a fully molten night-side by magma ocean circulation would require unrealistically low viscosities and therefore appears unlikely. Alternatively, the night-side may remain molten if the mush layer dissipates tidal energy at a rate of 8 × 10<sup>-4</sup> W kg<sup>-1</sup>, which is plausible for orbital eccentricities greater than 7 × 10<sup>-3</sup>. Night-side cooling, however, is a runaway process: increasing viscosity and mush solidification hinder both heat transport and tidal heating. Our results highlight the importance of measuring the night-side temperatures of lava planets, which would provide crucial insights into their thermal histories.

**Key words:** planets and satellites – planets and satellites: interiors – planets and satellites: physical evolution – planets and satellites: terrestrial planets.

## 1 INTRODUCTION

Lava planets are terrestrial planets with ultra-short orbital periods. Examples include CoRoT-7b (Léger et al. 2009), Kepler-10b (Batalha et al. 2011), TOI-561b (Patel et al. 2023), and K2-141b (Barragán et al. 2018; Malavolta et al. 2018). A lava planet's day-side is hot enough to melt and vapourize rock (Schaefer & Fegley 2009), while its night-side remains dark and cold if the planet is tidally locked into synchronous rotation. A lava planet therefore has two hemispheres experiencing different thermal evolutionary paths in parallel.

Thermal phase curves of lava planets constrain their day-side and night-side temperatures. With such data, we may be able to infer the history of lava planets by linking surface properties to an evolving interior structure and internal thermodynamics. To predict the characteristics of the magma ocean, we must create a thermal history model to simulate how lava planets evolve from their early days to the present.

The magma oceans on lava planets are analogous to those that existed in the early Solar system (Solomatov 2000; Elkins-Tanton 2012). Where they differ is in their boundary conditions: the magma ocean on a lava planet is heated from above, and may be in a steady-state. Depending on where the geotherm of the mantle meets the solidus and liquidus for silicates, the mantle could also have a mush layer (melt fraction between 30 and 60 per cent). Models of magma oceans have been developed for planets in the Solar system (Elkins-Tanton 2012) and applied to lava planets (Leger et al. 2011; Kite et al. 2016; Nguyen et al. 2020; Chao et al. 2021) predicts shallow day-side magma oceans. Boukaré, Cowan & Badro (2022), however, predicted deep magma oceans including extended mush zones.

In this paper, we develop a thermal model to simulate how a lava planet would evolve from a completely molten state to its present state. We use a pair of coupled 1D models for the two sides of a tidally locked planet to track how their interior structure and thermal profiles evolve. With this model, we simulated the thermal evolution of lava planets of different sizes, core mass fractions (CMFs), tidal heating, and day-to-night heat transport. We describe our model in Section 2, present our results in Section 3, and discuss them in Section 4.

<footnote n="1">¹Department of Earth & Planetary Science, McGill University, 3450 Rue University, Montreal, QC H3A 0E8, Canada</footnote>

<footnote n="2">²Trottier Space Institute, McGill University, 3550 Rue University, Montreal, QC H3A 2A7, Canada</footnote>

<footnote n="3">³Institut de Physique du Globe de Paris, 1 Rue Jussieu, 75005 Paris, France</footnote>

<footnote n="4">⁴Department of Physics & Astronomy, York University, 140 Campus Walk Room 128, North York, Toronto, ON M3J 1P3, Canada</footnote>

<footnote n="5">⁵Department of Physics, McGill University, 3600 Rue University, Montreal, QC H3A 2T8, Canada</footnote>

<footnote>* E-mail: mherath2@gmail.com (MH); nicolas.cowan@mcgill.ca (NC)</footnote>
## 2 MODEL

### 2.1 Single-column model

We considered a tidally locked planet with an Earth-like composition of peridotite (olivine and pyroxene). We additionally assumed that any existing atmosphere is optically thin (Castan & Menou 2011; Nguyen et al. 2020, 2022). We approximated the planet's day-side surface temperature $T_s$ with the local equilibrium temperature at the substellar point. Surface temperatures above 1800 K are necessary to sustain a molten surface (Kite et al. 2016). We used the substellar temperature and radius of K2-141b to test out our models for super-Earths, since K2-141b is the highest signal-to-noise ratio lava planet yet found (Barragán et al. 2018; Malavolta et al. 2018). It is also the subject of two approved JWST observing campaigns (Dang et al. 2021; Espinoza et al. 2021). $T_s$ was set to 3000 K on the day-side, the substellar temperature expected for K2-141b (Malavolta et al. 2018). We neglect the physics of solid-liquid phase separation and focus only on the thermal aspect of magma ocean evolution. For the night-side, the radiative equilibrium temperature is formally zero, so we allowed $T_s$ to vary in accordance with the flux escaping the surface. In the initial time-step, we started off with a value of $T_s = 3000$ on the night-side. In subsequent time steps, $T_s$ was found using the Stefan-Boltzmann law on flux radiated from the surface.

As shown in Fig. 1, we modelled the thermal evolution of this planet with a 1D parametrized model coupling silicate layers of different rheological properties with the core. Initially, the mantle was assumed to be fully molten, meaning that at the start of the simulation, we have the core and a layer of magma. The cooling of the magma ocean (MO) would lead to the formation of magma mush (MS) followed by the formation of solid rock (S) in the mantle. The formation of mush and solid rock was determined by the intersection of the geotherm with the liquidus and solidus curves for bulk silicate material. While the radial position of the core-mantle boundary (CMB) remains unchanged, the number of layers in the mantle as well as their respective thickness change with time. We applied a viscosity ($\eta$) of $\eta = 0.1$ Pa s (a melt fraction greater than 50 per cent) for the magma (Xie et al. 2021; Le Losq & Baldoni 2023). We assumed the mush is composed of 50 per cent of solid. The viscosity of the partially molten mushy layer was set to $10^{14}$ Pa s (a melt fraction between 10 and 40 per cent) (Xie et al. 2021; Le Losq & Baldoni 2023). For solid rock, we used $\eta = 10^{18}$ Pa s (Tyler, Henning & Hamilton 2015). As shown in Fig. 1, the parameters $T_1$ to $T_6$ give the temperatures at the boundaries of each layer in the mantle. When we have a completely molten mantle, we only consider $T_1$ and $T_2$. $T_1$ is the temperature just below the surface while $T_2$ is the temperature at the bottom of the magma ocean. Once it cools below the liquidus, we consider $T_3$ and $T_4$ in addition to $T_1$ and $T_2$. $T_3$ and $T_4$ are the temperatures at the top and bottom of the mush layer. When the temperature goes below the solidus, we include $T_5$ and $T_6$, the temperatures at the top and bottom of the solid layer, respectively. The average of $T_2$ and $T_3$ is $T_i$ while the average of $T_4$ and $T_5$ is $T_{i2}$. $T_c$ gives the temperature at the CMB. The differences between $T_3$ and $T_2$, $T_4$ and $T_5$, and $T_6$ and $T_c$, are due to the presence of the thermal boundary layers.

The parameters $r_1$ to $r_5$ give the radial dimensions of the mantle layers at a given time. The terms $r_2$ and $r_3$ denote the boundary lines just above and just below the interface at which the magma transitions to mush, and $r_4$ and $r_5$ are the boundary lines above and below the interface where the mush transitions to solid rock. The transition from magma to mush takes place at $r_i$ and is the average of $r_2$ and $r_3$. Similarly, $r_{i2}$ is the transition point from mush to rock and

is the average of $r_4$ and $r_5$. To find these points, we set the thermal profile to the liquidus and solidus. For the liquidus $T_{liq}$ and solidus $T_{sol}$, we used expressions derived from the high pressure experiments of Fiquet et al. (2010) and Zhang & Herzberg (1994).

$$T_{\text{liq}} = 2000 \left( 0.1169 \left( \frac{P}{\text{GPa}} \right) + 1 \right)^{0.32726} \quad (1)$$

$$T_{\text{sol}} = 1674 \left( 0.0971 \left( \frac{P}{\text{GPa}} \right) + 1 \right)^{0.35175} \quad (2)$$

We assumed the interior is fully convective, and therefore the layers have adiabatic temperature profiles. We start simulations from a fully molten mantle, so $T_c$ can be extrapolated from $T_s$ by following an adiabatic geotherm. The model begins its run at the point where the geotherm of the fully molten mantle begins to dip below the liquidus curve. Following the approach of Labrosse (2003, 2015), the relationship between $T_1$ and $T_2$ is

$$T_2 = T_1 \exp \left( \frac{R_p - R_c}{D_{\text{MO}}} \right), \quad (3)$$

where $D_{\text{MO}}$ is the temperature scale height of the magma ocean. Once mush and solid layer begin to form, we determine $T_4$ and $T_6$ with

$$T_4 = T_3 \exp \left( \frac{r_3 - r_4}{D_{\text{MS}}} \right) \quad (4)$$

$$T_6 = T_5 \exp \left( \frac{r_5 - r_6}{D_S} \right), \quad (5)$$

where $D_{\text{MS}}$ and $D_S$ are the temperature scale heights of the mush and rock, respectively.

Energy conservation in each layer gives the expressions describing the evolution of the temperature in each layer. These temperatures depend on the difference in heat flux entering and leaving the layer. For the magma ocean, we evolve the temperatures $T_s$ and $T_1$. We used the Stefan-Boltzmann law to calculate the flux (F) radiated from the surface using the difference in incident stellar flux ($T_{\text{eq}}$) and the surface temperature.

$$F = \sigma (T_s^4 - T_{\text{eq}}^4). \quad (6)$$

Then to evolve $T_s$, we used the equation

$$\frac{dT_s}{dt} = \frac{-F + Q_{\text{MO}}A_{\text{MO}}}{V_{\text{MO}}C\rho}. \quad (7)$$

The terms $Q$ and $A$ denote the heat flux and surface area for a given layer. $Q_{\text{MO}}A_{\text{MO}}$ is the heat lost by the magma ocean. $A_{\text{MO}} = 4\pi R_p^2$ and $A_{\text{MS}} = 4\pi r_i^2$. $V$, $C$, and $\rho$ denote the volume, specific heat capacity, and density of a given layer of the mantle. For simplicity, we use constant values for $C$ and $\rho$, which are given in Table A1. For the magma ocean, $V_{\text{MO}}$ is calculated as $\frac{4}{3}\pi(R_p^3 - r_i^3)$. For the evolution of $T_1$, we used the equation

$$\frac{dT_1}{dt} = \frac{-Q_{\text{MO}}A_{\text{MO}} + Q_{\text{MS}}A_{\text{MS}} + \frac{1}{2}L\rho A_{\text{MS}} \frac{d r_{i1}}{d t}}{V_{\text{MO}}C\rho}, \quad (8)$$

where $Q_{\text{MS}}A_{\text{MS}}$ is the heat flux entering the ocean from the mush below. $L$ is the latent heat released into the magma ocean due to the solidification from liquid to mushy state, i.e. 50 per cent of melt. $\frac{d r_{i1}}{d t}$ gives the rate at which the mush front grows. We use similar energy conservation equations for the magma mush, solid rock, and the CMB, respectively, denoted as

$$\frac{d T_3}{d t} = \frac{-Q_{\text{MS}}A_{\text{MS}} + Q_{\text{S}}A_{\text{S}} + \frac{1}{2} L \rho A_{\text{S}} \frac{d r_{i2}}{d t}}{V_{\text{MS}} C \rho} \quad (9)$$
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Schematic of the 1D thermal evolution model used in this work. $Q_R$ is the radiated flux, $Q_{MO}$ is the heat flux from the magma ocean, $Q_{MS}$ is the heat flux from the magma mush, $Q_S$ is the flux from solid rock, and $Q_c$ is the flux from the core. $Q_H$ is the heat removed from the day-side and added to the night-side when we include horizontal heat transport. $Q_T$ is the tidal heating when tidal dissipation is included in the model. The light blue line is the geotherm of the mantle while the line ‘L’ (dark blue line) is the mantle liquidus and ‘S’ (green line) is the mantle solidus. The pressure increases in the downward direction while the temperature increases to the right.</figcaption></figure>

$$ \frac{dT_5}{dt} = \frac{-Q_S A_S + Q_c A_c}{V_S C \rho} \qquad (10) $$

$$ \frac{dT_c}{dt} = \frac{-Q_c A_c}{V_S C \rho}. \qquad (11) $$

The terms $Q_{MS}$, $Q_S$, and $Q_c$ are the heat flux leaving the mush, solid rock, and CMB, respectively. The volumes of the mush and solid layers are calculated as $V_{MS} = \frac{4}{3}\pi(r_i^3 - r_{i2}^3)$ and $V_S = \frac{4}{3}\pi(r_{i2}^3 - R_c^3)$.

We next assume that the heat flux leaving each layer is proportional to the temperature gradient across that layer. The expression used for the magma ocean is

$$ Q_{MO} = 0.22\kappa \left( \frac{T_1 - T_S}{R_p - r_i} \right) (Ra_{MO})^{\frac{2}{7}} \text{Pr}^{-\frac{1}{7}}, \qquad (12) $$

where $\kappa$ is the thermal conductivity of the magma ocean, $Ra_{MO}$ the Rayleigh number, and Pr its Prandtl number. We look at the convective efficiency between $R_p$ and $r_i$ to compute the Rayleigh number of the magma ocean. We used the temperature difference between the surface and thermal interface ($T_i$ and $T_S$) corrected by the adiabatic gradient. We get the sum of this difference when calculating the Rayleigh number:

$$ Ra_{MO} = \frac{\alpha \rho g (T_i - T_S - (T_2 - T_1))(R_p - r_i)^3}{\eta_{MO} \kappa_d}, \qquad (13) $$

where $\alpha$ is the thermal expansivity and $\kappa_d$ is the thermal diffusivity. We found Pr with the equation

$$ Pr = \frac{\eta_k}{\kappa_d}, \qquad (14) $$

where $\eta_k$ is the kinematic viscosity. The heat flux through the mush layer, $Q_{MS}$, is controlled by the thermal properties of the mush. We used the temperature difference between the boundary layer at the upper portion of the mush and the ocean-mush interface to find the heat flux with

$$ Q_{MS} = 0.22\kappa \left( \frac{T_3 - T_i}{r_i - R_c} \right) (Ra_{MS})^{\frac{1}{3}} \qquad (15) $$

$$ Q_{MS} = 0.22\kappa \left( \frac{T_3 - T_i}{r_i - r_{i2}} \right) (Ra_{MS})^{\frac{1}{3}}. \qquad (16) $$

If we have the CMB right below the magma mush, we use equation (15) or if it is a solidification front right below, we use equation (16). Since we are looking at layers of mush and solid rock, in which inertial forces are negligible. The Rayleigh number for the mush is defined as

$$ Ra_{MS} = \frac{\alpha \rho g (T_c - T_i + T_3 - T_4)(r_i - R_c)^3}{\eta_{MS} \kappa_d} \qquad (17) $$

$$ Ra_{MS} = \frac{\alpha \rho g (T_{i2} - T_i + T_3 - T_4)(r_i - r_{i2})^3}{\eta_{MS} \kappa_d}, \qquad (18) $$

where we use equation (17) if there is no solid layer and equation (18) if the solid layer is below the mush. The equations of flux and Rayleigh number for the layer of solid rock are

$$ Q_S = 0.22\kappa \left( \frac{T_5 - T_{i2}}{r_{i2} - R_c} \right) (Ra_S)^{\frac{1}{3}} \qquad (19) $$

$$ Ra_S = \frac{\alpha \rho g (T_c - T_{i2} + T_5 - T_6)(r_{i2} - R_c)^3}{\eta_S \kappa_d}. \qquad (20) $$
Finally, for the CMB, we use the same Rayleigh number as what is given in equation (20), while the heat flux can be evaluated with the expression

$$Q_c = 0.22\kappa \frac{(T_c - T_i)}{r_i - R_c} (Ra_s)^{\frac{1}{3}} \quad (21)$$

$$Q_c = 0.22\kappa \frac{(T_c - T_6)}{r_{i2} - R_c} (Ra_s)^{\frac{1}{3}}, \quad (22)$$

where equations (21) and (22) are if we have mush or solid rock above the CMB, respectively.

## 2.2 Horizontal heat transfer

While the vertical convective currents are driven by the temperature differences between the surface and the CMB, we also get horizontal convective currents driven by the day-to-night temperature difference. We consider day-to-night heat transport through the deep magma ocean expected to persist when the planetary interior is hot (Boukare et al. 2022, 2023) rather than through the shallow magma ocean expected in a cold start scenario (Kite et al. 2016; Lai et al. 2024a, 2024b). In our model, the day-side surface temperature is fixed at $T_s = 3000$ K, while the night-side has a time varying $T_s$. We assume that most of the horizontal heat transfer happens in the layers with the lowest viscosities. We add a sink/source term to the energy conservation equations of the magma oceans on the day and night sides. We use a thermal Rayleigh number $Ra_H$ for horizontal convection from Hughes & Griffiths (2008). $Ra_H$ differs in that we use the difference in $T_s$ between the day-side and night-side for the temperature contrast, and we use half the planetary circumference as the length. We can then have the expression

$$Ra_H = \frac{\alpha \rho g (T_{sd} - T_{sn}) (\pi R_p)^3}{\eta_{ns} \kappa_d}, \quad (23)$$

where $T_{sd}$ and $T_{sn}$ are the day-side and night-side surface temperatures, respectively. The parameter $\eta_{ns}$ is the magma viscosity, while $\rho_{ns}$ is the density of that same layer. We then evaluate the horizontal heat flux $Q_H$ with the expression

$$Q_H = 0.22\kappa \frac{(T_{sd} - T_{sn})}{\pi R_p} (Ra_H)^\alpha. \quad (24)$$

We set the value of $\alpha$ to $\frac{2}{7}$ under the assumption we have turbulent convection in the magma (Hughes & Griffiths 2008). By adjusting $\eta_{ns}$, we tested different quantities of heat flux flowing to the night-side. When $\eta_{ns}$ increases, the magnitude of horizontal heat flow will get weaker. To account for horizontal heat flow from the day-side, we added $Q_H$ to equation (7) as applied to the day-side magma ocean. For the night-side, we reduced $Q_H$ from equation (7) as applied to the night-side magma ocean. This gave new expressions for the day-side and night-side evolution of the surface temperature.

$$\frac{dT_{i,d}}{dt} = \frac{-Q_{MO}A_{MO} + Q_{MS}A_{MS} + \frac{1}{2}L\rho A_{MS} \frac{dr_{i,d}}{dt} - Q_H}{V_{MO}C\rho} \quad (25)$$

$$\frac{dT_{i,n}}{dt} = \frac{-Q_{MO}A_{MO} + Q_{MS}A_{MS} + \frac{1}{2}L\rho A_{MS} \frac{dr_{i,n}}{dt} + Q_H}{V_{MO}C\rho}. \quad (26)$$

## 2.3 Tidal heating

Because of the proximity of lava planets to their host stars, it is likely that they experience strong tidal forces assuming they are not tidally locked (Leger et al. 2011). To include the effect of tidal heating $Q_T$, we assumed that most of the tidal dissipation happens

in the mush rather than in the magma or solid rock (Kervazo et al. 2021). Therefore, we scaled the magnitude of tidal dissipation with the mass of the mush at a given time. To find the amount of tidal energy dissipation per unit mass, we used the data from Tyler et al. (2015). Their work simulates the tidal heating taking place within the magma ocean of Io. They concluded that the solid part of Io could have tidal dissipation in the amount of 2.25 W m⁻². Converting the energy units from W m⁻² to total energy dissipation, we get ≈10¹⁴ W within Io. Dividing the total energy dissipation by the mass of Io and then scaling that value by the mass of mush in our model planets, we find $Q_T$.

$$Q_T = f \times \frac{10^{14}}{M_{Io}} \times M_{\text{mush}}. \quad (27)$$

The value of $Q_T$ can be increased or decreased by the scaling factor $f$. We add $Q_T$ to equation (7) along with $Q_H$.

$$\frac{dT_1}{dt} = \frac{-Q_{MO}A_{MO} + Q_{MS}A_{MS} + L\rho A_{MS} \frac{dr_1}{dt} + Q_H + Q_T}{V_{MO}C\rho}. \quad (28)$$

To investigate if the quantities of tidal heating needed to keep a molten night-side are possible, we calculated the relationship between orbital eccentricity ($e$) and tidal dissipation ($Q_{tidal}$). We used the tidal power equation from Driscoll & Barnes (2015).

$$Q_{\text{tidal}} = -\frac{21}{2} \text{Im}(k2) G^{\frac{3}{2}} M_s^{\frac{5}{2}} R_p^{\frac{5}{2}} \frac{e^2}{a^{15}}. \quad (29)$$

The parameters $G$, $M_s$, $a$, and $R_p$ are the gravitational constant, stellar mass, semimajor axis, and planetary radius, respectively. $\text{Im}(k2)$ is the imaginary component of the Love number $k2$, and was given a value of $3 \times 10^{-3}$ for an Earth-like composition. We used this value of $\text{Im}(k2)$ despite it being for a more viscous medium than mush, to investigate what quantities of $Q_T$ are possible with it.

The parameters $Q_H$ and $Q_T$ can be turned on or off to simulate different combinations of heating added to the night-side. The parameter $T_{eq}$ in equation (8) was fixed at 3000 K for the day-side, but given a value of zero for the night-side since there is no incoming stellar flux. The parameter $T_c$ was given the same value for both the day and night sides. The two hemispheres of the planet were evolved while exchanging heat from the day-side to the night-side only in the top-most layer. The simulations were run until no more changes were seen in the internal structure for each hemisphere of a given planet.

## 3 RESULTS

We applied our model to lava planets with radii of 1.0 and $1.5R_⊕$, the former for bench-marking with terrestrial models and the latter to simulate K2-141b. We evaluated each radius for a CMF of 32 and 70 per cent. The CMF of 32 per cent mimics an Earth-like interior, while the CMF of 70 per cent is for a Mercury-like interior. Our models showed that the latent heat released by the solidification of the magma ocean was smaller in magnitude by a factor of at least 100 at their highest value, compared to the radiated heat flux and the convective heat flux. This was because the rate of growth of mush and solid was too slow to have an effect on the latent heat quantity. The latent heat shows a steady decline in magnitude with a peak around $2 \times 10^{18}$ W for a $1.0R_⊕$ planet. This is when the growth rate is at its maximum value of $10^{-5}$ m s⁻¹. The peak energy was even lower for $1.5R_⊕$ on account of their mush growth rates being lower. It never reached a level within two orders of magnitude of the energy lost by radiation or the energy contribution via day-night heating and tidal
dissipation. Therefore, we will not take the effect of Latent heat into account in the rest of this study.

## 3.1 Earth-size planets

Fig. 2 shows the thermal and internal evolution of the day-side and night-side of a 1.0R⊕ planet with a CMF of 32 per cent. On the day-side, the surface reaches radiative equilibrium with the incoming stellar flux at a temperature of 3000 K. This is high enough to keep a stable layer of magma at the top of the mantle. The high level of flux leaving the magma ocean at the onset of cooling on the day-side reduces rapidly before reaching an equilibrium state with the flux from the magma mush. This only lasts for a brief time period before the flux from the mush reduces and itself reaches equilibrium with the flux from the newly formed solid rock at the bottom. On the night-side, the lower surface temperature and the absence of incoming flux leads to a higher magnitude of cooling. This results in the total solidification of the night-side after about 500 Myr. Fig. 2 shows the evolution of the two sides of the planet without horizontal heat transfer from the day-side to the night-side. The rapid cooling rate on the night-side would mean there is a large temperature contrast between the magma oceans of the two sides at a given time. Therefore, heat transfer from the hotter day-side to the night-side is highly probable.

We reran the simulations with the inclusion of horizontal heat transport $Q_H$ in the model. From equation (23), it can be shown that the viscosity of the magma controls $Ra_H$. The value of $\eta$ in equation (23) was the same for the day and night sides. The magma viscosity was varied from $1 \times 10^{-5}$ Pa s, and we recorded the values of $T_s$ resulting from each. A viscosity of $10^{-5}$ Pa s is an extreme scenario that we tested in order to explore the conditions under which a molten night-side surface could be maintained. In Fig. 3, we see how the state of the surface changes with $Q_H$. If we have a minimum $Q_H$ of $2.3 \times 10^{20}$ W, we observe a molten night-side surface. The energy transferred between hemispheres change with $\eta$. If the viscosity is between $1$ and $10^{-2}$ Pa s ($5 \times 10^{19} < Q_H < 8 \times 10^{19}$ W), the night-side completely solidifies after 2 Gyr. The surface temperature fluctuates between 1000 and 1500 K at these values. Reducing the viscosity to values between $10^{-2}$ and $10^{-4}$ Pa s ($8 \times 10^{19} < Q_H < 2.3 \times 10^{20}$ W) results in a mush surface as well as surface temperatures between 1500 and 2000 K. At $10^{-4}$ Pa s, we start to see a very thin (about 40 km) sustained magma ocean on the night-side. Theoretically there is enough energy at this particular viscosity to keep the night-side surface melted, which is consistent with the model result. However, the viscosity is expected to increase when temperature decreases and crystal fraction increases. Hence, the sustenance of a magma ocean could be impeded if we include a temperature-dependant viscosity. This means the existence of a magma ocean at $\eta = 10^{-4}$ Pa s may or may not be possible. Once the viscosity goes below $10^{-4}$ Pa s, $Q_H$ is high enough to maintain a night-side magma ocean. At $10^{-5}$ Pa s ($3.5 \times 10^{20}$ W), we could see a magma ocean about 120 km deep, meaning that a very low magma viscosity is necessary to have a melted night-side surface if we only take horizontal heat transfer into account. Given that a magma viscosity of $10^{-5}$ Pa s is an extreme scenario, it is unlikely that horizontal day-night convection alone could sustain a fully molten night-side magma ocean.

We increased the CMF of the 1.0R⊕ planet to 70 per cent for the next set of simulations. As before, the initial simulations did not have any horizontal heat transport before being added later on. The day-side and night-side did not show much variation compared to the models with a CMF of 32 per cent. The only noticeable difference was that there was less solid rock in the deep interior of the day-side

and the time-scale for total solidification on the night-side happens about 600 Myr earlier than it did for the 32 per cent CMF model. With the addition of horizontal heat transport, once again there was no discernible difference between the CMF = 32 and 70 per cent models. The value of $Q_H$ required to sustain a night-side magma ocean was the same as it was in the earlier models.

## 3.2 Super-Earths

For the next set of models, the radius of the planet was increased to 1.5R⊕ to approximate a Super-Earth. In Fig. 2, the right-side panels of 2(b), (d), (f), and (h) show the evolution of the internal structure of 1.5R⊕ model for the day-side and night-side. Compared to Earth-sized planets, there is no observable difference in the evolution of the night-side between them and 1.5R⊕ super-Earths. The day-side, however, shows a greater quantity of rock formation in the interior of the 1.5R⊕ model. The day-side magma ocean of the 1.5R⊕ model is also much shallower (≈190 km) than that of the 1.0R⊕ model (≈560 km). This is because of the higher mass of the 1.5R⊕ planet. With the increase in radius, we also had to increase the mass of the planet. For a 1.5R⊕ planet with a CMF of 32 per cent, the mass was found to be ≈7.0 M⊕, which contributes to a higher gravitational potential than in the 1.0R⊕ planet. Therefore, the pressure acting radially along the mantle would be higher. The higher pressure would necessitate a higher temperature to keep the mantle as a melt at a given point. In our model, the temperature profile falls below the solidus curve for a large portion of the mantle because of the aforementioned reasons, leading to more rock formation in the day-side interior.

As was done with the 1.0R⊕ model, we investigated the variation in $T_s$ of the night-side with respect to horizontal heat transport, where we changed the value of $Q_H$ by varying the magma viscosity. The model showed that the night-side surface solidifies in about 1.2 Gyr in the absence of $Q_H$ and $Q_T$. Once we activated $Q_H$, the values of night-side $T_s$ at the same viscosities used for the 1.0R⊕ model were slightly higher in the 1.5R⊕ model. The reason for this feature can be found in equation (24). The magnitude of $T_s$ on the night-side depends on the gravity and planetary radius, and scales approximately as $g^2/7 \times R_p^{-15}$. When we apply the difference in $g$ and $R_p$ between an Earth-sized and Super-Earth type planet into this scaling law, we find that the model derived $T_s$ for a Super-Earth is approximately 1.05 times the value of $T_s$ for an Earth-sized model.

The energy required to keep the night-side melted on a 1.5R⊕ planet is around $5.4 \times 10^{20}$ Wce again this makes sense as the energy radiated from the surface scales with the square of $R_p$. The panels in Fig. 5 demonstrate how the night-side internal structure varies for different magma viscosities. More time is needed for cooling, and as a result, we see the night-side takes approximately 4.8 billion years to fully solidify at $\eta = 1$ Pa s compared to the time taken for the 1.0R⊕ model at the same viscosity. The pattern follows for subsequent values of $\eta$ until we get to $\eta = 10^{-2}$ Pa s. The latter shows that it takes around 10 billion years to fully solidify. For Super-Earth's much younger than this time-scale, it is plausible that the night-side would be melted even at $10^{-2}$ Pa s. Between viscosities of $10^{-2}$ and $10^{-4}$ Pa s, we can see night-side surfaces of mush, until we get to $\eta = 10^{-5}$ Pa s, where we start to see a long-term night-side magma ocean. It is, however, a lot thinner (about 80 km) compared to the magma oceans possible on the surface of a 1.0R⊕ planet. This is once again due to the larger gravitational potential of a Super-Earth. It is possible that if $Q_H$ is the only source of
<figure><img src="image_3.png" /><figcaption>**Figure 2.** The thermal and internal evolution of a 1.0R⊕ planet (left) and a 1.5R⊕ planet (right) with a core mass fraction of 32 per cent in the absence of day-to-night heat transport or tidal heating. The top four panels show the evolution of the day-side while the bottom four panels show the night-side evolution. The coloured lines show the temperatures at the boundaries shown in Fig. 1. In panel 2(a), the temperatures $T_4$ and $T_6$ are overlapping $T_3$ and $T_5$, respectively. The temperature increases downwards in these plots to ease comparison with the structure model. The coloured shading shows how the depth of the magma, mush, and solid rock change in time. The sections where the colours are translucent are the regions where phase transitions are happening (i.e. mush and solid crystals). The day-side of the 1.0R⊕ planet maintains a 550 km deep magma ocean in a steady state, and the day-side of the 1.5R⊕ has a shallower magma ocean of about 200 km. The 1.5R⊕ planet has a shallow magma ocean and more solid rock in the day-side interior because of the higher gravity of the bigger planet. The night-side of both planets solidifies completely within approximately $10^9$ yr; the smaller planet does so faster due to its smaller thermal inertia.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>**Figure 3.** The interior evolution of the night-side of 1.0R⊕ (left) and 1.5R⊕ (right) planets with a core mass fraction of 32 per cent and the inclusion of horizontal heat transport $Q_H$. Panels (a) and (b) show the evolution at a magma viscosity of $\eta = 1$ Pa s, (c) and (d) for $\eta = 10^{-3}$ Pa s, and (e) and (f) for $\eta = 10^{-5}$ Pa s. In (a) and (b), we can see the night-side fully solidifying, while (c) and (d) show a scenario where the night-side surface remains mushy. Panels (e) and (f) show that the night-side surface can be kept molten for both Earth-sized and Super-Earth lava planets if the magma viscosity is sufficiently low.</figcaption></figure>

heat powering the night-side, we might see surfaces of mush but no magma oceans on the night-side of Super-Earth's because of the aforementioned reason. This is also assuming that the viscosity cannot be continuously reduced to smaller values.

Just as we observed in the 1.0R⊕ model, increasing the CMF to 70 per cent did not show a difference from the results with a CMF of 32 per cent. Hence, we conclude based on our current findings that the CMF does not affect the thermal evolution, with the parameters we are currently using in our model.

## 3.3 The inclusion of tidal heating

We added tidal heating $Q_T$ to both the 1.0 and 1.5R⊕ models (see Fig. 4). With the constraints we applied on the magnitude of tidal heating from Section 2.3, we initially included only the effect of $Q_T$ with $Q_H$ set to zero, and then ran a set of simulations with both $Q_T$ and $Q_H$ given non-zero values. Figs. 5(a) and (b) show the 2D parameter space of steady-state $Q_H$, $Q_T$, and the resulting night-side

surface conditions along with the results of numerical simulations conducted at different horizontal and tidal heating values.

For both 1.0 and 1.5R⊕ planets, at least $8 \times 10^{-4}$ W kg$^{-1}$ of tidal heating energy is required to keep the night-side melted if the model relies only on tidal heating without $Q_H$. This is approximately $8 \times 10^5$ times the amount of tidal heating observed in Io. Otherwise, tidal heating alone is not enough to keep the night-side molten. This is because $Q_T$ is tied to the mush mass, which increases initially as the magma ocean cools, and then decreases when the mush starts to solidify into rock. As this is not a fixed energy source, the reduction of mush leads to the reduction in tidal heating. The reduction of tidal heating in turn leads to the reduction of the surface temperature, which cycles around to the further reduction of mush. This phenomenon of runaway night-side cooling is why when we only have tidal heating as a heat source for the night-side, we would see high night-side surface temperatures between 100 million and 2 billion years, before it starts to rapidly reduce leading to the solidification of the night-side.
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Evolution of the interior with the inclusion of both day–night heat transport ($Q_H$) and tidal heating ($Q_T$). The wattage per kg of mush is $q_T = 7 \times 10^{-4}$ W kg$^{-1}$. The mush mass varies by a factor of ~2, peaking at $\le 10^8$ yr before dropping to approximately zero after a few billion years. Tidal heating is greater for larger planets because of the greater mass of mush.</figcaption></figure>

For $1.0R_⊕$ planets, even with the inclusion of $Q_H$, the amount of $q_T$ required for either melting or creating a sea of mush on the night-side remains at or above $7 \times 10^{-4}$ W kg$^{-1}$ for magma viscosities up to $10^{-3}$ Pa s. At $\eta = 10^{-3}$ Pa s and at lower magma viscosities, the combination of tidal heating and day-night heat transfer leads to a greater parameter space where a night-side that is either molten or mush is possible. Between viscosities of $10^{-4}$ and $10^{-5}$ Pa s, we see a molten night-side even at values of $q_T$ as low as $8 \times 10^{-5}$ W kg$^{-1}$. For $1.5R_⊕$ planets, the minimum wattage per unit mass at which melting occurs on the night-side surface is lower than what was observed for $1.0R_⊕$ planets. Consequently, there is a larger parameter space (see Fig. 5), in which a molten or mush night-side surface is possible. The reason for this is that a Super-Earth will have a higher mass of mush at a given point in its cooling history than an Earth-sized planet. Since $Q_T$ is proportional to the mass of mush, a higher mush mass means more tidal heating for a given amount of $q_T$. Because of the higher mass of Super-Earth's, the cooling rate is slower, leading to longer lasting quantities of mush on the night-side. Because of these reasons, the night-side of Super-Earth's can have a mushy night-side surface with $3 \times 10^{-4}$ W kg$^{-1}$ at a magma viscosity of $10^{-1}$ Pa s, which is

not possible for an Earth-sized planet. From Fig. 5(d), it can be seen that a sea of mush is possible for the night-side of a Super-Earth with energies below $1 \times 10^{-4}$ W kg$^{-1}$ at $\eta = 10^{-3}$ Pa s and a molten surface at $q_T < 5 \times 10^{-5}$ W kg$^{-1}$ at $\eta > 10^{-3}$ Pa s. A Super Earth with a magma viscosity of $\eta = 10^{-3}$ Pa s, which is comparable to the viscosity of heated water, could have a magma ocean or an ocean of mush if both tidal heating and horizontal heating is happening together.

Applying the equations from Driscoll & Barnes (2015) to K2-141b, we found what sort of eccentricity would enable the planet to generate the amount of tidal heating needed to have a molten night-side. If the eccentricity is between 0.007 and 0.01, it is possible to have the $Q_T$ needed to melt the night-side (see Fig. 6).

## 4 DISCUSSION AND CONCLUSIONS

Our thermal evolution model shows under what circumstances the night-side of a lava planets can be kept in a molten or mush state. We also found the time-scales at which the night-side can cool down from magma to solid rock and the physical parameters that determine
<figure><img src="image_6.png" /></figure>

**Figure 5.** Night-side surface end-state as a function of day-night heat transport ($Q_H$), and night-side tidal heating ($Q_T$). Panels 5(a) (1.0R⊕ planets) and (b) (1.5R⊕ planets) show the theoretical 2D parameter space based on the Stefan–Boltzmann formula for night-side cooling. At the top of each panel, we show the horizontal heat transport $Q_H$ as a percentage of the instellation. Panels 5(c) and (d) show a grid of 42 numerical simulations using our thermal evolution model. These figures show the viscosity $\eta$ and specific tidal heating, $q_T$, used in each simulation. The colour of each datum gives the state of the night-side surface after running the simulation for 10 Gyr. The upper and lower panels show qualitatively the same pattern.

<figure><img src="image_7.png" /></figure>

**Figure 6.** A comparison of the tidal heating energy generated in Earth (blue), Io (orange), and K2-141b (red) as a function of orbital eccentricity. For K2-141b, the solid line is for the planet's actual 1.5R⊕ radius, while the dashed line shows the tidal heating if the planet was Earth-sized. The calculations presume an Earth-like Im(k2) = 3 × 10⁻³ for an Earth-like composition. The black dots show the minimal tidal energy required for the night-side to remain melted in the absence of day-to-night transport (see top panels of Fig. 5). The black triangle and square show the average tidal energy dissipated inside Io and Earth, respectively.

these time-scales. In addition to considering the effect of planetary radius, we also increased the CMF in some of our models from 32 to 70 per cent to simulate a Super-Mercury. The day-side and night-side evolution showed no change with an increase in CMF for either 1.0 or 1.5R⊕ planets. The simulations showed that latent heat does not significantly affect the thermal evolution because the growth rates

of the mush and solid layers are too slow to boost the magnitude of the latent heat too slow. The latent heat is at most 100 times less than either $Q_R$, $Q_H$, or $q_T$. If the growth rate is fast enough, it could potentially have a slight effect on early evolution but not at billion year time-scales.

The solidus and liquidus curves used in this study (Zhang & Herzberg 1994; Fiquet et al. 2010) were extrapolated to pressures beyond 150 GPa. We conducted a sensitivity analysis of the extrapolations of these melt curves by generating synthetic data points above 150 GPa and refitting for the melt curves including the new data. We adjusted the solidus–liquidus equations with the new fitted parameters and reran the models. We found that the refitted melt curves do not alter time-scales for the thermal evolution of super-Earths. This study does, however, motivate the acquisition of new experimental data for silicate melts at extreme pressures.

For 1.0R⊕ planets, a minimum of $Q_H = 2.3 \times 10^{20}$ W of heat must be transferred from the day-side to the night-side to keep it molten. This is approximately 20 per cent of the incident stellar flux being transported across hemispheres. The viscosity of magma needed for this quantity of heat to flow in our model is $10^{-4}$ Pa s. Between viscosities of $10^{-1}$ and $10^{-3}$ Pa s, our model showed a night-side that was kept in a mush state for billions of years. At and below $10^{-1}$ Pa s, the energy transmitted horizontally is insufficient to keep the night-side from solidifying. To put these numbers in context, magma at thousands of K has a viscosity of $10^{-3}$ Pa s, similar to that of water (Malosso et al. 2022). In the case of a 1.5R⊕ Super-Earth, the steady-state energy requirement for a night-side magma ocean is $Q_H = 5.2 \times 10^{20}$ W, or 23 per cent of the stellar instellation flux. The viscosity of magma required to transport this quantity of energy
is about the same as for an Earth-sized planet. It should be noted,
however, that the magma oceans of the latter are far shallower than
that of the former because of the higher gravity of Super-Earth's.
Since the magma viscosity had to be lowered to extremely low
values (less than $10^{-4}$ Pa s) to sustain a night-side magma ocean
in our models, it is likely that horizontal day-night convection alone
is insufficient to maintain a molten night-side. Furthermore, the work
of McKenzie et al. (2000) suggests that a gradient in the solid
fraction between the day-side and night-side mantles could cause
cold material to flow from the night-side to the day-side, impacting
the transfer of heat across hemispheres.

For a tidally heated lava planet, the energy requirement for night-
side melting is the same as for day-night heat transfer (see Fig. 5(a)
and (b)). A $1.0R_⊕$ planet would require $Q_T = 1.2 \times 10^{-4} \text{ W kg}^{-1}$
of tidal heating while a $1.5R_⊕$ planet would need $Q_T = 4.8 \times$
$10^{-5} \text{ W kg}^{-1}$ of energy. But because tidal heating depends on the
mush mass, the simulated quantities deviate from the idealized
steady-state values. The tidal energy required for a molten night-
side is much higher because the mass of mush decreases as the
interior cools down. Specific tidal heating of $q_T = 8 \times 10^{-4} \text{ W kg}^{-1}$
is necessary for a sustained night-side magma ocean that would last
10 billion years. Values of $q_T$ down to $q_T = 5 \times 10^{-4} \text{ W kg}^{-1}$ could
keep the night-side molten for several billion years before it cools
to either mush or solid. Our simulations show that lava planets can
experience a brief period of internal warming when mush mass is at
its maximum. The night-side magma ocean is sustained for periods
between 800 million and 3 billion years in such instances before they
experience runaway cooling as the mush becomes more viscous and
solidifies. In super-Earths, it is more likely that the night-side would
either stay molten or mushy with a specific tidal heating $q_T$. But
their magma oceans and mush oceans will be shallower due to the
higher gravity of the Super-Earth. It should be noted, however, that a
sustained partially molten interior due to tidal heating is only possible
if the eccentricity required to do so is maintained over billions of
years. Studying the orbital evolution and concurrent interior thermal
evolution would shed light on this aspect.

By the same token, the day-side of $1.5R_⊕$ planets had a thicker
solid layer compared to Earth-sized planets, which had more mush
(see Figs 2c and d), again because of the higher gravitational potential
of Super-Earth's leads to greater compression of the interior layers.
We see this feature on the night-side as well. Magma and mush
oceans, if they can be maintained, will be shallower on super-Earths.

Applying the equations from Driscoll & Barnes (2015), we find
that tidal heating can maintain a molten night-side on a $1.5R_⊕$ planet
for eccentricity of approximately 0.004, within the uncertainty of
current orbital constraints (Barragán et al. 2018; Malavolta et al.
2018; Zieba et al. 2022). Precise night-side measurements of K2-
141b, and ideally lava planets of different ages, could provide
sensitive constraints on the thermal evolution of their deep interior.

ACKNOWLEDGEMENTS

NBC acknowledges support from an NSERC Discovery Grant, a Tier 2 Canada Research Chair, and an Arthur B. McDonald Fellowship. MH would like to thank the Fonds de recherche du Québec for a doctoral fellowship. The authors also thank the Trottier Space Institute and L'Institut Trottier de recherche sur les exoplanétes for their financial support and dynamic intellectual environment. The authors thank Edwin Kite for thoughtful feedback on the manuscript.

DATA AVAILABILITY

The data underlying this article are available in the article itself and
in any online supplementary material that may be made available.

REFERENCES

Barragán O. et al., 2018, A&A, 612, A95

Batalha N. M. et al., 2011, ApJ, 729, 27

Boukaré C.-E., Cowan N., Badro J., 2022, ApJ, 936, 148

Boukaré C.É., Lemasquerier D., Cowan N., Samuel H., Badro J., 2023,
preprint (arXiv:2308.13614)

Castan T., Menou K., 2011, ApJ, 743, L36

Chao K.-H., deGraffenried R., Lach M., Nelson W., Truax K., Gaidos E.,
2021, Geochemistry, 81, 125735

Dang L. et al., 2021, A Hell of a Phase Curve: Mapping the Surface and Atmosphere of a Lava Planet K2-141b, JWST Proposal, Cycle 1, ID., p. 2347

Driscoll P. E., Barnes R., 2015, Astrobiology, 15, 739

Elkins-Tanton L. T., 2012, Annu. Rev. Earth Planet. Sci., 40, 113

Espinoza N. et al., 2021, The first near-infrared spectroscopic phase-curve of a super-Earth. JWST Proposal, Cycle 1, ID., p. 2159

Fiquet G., Auzende A., Siebert J., Corgne A., Bureau H., Tateno H., Garbarino G., 2010, Science, 329, 1516

Hughes G. O., Griffiths R. W., 2008, Annu. Rev. Fluid Mech., 40, 185

Kervazo M., Tobie G., Choblet G., Dumoulin C., Běhounková M., 2021,
A&A, 650, A72

Kite E. S., Fogley Bruce J., Schaefer L., Gaidos E., 2016, ApJ, 828, 80

Labrosse S., 2003, Phys. Earth Planet. Inter., 140, 127

Labrosse S., 2015, Phys. Earth Planet. Inter., 247, 36

Lai Y., Yang J., Kang W., 2024a, Planet. sci. j., 5, 204

Lai Y., Kang W., Yang J., 2024b, Planet. sci. j., 5, 205

Le Losq C., Baldoni B., 2023, J. Non-Cryst. Solids, 617, 122481

Lebrun T., Massol H., ChassefièRe E., Davaille A., Marcq E., Sarda P.,
Leblanc F., Brandeis G., 2013, J. Geophys. Res. (Planets), 118, 1155

Léger A. et al., 2009, A&A, 506, 287

Leger A. et al., 2011, Icarus, 213, 1

Malavolta L. et al., 2018, AJ, 155, 107

Malosso C., Zhang L., Car R., Baroni S., Tisi D., 2022, Npj Comput. Mater.,
8, 139

McKenzie D., Nimmo F., Jackson J. A., Gans P. B., Miller E. L., 2000, J.
Geophys. Res., 105, 11029

Nguyen T. G., Cowan N. B., Banerjee A., Moores J. E., 2020, MNRAS, 499,
4605

Nguyen T. G., Cowan N. B., Pierrehumbert R. T., Lupu R. E., Moores J. E.,
2022, MNRAS, 513, 6125

Ohta K., Yagi T., Hirose K., Ohishi Y., 2017, Earth Planet. Sci. Lett., 465, 29

Patel J. A. et al., 2023, A&A, 679, A92

Schaefer L., Fogley B., 2009, ApJ, 703, L113

Schaefer L., Wordsworth R. D., Berta-Thompson Z., Sasselov D., 2016, ApJ,
829, 63

Solomatov V.S., 2000, Philos. Trans. A Math. Phys. Eng. Sci., 366, 4105

Tachinami C., Senshu H., Ida S., 2011, ApJ, 726, 70

Tyler R. H., Henning W. G., Hamilton C. W., 2015, ApJS, 218, 22

Xie L., Yoneda A., Katsura T., Andrault D., Tange Y., Higo Y., 2021, Geophys.
Res. Lett., 48, 105

Zhang J., Herzberg C., 1994, J. Geophys. Res., 99, 17

Zieba S. et al., 2022, A&A, 664, A79

APPENDIX A: TABLE OF DATA

This paper has been typeset from a TEX/LATEX file prepared by the author.
<table><caption>**Table A1.** Table of constants used in this study.</caption><thead><tr><th>Parameter</th><th>Description</th><th>Value</th><th>Units</th><th>Reference</th></tr></thead><tbody><tr><td>C</td><td>specific heat capacity of silicates</td><td>1200</td><td>J Kg<sup>-1</sup> K<sup>-1</sup></td><td>Schaefer et al. (2016)</td></tr><tr><td>ρ</td><td>Bulk density of silicates</td><td>4000</td><td>Kg m<sup>-3</sup></td><td>Lebrun et al. (2013)</td></tr><tr><td>σ</td><td>Stefan-Boltzmann constant</td><td>5.67 × 10<sup>-8</sup></td><td>W m<sup>-2</sup> K<sup>-4</sup></td><td>–</td></tr><tr><td>κ</td><td>thermal conductivity of magma</td><td>13</td><td>W m<sup>-1</sup> K<sup>-1</sup></td><td>Ohta et al. (2017)</td></tr><tr><td>α</td><td>thermal expansivity</td><td>3 × 10<sup>-5</sup></td><td>K<sup>-1</sup></td><td>Tachinami, Senshu & Ida (2011)</td></tr><tr><td>κ<sub>d</sub></td><td>thermal diffusivity</td><td>7.5 × 10<sup>-7</sup></td><td>m<sup>2</sup> s<sup>-1</sup></td><td>Schaefer et al. (2016)</td></tr><tr><td>η<sub>m</sub></td><td>viscosity of magma</td><td>0.1</td><td>Pa s</td><td>Boukaré et al. (2023)</td></tr></tbody></table>