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 SiO_{2} atmospheres would impact future observations. The more volatile Na atmosphere is thickest followed by SiO and SiO_{2}, as expected. Despite its low vapour pressure, we find that a SiO_{2} 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 (SiO_{2}) 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 THEORY

### 2.1 The surface and below

*T*

_{s}) by Castan & Menou (2011) and Kite et al. (2016) accounted for the zero-albedo local radiative equilibrium temperature and geothermal heating:

*T*

_{ss}and

*T*

_{as}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.

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 (SiO_{2}). 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 SiO_{2} (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 SiO_{2} 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 SiO_{2} 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 SiO_{2} and may recombine with oxygen even before condensing, we finally use a pure SiO_{2} atmosphere.

*h*is the fluid’s column thickness, and |$v$| is the horizontal flow. Applying this equation to our case requires time invariance, |$\partial$|

_{t}

*h*= 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:

*r*is the planet’s radius (9.62 × 10

^{6}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:

_{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:

*C*

_{p}is the heat capacity, β is the ratio

*R*/(

*R*+

*C*

_{p}), 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:

*C*

_{p}

*T*. Considering 1D steady-state flow in spherical coordinates yields the equation for energy conservation:

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 RESULTS

### 3.1 The atmosphere

We used the model described above to determine the steady-state flow of a pure Na, SiO, and SiO_{2} 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.

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 SiO_{2} 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 SiO_{2} atmosphere is much warmer as the temperature remains above 1500 K.

Evaporation and condensation rates vary significantly between the Na, SiO, and SiO_{2} 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.

### 3.2 The magma ocean

*M*

_{T}) at some θ is given by

_{o}, that balances the atmospheric mass transport is therefore calculated via

^{−3}for SiO

_{2}or 2270 kg m

^{−3}for Na

_{2}O (Na likely exists as Na

_{2}O; 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 SiO

_{2}atmosphere, but not for the SiO and Na scenarios. To account for the oxygen dissociated from SiO

_{2}or Na

_{2}O, we simply scale the

*m*value in equation (7) to conserve the mass of oxygen during the evaporation and condensation processes: scaling factor = (1 +

*m*

_{O}/

*m*

_{SiO}) for SiO and (1 + 0.5 ×

*m*

_{O}/

*m*

_{Na}) for Na.

As the magma is predominantly SiO_{2}, the ocean circulation is dictated by our SiO atmosphere (recall that SiO has a much greater vapour pressure than SiO_{2}). 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 SiO_{2} to Na ratio similar to that of a bulk silicate Earth (BSE; where Na_{2}O 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 SiO_{2} 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 DISCUSSION

### 4.1 Implications for observation

We presented three scenarios where the atmosphere of K2-141b is either purely Na, SiO, or SiO_{2} 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 SiO_{2}, 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 SiO_{2} 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.

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 SiO_{2} atmosphere produces the highest radiative flux as it is significantly warmer than the other two. However, the SiO_{2} atmosphere will have fully condensed out before θ = 90°. This suggests that transit spectroscopy is impossible for a SiO_{2} 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 SiO_{2} 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°).

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.

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 SiO_{2} 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 SiO_{2} 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 10^{2} 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 (∼10^{2} W m^{−2}), silicon monoxide, and sodium (both ∼10^{5} W m^{−2}) are at least an order of magnitude smaller than the instellation (∼10^{6} 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. SiO_{2} −> 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 SiO_{2} or Na_{2}O dissociate to form O and O_{2}, both of which are more volatile than SiO_{2} 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.

## 5 CONCLUSIONS

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 SiO_{2} 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 SiO_{2} 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.

## ACKNOWLEDGEMENTS

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).

## DATA AVAILABILITY

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

## REFERENCES

Jr, ,

, , , , ,

, , , , , , , , ,

et al. , , , ,

, , , , , , ,

, , , , ,

### APPENDIX A: SiO_{2} PHASE DIAGRAM

See Fig. A1.

### APPENDIX B: GAS-SPECIFIC PARAMETERS

This subsection includes molecule-specific parameters of Na, SiO, and SiO_{2}. The mass per molecule for the three constituents are *m*_{Na} = 3.82 × 10^{−26} kg molecule^{−1}, *m*_{SiO} = 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 *C*_{p} of Na, SiO, and SiO_{2} 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* + *C*_{p}).

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}}$|.

### APPENDIX C: CALCULATING THE FLUXES

#### 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:

_{c}, which for any species sublimating into a near vacuum is |$\sqrt{kT_\mathrm{ s}/m}$|, where

*k*is the Boltzmann constant.

#### C2 Momentum

*V*

_{e}and

*V*

_{d}. 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.

*V*

_{e}represents the mean flow and is calculated by

*V*

_{e}=

*m*

*E*/ρ

_{s}.

*V*

_{d}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

*H*is the scaled height, η is the dynamic viscosity calculated by η(

*T*

_{s}) [kg m

^{−1}s

^{−1}] = 1.8 × 10

^{−5}(

*T*

_{s}/291)

^{3/2}(411/(

*T*

_{s}+ 120)) (Castan & Menou 2011).

#### C3 Energy

*Q*, a combination of enthalpy at surface and top of boundary layer and kinetic energy, is calculated by the following:

_{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):

### APPENDIX D: NUMERICALLY SOLVING THE SYSTEM OF EQUATIONS

*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

We solve for a solution through finite differences using a Runge–Kutta scheme to the second-order accuracy. After acquiring *f*_{state}, 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) = *T*_{ss}. 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.

### APPENDIX E: OBSERVATIONAL INFERENCE CALCULATIONS

*T*, to blackbody radiative intensity,

*I*[W m

^{−2}m

^{−1}sr

^{−1}], at wavelength λ:

*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:

*F*, is calculated by

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 (*r*_{star}/*r*_{planet})^{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.

*T*

_{b}is calculated from (Cowan & Agol 2011)