Transit searches have uncovered Earth-size planets orbiting so close to their host star that their surface should be molten, so-called lava planets. We present idealized simulations of the atmosphere of lava planet K2-141b and calculate the return flow of material via circulation in the magma ocean. We then compare how pure Na, SiO, or SiO2 atmospheres would impact future observations. The more volatile Na atmosphere is thickest followed by SiO and SiO2, as expected. Despite its low vapour pressure, we find that a SiO2 atmosphere is easier to observe via transit spectroscopy due to its greater scale height near the day–night terminator and the planetary radial velocity and acceleration are very high, facilitating high dispersion spectroscopy. The special geometry that arises from very small orbits allows for a wide range of limb observations for K2-141b. After determining the magma ocean depth, we infer that the ocean circulation required for SiO steady-state flow is only 10−4 m s−1, while the equivalent return flow for Na is several orders of magnitude greater. This suggests that a steady-state Na atmosphere cannot be sustained and that the surface will evolve over time.

Lava planets are a class of rocky exoplanets that orbit so close to their star that parts of their surface are molten (for a review of ultrashort-period planets, see Winn, Sanchis-Ojeda & Rappaport 2018). Indeed, dayside temperatures can be hot enough to maintain a rock vapour atmosphere detectable through transit and eclipse spectroscopy, as well as phase curves (Schaefer & Fegley 2009). Theoretical studies of lava planets and their relatively volatile sodium atmospheres have been published for CoRoT-7b (Léger et al. 2011), Kepler-10b, 55 Cnc e (Castan & Menou 2011), KIC 12556548, and HD 219134b (Kite et al. 2016). These studies utilized a 1D model that solves for the flow of mass, momentum, and energy from the hot dayside to the cold nightside.

We study K2-141b, the highest signal-to-noise ratio lava planet discovered to date, as its phase variations have been measured by K2 (Barragán et al. 2018; Malavolta et al. 2018), Spitzer (Kreidberg et al. 2019), and possibly in the future with the James Webb Space Telescope (JWST; e.g. Stevenson et al. 2016). We use essentially the same model as the above studies, but we employ more accurate surface temperature calculations. We consider three compositional end-members for the atmosphere: (1) pure sodium (Na), due to its high volatility, (2) silicon monoxide (SiO), or (3) pure silica (SiO2) due to their abundance in the crust of rocky planets. These three cases bracket the range of plausible atmospheric compositions (Schaefer & Fegley 2009). We focus on the mass exchange between surface and atmosphere to provide insights on the magma ocean return flow. We also use the results to anticipate future observations with space telescopes and ground-based high-dispersion spectroscopy.


2.1 The surface and below

We adopt planetary parameters from Barragán et al. (2018) and Malavolta et al. (2018). The subsolar equilibrium temperature (the planet’s irradiation temperature) is 3038 ± 64 K; we adopt the round number of 3000 K for the remainder of this study. Previous estimates of the surface temperature (Ts) by Castan & Menou (2011) and Kite et al. (2016) accounted for the zero-albedo local radiative equilibrium temperature and geothermal heating:
$$\begin{eqnarray*} T_\mathrm{ s} = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}(T_{\rm ss} - T_{\rm as}) \cos ^{1/4}(\theta) + T_{\rm as} & \mathrm{for}\ \theta \le 90{^\circ }, \\ T_{\rm as} & \text{for}\ \theta \gt 90{^\circ }, \end{array}\right. \end{eqnarray*}$$
where Tss and Tas are the temperature at the subsolar and antisolar point, respectively, and θ is the angular distance from the subsolar point. The cos 1/4θ expression, however, is only valid in the limit of a distant star, R*/a << 1. K2-141b has such a tight orbit that the star has an angular diameter of 50° as seen from the planet and the surface can still be illuminated well past θ = 90°. Here, we employ a more accurate analytic expression for the surface temperature originally developed for binary star systems (Kopal 1954). The surface temperature of K2-141b calculated using these two schemes is shown in Fig. 1. Beyond θ ≃ 115°, the illumination flux is so low that the temperature is mainly dependent on a geothermal flux, which we neglect. Our model will show that the atmosphere fully condenses out well before this point.

Top: the illumination received by K2-141b as a function of angular distance from the substellar point. Bottom: equilibrium temperature. In both panels, the black line shows the cos1/4θ flux dependence adopted by previous researchers, while the purple line shows more accurate formulation that accounts for the large angular size of the star as seen from the planet.

With the surface temperature now calculated, we can infer physical properties of the magma ocean by making a few assumptions and simplifications. As silicate materials dominate the composition of rocky planets (Schaefer & Fegley 2009) and K2-141b’s bulk density is comparable to that of the Earth (Malavolta et al. 2018), we assume that K2-141b’s crust is entirely silica (SiO2). We further assume that the temperature is vertically constant below the surface for the top 150 km. We neglect geothermal heating, which is reasonable if the ocean is relatively shallow and vertically well mixed.

To determine the magma ocean depth, we combine the planet’s known surface gravity (21.8 m s−2) with the density of SiO2 (2650 kg m−3) to obtain the pressure as a function of depth. These calculated pressures and temperatures then allow us to predict the phase (liquid versus solid) of the SiO2 as a function of depth using the phase diagram shown in Fig. A1 (Swamy et al. 1994). Using these data, we calculate the depth of magma ocean as a function of the surface temperature at the top of the column and hence as a function of θ, also shown in Fig. A1. In Section 3.2, we will close the loop between the atmosphere and the surface: mass must be redistributed via ocean currents to account for the constant evaporation and condensation.

2.2 The atmosphere

As a short-period planet on a circular orbit, K2-141b is expected to be tidally locked into synchronous rotation, leading to a permanent dayside and nightside. Not only does this introduce symmetry in how the planet is illuminated, it also leads to equilibrium dynamics (Schaefer & Fegley 2009). The atmosphere of a lava planet originates from the sublimation/evaporation of the surface and the flow is maintained by the pressure gradient between the nightside and dayside. The solutions are therefore dependent on the surface temperature and solely a function of the angular distance from the subsolar point (θss).

We assume a hydrostatically bound ‘thin’ atmosphere with negligible absorption of radiation. Although radiative transfer simulations of mineral-vapour atmospheres suggest that this is not the case (Ito et al. 2015), it is a useful first approximation and makes the problem tractable. We also neglect Coriolis forces to keep the problem one-dimensional although our estimated Rossby number of 10−1 suggests that Coriolis forces may be stronger than on other lava planets (Castan & Menou 2011). Lastly, we assume that the atmosphere is turbulent and well mixed, leading the vertically uniform wind velocity, V. Other state variables, pressure P and temperature T, are values calculated at the top of the boundary layer (Ingersoll, Summers & Schlipf 1985). We therefore use the shallow-water equations, commonly used for hydrodynamical calculations of non-linear evaporation/sublimation processes (e.g. Ingersoll et al. 1985; Ding & Pierrehumbert 2018).

We compute the surface temperature ahead of time because, as stated above, we assume it is unaffected by the atmosphere. The atmosphere is idealized to be pure and condensible because of the model’s limits. We first consider a pure Na atmosphere, neglecting dissociation processes that unbind Na from other molecules as is done by previous studies of lava planets (e.g. Castan & Menou 2011; Kite et al. 2016). Na is relatively scarce in a rocky planet’s crust (Schaefer & Fegley 2009), whereas SiO2 is the dominant constituent. Therefore, we also run the model using SiO in the same spirit as Na. Since gas-phase SiO will likely condense as SiO2 and may recombine with oxygen even before condensing, we finally use a pure SiO2 atmosphere.

The model was derived from the shallow-water equations by Ingersoll et al. (1985) and adapted by Castan & Menou (2011) for exoplanets. It describes the conservation of mass (equation 3), momentum (equation 5), and energy (equation 6). Mass conservation is as follows:
$$\begin{eqnarray*} \partial _t h + \nabla \, (\rho hv) = 0, \end{eqnarray*}$$
where ρ is the fluid density, h is the fluid’s column thickness, and |$v$| is the horizontal flow. Applying this equation to our case requires time invariance, |$\partial$|th = 0, hydrostatic equilibrium, h = P/(ρg), and a source/sink to account for evaporation/condensation such that the right-hand side of equation (2) is no longer 0 but some mass flux rate, E [molecules s−1 m−2]. The final step is to change to spherical coordinates that gives our mass conservation equation:
$$\begin{eqnarray*} \frac{1}{r \sin (\theta)}\frac{\mathrm{ d}}{\mathrm{ d}\theta }\left (\frac{V \, P \sin (\theta)}{g}\right) = m \, E, \end{eqnarray*}$$
where r is the planet’s radius (9.62 × 106 m), m is the mass per molecule, and g is the surface gravity (21.8 m s−2). We can do the same for momentum, which is affected by pressure gradient forces and surface drag τ [kg m−1 s−2]. This relation has the general form:
$$\begin{eqnarray*} \partial _t (\rho v) + \nabla \, (\rho v^2) = -\nabla P + \tau . \end{eqnarray*}$$
A steady-state flow means that |$\partial$|t|$v$|⁠) is 0. Substituting the mass for a product of V, P, and g, integrating ∇P vertically, collecting like terms, and changing to spherical coordinates yield the momentum conservation equation:
$$\begin{eqnarray*} \frac{1}{r \sin (\theta)}\frac{\mathrm{ d}}{\mathrm{ d}\theta }\left (\frac{(V^2 + \beta \, C_p \, T)P \sin (\theta)}{g}\right) = \frac{\beta \, C_p \, T \, P}{g \, r \tan (\theta)} + \tau , \end{eqnarray*}$$
where Cp is the heat capacity, β is the ratio R/(R + Cp), and R is the gas constant. Parameters specific to each atmospheric constituent are included in Appendix  B. Energy is conserved if the vertically integrated kinetic energy flux, internal heat energy, gravitational energy, and work done via adiabatic expansion are balanced out by energy exchanged with the surface, Q [W m−2]. The sum of the internal heat, gravity, and work terms becomes the enthalpy per unit mass at the top of the boundary layer: CpT. Considering 1D steady-state flow in spherical coordinates yields the equation for energy conservation:
$$\begin{eqnarray*} \frac{1}{r \sin (\theta)}\frac{\mathrm{ d}}{\mathrm{ d}\theta }\left (\frac{(V^2/2 + \beta \, C_p \, T)V \, P \sin (\theta)}{g}\right) = Q. \end{eqnarray*}$$

Calculations of fluxes E, τ, and Q are included in Appendix  C. Detailed derivation and steps shown between the generalized equations and our conservation equations can be found in Ingersoll et al. (1985). As the system only depends on θ through the state variables P [Pa], V [m s−1], and T [K], the equilibrium solution is obtained by iteratively integrating the fluxes. This computationally solves the ordinary differential equations as an initial boundary condition problem (details in Appendix  D).


3.1 The atmosphere

We used the model described above to determine the steady-state flow of a pure Na, SiO, and SiO2 atmosphere for lava planet K2-141b. Atmospheric pressure, wind velocity, and temperature for the three cases are shown in Fig. 2. Na is the most volatile molecule, making the Na atmosphere the thickest of the three with the maximum pressure at the subsolar point of 13.8 kPa. The Na atmosphere fully condenses out at θ = 102° where wind speed can reach up to 2.3 km s−1 (Mach 30). The atmospheric temperature of Na can drop to near zero before fully collapsing.

Top panel: atmospheric boundary-layer pressure for a pure sodium (Na), silicon monoxide (SiO), and silica (SiO2) atmosphere. Middle panel: wind velocity (left/black) and Mach number (right/purple). Bottom panel: atmospheric and surface temperatures (solid purple line).

SiO is not as volatile as Na making the SiO atmosphere thinner with a maximum pressure of 5.25 kPa. The wind, however, accelerates much faster than the Na case with velocities reaching up to 2.2 km s−1 and the atmosphere condenses out at θ = 92°. The temperature also drops much faster. Running the model for SiO2 produces the thinnest atmosphere with a maximum pressure of 240 Pa and collapses at θ = 89°. The wind is slower with a maximum speed of 1.75 km s−1 but the SiO2 atmosphere is much warmer as the temperature remains above 1500 K.

Evaporation and condensation rates vary significantly between the Na, SiO, and SiO2 cases. Although Na and SiO have similar evaporation rates near the subsolar point, the SiO atmosphere starts to condense sooner than its Na counterpart. The evaporation rates for the three scenarios are shown in Fig. 3.

Left/black: evaporation rate per square metre (negative numbers denote condensation). Right/purple: velocity of magma ocean return flow. Negative values denote mass flow towards the subsolar point (in the negative θ direction). Solid purple lines denote calculations neglecting the transport of oxygen, while dashed lines include oxygen to maintain mass balance.

3.2 The magma ocean

To maintain steady-state atmospheric flow, material must be transported back to the substellar region along the surface, most likely through currents in the magma ocean. Horizontal mass transport is determined by calculating the cumulative mass that has evaporated and subsequently condensed. Therefore, the atmospheric mass transport (MT) at some θ is given by
$$\begin{eqnarray*} M_\mathrm{ T}(\theta) = \int _{0}^{\theta } m E(\theta) \, \mathrm{ d}\theta . \end{eqnarray*}$$
The ocean return flow, |$v$|o, that balances the atmospheric mass transport is therefore calculated via
$$\begin{eqnarray*} v_{\mathrm{ o}}(\theta) = M_\mathrm{ T}(\theta)/(\rho \, A(\theta)), \end{eqnarray*}$$
where ρ is the density of the magma ocean: 2650 kg m−3 for SiO2 or 2270 kg m−3 for Na2O (Na likely exists as Na2O; Miguel et al. 2011). A(θ) is the cross-sectional area of the magma ocean given the magma ocean depth (estimated in Section 2.1 and shown in Fig. A1). Equations (7) and (8) predict the ocean velocity without accounting for the dissociation of oxygen. This approximation may be fine for the SiO2 atmosphere, but not for the SiO and Na scenarios. To account for the oxygen dissociated from SiO2 or Na2O, we simply scale the m value in equation (7) to conserve the mass of oxygen during the evaporation and condensation processes: scaling factor = (1 + mO/mSiO) for SiO and (1 + 0.5 × mO/mNa) for Na.

As the magma is predominantly SiO2, the ocean circulation is dictated by our SiO atmosphere (recall that SiO has a much greater vapour pressure than SiO2). These calculations yield a maximum magma ocean return flow of 2–3 × 10−4 m s−1, very slow compared to ocean circulation on the Earth. The evaporation rate and inferred ocean velocities are shown as in Fig. 3.

If we assume a SiO2 to Na ratio similar to that of a bulk silicate Earth (BSE; where Na2O accounts for 0.7 per cent of the BSE), then Na accounts for 1.52 per cent by mass of our idealized crust (Miguel et al. 2011). Applying the same ratio to the magma ocean, we calculate the ocean circulation required to maintain mass balance for Na that resulted in a maximum velocity of 2–3 cm s−1, two orders of magnitude faster than the expected circulation in the magma ocean. While not much mass condenses beyond the magma ocean for our SiO and SiO2 cases, a significant amount of Na is deposited on to solid ground. Here, material can no longer be transported by ocean circulation, but rather by solid flow akin to glaciers flowing into the ocean on the Earth.


4.1 Implications for observation

We presented three scenarios where the atmosphere of K2-141b is either purely Na, SiO, or SiO2 and modelled the pressure, flow, and temperature for each. While qualitatively similar, the results differ significantly depending on the atmospheric constituent, suggesting that observations could distinguish between the three scenarios. For example, the JWST (Gardner et al. 2006) could observe the near-infrared (NIR) phase curve of K2-141b, while the Hubble Space Telescope could observe it in the visible spectrum. Spectroscopy spanning the NIR SiO2, infrared (IR) SiO, or visible Na features could simultaneously probe the planetary surface and atmospheric temperature.

In Fig. 4, we show phase variations at wavelengths in and out of the 589 nm Na feature; we likewise show phase variations at wavelengths in and out of the 3.4 |$\mu$|m SiO2 feature (Coblentz Society & Craver 1977) and 10.1 |$\mu$|m SiO feature (Rawlings 1968). Phase variations are determined by calculating the expected wavelength-dependent blackbody radiative flux of the planet throughout an orbit. This process, detailed in Appendix  E, involves transforming our 1D results to a 2D map, which we then integrate over the planet’s visible hemisphere to get the expected flux measurements.

Top: predicted blackbody radiation of K2-141b at 589 nm Na spectral feature (left/black) and phase-dependent disc-integrated planetary brightness temperature (right/purple). Middle: the same for the 3.4 |$\mu \mathrm{ m}$| SiO2 spectral feature. Bottom: the same for the 10.1 |$\mu \mathrm{ m}$| SiO spectral feature. Note that eclipses and transit have been omitted but would occur at 0 and 0.5, respectively. The phase curve is shown as the flux ratio between the planet and star. Solid lines refer to the surface (a wavelength just outside the spectral feature), and dashed lines refer to the atmosphere (a wavelength in the spectral feature).

As seen in Fig. 4, the nightside flux and brightness temperature of both the surface and the Na and SiO atmospheres do not go to 0 at inferior conjunction (orbital phase of 0.5). This is because the visible hemisphere is still somewhat illuminated, and the Na and SiO atmospheres have not fully condensed. Despite being the thinnest, the SiO2 atmosphere produces the highest radiative flux as it is significantly warmer than the other two. However, the SiO2 atmosphere will have fully condensed out before θ = 90°. This suggests that transit spectroscopy is impossible for a SiO2 atmosphere, unlike Na and SiO, as the atmosphere is entirely hidden from view on the planet’s dayside. We show below that this is not necessarily the case.

In addition to its effect on illumination patterns and surface temperature, the especially close orbit of K2-141b has implications for transit geometry. It is generally assumed that transit spectroscopy only probes the atmosphere 90° from the subsolar point. As illustrated in Fig. 5, regions probed during transit can be quite far from θ = 90° for ultrashort-period planets. The start and end of transit therefore probe regions with a significant atmosphere even in the SiO2 case. For the specific case of K2-141b and assuming an impact parameter of zero (consistent with observational constraints), ingress and egress probe equatorial regions at an angle tan −1(R*/a) ≃ 26° away from the usually assumed θ = 90° (for sightlines offset from the planetary equator, the observed θ is closer to 90°).

Illustration of the observed longitude during a transit (the observer is at the bottom of the page). The star and planetary orbit are drawn to scale.

This opens up an avenue for high-dispersion spectroscopy where Doppler shifts from the changing radial velocity of the planet can help disentangle planetary and stellar absorption features and potentially detect atmospheric winds (Snellen 2014). Note that the extreme orbit of K2-141b leads to significant changes in radial velocity over the course of the transit. The radial velocity and the observed scale height along the planetary equator are shown in Fig. 6.

Anatomy of a transit of K2-141b, in degrees from transit. Left/black: observed scale height for the Na, SiO, and SiO2 cases for regions along the equator. Even though the SiO2 atmosphere collapses before θ = 90° and hence is invisible in the centre of eclipse, the scale height of the SiO2 atmosphere is bigger and observable for most of the planet’s transit. Right/purple: the planet’s velocity along the observer’s line of sight (negative velocities denote motion towards the observer).

If we assume the composition of K2-141b to be similar to the BSE, then there are 100× more silicate material than Na, by mass (Miguel et al. 2011). However, our results show that the mass transport required for a steady-state flow of Na is much greater than the ≃10−4 m s−1 provided by the SiO2 magma ocean. This implies that a more realistic Na atmosphere is much thinner than what we predicted above because the steady-state flow cannot be sustained, leading to the ‘evolved’ state anticipated by Kite et al. (2016). Also, any precipitation that falls beyond the magma ocean shore (θ > 79°) must be brought back to maintain mass conservation. This may occur via solid state flow analogous to glaciers on the Earth, or via isostatic adjustment. If mass is transported back too slowly, the resulting mass imbalance could lead to reorientation of the planetary spin (Leconte 2018).

4.2 Model limitations

Our model makes promising predictions for observation of the lava planet K2-141b. However, our idealized approach to characterize K2-141b’s atmosphere and magma ocean has limits that must be addressed. This subsection will discuss these limitations and list the caveats that come with our expected observation measurements as seen in Figs 4 and 6.

Our proposed pure or nearly pure SiO2 magma ocean is likely affected by interior dynamics that we neglect. The calculated ocean flow only accounts for the bulk mass transport along θ to balance the evaporation and condensation rates. We also neglect Coriolis forces despite the fast rotation of K2-141b. Qualitatively, Coriolis force would deflect the wind en route to the cold nightside and the mass transport of both the atmosphere and the ocean would be constrained to a smaller area, presumably near the equator. Despite this, there should still be net mass transport in the atmosphere as the evaporated volatiles would still move from an area of high vapour pressure to areas of low vapour pressure. This also applies to the magma ocean as mass would need to recirculate to maintain net mass balance between the surface and atmosphere.

If atmospheric dynamics were dictated by both vapour pressure gradients and Coriolis forces, then it would lead to the formation of Rossby waves and possibly weather systems. The ocean flow is slower than the atmosphere, so it would be even more affected by Coriolis forces. Ultimately, a global circulation model (GCM) is required to accurately infer atmospheric dynamics that arise from the planet’s rotation. However, this is difficult for a lava planet because the GCM would need to handle supersonic winds in an atmosphere that fully collapses far away from the subsolar point.

Another approximation is the fixed surface temperature that allowed for the energy of the surface and atmosphere to be decoupled. This adds significant stability to the model especially as supersonic flow is guaranteed. Coupling the energy budgets of the surface and atmosphere in post-processing can help to estimate the actual surface temperature. We calculated how much energy it would take to heat up the advected ocean mass to the equilibrium surface temperature from the model’s results. This is done by taking the bulk transported mass per θ and multiplying it with the heat capacity and the temperature gradient. This energy peaks at  102 W m−2 that is minuscule compared to other energy contributions.

Latent heat, may be a more significant source of energy for the surface. But the maximum latent heat produced by the evaporation of silica (∼102 W m−2), silicon monoxide, and sodium (both ∼105 W m−2) are at least an order of magnitude smaller than the instellation (∼106 W m−2). This suggests that latent heat would barely affect the surface temperature, in agreement with Castan & Menou (2011) and Kite et al. (2016). It remains to be seen whether the energy of dissociation (e.g. SiO2 −> SiO + O) is as important for lava planets as it has proven to be for ultrahot Jupiters (Bell & Cowan 2018; Tan & Komacek 2019; Mansfield et al. 2020).

Another important limitation if our model is the lack of consideration for radiative processes in the atmosphere. Although the work of Ito et al. (2015) suggests that there should be strong radiation absorption in the upper atmosphere, their results are sensitive to specific atmospheric compositions, stellar spectra in the ultraviolet, and atomic/molecular line lists, all of which are difficult to know for exoplanets.

The composition of the atmosphere is also much idealized as we ignore chemical processes where SiO2 or Na2O dissociate to form O and O2, both of which are more volatile than SiO2 but less volatile than Na. This has huge implications on transit spectroscopy as each molecule has its own distinct spectral features. Even in our idealized case where we focused on evaporation and condensation at great length, we neglected macroscale meteorological processes, chiefly cloud formation. This also greatly affects observations as clouds can trap heat via longwave radiation absorption but reflect shortwave radiation.


K2-141b is the best target to study an exotic situation where a planet’s interior, ocean, and atmosphere are compositionally similar. It will soon be observable with the JWST and with ground-based telescopes equipped with high-resolution spectrographs. Although previous studies focused on more volatile species such as Na or K (Castan & Menou 2011; Kite et al. 2016), we have shown that the significantly less volatile SiO2 atmosphere may be easier to observe due to the extreme geometry of this system. Transit spectroscopy of K2-141b is possible due to its proximity to its star that makes dayside regions accessible at the start and end of transit, and leads to large acceleration throughout transit. Doppler shifts in transmission spectra could confirm the km s−1 wind speed of the atmosphere. Spectral phase measurements should be able to distinguish between the planetary surface and the somewhat cooler atmosphere. Since radiative transfer simulations predict a dayside temperature inversion (Ito et al. 2015), phase curves may show a transition from inverted to non-inverted spectra.

Our calculated return flow of SiO2 through the magma ocean is slow (10−4 m s−1), so mass balance can be easily achieved. This supports a steady-state silicate flow, whereas the volatile Na atmosphere would be limited by the return flow through the magma ocean. Although K2-141b is an especially good target for atmospheric observation, our results regarding atmospheric dynamics, replenishment via the magma ocean, and transit geometry may apply to other lava planets.


We would like to thank Raymond Pierrehumbert for his comments on the pre-submission draft and mentorship in hydrodynamical modelling. This work was made possible by the Natural Science and Engineering Research Council (NSERC) of Canada’s Collaborative Research and Training Experience Program (CREATE) for Technology for Exo-Planetary Sciences (TEPS). This work is also made possible by Mathematics of Information Technology and Complex Systems (Mitacs) Globalink Internship program. NBC acknowledges support from the McGill Space Institute (MSI) and l’Institut de recherche sur les exoplanètes (iREx).


The data underlying this paper will be shared on reasonable request to the corresponding author.


   Jr, ,

J. Phys. Chem. Ref. Data Monograph
, ,

  ,   , , , ,  

  ,   ,   ,   ,   ,   , , , ,

   et al. , , , ,

  ,   ,   ,   ,   ,   , ,

J. Analytical At. Spectrometry
, ,

  ,   ,   ,   , ,

J. Geophys. Res.: Solid Earth
, ,


See Fig. A1.

Top: phase diagram (black) and depth (purple) of a SiO2 magma ocean. Left axis shows the pressure while the right axis shows the depth at which the phase will change between solid and liquid given the temperature and pressure. The inflection points reflect transitions between different SiO2 crystal structures that are, in increasing pressure: cristobalite, β-quartz, coesite, and stishovite (Swamy et al. 1994). Bottom: surface temperature (left/black) and magma ocean depth (right/purple) as a function of angular distance from subsolar point.


This subsection includes molecule-specific parameters of Na, SiO, and SiO2. The mass per molecule for the three constituents are mNa = 3.82 × 10−26 kg molecule−1, mSiO = 7.32 × 10−26 kg molecule−1, and |$m_{\rm SiO_2} = 9.98\times 10^{-26}$| kg molecule−1. With this, one can compute the respective gas constant R with R = k/m, where k is the Boltzmann constant.

The heat capacity Cp of Na, SiO, and SiO2 is 904 J K−1 kg−1 (Castan & Menou 2011), 851 J K−1 kg−1 (evaluated from the Shomate equation at 2400 K provided by Chase 1998), and 1428 J K−1 kg−1 (Lim et al. 2002), respectively. With this, β can be calculated with the gas constant R, where β = R/(R + Cp).

The final component is the saturated vapour pressure approximated from the model of Schaefer & Fegley (2009) where |$P_\mathrm{ v}^{\rm Na} \, [\mathrm{ Pa}] = 10^{9.6}\,\mathrm{ e}^{-38000/T_\mathrm{ s}}$|⁠, |$P_\mathrm{ v}^{\rm SiO} \, [\mathrm{ Pa}] = 6.16\times 10^{13}\, \mathrm{ e}^{-69400/T_\mathrm{ s}}$|⁠, and |$P_\mathrm{ v}^{\rm SiO_2} \, [\mathrm{ Pa}] = 1.07\times 10^{12} \, \mathrm{ e}^{-66600/T_\mathrm{ s}}$|⁠.


C1 Mass

E is the mass flux from equation (3), which is essentially calculated by the evaporation/sublimation of the surface through the difference between the saturated vapour pressure and the atmospheric pressure:
$$\begin{eqnarray*} E = \frac{(P_\mathrm{ v}-P)v_\mathrm{ c}}{\sqrt{2\pi } \, m \, R \, T_\mathrm{ s}}. \end{eqnarray*}$$
The kinetic speed of sublimation/evaporation is denoted by |$v$|c, which for any species sublimating into a near vacuum is |$\sqrt{kT_\mathrm{ s}/m}$|⁠, where k is the Boltzmann constant.

C2 Momentum

The momentum flux τ is calculated by the following:
$$\begin{eqnarray*} \tau = \rho _\mathrm{ s} (-V \, w_\mathrm{ a}), \end{eqnarray*}$$
where |$\omega _\mathrm{ a} \, [\mathrm{m\, s}^{-1}]$|⁠, the transfer coefficient as defined by Ingersoll et al. (1985), is dependent on two velocity values: Ve and Vd. There is a negative sign to denote that this force always acts against the flow. The vapour density, ρ, is determined via the ideal gas law. Ve represents the mean flow and is calculated by Ve = mEs. Vd represents the effects of eddies that is calculated by |$V_\mathrm{ d} = V_*^2/V$|⁠, where V* is the friction velocity; this must be calculated iteratively from
$$\begin{eqnarray*} V_{*,n+1} = \frac{V}{2.5 \log(9 \, V_{*,n} \, H \, \rho _\mathrm{ s}/(2\eta))}, \end{eqnarray*}$$
where H is the scaled height, η is the dynamic viscosity calculated by η(Ts) [kg m−1 s−1] = 1.8 × 10−5(Ts/291)3/2(411/(Ts + 120)) (Castan & Menou 2011).

C3 Energy

The energy flux Q, a combination of enthalpy at surface and top of boundary layer and kinetic energy, is calculated by the following:
$$\begin{eqnarray*} Q = \rho _\mathrm{ s}(w_\mathrm{ s} \, C_p \, T_\mathrm{ s} - w_\mathrm{ a} \, (V^2/2 + C_p \, T)), \end{eqnarray*}$$
where |$w$|s is the other transfer coefficient defined by Ingersoll et al. (1985). Both |$w$|a and |$w$|s are evaluated from equations (C5) to (C8):
$$\begin{eqnarray*} w_\mathrm{ a} = \frac{V_\mathrm{ e}^2 - 2V_\mathrm{ d}\mathrm{ } \, V_\mathrm{ e} + 2V_\mathrm{ d}^2}{-V_\mathrm{ e} + 2 V_\mathrm{ d}}, \ V_\mathrm{ e} \lt 0, \end{eqnarray*}$$
$$\begin{eqnarray*} w_\mathrm{ a} = \frac{2V_\mathrm{ d}^2}{V_\mathrm{ e} + 2 V_\mathrm{ d}\mathrm{ }}, \ V_\mathrm{ e} \ge 0, \end{eqnarray*}$$
$$\begin{eqnarray*} w_\mathrm{ s} = \frac{2V_\mathrm{ d}^2}{-V_\mathrm{ e} + 2 V_\mathrm{ d}}, \ V_\mathrm{ e} \le 0, \end{eqnarray*}$$
$$\begin{eqnarray*} w_\mathrm{ s} = \frac{V_\mathrm{ e}^2 - 2V_\mathrm{ d} \, V_\mathrm{ e} + 2V_\mathrm{ d}^2}{V_\mathrm{ e} + 2 V_\mathrm{ d}}, \ V_\mathrm{ e} \gt 0. \end{eqnarray*}$$


The state variables (P, V, T) can be calculated given initial values at the boundary (subsolar and antisolar point) by integration of fluxes in the right-hand side of equations (3)–(6). If we treat these equations as a single equation, we get
$$\begin{eqnarray*} \frac{\mathrm{ d}}{\mathrm{ d}\theta }(f_{\rm state}(\theta)) = f_{\rm fluxes}(\theta). \end{eqnarray*}$$
The state variable can be calculated at the next spatial step (the change in angular distance Δθ) through finite differences,
$$\begin{eqnarray*} f_{\rm state}(\theta + \Delta \theta) = f_{\rm fluxes}(\theta) \, \Delta \theta + f_{\rm state}(\theta). \end{eqnarray*}$$

We solve for a solution through finite differences using a Runge–Kutta scheme to the second-order accuracy. After acquiring fstate, the state variables can be solved arithmetically but due to the quadratic nature of V in the equations, two branches can be extracted corresponding to the flow being either subsonic or supersonic. Although, a higher order scheme could remedy the sensitivity of modelling fluid transitioning from subsonic to supersonic flow, the calculations to extract the state variables from equations (3)–(6) would take longer. Therefore, we elect to employ other schemes to control the instability inherent in the model, as described below.

We use the shooting method where we guess the initial P(θ = 0) and find a solution that fits the boundary conditions: V is 0 at θ = 0° and 180°, T(θ = 0) = Tss. The instability of the non-linear system of equations can lead to solutions either exponentially expanding or decaying. The goal, therefore, is to solve for a solution that exists in between the bounds of solutions that blew up (upper) or decayed to 0 (lower) exponentially. Conditions that define whenever a solution expands or decays exponentially is explicitly described by Ingersoll et al. (1985).

Solving for velocity will yield two solutions, one for subsonic flow and the other for supersonic flow. Getting to the critical point where the flow transition between different subsonic and supersonic regimes can be difficult. Therefore, we dynamically increase the spatial resolution whenever the upper bound and lower bound solution start to diverge past a set parameter. A solution close to the ‘true’ solution (which lies in between the upper and lower bound) will approach Mach 1 more reliably, after which we extrapolate the flow into the supersonic regime. This extrapolation process is done by linearizing the state variables with respect to θ until Mach 1 is achieved. The flow is then calculated using the supersonic branch of the solutions.


The predicted blackbody radiation of K2-141b required to produce Fig. 4 is calculated using the methods from Cowan & Agol (2008, 2011). The first step is to convert temperature, T, to blackbody radiative intensity, I [W m−2 m−1 sr−1], at wavelength λ:
$$\begin{eqnarray*} I = \frac{2 h c^2}{\lambda ^5} \frac{1}{\mathrm{ e}^{h c/(\lambda k T)} - 1}, \end{eqnarray*}$$
where h is the Planck constant, c is the speed of light, and k is the Boltzmann constant. The next step is to create an intensity map with latitude (α) and longitude (ϕ) from angular distance (θ) from the subsolar point using the spherical Pythagorean theorem:
$$\begin{eqnarray*} I(\alpha ,\phi) = I(\theta = \cos ^{-1}[\cos (\alpha)\cos (\phi)]). \end{eqnarray*}$$
Now that each latitude (α) and longitude (ϕ) coordinate has a radiative intensity, we then integrate the intensity across the hemisphere visible to the observer. Assuming an edge-on viewing geometry where the observer is on the planet’s orbital plane, the observed radiative flux, F, is calculated by
$$\begin{eqnarray*} F(\xi) = \int _{0}^{\pi } \int _{\xi -\pi /2}^{\xi +\pi /2} I(\alpha ,\phi) \sin ^2 \alpha \cos (\phi -\xi) \, \mathrm{ d}\alpha \, \mathrm{ d}\phi , \end{eqnarray*}$$

 where ξ is the planet’s position around the orbit. The stellar flux at the corresponding wavelength is calculated in the exact same way but with a uniform temperature of 4599 K (Malavolta et al. 2018) and scaled to the ratio of (rstar/rplanet)2. This therefore yields the ratio of planetary to stellar flux at that point as a function of the planet’s orbital motion about its star, the so-called phase curve.

Brightness temperature Tb is calculated from (Cowan & Agol 2011)
$$\begin{eqnarray*} T_\mathrm{ b} = \frac{h c}{\lambda k}\left [ \log \left (1+\frac{\mathrm{ e}^{hc/\lambda k T_\mathrm{ b}^*} -1}{\psi } \right) \right ]^{-1}, \end{eqnarray*}$$
where |$T_\mathrm{ b}^*$| is the star’s effective temperature (4599 K) and ψ is the radiative flux ratio between the planet and the star.