One-Dimensional Climate Models: The Basics of Planets¶
In this tutorial you will learn the very basics of running 1D climate runs. For a more in depth look at the climate code check out Mukherjee et al. 2022 (note this should also be cited if using this code/tutorial).
What you should already be familiar with:
What you will need to download to use this tutorial:
Download 196 Correlated-K Tables from Roxana Lupu to be used by the climate code for opacity (if you completed the brown dwarf tutorial you should already have this)
Note: the two files above are dependent on metallicity and C/O.
[1]:
import warnings
warnings.filterwarnings('ignore')
import picaso.justdoit as jdi
import picaso.justplotit as jpi
import astropy.units as u
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
[2]:
#1 ck tables from roxana
mh = '+100'#'+1.0' #log metallicity, 10xSolar
CtoO = '100'#'1.0' # CtoO ratio, Solar C/O
ck_db = f'/data/kcoeff_2020_v3/sonora_2020_feh{mh}_co_{CtoO}.data.196'
[3]:
# Notice The keyword ck is set to True because you want to use the correlated-k opacities for your calculation
# and not the line by line opacities
opacity_ck = jdi.opannection(ck_db=ck_db) # grab your opacities
Starting up the Run¶
You will notice that starting a run is nearly identical as running a spectrum and brown dwarf climate model. However, how we will add climate=True to our inputs flag. We will also specify planet in this case, which will turn on the irradiation the object is receiving.
New Parameter (was also used in the Brown Dwarf tutorial): Effective Temperature. This excerpt from Modeling Exoplanetary Atmospheres (Fortney et al) provides a thorough description and more reading, if you are interested.
If the effective temperature, \(T_{eff}\), is defined as the temperature of a blackbody of the same radius that would emit the equivalent flux as the real planet, \(T_{eff}\) and \(T_{eq}\) can be simply related. This relation requires the inclusion of a third temperature, \(T_{int}\), the intrinsic effective temperature, that describes the flux from the planet’s interior. These temperatures are related by:
\(T_{eff}^4 = T_{int}^4 + T_{eq}^4\)
We then recover our limiting cases: if a planet is self-luminous (like a young giant planet) and far from its parent star, \(T_{eff} \approx T_{int}\); for most rocky planets, or any planets under extreme stellar irradiation, \(T_{eff} \approx T_{eq}\).
[4]:
cl_run = jdi.inputs(calculation="planet", climate = True) # start a calculation
#note you need to put the climate keyword to be True in order to do so
# now you need to add these parameters to your calculation
tint= 200 # Intrinsic Temperature of your Planet in K
grav = 4.5 # Gravity of your Planet in m/s/s
cl_run.gravity(gravity=grav, gravity_unit=u.Unit('m/(s**2)')) # input gravity
cl_run.effective_temp(tint) # input effective temperature
Let’s now input the host-star properties
[5]:
T_star =5326.6 # K, star effective temperature
logg =4.38933 #logg , cgs
metal =-0.03 # metallicity of star
r_star = 0.932 # solar radius
semi_major = 0.0486 # star planet distance, AU
cl_run.star(opacity_ck, temp =T_star,metal =metal, logg =logg, radius = r_star,
radius_unit=u.R_sun,semi_major= semi_major , semi_major_unit = u.AU)#opacity db, pysynphot database, temp, metallicity, logg
Initial T(P) Guess¶
Every calculation requires an initial guess of the pressure temperature profile. The code will iterate from there to find the correct solution. A few tips:
We recommend using typically 51-91 atmospheric pressure levels. Too many pressure layers increases the computational time required for convergence. Too little layers makes the atmospheric grid too coarse for an accurate calculation.
Start with a guess that is close to your expected solution. One easy way to get fairly close is by using the Guillot et al 2010 temperature-pressure profile approximation
[6]:
nlevel = 91 # number of plane-parallel levels in your code
#Lets set the max and min at 1e-4 bars and 500 bars
Teq=1000 # planet equilibrium temperature
pt = cl_run.guillot_pt(Teq, nlevel=nlevel, T_int = tint, p_bottom=2, p_top=-6)
temp_guess = pt['temperature'].values
pressure = pt['pressure'].values
Initial Convective Zone Guess¶
You also need to have a crude guess of the convective zone of your atmosphere. Generally the deeper atmosphere is always convective. Again a good guess is always the published SONORA grid of models for this. But lets assume that the bottom 7 levels of the atmosphere is convective.
New Parameters:
nofczns: Number of convective zones. Though the code has functionality to solve for more than one. In this basic tutorial, let’s stick to 1 for now.rfacv: (See Mukherjee et al Eqn. 20r_st) https://arxiv.org/pdf/2208.07836.pdf
Non-zero values of rst (aka “rfacv” legacy terminology) is only relevant when the external irradiation on the atmosphere is non-zero. In the scenario when a user is computing a planet-wide average T(P) profile, the stellar irradiation is contributing to 50% (one hemisphere) of the planet and as a result rst = 0.5. If instead the goal is to compute a night-side average atmospheric state, rst is set to be 0. On the other extreme, to compute the day-side atmospheric state of a tidally locked planet rst should be set at 1.
[7]:
nofczns = 1 # number of convective zones initially. Let's not play with this for now.
nstr_upper = 85 # top most level of guessed convective zone
nstr_deep = nlevel -2 # this is always the case. Dont change this
nstr = np.array([0,nstr_upper,nstr_deep,0,0,0]) # initial guess of convective zones
# Here are some other parameters needed for the code.
rfacv = 0.5 #we are focused on a brown dwarf so let's keep this as is
Now we would use the inputs_climate function to input everything together to our cl_run we started.
[8]:
cl_run.inputs_climate(temp_guess= temp_guess, pressure= pressure,
nstr = nstr, nofczns = nofczns , rfacv = rfacv)
Run the Climate Code¶
The actual climate code can be run with the cl_run.run command. The save_all_profiles is set to True to save the T(P) profile at all steps. The code will now iterate from your guess to reach the correct atmospheric solution for your brown dwarf of interest.
[9]:
out = cl_run.climate(opacity_ck, save_all_profiles=True,with_spec=True)
Iteration number 0 , min , max temp 859.054155942838 2289.478427903662 , flux balance -46.15847181591005
Iteration number 1 , min , max temp 835.3574006260803 2777.6878509379635 , flux balance 17.82965190850441
Iteration number 2 , min , max temp 834.2722222145496 2653.4220211996176 , flux balance 0.2594411932152453
Iteration number 3 , min , max temp 834.2651124752114 2639.2957580270786 , flux balance 0.0012104268509236513
Iteration number 4 , min , max temp 834.2650775441722 2639.0924332789596 , flux balance 5.499585765882438e-06
In t_start: Converged Solution in iterations 4
Big iteration is 834.2650775441722 0
Iteration number 0 , min , max temp 834.0393684187401 2734.9589459834865 , flux balance -3.660621966940465
Iteration number 1 , min , max temp 833.7618986334929 2823.879709982079 , flux balance 0.12744719833424106
Iteration number 2 , min , max temp 833.7602101385852 2819.106469079621 , flux balance 0.0007067647384675832
In t_start: Converged Solution in iterations 2
Big iteration is 833.7602101385852 1
Iteration number 0 , min , max temp 834.6507021802025 2930.3758747219267 , flux balance -0.14424565099125855
Iteration number 1 , min , max temp 834.7484053262953 2933.850480135399 , flux balance -0.0002867399420656116
In t_start: Converged Solution in iterations 1
Big iteration is 834.7484053262953 2
Iteration number 0 , min , max temp 835.3849918816801 3018.1228722038486 , flux balance -0.0001937036980519201
Iteration number 1 , min , max temp 835.3874378060507 3014.557263577295 , flux balance 2.4716598291939565e-07
In t_start: Converged Solution in iterations 1
Profile converged before itmx
Iteration number 0 , min , max temp 835.7816725005212 3070.631061776754 , flux balance 8.60778296513762e-05
Iteration number 1 , min , max temp 835.7834656038983 3069.1460884516528 , flux balance 3.740806372413759e-07
In t_start: Converged Solution in iterations 1
Big iteration is 835.7834656038983 0
Iteration number 0 , min , max temp 836.0058826361692 3105.051279473138 , flux balance 4.2639690072613516e-05
Iteration number 1 , min , max temp 836.0069497529855 3104.474012120312 , flux balance 1.8151256338186683e-07
In t_start: Converged Solution in iterations 1
Profile converged before itmx
Move up two levels
Iteration number 0 , min , max temp 836.1272490822133 2814.4533222676664 , flux balance 5.112756588168638e-06
In t_start: Converged Solution in iterations 0
Big iteration is 836.1272490822133 0
Iteration number 0 , min , max temp 836.1924392218433 2814.4531485178004 , flux balance -8.640664528425539e-07
In t_start: Converged Solution in iterations 0
Profile converged before itmx
final [ 0 83 89 0 0 0]
Iteration number 0 , min , max temp 836.2271503391952 2814.453648097963 , flux balance -1.987533884891869e-06
In t_start: Converged Solution in iterations 0
Big iteration is 836.2271503391952 0
Iteration number 0 , min , max temp 836.2456450333358 2814.4538672423496 , flux balance -1.4208921069493005e-06
In t_start: Converged Solution in iterations 0
Profile converged before itmx
YAY ! ENDING WITH CONVERGENCE
[10]:
base_case = jdi.pd.read_csv(jdi.HJ_pt(), delim_whitespace=True)
plt.figure(figsize=(10,10))
plt.ylabel("Pressure [Bars]", fontsize=25)
plt.xlabel('Temperature [K]', fontsize=25)
plt.ylim(500,1e-4)
plt.xlim(200,3000)
plt.semilogy(out['temperature'],out['pressure'],color="r",linewidth=3,label="Our Run")
plt.semilogy(base_case['temperature'],base_case['pressure'],color="k",linestyle="--",linewidth=3,label="WASP-39 b ERS Run")
plt.minorticks_on()
plt.tick_params(axis='both',which='major',length =30, width=2,direction='in',labelsize=23)
plt.tick_params(axis='both',which='minor',length =10, width=2,direction='in',labelsize=23)
plt.legend(fontsize=15)
plt.title(r"T$_{\rm eff}$= 1000 K, log(g)=5.0",fontsize=25)
[10]:
Text(0.5, 1.0, 'T$_{\\rm eff}$= 1000 K, log(g)=5.0')
Brightness Temperature¶
Checking the brightness temperature serves many useful purposes:
Intuition building. Allows you to see what corresponding temperature are you sensitive to at each wavelength
Note that this temperature doesn’t need to be the physical temperature of your atmosphere but if you can find the physical converged atmospheric temperature closest to this brightness temperature you can also get an idea of the atmospheric pressure from where the flux you are seeing is originating from.
Determining if your choice in bottom boundary pressure grid was correct.
If your brightness temperature is such that you bottom out at the temperature corresponding to the highest pressure, you have not extended your grid to high enough pressures.
Brightness Temperature Equation:
\(T_{\rm bright}=\dfrac{a}{{\lambda}log\left(\dfrac{{b}}{F(\lambda){\lambda}^5}+1\right)}\)
where a = 1.43877735$:nbsphinx-math:times\(10\)^{-2}$ m.K and b = 11.91042952$:nbsphinx-math:times\(10\)^{-17}$ m\(^4\)kg/s\(^3\)
Let’s calculate the brigthness temperature of our current run and check if our pressure grid was okay.
[11]:
brightness_temp, figure= jpi.brightness_temperature(out['spectrum_output'])
Check Adiabat¶
This plot and datareturn is helpful to check that there have been no issues with where the code has found the location of the convective zone(s). See below, dTdP never exceeds the adiabat.
[12]:
adiabat, dtdp, pressure = jpi.pt_adiabat(out, cl_run)
