Author: Dr. Megan Tannock

For a quick introduction to the SpaceEngine climate model and its main features, see our blog post here. In this article we provide a closer look at the physics and calculations behind planetary temperatures.

The temperatures of the SpaceEngine climate model are based on energy transport calculations and account for planetary albedo, presence of an atmosphere, atmosphere properties (including wind speeds, radiation and advection, and greenhouse effects), internal planetary heating, day sides, night sides, planet obliquity (for seasons and varying daylight hours, polar days and polar nights), eccentric orbits, tidal locking, and incident light of all stars in the system. Additionally, the altitude dependence had a major overhaul from its previous implementation. Now, we use real temperature-pressure profile data for different types of atmospheres, allowing for exciting behavior like temperature inversions (like we see at the tropopause, stratopause, and mesopause for Earth).

All planets are born hot, but without a sustained energy source like the fusion taking place inside of stars, they cool quickly. So unless planets are very young, almost all of their energy comes from stellar irradiation. Therefore, the main focus of SpaceEngine’s climate model is on stellar irradiation. We also account for internal planetary heating and greenhouse effects when an atmosphere is present.

The stellar part of the model is implemented in three main parts: 1) a starting temperature, based on the distance to and properties of the host star(s); 2) a longitudinal dependence based on the type of planet; and 3) a latitudinal dependence that also accounts for obliquity (axial tilt).

#### Where do we start? Calculating a Planetary Temperature

We start with the Planetary Equilibrium Temperature (the theoretical temperature of a planet if it were a blackbody heated only by its parent star):

T_{\mathrm{equilibrium}} = \left ( F_*\frac{(1-A)}{\sigma_{\mathrm{SB}} f} \right )^{1/4}= \left ( \frac{L_*}{4\pi d^2}\frac{(1-A)}{\sigma_{\mathrm{SB}} f} \right )^{1/4} ,

where $$F_*$$ and $$L_*$$ are the flux and luminosity of the host star, respectively, $$d$$ is the distance between the star and planet, $$A$$ is the bond albedo of the planet, $$\sigma_{SB}$$ is the Stefan-Boltzmann constant, and $$f$$ is a redistribution factor. Usually, $$f=4$$ for planets where all incident energy is uniformly distributed, and $$f=2$$ for planets where energy is distributed over the starlit hemisphere only, like a tidally locked planet (see, for example, Selsis et al. 2007 A&A 476).

This value is computed for each star in the system at the current point in the orbit ($$d$$) changes as time progresses for elliptical orbits) and is used to calculate a fiducial temperature ($$T_{\rm fiducial}$$) that a surface temperature map is built out from.

#### Day/Night Cycles and Longitude

We split planets into three categories: terrestrial planets, gas giant planets, and tidally locked planets (that can be terrestrial or gas giant). Each category has a different base temperature model in SpaceEngine. Our first step in all three categories was to model how the temperature changes around the planet with longitude, i.e.: how should the temperature change from day to night?

Starting with terrestrial planets with atmospheres, we used a thermal transport model to calculate temperatures across planet surfaces to account for day and night. We considered the temperature of a gas parcel as it travels over the planet surface and the planet rotates away from the host star and then back towards it. This gas parcel absorbs incident flux from the host star on the day side and radiates heat away during both day and night. We made the simple assumption that this gas parcel radiates as a blackbody. To determine the temperature at any time during the day or night, we set up and solved a thermodynamic energy equation (e.g.: Showman et al. 2009, Crossfield course notes 2019):

\frac{\mathrm{d}T}{\mathrm{d}t} = \frac{1}{c_p \rho H} \Delta F = \frac{1}{c_p \rho H} \left( \sigma_{\mathrm{SB}} {T_{\rm fiducial}}^4 T_{\rm latitude} ~{\rm max}(\cos \theta, 0) - \sigma_{\mathrm{SB}} T^4 \right) ,

where $$T$$ is the temperature we are solving for, $$t$$ is time, $$c_p$$ is the heat capacity (under constant pressure, since this gas parcel moves along a constant surface on the planet), $$\rho$$ is the atmospheric density, $$H$$ is the thickness of the atmosphere, and $$\Delta F$$ is the net flux on the gas parcel. $$T_{\rm fiducial}$$ (in Kelvin) is a function of starting planetary equilibrium temperature, $$T_{\rm latitude}$$ is a unitless latitude scaling (described below), $$\theta$$ is the longitude measured relative to the subsolar point (where the host star is directly overhead), and "max" is a function that gives the maximum of the two arguments (to account for the fact that when the parcel is on the night side of the planet is absorbs zero flux, not negative flux).

We rewrote this equation and substituted a longitude dependence for the time dependence so solving the equation resulted in a longitudinal temperature map for the planet, where $$\theta$$ is a proxy for time of day. We used a radiative time scale ($$\tau_{\rm rad}$$) and an advective time scale ($$\tau_{\rm adv}$$) to convert time to longitude:

\tau_{\rm rad} = \frac{c_p \rho H}{\sigma_{\mathrm{SB}} {T_0}^3},

and

\tau_{\rm adv} = \frac{2 \pi R_p}{v_{\rm wind}} ,

where $$T_0$$ is the fiducial temperature combined with the latitude scaling, $$R_p$$ is the planet’s radius, and $$v_{\rm wind}$$ is a global wind speed.

The ratio of these two timescales tells us how heat is carried around the planet. When $$\tau_{\rm rad} / \tau_{\rm adv} \gg 1$$, winds carry the heat around the planet before it is radiated away, the planet stays warmer overall. If $$\tau_{\rm rad} / \tau_{\rm adv} \ll 1$$, radiation is more efficient than advection, so energy is radiated away before it is carried around the planet; the planet stays cooler overall, and we get more “averaged” temperatures over the surface (the amplitude between the daytime maximum and nighttime minimum temperatures is smaller). This also produces a phase offset of the maximum temperature in our model. This means the hottest point on a planet’s surface is not necessarily at noon, when the host star is directly overhead (note that on Earth, winds are not the reason for this offset — on Earth it is warmer shortly after noon than directly at noon due to “thermal inertia.” We have not modeled that effect here, but our phase-offset due to winds makes a nice approximation for this in SpaceEngine).

The resulting longitudinal temperature map is composed of two solutions of the thermodynamic energy equation: a daytime function where the gas parcel absorbs incident flux and radiates away heat, and a night time function where the gas parcel only radiates away heat. The longitudinal temperature map of an arbitrary planet with varying wind speeds is shown in the following figure:

A longitudinal temperature map for a random terrestrial planet with an atmosphere. The temperature along the equator of a planet with 0 obliquity (axial tilt) for different global wind speeds is shown. There are three curves for a planet with identical atmospheric and orbital parameters in SpaceEngine and they only differ in global wind speed. $$\theta = 0^\circ, 360^\circ$$ is the subsolar point, or noon. $$\theta = 90^\circ$$ and $$\theta = 270^\circ$$ are the relative longitudes of sunset and sunrise, respectively (for a planet with 0 obliquity), and $$\theta = 180^\circ$$ is the relative longitude of midnight, when the host star is on the opposite side of the planet. Higher wind speeds means heat is carried around the planet, before it radiates away. This results in more average temperatures over the surface, and an offset of the maximum temperature from the position where the host-star is directly overhead.

For terrestrial planets without atmospheres, heat is not carried around the planet by moving air. For such planets, we removed the atmospheric component of the model and the resulting longitudinal temperature maps are similar to the “terrestrial with atmospheres” case, but with the hottest point being always at the subsolar point, and without the “averaging” effect over the surface of the planet. This results in a more extreme temperature variation from day to night.

Gas giant planets have massive atmospheres with significant convection, resulting in longitudinal temperature functions that are smoother, with a less sharp boundary at sunrise and sunset points. We adopted a sinusoidal function for the temperatures of gas giants, where the wind speed corresponds to a phase offset between the subsolar point and the hottest point on the surface of the planet (where, as per convention, the “surface” is the level in the atmosphere where the pressure is equal to 1 atm).

A longitudinal temperature map for the northern hemisphere of a gas giant planet on its northern summer solstice (when the host star’s subsolar point is at its highest latitude —see the next section for details on accounting for latitude!). As with the terrestrial temperature map, $$\theta = 0^\circ, 360^\circ$$ is the subsolar point, or noon. $$\theta = 90^\circ$$ and $$\theta = 270^\circ$$ are the relative longitudes of sunset and sunrise, respectively, at the equator, and $$\theta = 180^\circ$$ is the relative longitude of midnight, when the host star is on the opposite side of the planet. Following convention, the “surface” of a gas giant planet is the level in the atmosphere where the pressure is equal to 1 atm.

A longitudinal temperature map for a tidally locked planet. As with the other temperature maps, $$\theta = 0^\circ, 360^\circ$$ is the subsolar point, or noon and $$\theta = 180^\circ$$ is the relative longitude of midnight, where the host star is on the opposite side of the planet. For tidally locked planets, $$\theta = 90^\circ$$ and $$\theta = 270^\circ$$ are the relative longitudes of sunset and sunrise, respectively, for all latitudes.

Since tidally locked planets, as the name suggests, always have the same side facing the host star, our longitudinal temperature map for tidally locked planets is based on the Lambert cosine law. This law states that the relative intensity of light on a surface is proportional to the cosine of the angle of incidence. This results in a hot spot at the subsolar point, with temperatures dropping as angular distance from the subsolar point increases, with a fairly constant temperature across the (permanent) night side of the planet. Our model also incorporates a phase offset calculated from the planet’s wind speed. For tidally locked planets, the temperature, on the night side of the planet, is typically 1/3 the day time high (e.g.: Lunine & Lorenz 2002).

For a tidally locked planet on an eccentric orbit, the exact position of the subsolar point will oscillate over the planet’s surface with time. The subsolar point in SpaceEngine is not fixed to a particular latitude/longitude position on the planet’s surface, and this oscillation is included in our model.

#### Obliquity, Seasons, and Latitude

Our model includes a latitudinal dependence that fixes the maximum temperature at the subsolar latitude (for a planet with zero obliquity/no axial tilt, this is the equator) and a minimum temperature at the maximum relative latitude from the subsolar point (for a planet with zero obliquity/no axial tilt, this would be the poles).

We started with Lambert’s cosine law, which states that the relative intensity of light on a surface is proportional to the cosine of the angle of incidence. The angle in our case is the latitude ($$\phi$$) on the planet relative to the host star’s subsolar latitude. So, if a planet has zero obliquity, at the equator ($$\phi = 0^\circ$$) where the surface normal is parallel to the incident light, this is $$\cos 0^\circ = 1$$, and at the poles ($$\phi = \pm 90^\circ$$) we have $$\cos \pm 90^\circ = 0$$.

We multiplied our latitudinal dependence with our longitudinal dependence to scale the temperature appropriately, but wanted to modify this slightly. We didn’t want the temperature to go to 0 at the poles, and we wanted to be able to account for obliquities different from zero.

The first modification we made was to re-scale the cosine so we could accommodate obliquities that are not zero, and for latitude we used the absolute value of the latitude relative to the subsolar latitude. This means that on the solstices, when the subsolar latitude is equal to the obliquity (for Earth on the June solstice this would be $$\phi = 23.44^\circ$$), the winter pole is at a relative latitude of $$90^\circ + obliquity$$ and the summer pole is at $$90^\circ - obliquity$$. So our latitude dependence needed to go to $$90^\circ + obliquity$$, not just $$90^\circ$$ . Our next modification was to prevent the minimum value from reaching zero. We defined a unitless minimum temperature fraction for the winter pole ($$T_{\rm pole}$$) to scale our latitude dependence for temperature. This minimum temperature is based on the relationship between surface pressure, rotation period, and equator-to-pole temperature difference (e.g.: Kaspi & Showman 2015 ApJ 804), and is calculated for each planet/moon.

We applied this scaling to our longitudinal temperature map to account for the latitude in our temperature calculations.

Left panel: the cosine of $$\phi$$ goes to 0 at $$\phi = 90^\circ$$. Right panel: Our latitude dependence as a function of latitude relative to the subsolar latitude for a planet with obliquity = $$23.44^\circ$$. This function goes to a minimum value of $$T_{pole}$$ at $$\phi = 90^\circ + 23.44^\circ$$.

The next aspect related to latitude we accounted for is the varying daylight hours throughout the year for planets with non-zero obliquity/axial tilt. For example, on Earth’s June solstice, we know that the northern hemisphere experiences their longest day of the year, and the southern hemisphere experiences their shortest day of the year. At this time, the sun never sets inside of the northern polar circle (polar day), and the sun never rises inside of the southern polar circle (polar night).

The amount of time a point on the planet surface experiences incident sunlight affects the temperature, and is needed to set the boundaries between our day and night longitudinal temperature maps. To include this in our climate model, we needed to compute the fraction of the day a surface point is in sunlight at a given latitude. This is a straightforward calculation that can be done a few different ways, typically with the true anomaly (an orbital parameter). It can also be computed using the subsolar latitude:

fraction~of~day~receiving~sunlight = \frac{\arccos \left ( -\tan \phi \tan \phi_{SS} \right )}{\pi} ,

where $$\phi$$ is the latitude of interest and $$\phi_{SS}$$ is the subsolar latitude (e.g.: Rauscher 2017 ApJ 846).

The number of hours a given latitude experiences incident sunlight on Earth. Ignoring atmospheric effects and the angular size of the sun, at the equator there are equal hours of daylight and night on every day of the year. Other latitudes experience different hours of daylight throughout the year. On the vernal and autumnal equinoxes, all points on the Earth’s surface experience equal hours of daylight and night (ignoring the angular size of the Sun).

#### A Complete Surface Temperature Map

By combining the fiducial temperature, the longitudinal and latitudinal dependencies, and the varying length of day, we obtained a complete surface temperature map for a planet. Our model accounts for the properties of the host star, the planet’s orbit, obliquity/axial tilt (resulting in seasons!), rotation, and atmospheric properties (if an atmosphere is present). At this point, we also accounted for any internal heating from the planet itself, and greenhouse effects in the atmosphere.

We added a $$20^\circ$$ obliquity to the planet from the terrestrial longitudinal temperature map figure and applying the latitude dependence for a few sample latitudes, we obtain a more complete picture of the surface temperatures on a planet. The time of local sunrises and sunsets are marked for each latitude curve. For this planet, the atmospheric properties result in a small offset between noon and the hottest temperature of the day.

The surface temperatures for a planet with a $$20^\circ$$ obliquity on the northern summer solstice in SpaceEngine, as determined from our model, plotted on a sphere.

#### Tatooine and Beyond: Systems with Multiple Stars

The SpaceEngine climate model also accounts for systems with multiple stars! Whether your favorite planet is orbiting a star orbiting a star (an S-type orbit), orbiting two stars (a P-type orbit), or some other configuration with even more stars, it is accounted for in the model. Temperatures are computed as time progresses and as stars and planets move along their orbits. Binary planets and moons (and moons around moons!) are properly accounted for as well.

We treated each star in the system independently, calculating a longitudinal temperature map of the planet for each star. We then combined the effects of all stars in the system. Energies are additive, and energy is proportional to temperature to the power of 4, so our final temperature at a point on the surface will be:

T_{\mathrm{final,~current~position}} = [~( T_{\mathrm{current~position,~star~1}})^4 + ( T_{\mathrm{current~position,~star~2}})^4 + ( T_{\mathrm{current~position,~star~3}})^4 + \dots~]^{1/4}.

Any internal heating from the planet itself, and greenhouse effects in the atmosphere are only added to the model after all stars are accounted for.

#### Accounting for Altitude and Layered Atmospheres

The final part of calculating temperature was to account for the altitude. Once the surface temperature calculation included the incident light from all stars in the system, the longitude and latitude, the internal planetary heating, and any greenhouse effects due to an atmosphere, we calculated the temperature as a function of altitude.

The main challenge in accounting for altitude is that temperature depends on both altitude and pressure, but pressure also depends on altitude and temperature! In addition, vertical temperature profiles have varying slopes and inversions with altitude, because pressure and atmospheric chemical composition change with altitude.

To avoid the pressure-temperature mutual-dependence, in the SpaceEngine climate model we assumed pressure decreases exponentially with altitude (e.g.: Wikipedia - Scale Height).

Note: Following convention, zero altitude on gas giant planets is set to pressure=1 atm in SpaceEngine.

Then, to allow for temperature inversions and to account for different layers of atmospheres, we used real and simulated pressure-temperature profiles to determine the temperature as a function of pressure. The majority of the pressure-temperature profiles we use in SpaceEngine were computed using the NASA Planetary Spectrum Generator (Villanueva et al. 2018). Other profiles and simulated profiles were sourced from assorted publications (Thorngren et al. 2019, Piette & Madhusudhan et al. 2020, Zhang 2020, Ohno & Fortney 2023).

The suite of pressure-temperature profiles included in SpaceEngine. Surfaces of terrestrial solar system planets are marked with dots. When these profiles are selected for planets with higher surface pressures, or differing surface temperatures, the profiles are scaled appropriately in SpaceEngine.

For each planet, a suitable profile is selected, scaled for the surface pressure and temperature (calculated as described above), and then interpolated at the user’s current altitude/pressure position to determine the local temperature.

And that’s the final step in the temperature component of the new SpaceEngine Climate Model! We hope you enjoyed this peek inside the model.

For a look at how these quantities are displayed in-game, be sure to check out our related blog post here.

Note: None of the figures shown in this article are displayed "in-program": this is simply a look “under the hood” of the SpaceEngine climate model.