Atmosphere Tutorial
This run mode generates a 1D atmospheric model (pressure, temperature,
abundances, and altitude profiles), and saves it to a file. At
minimum, the user needs to provide the arguments required to compute
the pressure and temperature profiles. Further, Pyrat Bay
will
compute volume-mixing ratio (abundance) profiles only if the user sets
the chemistry
and respective arguments. Likewise, the code will
compute the altitude profiles only if the user provides the
radmodel
and respective arguments (which also require that the
abundance profiles to be defined).
Regardless of which profiles are computed, in an interactive run the
code returns a five-element tuple containing the pressure profile
(barye), the temperature profile (Kelvin), the abundance profiles
(volume mixing fraction), the species names, and the altitude profile
(cm). The outputs that were not calculated are set to None
.
Also, regardless of the input units, the output variables will always
be in CGS units.
Pressure Profile
The pressure profile model is an equi-spaced array in log-pressure,
determined by the pressure at the top of the atmosphere ptop
, at
the bottom of the atmosphere pbottom
, and the number of layers
nlayers
.
The units for the ptop
and pbottom
pressures may be defined
in-place along with the variable value (e.g., pbottom = 100 bar
)
or through the punits
key (in-place units take precedence over
punits
). See Units for a list of available units.
Temperature Profiles
Currently, there are three temperature models (tmodel
argument) to
compute temperature profiles. Each one requires a different set of
parameters (tpars
), which is described in the table and sections below:
Models ( |
Parameters ( |
References |
---|---|---|
isothermal |
\(T_0\) |
— |
guillot |
\(\log_{10}(\kappa')\), \(\log_{10}(\gamma_1)\), \(\log_{10}(\gamma_2)\), \(\alpha\), \(T_{\rm irr}\), \(T_{\rm int}\) |
|
madhu |
\(\log_{10}(p_1)\), \(\log_{10}(p_2)\), \(\log_{10}(p_3)\), \(a_1\), \(a_2\), \(T_0\) |
Isothemal profile
The isothermal model is the simplest model, having just one free
parameter (tpars
): the temperature (\(T_0\)) at all layers.
Here is an example of an isothermal atmosphere configuration file:
[pyrat]
# Pyrat Bay run mode, select from: [tli atmosphere spectrum opacity mcmc]
runmode = atmosphere
# Output atmospheric file:
atmfile = pt_isothermal_tutorial.atm
# Pressure at the top and bottom of the atmosphere, and number of layers:
ptop = 1.0e-6 bar
pbottom = 100.0 bar
nlayers = 100
# Temperature-profile model, select from [isothermal guillot madhu]
tmodel = isothermal
tpars = 990.0
# Verbosity level (<0:errors, 0:warnings, 1:headlines, 2:details, 3:debug):
verb = 2
Guillot profile
The guillot model has six parameters as defined in [Line2013]: \(\log_{10}(\kappa')\), \(\log_{10}(\gamma_1)\), \(\log_{10}(\gamma_2)\), \(\alpha\), \(T_{\rm irr}\) and \(T_{\rm int}\). The temperature profile is given as:
with
where \(E_{2}(\gamma_{i}\tau)\) is the second-order exponential integral; \(T_{\rm int}\) is the internal heat temperature; and \(\tau(p) = \kappa' p\) (note that this parameterization differs from that of [Line2013], which are related as \(\kappa' \equiv \kappa/g\)). \(T_{\rm irr}\) parametrizes the stellar irradiation absorbed by the planet:
Here is an example of a guillot
atmosphere configuration file
(pt_guillot.cfg):
[pyrat]
# Pyrat Bay run mode, select from: [tli atmosphere spectrum opacity mcmc]
runmode = atmosphere
# Output atmospheric file:
atmfile = pt_guillot_tutorial.atm
# Pressure at the top and bottom of the atmosphere, and number of layers:
ptop = 1.0e-6 bar
pbottom = 100.0 bar
nlayers = 100
# Temperature-profile model, select from [isothermal guillot madhu]
tmodel = guillot
# log(kappa) log(g1) log(g2) alpha T_irr T_int
tpars = -6.0 -0.25 0.0 0.0 950.0 100.0
# Verbosity level (<0:errors, 0:warnings, 1:headlines, 2:details, 3:debug):
verb = 2
Madhu profile
The madhu model has six parameters: \(\log_{10}(p_1)\), \(\log_{10}(p_2)\), \(\log_{10}(p_3)\), \(a_1\), \(a_2\), and \(T_0\), as described in [Madhusudhan2009], where the pressure values must be given in bars. The temperature profile is given as:
A thermally inverted profile will result when \(p_1 < p_2\); a non-inverted profile will result when \(p_2 < p_1\). The pressure parameters must also satisfy: \(p_1 < p_3\).
Here is an example of a madhu
atmosphere configuration file
(pt_madhu.cfg):
[pyrat]
# Pyrat Bay run mode, select from: [tli atmosphere spectrum opacity mcmc]
runmode = atmosphere
# Output atmospheric file:
atmfile = pt_madhu_tutorial.atm
# Pressure at the top and bottom of the atmosphere, and number of layers:
ptop = 1.0e-6 bar
pbottom = 100.0 bar
nlayers = 100
# Temperature-profile model, select from [isothermal guillot madhu]
tmodel = madhu
# log(p1) log(p2) log(p3) a1 a2 T0
tpars = -1.0 -3.6 0.50 0.95 0.67 850.0
# Verbosity level (<0:errors, 0:warnings, 1:headlines, 2:details, 3:debug):
verb = 2
Temperature-profile Examples
Note
Before running this example, download the configuration files shown above, e.g., with these shell commands:
tutorial_path=https://raw.githubusercontent.com/pcubillos/pyratbay/master/examples/tutorial
wget $tutorial_path/pt_isothermal.cfg
wget $tutorial_path/pt_guillot.cfg
wget $tutorial_path/pt_madhu.cfg
The following Python script creates and plots the pressure-temperature profiles for the configuration files shown above:
import matplotlib.pyplot as plt
plt.ion()
import pyratbay as pb
import pyratbay.constants as pc
# Generate PT profiles:
press, t_iso = pb.run("pt_isothermal.cfg")[0:2]
press, t_guillot = pb.run("pt_guillot.cfg")[0:2]
press, t_madhu = pb.run("pt_madhu.cfg")[0:2]
# Plot the PT profiles:
plt.figure(11)
plt.clf()
plt.semilogy(t_iso, press/pc.bar, color='mediumblue', lw=2, label='isothermal')
plt.semilogy(t_guillot, press/pc.bar, color='orange', lw=2, label='guillot')
plt.semilogy(t_madhu, press/pc.bar, color='tab:green', lw=2, label='madhu')
plt.ylim(100, 1e-5)
plt.xlim(800, 1200)
plt.legend(loc="best", fontsize=12)
plt.xlabel("Temperature (K)", fontsize=12)
plt.ylabel("Pressure (bar)", fontsize=12)
plt.tight_layout()
And the results should look like this:
Note
If any of the required variables is missing form the
configuration file, Pyrat Bay
will throw an error
indicating the missing value, and stop executing the
run. Similarly, Pyrat Bay
will throw a warning for a
missing variable that was defaulted, and continue
executing the run.
Abundance Profiles
Currently, there are two chemistry models (chemistry
argument) to
compute volume-mixing-ratio abundances: uniform or thermochemical
equilibrium. Each one requires a different set of arguments, which is
described in the table and sections below:
Models ( |
Required arguments [optional arguments] |
References |
---|---|---|
uniform |
|
— |
tea |
|
Uniform abundances
To produce a uniform-abundance model, the configuration file must
contain the species
key specifying a list of the name of the
species to include in the atmosphere, and the uniform
key
specifying the mole mixing fraction for each of the species listed in
species
. An atmosphere config file must also set the atmfile
key specifying the output atmospheric file name.
Here is an example of a uniform atmosphere configuration file (atmosphere_uniform.cfg):
[pyrat]
# Pyrat Bay run mode, select from: [tli atmosphere spectrum opacity mcmc]
runmode = atmosphere
# Pressure at the top and bottom of the atmosphere, and number of layers:
ptop = 1.0e-6 bar
pbottom = 100.0 bar
nlayers = 100
# Temperature-profile model, select from [isothermal guillot madhu]
tmodel = guillot
# log(kappa') log(g1) log(g2) alpha T_irr T_int
tpars = -6.0 -0.25 0.0 0.0 950.0 100.0
# Output atmospheric model:
atmfile = WASP-00c.atm
# Chemistry model, select from [uniform tea]
chemistry = uniform
# Atmospheric composition and abundances:
species = H2 He Na H2O CH4 CO CO2 NH3 HCN N2
uniform = 0.85 0.149 5e-6 2e-3 1e-4 5e-4 1e-6 1e-5 1e-9 2e-4
# Verbosity level (<0:errors, 0:warnings, 1:headlines, 2:details, 3:debug):
verb = 2
Thermochemical-equilibrium Abundances
Pyrat Bay
computes abundances in thermochemical equilibrium via
the TEA package [Blecic2016], by minimizing the Gibbs free energy at
each layer. To produce a TEA model, the configuration file must
contain the species
key specifying the species to include in the
atmosphere, the elements
key specifying the elemental composition.
An atmosphere config file must also set the atmfile
key specifying
the output atmospheric file name.
The TEA run assumes a solar elemental composition from [Asplund2009];
however, the user can scale the abundance of metals by setting the
xsolar
key, or can scale individual elemental abundances by setting
the escale
key element and scale factor pairs:
Here is an example of a thermochemical-equilibrium atmosphere configuration file (atmosphere_tea.cfg):
[pyrat]
# Pyrat Bay run mode, select from: [tli atmosphere spectrum opacity mcmc]
runmode = atmosphere
# Pressure at the top and bottom of the atmosphere, and number of layers:
ptop = 1.0e-6 bar
pbottom = 100.0 bar
nlayers = 100
# Temperature-profile model, select from [isothermal guillot madhu]
tmodel = guillot
# log(kappa') log(g1) log(g2) alpha T_irr T_int
tpars = -6.0 -0.25 0.0 0.0 950.0 100.0
# Output atmospheric model:
atmfile = WASP-00b.atm
# Chemistry model, select from [uniform tea]
chemistry = tea
# Output atmospheric composition:
species = H He C N O Na H2 H2O CH4 CO CO2 NH3 HCN N2
# Set elemental abundance of all metals (dex units, relative to solar):
# (2x solar metallicity for everything except H, He)
metallicity = 0.3
# Scale abundances of specific elements (dex units):
# (Further enhance carbon by 10x and nitrogen by 2x)
e_scale =
C 1.0
N 0.3
# Verbosity level (<0:errors, 0:warnings, 1:headlines, 2:details, 3:debug):
verb = 2
Abundance-profile Examples
Note
Before running this example, download the configuration files shown above, e.g., with these shell commands:
tutorial_path=https://raw.githubusercontent.com/pcubillos/pyratbay/master/examples/tutorial
wget $tutorial_path/atmosphere_tea.cfg
wget $tutorial_path/atmosphere_uniform.cfg
The following Python script creates and plots the abundance Aprofiles for the configuration files shown above:
import matplotlib.pyplot as plt
plt.ion()
import pyratbay as pb
import pyratbay.plots as pp
# Generate a uniform and a thermochemical-equilibrium atmospheric model:
pressure, temp, vmr_tea, species_tea, radius = pb.run("atmosphere_tea.cfg")
pressure, temp, vmr_uni, species_uni, radius = pb.run("atmosphere_uniform.cfg")
# Plot the results:
plt.figure(12, (6,5))
plt.clf()
ax = plt.subplot(211)
ax1 = pp.abundance(
vmr_tea, pressure, species_tea,
colors='default', xlim=[1e-11, 10.0], legend_fs=0, ax=ax,
)
ax = plt.subplot(212)
ax2 = pp.abundance(
vmr_uni, pressure, species_uni,
colors='default', xlim=[1e-11, 10.0], legend_fs=8, ax=ax,
)
plt.tight_layout()
And the results should look like this:
Altitude Profiles
If the user sets the radmodel
key, the code will to compute the
atmospheric altitude profile (radius profile). The currently
available models solve for the hydrostatic-equilibrium equation,
combined with the ideal gas law with a pressure-dependent gravity
(radmodel=hydro_m
, recommended):
or a constant surface gravity (radmodel=hydro_g
):
where \(M_{\rm p}\) is the mass of the planet, \(T(p)\) is the atmospheric temperature, \(\mu(p)\) is the atmospheric mean molecular mass, \(k_{\rm B}\) is the Boltzmann constant, and \(G\) is the gravitational constant. Note that \(T(p)\) and \(\mu(p)\) are computed from the models of the Temperature Profiles and Abundance Profiles, respectively.
To obtain the particular solution of these differential equations,
the user needs to supply a pair of radius–pressure reference values
to define the boundary condition \(r(p_0) = R_0\). The
rplanet
and refpressure
keys set \(R_0\) and \(p_0\),
respectively.
The mplanet
and gplanet
keys set the
planetary mass (\(M_p\)) and surface gravity (\(g\)),
respectively.
Note
Note that the user needs only to define two out of the three \(\{R_0, M_p, g\}\) variables, since they are related through the equation: \(g(R_0) = G M_p / R_0^2\).
Note that the selection of the \(\{p_0,R_0\}\) pair is arbitrary. A good practice is to choose values close to the transit radius of the planet. Although the pressure at the transit radius is a priori unknown for a give particular case [Griffith2014], its value lies at around 0.1 bar.
Here is an example of a hydrostatic-equilibrium atmosphere configuration file (atmosphere_hydro_m.cfg):
[pyrat]
# Pyrat Bay run mode, select from: [tli atmosphere spectrum opacity mcmc]
runmode = atmosphere
# Pressure at the top and bottom of the atmosphere, and number of layers:
ptop = 1.0e-6 bar
pbottom = 100.0 bar
nlayers = 100
# Temperature-profile model, select from [isothermal guillot madhu]
tmodel = guillot
# log(kappa') log(g1) log(g2) alpha T_irr T_int
tpars = -6.0 -0.25 0.0 0.0 950.0 100.0
# Output atmospheric model:
atmfile = WASP-00c.atm
# Chemistry model, select from [uniform tea]
chemistry = uniform
# Atmospheric composition and abundances:
species = H2 He Na H2O CH4 CO CO2 NH3 HCN N2
uniform = 0.85 0.149 5e-6 2e-3 1e-4 5e-4 1e-6 1e-5 1e-9 2e-4
# Altitude/radius profile model, select from [hydro_m, hydro_g]:
radmodel = hydro_m
rplanet = 2.87 rearth
mplanet = 2.86 mearth
refpressure = 0.1 bar
# Verbosity level (<0:errors, 0:warnings, 1:headlines, 2:details, 3:debug):
verb = 2
Altitude-profile Examples
Note
Before running this example, download the configuration files shown above, e.g., with these shell commands:
tutorial_path=https://raw.githubusercontent.com/pcubillos/pyratbay/master/examples/tutorial
wget $tutorial_path/atmosphere_hydro_m.cfg
wget $tutorial_path/atmosphere_hydro_g.cfg
The following Python script creates and plots the profiles for the configuration file shown above:
import matplotlib.pyplot as plt
plt.ion()
import pyratbay as pb
import pyratbay.constants as pc
# Kepler-11c mass and radius:
pressure, temp, q, species, radius = pb.run("atmosphere_hydro_m.cfg")
pressure, temp, q, species, radius_g = pb.run("atmosphere_hydro_g.cfg")
# Plot the results:
plt.figure(12, (6,5))
plt.clf()
ax = plt.subplot(111)
ax.semilogy(radius_g/pc.rearth, pressure/pc.bar, lw=2, c='navy', label='constant g')
ax.semilogy(radius/pc.rearth, pressure/pc.bar, lw=2, c='orange', label='g = g(p)')
ax.set_ylim(1e2, 1e-6)
ax.set_xlabel(r'Radius $(R_{\oplus})$', fontsize=12)
ax.set_ylabel('Pressure (bar)', fontsize=12)
ax.tick_params(labelsize=11)
ax.legend(loc='upper left', fontsize=12)
plt.tight_layout()
And the results should look like this: