A Tutorial for the Bell Energy Balance Model (Bell_EBM)

Author: Taylor James Bell (McGill University)¶

Last Modified: 2019-01-24¶

import numpy as np
import matplotlib.pyplot as plt
import astropy.constants as const

import Bell_EBM as ebm

from Bell_EBM import pretty_plotting # Makes the plots look much nicer!

First, let's setup a planet, star, and system

We'll use the WASP-12b system and include the effects of H₂ dissociation (Bell & Cowan 2018)

# WASP-12b
# Periods are in units of days
# All other units are in SI (m, kg, s, K, J, Pa, etc.)
# All angles are in degrees
planet = ebm.Planet('bell2018', rad=1.900*const.R_jup.value, mass=1.470*const.M_jup.value,
                    Porb=1.09142030, a=0.02340*const.au.value, inc=83.37, vWind=5e3)

# WASP-12
# Radius and mass are in solar units
# Temperatures are always in K
star = ebm.Star(teff=6300., rad=1.59, mass=1.20)

# Load the star and planet into a system
system = ebm.System(star, planet)

Let's show the star-planet system¶

fig = planet.orbit.plot_orbit()
plt.show()

Now, let's model the planet's atmoshere¶

# Run initial burn-in - try guessing the equilibrium temperature to begin with
Teq = system.get_teq()
T0 = np.ones_like(system.planet.map.values)*Teq
t0 = 0.
t1 = t0+system.planet.Porb*2
dt = system.planet.Porb/1000.
times, maps = system.run_model(T0, t0, t1, dt, verbose=False)

Then we can visualize the planet's temperature and H₂ dissociation map¶

Even with only moderate winds, there is a strong hotspot offset due to the effects of H₂ dissociation/recombination¶

fig = system.planet.plot_map()
plt.show()

fig = system.planet.plot_H2_dissociation()
plt.show()

Finally, we can look at the orbital phase curve¶

# Plot the lightcurves - by default it will use one planetary orbit
fig = system.plot_lightcurve()
plt.show()

# Plot the temperature curves
fig = system.plot_tempcurve()
plt.show()

Now let's look at an eccentric rocky planet in a 3:2 spin-orbit resonance

# We can also specify a rotation rate for the planet: we'll say it somehow isn't
# tidally locked yet. We also don't have to specify the planet's orbital period:
# it will be calculated using Kepler's equations
planet = ebm.Planet('rock', rad=const.R_earth.value, mass=const.M_earth.value,
                    a=0.05*const.au.value, e=0.4, argp=45.)

# The default values are the sun's values already
star = ebm.Star()

# Load everything into the system: this is when the orbital period is calculated
system = ebm.System(star, planet)

# Now that we know the orbital period, we can set the planet's rotation period
system.planet.Prot = system.planet.Porb*(2./3.)

# The plotted dots have equal spacing in time, so the get scruntched up near apastron
fig = planet.orbit.plot_orbit()
plt.show()

# Run initial burn-in - try guessing the median irradiation temperature to begin with
Teq = np.median(system.get_teq(np.linspace(0.,system.planet.Porb,1000, endpoint=False)))
T0 = Teq*np.ones_like(system.planet.map.values)
t0 = 0.
t1 = t0+system.planet.Porb*10
dt = system.planet.Porb/1000.
times, maps = system.run_model(T0, t0, t1, dt, verbose=False)

Because this is an eccentric planet and it isn't synchronously rotating, it's map isn't constant with time.
As a result, we need to grab the map and time stamps for one orbit to do our plotting¶

# Run a production run with intermediates=True to get the times and maps at each time step
T0 = maps[-1]
t0 = times[-1]
t1 = t0+system.planet.Porb
dt = system.planet.Porb/1000.
times, maps = system.run_model(T0, t0, t1, dt, verbose=False, intermediates=True)

Let's show the map during periastron¶

For planet's which are not synchronously rotating, the longitude of zero in all plots refers to the sub-stellar point. This requires that we pass in the time of the map though, so it can rotate it as necessary¶

phases = system.get_phase(times)
phasePeri = system.get_phase_periastron()
indexPeri = np.argmin(np.abs(phases-phasePeri))

fig = system.planet.plot_map(maps[indexPeri], times[indexPeri])
plt.show()

This time, let's look at the lightcurve at the Spitzer/IRAC 4.5 micron band¶

# Plot the lightcurves
fig = system.plot_lightcurve(times, maps, bolo=False, wav=4.5e-6)
plt.show()

# Plot the temperature curves
fig = system.plot_tempcurve(times, maps, bolo=False, wav=4.5e-6)
plt.show()

It might seem odd that the coldest temperature on the plot above is near periastron when the planet is receiving the most flux, but this is because we are looking at the planet's cold nightside at that point in time¶

Now let's try to model the Earth (a pretty challenging task...) to demonstrate the flexibility, strengths, and weaknesses of this EBM.¶

We'll need to make some simplifying assumptions. We'll pretend that we only probe the atmosphere (pretending the atmosphere is opaque) and assume the atmosphere is made of pure N2. Feel free to instead model the Earth as pure rock or water, and see how the results change!¶

If you have healpy installed on your computer, you can also use healpix maps by setting "useHealpix=True" when initializing the planet. If the following cell fails for you, either install healpy or remove "useHealpix=True".¶

# The Earth's atmosphere is mostly N2, so we'll change the heat capacity
# of the gas from H2 (the default) to N2
cp_N2 = 1.039e3 # J/(kg K)

# We'll pretend the whole atmosphere absorbs and radiates
P0 = const.atm.value

a = const.au.value # AU

Rp = const.R_earth.value # m

Mp = const.M_earth.value # kg

Prot = 1. # days

albedo = 0.30

e = 0.0167

obliquity = 23.5 # degrees

# We'll assume that the Earth transits during the northern hemisphere's summer solstice
argp = -90.
t0 = 0.

planet = ebm.Planet('gas', rad=Rp, mass=Mp, Prot=Prot, a=a, e=e, obliq=obliquity,
                    argp=argp, t0=t0, cp=cp_N2, mlDepth=P0, albedo=albedo,
                    useHealpix=True)

# The default values are the sun's values already
star = ebm.Star()

# Load everything into the system
system = ebm.System(star, planet)

Given the extreme difference between the spin and orbital periods, burning in this planet is a little trickier and takes longer¶

# Run initial burn-in - try guessing the median irradiation temperature to begin with
# Run for two full-orbits to work out seasonal variations
T0 = system.get_teq()*np.ones_like(system.planet.map.values)
t0 = planet.t0
t1 = t0+system.planet.Porb*2.-system.planet.Prot
dt = system.planet.Prot/10.
times, maps = system.run_model(T0, t0, t1, dt, verbose=False)

# Run one day at high-rez to get good temperatures around the planet
T0 = maps[-1]
t0 = times[-1]
t1 = t0+system.planet.Prot
dt = system.planet.Prot/100.
times, maps = system.run_model(T0, t0, t1, dt, verbose=False)

Run a production run¶

T0 = maps[-1]
t0 = times[-1]
t1 = t0+system.planet.Porb
# We'll use a small time step to make sure we're doing a good job
dt = system.planet.Prot/10.
times, maps = system.run_model(T0, t0, t1, dt, verbose=False, intermediates=True)

# Then we'll cut down on the number of points, since our full-orbit phasecurves
# don't need that high of a time resolution. Just daily outputs are fine
times = times[::10]
maps = maps[::10]

Let's get the lightcurves ignoring and accounting for reflected light to demonstrate the difference at blue wavelengths¶

lc_reflected = system.lightcurve(times, maps, bolo=False, wav=0.5e-6,
                                 allowThermal=False, allowReflect=True)
lc_thermal = system.lightcurve(times, maps, bolo=False, wav=0.5e-6,
                               allowThermal=True, allowReflect=False)

x = times.flatten()
x = (x-planet.Porb*int(np.median(x/planet.Porb)))

fig = plt.figure(figsize=(12,4))
plt.plot(x, lc_reflected*1e9, c='b', label=r'$\rm Reflected~Flux$')
plt.plot(x, lc_thermal*1e9, c='r', label=r'$\rm Thermal~Flux$')
plt.ylabel(r'$F_p/F_* \rm ~(ppb)$')
plt.xlabel(r'$\rm Time~from~Transit~(days)$')
plt.xlim(0, 365)
plt.legend(loc='best')
plt.show()
plt.close(fig)

We're seeing here that at short wavelengths nearly all of the flux comes from reflected light, and the amount of light we see varies due to the Earth going through different phases (like the Moon does from our perspective). The blue curve just ends up being the Lambertian phase function¶

This time, let's look at the peak in Earth's thermal emission (~10 microns) where reflected light shouldn't matter near as much¶

lc_reflected = system.lightcurve(times, maps, bolo=False, wav=10e-6,
                                 allowThermal=False, allowReflect=True)
lc_thermal = system.lightcurve(times, maps, bolo=False, wav=10e-6,
                               allowThermal=True, allowReflect=False)

x = times.flatten()
x = (x-planet.Porb*int(np.median(x/planet.Porb)))

fig = plt.figure(figsize=(12,4))
plt.plot(x, lc_reflected*1e9, c='b', label=r'$\rm Reflected~Flux$')
plt.plot(x, lc_thermal*1e9, c='r', label=r'$\rm Thermal~Flux$')
plt.ylabel(r'$F_p/F_* \rm ~(ppb)$')
plt.xlabel(r'$\rm Time~from~Transit~(days)$')
plt.xlim(0, 365)
plt.legend(loc='best')
plt.show()
plt.close(fig)

We see that there is a significant offset in the peak of the phasecurve due to the rapid rotation of the planet, and we also see that thermal emission indeed dominates¶

While this EBM is quite flexible, it is by no means perfect, and you should evaluate your findings with caution and scepticism (as with all things in science).¶

Let's show the planet's map at the time of the last timestep in the production run (which is during the Northern Hemisphere's Summer Solstice)¶

fig = system.planet.plot_map()
plt.show()

As you can see, the map doesn't appear to make a lot of sense - the south pole is far too cold (-100°C), the tropics are too cold, and the maximum temperature is only ~6°C. This is because we haven't taken into account many important features like any north-south heat transport (e.g. Hadley circulation), ocean currents, spatially varying albedo, the contribution of flux from the surface (the Earth's atmosphere is not opaque), latent heat, greenhouse gasses...: the list goes on. This EBM is pretty flexible, but it is definitely idealized!¶

We do get the year-averaged equator-to-pole temperature difference (ETPD) in roughly the correct ballpark though. For reference, on earth the ETPD is ~40°C for the north pole, and ~70°C for the south pole.¶

meanMap = np.mean(maps, axis=0)
EPTD = np.max(meanMap)-np.min(meanMap) # Average Equator-to-pole temperature difference
print('Average equator-to-pole temperature difference is ~'+str(int(np.rint(EPTD)))+'°C.')

Average equator-to-pole temperature difference is ~66°C.

A Tutorial for the Bell Energy Balance Model (Bell_EBM)

Author: Taylor James Bell (McGill University)¶

Last Modified: 2019-01-24¶

First, let's setup a planet, star, and system

We'll use the WASP-12b system and include the effects of H₂ dissociation (Bell & Cowan 2018)

Let's show the star-planet system¶

Now, let's model the planet's atmoshere¶

Then we can visualize the planet's temperature and H₂ dissociation map¶

Even with only moderate winds, there is a strong hotspot offset due to the effects of H₂ dissociation/recombination¶

Finally, we can look at the orbital phase curve¶

Now let's look at an eccentric rocky planet in a 3:2 spin-orbit resonance

Because this is an eccentric planet and it isn't synchronously rotating, it's map isn't constant with time.
As a result, we need to grab the map and time stamps for one orbit to do our plotting¶

Let's show the map during periastron¶

For planet's which are not synchronously rotating, the longitude of zero in all plots refers to the sub-stellar point. This requires that we pass in the time of the map though, so it can rotate it as necessary¶

This time, let's look at the lightcurve at the Spitzer/IRAC 4.5 micron band¶

It might seem odd that the coldest temperature on the plot above is near periastron when the planet is receiving the most flux, but this is because we are looking at the planet's cold nightside at that point in time¶

Now let's try to model the Earth (a pretty challenging task...) to demonstrate the flexibility, strengths, and weaknesses of this EBM.¶

We'll need to make some simplifying assumptions. We'll pretend that we only probe the atmosphere (pretending the atmosphere is opaque) and assume the atmosphere is made of pure N2. Feel free to instead model the Earth as pure rock or water, and see how the results change!¶

If you have healpy installed on your computer, you can also use healpix maps by setting "useHealpix=True" when initializing the planet. If the following cell fails for you, either install healpy or remove "useHealpix=True".¶

Given the extreme difference between the spin and orbital periods, burning in this planet is a little trickier and takes longer¶

Run a production run¶

Let's get the lightcurves ignoring and accounting for reflected light to demonstrate the difference at blue wavelengths¶

We're seeing here that at short wavelengths nearly all of the flux comes from reflected light, and the amount of light we see varies due to the Earth going through different phases (like the Moon does from our perspective). The blue curve just ends up being the Lambertian phase function¶

This time, let's look at the peak in Earth's thermal emission (~10 microns) where reflected light shouldn't matter near as much¶

We see that there is a significant offset in the peak of the phasecurve due to the rapid rotation of the planet, and we also see that thermal emission indeed dominates¶

While this EBM is quite flexible, it is by no means perfect, and you should evaluate your findings with caution and scepticism (as with all things in science).¶

Let's show the planet's map at the time of the last timestep in the production run (which is during the Northern Hemisphere's Summer Solstice)¶

We do get the year-averaged equator-to-pole temperature difference (ETPD) in roughly the correct ballpark though. For reference, on earth the ETPD is ~40°C for the north pole, and ~70°C for the south pole.¶

Have fun, explore, learn!

The full API can be found at https://bell-ebm.readthedocs.io

The most up-to-date version of the package can be found at https://github.com/taylorbell57/Bell_EBM

Taylor Bell may be contacted (taylor DOT bell AT mail DOT mcgill DOT ca) to report any issues discovered with the code

A Tutorial for the Bell Energy Balance Model (Bell_EBM)

Author: Taylor James Bell (McGill University)¶

Last Modified: 2019-01-24¶

First, let's setup a planet, star, and system

We'll use the WASP-12b system and include the effects of H2 dissociation (Bell & Cowan 2018)

Let's show the star-planet system¶

Now, let's model the planet's atmoshere¶

Then we can visualize the planet's temperature and H2 dissociation map¶

Even with only moderate winds, there is a strong hotspot offset due to the effects of H2 dissociation/recombination¶

Finally, we can look at the orbital phase curve¶

Now let's look at an eccentric rocky planet in a 3:2 spin-orbit resonance

Because this is an eccentric planet and it isn't synchronously rotating, it's map isn't constant with time. As a result, we need to grab the map and time stamps for one orbit to do our plotting¶

Let's show the map during periastron¶

For planet's which are not synchronously rotating, the longitude of zero in all plots refers to the sub-stellar point. This requires that we pass in the time of the map though, so it can rotate it as necessary¶

This time, let's look at the lightcurve at the Spitzer/IRAC 4.5 micron band¶

It might seem odd that the coldest temperature on the plot above is near periastron when the planet is receiving the most flux, but this is because we are looking at the planet's cold nightside at that point in time¶

Now let's try to model the Earth (a pretty challenging task...) to demonstrate the flexibility, strengths, and weaknesses of this EBM.¶

We'll need to make some simplifying assumptions. We'll pretend that we only probe the atmosphere (pretending the atmosphere is opaque) and assume the atmosphere is made of pure N2. Feel free to instead model the Earth as pure rock or water, and see how the results change!¶

If you have healpy installed on your computer, you can also use healpix maps by setting "useHealpix=True" when initializing the planet. If the following cell fails for you, either install healpy or remove "useHealpix=True".¶

Given the extreme difference between the spin and orbital periods, burning in this planet is a little trickier and takes longer¶

Run a production run¶

Let's get the lightcurves ignoring and accounting for reflected light to demonstrate the difference at blue wavelengths¶

We're seeing here that at short wavelengths nearly all of the flux comes from reflected light, and the amount of light we see varies due to the Earth going through different phases (like the Moon does from our perspective). The blue curve just ends up being the Lambertian phase function¶

This time, let's look at the peak in Earth's thermal emission (~10 microns) where reflected light shouldn't matter near as much¶

We see that there is a significant offset in the peak of the phasecurve due to the rapid rotation of the planet, and we also see that thermal emission indeed dominates¶

While this EBM is quite flexible, it is by no means perfect, and you should evaluate your findings with caution and scepticism (as with all things in science).¶

Let's show the planet's map at the time of the last timestep in the production run (which is during the Northern Hemisphere's Summer Solstice)¶

We do get the year-averaged equator-to-pole temperature difference (ETPD) in roughly the correct ballpark though. For reference, on earth the ETPD is ~40°C for the north pole, and ~70°C for the south pole.¶

Have fun, explore, learn!

The full API can be found at https://bell-ebm.readthedocs.io

The most up-to-date version of the package can be found at https://github.com/taylorbell57/Bell_EBM

Taylor Bell may be contacted (taylor DOT bell AT mail DOT mcgill DOT ca) to report any issues discovered with the code

We'll use the WASP-12b system and include the effects of H₂ dissociation (Bell & Cowan 2018)

Then we can visualize the planet's temperature and H₂ dissociation map¶

Even with only moderate winds, there is a strong hotspot offset due to the effects of H₂ dissociation/recombination¶

Because this is an eccentric planet and it isn't synchronously rotating, it's map isn't constant with time.
As a result, we need to grab the map and time stamps for one orbit to do our plotting¶