Quickstart demonstration for how to use the easychem package
You can find an executable version of this notebook on Google Colabhere.
We begin by loading the relevant packages: numpy, matplotlib, and easychem itself.
[1]:
import matplotlib.pyplot as plt
import numpy as np
import easychem as ec
plt.rcParams['figure.figsize'] = (10, 6)
Chemistry is calculated by creating objects from easychem’s ExoAtmos class.
[2]:
exo = ec.ExoAtmos()
By default, easychem assumes solar atmospheric abundances, following Asplund et al. (2009). You can switch to any other composition, however, and then update the abundances using exo.updateAtomAbund(), see an example further below.
The reactant species which are added by default are
'H', 'H2', 'He', 'O', 'C', 'N', 'Mg', 'Si', 'Fe', 'S', 'Al', 'Ca', 'Na', 'Ni', 'P', 'K', 'Ti', 'CO', 'OH', 'SH', 'N2', 'O2', 'SiO', 'TiO', 'SiS', 'H2O', 'C2', 'CH', 'CN', 'CS', 'SiC', 'NH', 'SiH', 'NO', 'SN', 'SiN', 'SO', 'S2', 'C2H', 'HCN', 'C2H2,acetylene', 'CH4', 'AlH', 'AlOH', 'Al2O', 'CaOH', 'MgH', 'MgOH', 'PH3', 'CO2', 'TiO2', 'Si2C', 'SiO2', 'FeO', 'NH2', 'NH3', 'CH2', 'CH3', 'H2S', 'VO', 'VO2', 'NaCl', 'KCl', 'e-', 'H+', 'H-', 'Na+', 'K+', 'PH2', 'P2', 'PS', 'PO', 'P4O6', 'PH', 'V', 'VO(c)','VO(L)','MgSiO3(c)','SiC(c)', 'Fe(c)','Al2O3(c)','Na2S(c)','KCl(c)','Fe(L)','SiC(L)', 'MgSiO3(L)','H2O(L)','H2O(c)','TiO(c)','TiO(L)', 'TiO2(c)','TiO2(L)','H3PO4(c)','H3PO4(L)'.
More species (molecular, ions, condensates, etc.) can be added if they have a corresponding entry in easychem’s thermo_easy_chem_simp_own.inp file. This file follows the typical CEA input data file structure (NASA polynomials), so you can add more species if you have the corresponding thermodynamic data.
This is the how to construct the path of thermo_easy_chem_simp_own.inp, so you can modify it. You can also create a local copy elsewhere, and then store its absolute path in a hidden file called .easychem in your home directory.
[3]:
thermo_file_path = ec.__file__.replace('__init__.py','thermo_easy_chem_simp_own.inp')
This is what the first few lines of the default file look like, copied from the CEA input data file. This is for the species H, so atomic hydrogen.
[4]:
# show the first few lines:
f = open(thermo_file_path, 'r')
lines = f.readlines()
f.close()
for i in range(14):
print(lines[i][:-1])
H D0(H2):Herzberg,1970. Moore,1972. Gordon,1999.
4 g 6/97 H 1.00 0.00 0.00 0.00 0.00 0 1.0079400 217998.828
60.000 300.0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 6197.428
0.00000000D+00 0.00000000D+00 0.25000000D+01 0.00000000D+00 0.00000000D+00
0.00000000D+00 0.00000000D+00 0.00000000D+00 0.25473708D+05 -0.44668285D+00
300.000 1000.0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 6197.428
0.000000000D+00 0.000000000D+00 2.500000000D+00 0.000000000D+00 0.000000000D+00
0.000000000D+00 0.000000000D+00 2.547370801D+04-4.466828530D-01
1000.000 6000.0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 6197.428
6.078774250D+01-1.819354417D-01 2.500211817D+00-1.226512864D-07 3.732876330D-11
-5.687744560D-15 3.410210197D-19 2.547486398D+04-4.481917770D-01
6000.000 20000.0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 6197.428
2.173757694D+08-1.312035403D+05 3.399174200D+01-3.813999680D-03 2.432854837D-07
-7.694275540D-12 9.644105630D-17 1.067638086D+06-2.742301051D+02
Calculating equilibrium abundances at single pressure-temperature point:
This is how you would calculate the atmospheric composition at a single point in the atmosphere, in this case at 1 bar and 1000 K:
[5]:
%time exo.solve(1, 1000) # pressure = 1 bar; temperature = 1000 K
CPU times: user 13.7 ms, sys: 1.01 ms, total: 14.7 ms
Wall time: 14.6 ms
You can either get the result in the form of a numpy.ndarray, or a dictionary. Here first the arrays:
[6]:
print(exo.reacMass) # mass fractions
print(exo.reacMols) # number fractions (aka volume mixing ratios, VMRs)
[9.06021531e-10 7.37035941e-01 2.49582309e-01 7.14991129e-23
1.71120160e-31 1.32997187e-23 1.32552621e-04 1.69366974e-38
3.59250117e-13 1.75854880e-14 9.25083891e-26 1.02065921e-05
2.49918731e-05 7.13609915e-05 8.04922751e-14 1.28500554e-06
7.75264035e-30 1.47737548e-04 3.29252917e-14 1.44426842e-09
6.73058823e-04 6.02815405e-26 5.03657374e-26 1.75359098e-23
7.51799992e-27 5.00502878e-03 1.99056606e-37 1.58180608e-27
1.36952560e-21 8.43062896e-15 0.00000000e+00 1.11254028e-19
2.96316322e-36 4.90846769e-19 3.12107218e-18 1.55689884e-40
6.04280958e-16 9.55662579e-13 1.79368225e-26 3.36098596e-09
4.47393261e-13 3.07857565e-03 6.95996723e-24 2.84105592e-17
9.34852915e-32 7.69568806e-05 3.94895382e-10 2.14643882e-07
6.30982696e-06 2.52917564e-07 5.60594704e-22 0.00000000e+00
8.42650905e-32 5.61152874e-20 7.99163549e-14 2.51643916e-05
1.05201637e-19 1.25782953e-10 3.29138971e-04 1.11000627e-19
1.28934651e-16 1.08724834e-05 3.40209843e-06 2.44241088e-18
0.00000000e+00 1.41372724e-21 3.19478213e-16 1.73529394e-13
3.81727878e-08 4.80369239e-08 1.62565617e-10 3.64336846e-11
4.15672352e-32 3.43966910e-11 7.20735796e-23 4.17410293e-07
0.00000000e+00 2.37988650e-03 0.00000000e+00 1.29376650e-03
1.05261688e-04 0.00000000e+00 0.00000000e+00 0.00000000e+00
0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00
0.00000000e+00 0.00000000e+00 5.21473752e-06 0.00000000e+00
0.00000000e+00 0.00000000e+00]
[2.09757275e-09 8.53173161e-01 1.45501532e-01 1.04278054e-23
3.32451899e-32 2.21565731e-24 1.27259043e-05 1.40715691e-39
1.50109640e-14 1.27973232e-15 8.00038686e-27 5.94252036e-07
2.53665534e-06 2.83705377e-06 6.06395085e-15 7.66907676e-08
3.77928277e-31 1.23075664e-05 4.51741738e-15 1.01899393e-10
5.60638811e-05 4.39588791e-27 2.66588430e-27 6.40695723e-25
2.91648021e-28 6.48281235e-04 1.93363385e-38 2.83520794e-28
1.22829369e-22 4.46330656e-16 0.00000000e+00 1.72901172e-20
2.37660527e-37 3.81708951e-20 1.58075853e-19 8.63087307e-42
2.93367056e-17 3.47727710e-14 1.67221750e-27 2.90196476e-10
4.00950831e-14 4.47794912e-04 5.80242540e-25 1.50707070e-18
3.11798877e-33 3.14571322e-06 3.64029108e-11 1.21236977e-08
4.33078303e-07 1.34099919e-08 1.63788763e-23 0.00000000e+00
3.27252253e-33 1.82257005e-21 1.16386162e-14 3.44792290e-06
1.75012529e-20 1.95223298e-11 2.25354111e-05 3.86927721e-21
3.62743806e-18 4.34104653e-07 1.06484798e-07 1.03890249e-14
0.00000000e+00 3.27120550e-21 3.24275595e-17 1.03566006e-14
2.70005610e-09 1.80945032e-09 6.01750859e-12 1.80987636e-12
4.41101423e-34 2.50963799e-12 3.30141654e-24 1.45501532e-08
0.00000000e+00 5.53181153e-05 0.00000000e+00 5.40589451e-05
2.40896978e-06 0.00000000e+00 0.00000000e+00 0.00000000e+00
0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00
0.00000000e+00 0.00000000e+00 1.52358808e-07 0.00000000e+00
0.00000000e+00 0.00000000e+00]
And here as dictionaries:
[7]:
print(exo.result_mass()) # mass fractions
print(exo.result_mol()) # number fractions (VMRs)
{'H': 9.060215308072591e-10, 'H2': 0.737035941363717, 'He': 0.24958230902520975, 'O': 7.149911286224232e-23, 'C': 1.7112016027774863e-31, 'N': 1.329971874554925e-23, 'Mg': 0.00013255262131084644, 'Si': 1.693669736736222e-38, 'Fe': 3.5925011694060837e-13, 'S': 1.758548798129851e-14, 'Al': 9.250838911088532e-26, 'Ca': 1.0206592126054387e-05, 'Na': 2.499187313428508e-05, 'Ni': 7.136099147838817e-05, 'P': 8.049227508563026e-14, 'K': 1.2850055396102771e-06, 'Ti': 7.752640352282779e-30, 'CO': 0.00014773754762171913, 'OH': 3.292529169406808e-14, 'SH': 1.4442684247517603e-09, 'N2': 0.0006730588233870282, 'O2': 6.028154046711741e-26, 'SiO': 5.0365737409058804e-26, 'TiO': 1.753590975315518e-23, 'SiS': 7.517999920493922e-27, 'H2O': 0.005005028779706047, 'C2': 1.9905660636964318e-37, 'CH': 1.5818060797453722e-27, 'CN': 1.3695256041318704e-21, 'CS': 8.430628961114065e-15, 'SiC': 0.0, 'NH': 1.1125402847591055e-19, 'SiH': 2.963163218353625e-36, 'NO': 4.908467691460504e-19, 'SN': 3.1210721777243077e-18, 'SiN': 1.5568988407771779e-40, 'SO': 6.04280957693754e-16, 'S2': 9.55662578551513e-13, 'C2H': 1.7936822458495732e-26, 'HCN': 3.3609859618708623e-09, 'C2H2,acetylene': 4.4739326061170604e-13, 'CH4': 0.0030785756451886057, 'AlH': 6.959967231353941e-24, 'AlOH': 2.841055916473596e-17, 'Al2O': 9.348529153791191e-32, 'CaOH': 7.695688062248463e-05, 'MgH': 3.94895382464083e-10, 'MgOH': 2.1464388230468712e-07, 'PH3': 6.3098269561914736e-06, 'CO2': 2.529175638081966e-07, 'TiO2': 5.605947041447451e-22, 'Si2C': 0.0, 'SiO2': 8.426509045047038e-32, 'FeO': 5.611528735219297e-20, 'NH2': 7.991635489035828e-14, 'NH3': 2.51643916178121e-05, 'CH2': 1.0520163693196814e-19, 'CH3': 1.2578295314971457e-10, 'H2S': 0.0003291389708739827, 'VO': 1.1100062736219158e-19, 'VO2': 1.2893465080550783e-16, 'NaCl': 1.0872483390650977e-05, 'KCl': 3.4020984300377294e-06, 'e-': 2.442410880591701e-18, 'H+': 0.0, 'H-': 1.4137272407676952e-21, 'Na+': 3.194782125264179e-16, 'K+': 1.735293941070376e-13, 'PH2': 3.8172787845697383e-08, 'P2': 4.8036923852053564e-08, 'PS': 1.6256561708343488e-10, 'PO': 3.643368458483637e-11, 'P4O6': 4.156723519087943e-32, 'PH': 3.4396690988501726e-11, 'V': 7.207357959439425e-23, 'VO(c)': 4.174102930089093e-07, 'VO(L)': 0.0, 'MgSiO3(c)': 0.0023798865025368332, 'SiC(c)': 0.0, 'Fe(c)': 0.0012937665025954154, 'Al2O3(c)': 0.00010526168807752692, 'Na2S(c)': 0.0, 'KCl(c)': 0.0, 'Fe(L)': 0.0, 'SiC(L)': 0.0, 'MgSiO3(L)': 0.0, 'H2O(L)': 0.0, 'H2O(c)': 0.0, 'TiO(c)': 0.0, 'TiO(L)': 0.0, 'TiO2(c)': 5.214737521858054e-06, 'TiO2(L)': 0.0, 'H3PO4(c)': 0.0, 'H3PO4(L)': 0.0}
{'H': 2.097572751445111e-09, 'H2': 0.8531731613832432, 'He': 0.14550153221245757, 'O': 1.0427805363204739e-23, 'C': 3.3245189926217705e-32, 'N': 2.215657306831546e-24, 'Mg': 1.2725904348609038e-05, 'Si': 1.4071569076282964e-39, 'Fe': 1.5010963963258873e-14, 'S': 1.279732322873811e-15, 'Al': 8.000386860162457e-27, 'Ca': 5.942520360950701e-07, 'Na': 2.53665534086733e-06, 'Ni': 2.8370537687356277e-06, 'P': 6.063950848282071e-15, 'K': 7.669076761288391e-08, 'Ti': 3.7792827681465778e-31, 'CO': 1.2307566351336626e-05, 'OH': 4.517417377088386e-15, 'SH': 1.0189939331478205e-10, 'N2': 5.606388106755884e-05, 'O2': 4.395887905605103e-27, 'SiO': 2.665884295834736e-27, 'TiO': 6.4069572252552075e-25, 'SiS': 2.916480214291209e-28, 'H2O': 0.0006482812347429341, 'C2': 1.9336338494791793e-38, 'CH': 2.8352079403065337e-28, 'CN': 1.2282936920805474e-22, 'CS': 4.463306559567085e-16, 'SiC': 0.0, 'NH': 1.7290117246791445e-20, 'SiH': 2.376605268636456e-37, 'NO': 3.817089507329063e-20, 'SN': 1.5807585303872394e-19, 'SiN': 8.630873065122164e-42, 'SO': 2.933670562491673e-17, 'S2': 3.4772770958471805e-14, 'C2H': 1.6722175040431014e-27, 'HCN': 2.9019647603728134e-10, 'C2H2,acetylene': 4.009508311948472e-14, 'CH4': 0.00044779491248561057, 'AlH': 5.8024254004087625e-25, 'AlOH': 1.5070706987087033e-18, 'Al2O': 3.117988772592058e-33, 'CaOH': 3.145713219052891e-06, 'MgH': 3.640291080425702e-11, 'MgOH': 1.2123697732222612e-08, 'PH3': 4.3307830275452847e-07, 'CO2': 1.3409991878781424e-08, 'TiO2': 1.6378876336305845e-23, 'Si2C': 0.0, 'SiO2': 3.2725225276537845e-33, 'FeO': 1.8225700463999145e-21, 'NH2': 1.1638616175047716e-14, 'NH3': 3.4479229020391146e-06, 'CH2': 1.750125291164255e-20, 'CH3': 1.9522329817687155e-11, 'H2S': 2.2535411131173248e-05, 'VO': 3.869277213828102e-21, 'VO2': 3.627438064406813e-18, 'NaCl': 4.3410465269872033e-07, 'KCl': 1.0648479840880163e-07, 'e-': 1.0389024904381006e-14, 'H+': 0.0, 'H-': 3.2712054951772588e-21, 'Na+': 3.242755947513347e-17, 'K+': 1.035660061611145e-14, 'PH2': 2.700056103664487e-09, 'P2': 1.809450315770265e-09, 'PS': 6.01750859118075e-12, 'PO': 1.8098763564695e-12, 'P4O6': 4.4110142301718365e-34, 'PH': 2.5096379854763146e-12, 'V': 3.301416540231318e-24, 'VO(c)': 1.455015321928534e-08, 'VO(L)': 0.0, 'MgSiO3(c)': 5.531811528558553e-05, 'SiC(c)': 0.0, 'Fe(c)': 5.405894509574202e-05, 'Al2O3(c)': 2.4089697795947237e-06, 'Na2S(c)': 0.0, 'KCl(c)': 0.0, 'Fe(L)': 0.0, 'SiC(L)': 0.0, 'MgSiO3(L)': 0.0, 'H2O(L)': 0.0, 'H2O(c)': 0.0, 'TiO(c)': 0.0, 'TiO(L)': 0.0, 'TiO2(c)': 1.523588081823083e-07, 'TiO2(L)': 0.0, 'H3PO4(c)': 0.0, 'H3PO4(L)': 0.0}
Calculating equilibrium abundances for a whole atmospheric P-T structure:
Let’s assume an isothermal (1000 K) atmospheric temperature structure (100 layers) going from 1 µbar to 100 bar (equidistantly log-spaced). We can then calculate the atmospheric structure as follows:
[8]:
press = np.logspace(-6,2,100)
temp = np.ones(100) * 1000.
%time exo.solve(press, temp)
CPU times: user 478 ms, sys: 48 ms, total: 526 ms
Wall time: 529 ms
We can plot the results as follows:
[9]:
mass_fractions = exo.result_mass()
for spec in exo.reactants:
plt.loglog(mass_fractions[spec], press, label = spec)
plt.ylim([1e2,1e-6])
plt.xlim([1e-10, 1])
plt.xlabel('Mass fractions')
plt.ylabel('Pressure (bar)')
[9]:
Text(0, 0.5, 'Pressure (bar)')
Changing metallicity
In the following example, we are changing the considered atmosphere’s metallicity. The elemental composition exo.atomAbunds will then automatically be updated with the new value of exo.metallicity taken into account. exo.metallicity is defined with respect to solar, so if you change it all elemental abundances (except H and He) are scaled by 10**exo.metallicity (a value of 0 is thus solar).
[10]:
exo.metallicity = 0.5
%time exo.solve(press, temp)
mass_fractions = exo.result_mass()
CPU times: user 472 ms, sys: 44.8 ms, total: 517 ms
Wall time: 516 ms
Sanity check: let’s compare to the chemical equilibrium interpolation of petitRADTRANS:
[11]:
from petitRADTRANS.poor_mans_nonequ_chem import interpol_abundances
COs = 0.55 * np.ones_like(press) # solar
FeHs = 0.5 * np.ones_like(press) # metals are scaled by 10^0.5 * solar
pRT_mass_fractions = interpol_abundances(COs,
FeHs,
temp,
press)
Now we plot and compare the two resulting compositions:
[12]:
for i_spec, spec in enumerate(pRT_mass_fractions.keys()):
if 'nabla' in spec or 'MMW' in spec:
continue
plt.loglog(pRT_mass_fractions[spec],
press,
label = spec,
linestyle = '-',
color = 'C'+str(i_spec))
try:
plt.loglog(mass_fractions[spec],
press,
linestyle = ':',
color = 'C'+str(i_spec))
except KeyError:
pass
plt.ylim([1e2,1e-6])
plt.xlim([1e-10, 1])
plt.xlabel('Mass fractions')
plt.ylabel('Pressure (bar)')
plt.title('Solid = pRT chemistry interpolation, dashed = easyCHEM')
plt.legend()
[12]:
<matplotlib.legend.Legend at 0x7f85635d4af0>
Note that there are slight differences, because pRT interpolates in an equilibrium abundance table, while easychem calculates the abundances at every P-T point from scratch.
Changing C/O
Next we change the carbon-to-oxygen number ratio. In easyCHEM this is done by changing the oxygen content of the atmosphere. You can also change carbon instead, however, see below.
[13]:
exo.co = 1.2
%time exo.solve(press, temp)
mass_fractions = exo.result_mass()
CPU times: user 500 ms, sys: 46 ms, total: 546 ms
Wall time: 545 ms
Again we compare to pRT here:
[14]:
COs = 1.2 * np.ones_like(press)
FeHs = 0.5 * np.ones_like(press)
pRT_mass_fractions = interpol_abundances(COs,
FeHs,
temp,
press)
for i_spec, spec in enumerate(pRT_mass_fractions.keys()):
if 'nabla' in spec or 'MMW' in spec:
continue
plt.loglog(pRT_mass_fractions[spec],
press,
label = spec,
linestyle = '-',
color = 'C'+str(i_spec))
try:
plt.loglog(mass_fractions[spec],
press,
linestyle = ':',
color = 'C'+str(i_spec))
except KeyError:
pass
plt.ylim([1e2,1e-6])
plt.xlim([1e-10, 1])
plt.xlabel('Mass fractions')
plt.ylabel('Pressure (bar)')
plt.title('Solid = pRT chemistry interpolation, dashed = easyCHEM')
plt.legend()
[14]:
<matplotlib.legend.Legend at 0x7f856542fee0>
Again note the differences due to pRT’s interpolation.
Varying atmospheric abundances freely
This is an example for how to modify atmospheric elemental abundances freely. Here we chose to allow the user to scale [C/H], [O/H] and [Fe/H] independently, where [X/H] is defined as \(\rm log_{10}(X/H)-log_{10}(X/H)_\odot\). “Fe” in this case stands for all elements except for C, H, O, He. Similar setups could also be used to implement atmospheric condensate rainout.
[15]:
# Start with a clean ExoAtmos object (so solar composition)
exo = ec.ExoAtmos()
# Copy the solar abundance setup
atom_abundances_solar = exo._atomAbunds.copy()
# Get atom names
atom_names = exo.atoms
# Implement function to return the modified abundances.
def update_atom_abundances(C_H, O_H, Fe_H):
modif_abundances = atom_abundances_solar.copy()
for i_spec, name in enumerate(atom_names):
if name == 'C':
modif_abundances[i_spec] = modif_abundances[i_spec]*10**C_H
elif name == 'O':
modif_abundances[i_spec] = modif_abundances[i_spec]*10**O_H
elif name not in ['H', 'He', 'C', 'O']:
modif_abundances[i_spec] = modif_abundances[i_spec]*10**Fe_H
return modif_abundances
Let’s test how this setup behaves now! We start by changing [C/H] from -1 to 1.
[16]:
C_Hs = np.linspace(-1, 1, 10) # from 0.1 to 10 x solar
H2Os = np.zeros(10) # to store H2O mass fractions
COs = np.zeros(10) # to store CO mass fractions
CH4s = np.zeros(10) # to store CH4 mass fractions
CO2s = np.zeros(10) # to store CO2 mass fractions
for i, C_H in enumerate(C_Hs):
update_abunds = update_atom_abundances(C_H,0,0) # make update abundance list
exo.updateAtomAbunds(update_abunds) # update abundances
exo.solve(0.001, 1200) # calculate composition at 1 mbar, 1200 K.
mass_fractions = exo.result_mass()
H2Os[i] = mass_fractions['H2O']
COs[i] = mass_fractions['CO']
CH4s[i] = mass_fractions['CH4']
CO2s[i] = mass_fractions['CO2']
plt.plot(C_Hs, H2Os, label = 'H2O')
plt.plot(C_Hs, COs, label = 'CO')
plt.plot(C_Hs, CH4s, label = 'CH4')
plt.plot(C_Hs, CO2s, label = 'CO2')
plt.yscale('log')
plt.xlabel('[C/H]')
plt.ylabel('Mass fraction')
plt.legend()
plt.show()
Now we change [O/H] from -1 to 1, [Fe/H] from -1 to 1, then all three ([C/H], [O/H], [Fe/H]) at the same time.
[17]:
O_Hs = np.linspace(-1, 1, 20) # from 0.1 to 10 x solar
H2Os = np.zeros(20) # to store H2O mass fractions
COs = np.zeros(20) # to store CO mass fractions
CH4s = np.zeros(20) # to store CH4 mass fractions
CO2s = np.zeros(20) # to store CO2 mass fractions
for i, O_H in enumerate(O_Hs):
update_abunds = update_atom_abundances(0,O_H,0)
exo.updateAtomAbunds(update_abunds)
exo.solve(0.001, 1200)
mass_fractions = exo.result_mass()
H2Os[i] = mass_fractions['H2O']
COs[i] = mass_fractions['CO']
CH4s[i] = mass_fractions['CH4']
CO2s[i] = mass_fractions['CO2']
plt.plot(O_Hs, H2Os, label = 'H2O')
plt.plot(O_Hs, COs, label = 'CO')
plt.plot(O_Hs, CH4s, label = 'CH4')
plt.plot(O_Hs, CO2s, label = 'CO2')
plt.yscale('log')
plt.xlabel('[O/H]')
plt.ylabel('Mass fraction')
plt.legend()
plt.show()
Fe_Hs = np.linspace(-1, 1, 20) # from 0.1 to 10 x solar
H2Os = np.zeros(20)
COs = np.zeros(20)
CH4s = np.zeros(20)
CO2s = np.zeros(20)
for i, Fe_H in enumerate(Fe_Hs):
update_abunds = update_atom_abundances(0,0,Fe_H)
exo.updateAtomAbunds(update_abunds)
exo.solve(0.001, 1200)
mass_fractions = exo.result_mass()
H2Os[i] = mass_fractions['H2O']
COs[i] = mass_fractions['CO']
CH4s[i] = mass_fractions['CH4']
CO2s[i] = mass_fractions['CO2']
plt.plot(Fe_Hs, H2Os, label = 'H2O')
plt.plot(Fe_Hs, COs, label = 'CO')
plt.plot(Fe_Hs, CH4s, label = 'CH4')
plt.plot(Fe_Hs, CO2s, label = 'CO2')
plt.yscale('log')
plt.xlabel('[Fe/H] (increasing all metals except C, O)')
plt.ylabel('Mass fraction')
plt.legend()
plt.show()
H2Os = np.zeros(20)
COs = np.zeros(20)
CH4s = np.zeros(20)
CO2s = np.zeros(20)
for i, Fe_H in enumerate(Fe_Hs):
update_abunds = update_atom_abundances(Fe_H,Fe_H,Fe_H)
exo.updateAtomAbunds(update_abunds)
exo.solve(0.001, 1200)
mass_fractions = exo.result_mass()
H2Os[i] = mass_fractions['H2O']
COs[i] = mass_fractions['CO']
CH4s[i] = mass_fractions['CH4']
CO2s[i] = mass_fractions['CO2']
plt.plot(Fe_Hs, H2Os, label = 'H2O')
plt.plot(Fe_Hs, COs, label = 'CO')
plt.plot(Fe_Hs, CH4s, label = 'CH4')
plt.plot(Fe_Hs, CO2s, label = 'CO2')
plt.yscale('log')
plt.xlabel('[Fe/H] (increasing all metals including C, O)')
plt.ylabel('Mass fraction')
plt.legend()
plt.show()
Ancillary outputs
Adiabatic temperature gradient \(\nabla_{\rm ad}\)
easyCHEM also outputs the atmospheric’s adiabatic temperature gradient, that is, \(\nabla_{\rm ad} = (d {\rm log T}/d {\rm log P})_{\rm ad}\). This is the gradient an atmosphere follows if it is convectively unstable. Note that easyCHEM is a CEA clone, so \(\nabla_{\rm ad}\) includes changes in chemical composition as the temperature is changed at a given pressure. \(\nabla_{\rm ad}\) therefore corresponds to the moist adiabatic lapse rate. Here an example for how to access \(\nabla_{\rm ad}\) in a solar composition atmosphere:
[18]:
exo = ec.ExoAtmos()
press = np.logspace(-6,2,100)
temp = np.ones(100) * 1000.
%time exo.solve(press, temp)
CPU times: user 477 ms, sys: 45.5 ms, total: 522 ms
Wall time: 526 ms
[19]:
nabla_ad = (exo.gamma2-1.)/exo.gamma2
plt.semilogy(nabla_ad, press)
plt.xlabel(r'$\nabla_{\rm ad}$')
plt.xlim([0.2, 0.4])
plt.ylim([100, 1e-6])
plt.ylabel('Pressure (bar)')
[19]:
Text(0, 0.5, 'Pressure (bar)')
Mean molecular weight (aka mean molar mass)
Can be accessed like this:
[20]:
plt.semilogy(exo.mmw, press)
plt.xlabel(r'Mean molar mass')
plt.ylim([100, 1e-6])
plt.ylabel('Pressure (bar)')
[20]:
Text(0, 0.5, 'Pressure (bar)')
Condensates: a word on stability
Condensates are tricky to handle in any equilibrium chemistry code. This is because condensates cannot occur together if one can write down hypothetical reaction equations that transform them into each other while conserving the atom number. This only occurs if the formation of the involved condensates is thermodynamically favored (i.e., the temperature is cold enough and the pressure is high enough). An example would be
For the Gibbs minimizer that is the backbone of easyCHEM the actual problem is that the matrix to be inverted becomes rank-deficient. easyCHEM’s standard selection is pretty stable against this problem, but things can still go bad in a retrieval when many different temperature - pressure combinations are tried. In this case easyCHEM will become very slow and will give non-sensical solutions. See an example here when we add species that trigger the described problem:
[34]:
exo = ec.ExoAtmos()
%time exo.solve(1, 800)
print()
print(exo.result_mass())
print()
reactant_list = exo.reactants.copy()
reactant_list.append('Mg2SiO4(c)')
reactant_list.append('Mg2SiO4(L)') # L = liquid
reactant_list.append('MgAl2O4(c)')
reactant_list.append('FeO(c)')
reactant_list.append('Fe2SiO4(c)')
exo.updateReactants(reactant_list)
%time exo.solve(1, 800)
print()
print(exo.result_mass())
CPU times: user 17.2 ms, sys: 884 µs, total: 18 ms
Wall time: 18.2 ms
{'H': 1.1562719305509502e-12, 'H2': 0.7369883299709381, 'He': 0.24958230894433125, 'O': 2.1911616919487047e-29, 'C': 5.024901328789189e-42, 'N': 7.592855912936709e-30, 'Mg': 0.0001323378393782409, 'Si': 0.0, 'Fe': 1.5646190348006107e-18, 'S': 1.7064638276791544e-18, 'Al': 4.159314806444928e-36, 'Ca': 5.185032613251398e-07, 'Na': 1.3697489677399486e-08, 'Ni': 7.136099145526319e-05, 'P': 5.454813210872906e-18, 'K': 9.554110876025068e-10, 'Ti': 8.295820646107013e-42, 'CO': 1.855585229495953e-07, 'OH': 6.59223202905385e-18, 'SH': 8.960064009518709e-12, 'N2': 0.0005939228960091637, 'O2': 2.1970243045110852e-32, 'SiO': 2.745268655737504e-40, 'TiO': 4.0968179035921114e-33, 'SiS': 5.860437258622254e-41, 'H2O': 0.005095621456313255, 'C2': 0.0, 'CH': 1.8529902314158584e-36, 'CN': 1.7689418087324437e-28, 'CS': 6.333792360164624e-20, 'SiC': 0.0, 'NH': 2.281467331064815e-24, 'SiH': 0.0, 'NO': 1.7827722830099404e-23, 'SN': 4.0793997788504437e-22, 'SiN': 0.0, 'SO': 1.3956794625276365e-19, 'S2': 3.86500924422296e-15, 'C2H': 3.8474057525607355e-36, 'HCN': 4.411023338101566e-12, 'C2H2,acetylene': 2.7028003940500914e-18, 'CH4': 0.0031631774651549164, 'AlH': 3.7926518279348396e-33, 'AlOH': 5.340820435212119e-24, 'Al2O': 4.1649824347375293e-45, 'CaOH': 9.075615831473653e-05, 'MgH': 3.451933972151145e-11, 'MgOH': 5.803064605426587e-07, 'PH3': 6.4001887614241945e-06, 'CO2': 9.506929033569126e-10, 'TiO2': 3.3640168271827565e-30, 'Si2C': 0.0, 'SiO2': 2.133581511848544e-46, 'FeO': 2.2877498941823566e-26, 'NH2': 2.8849305984407644e-16, 'NH3': 0.00012138592287730278, 'CH2': 6.832509858803456e-26, 'CH3': 1.3908252771067763e-13, 'H2S': 0.00031006661151307397, 'VO': 4.354175860294414e-27, 'VO2': 2.9187184480188256e-22, 'NaCl': 8.95313300336674e-06, 'KCl': 5.850855282771632e-06, 'e-': 9.272680238001707e-23, 'H+': 0.0, 'H-': 1.0671415331114668e-27, 'Na+': 8.844833751496268e-22, 'K+': 6.6072131839109835e-18, 'PH2': 1.1581055358965007e-09, 'P2': 6.015938984252903e-10, 'PS': 7.388487764174484e-13, 'PO': 4.9301528394670126e-14, 'P4O6': 9.925323678139046e-31, 'PH': 3.358538146617307e-14, 'V': 6.386665499365991e-32, 'VO(c)': 4.1741029297781614e-07, 'VO(L)': 0.0, 'MgSiO3(c)': 0.002379886501765602, 'SiC(c)': 0.0, 'Fe(c)': 0.0012937665025354033, 'Al2O3(c)': 0.0001052616880434485, 'Na2S(c)': 4.367894495983109e-05, 'KCl(c)': 0.0, 'Fe(L)': 0.0, 'SiC(L)': 0.0, 'MgSiO3(L)': 0.0, 'H2O(L)': 0.0, 'H2O(c)': 0.0, 'TiO(c)': 0.0, 'TiO(L)': 0.0, 'TiO2(c)': 5.2147375201681555e-06, 'TiO2(L)': 0.0, 'H3PO4(c)': 0.0, 'H3PO4(L)': 0.0}
SLOW ! Press, Temp 1.0000000000000000 800.00000000000000
CPU times: user 32.4 s, sys: 893 ms, total: 33.2 s
Wall time: 33.3 s
{'H': 1.9684336459490452e-14, 'H2': 0.012546263651379166, 'He': 0.004249065446925054, 'O': 4.652509731338093e-31, 'C': 8.556937155115721e-44, 'N': 1.2927245443721003e-31, 'Mg': 1.2029841252362701e-05, 'Si': 1.1296859002355934e-07, 'Fe': 2.2025994102299282e-05, 'S': 1.3372059813014334e-20, 'Al': 9.531580524087138e-19, 'Ca': 7.090381371542072e-09, 'Na': 4.282463349545213e-07, 'Ni': 1.2148999034963767e-06, 'P': 9.288498932810395e-20, 'K': 1.7170971690720256e-08, 'Ti': 9.081698717069078e-44, 'CO': 3.94070042095599e-09, 'OH': 4.5457753212419663e-08, 'SH': 6.860366895829447e-14, 'N2': 1.011161640185868e-05, 'O2': 5.817674706480911e-34, 'SiO': 0.00010821434893440605, 'TiO': 5.593145355851912e-35, 'SiS': 8.333854699237142e-06, 'H2O': 0.00010818141696994615, 'C2': 0.0, 'CH': 3.1552607582051887e-38, 'CN': 3.0122709426632844e-30, 'CS': 4.852574809416103e-22, 'SiC': 2.6384502911981266e-40, 'NH': 3.8840643151763735e-26, 'SiH': 6.703298794326528e-18, 'NO': 3.785283561785818e-25, 'SN': 3.1247459760398593e-24, 'SiN': 1.0534508135812108e-22, 'SO': 1.3332640372165846e-21, 'S2': 1.3319775899616303e-17, 'C2H': 6.552541288034581e-38, 'HCN': 7.510885884397768e-14, 'C2H2,acetylene': 4.6028526469232625e-20, 'CH4': 5.385172890209732e-05, 'AlH': 8.690751833930731e-16, 'AlOH': 1.526244314474471e-06, 'Al2O': 1.6021007771131828e-08, 'CaOH': 1.5476303511645028e-06, 'MgH': 3.137687291484777e-12, 'MgOH': 6.578188186201669e-08, 'PH3': 1.0896133138637656e-07, 'CO2': 2.5178767909242015e-11, 'TiO2': 5.727555602565779e-32, 'Si2C': 1.0078687104036557e-25, 'SiO2': 1.0488441203267311e-10, 'FeO': 4.016401885215563e-13, 'NH2': 4.9111008262455304e-18, 'NH3': 2.066251239540588e-06, 'CH2': 1.1633588969939559e-27, 'CH3': 2.3679762687697728e-15, 'H2S': 2.374796582362191e-06, 'VO': 1.0532282174706015e-13, 'VO2': 8.804616261755816e-09, 'NaCl': 1.7806666822790585e-07, 'KCl': 6.689272298115187e-08, 'e-': 5.129791545763748e-23, 'H+': 0.0, 'H-': 5.9032038270103885e-28, 'Na+': 8.510578226665732e-22, 'K+': 3.654598196619834e-18, 'PH2': 1.971770813361958e-11, 'P2': 1.0245273870451392e-11, 'PS': 5.660291840131858e-15, 'PO': 1.0469543462056104e-15, 'P4O6': 6.36039486045353e-32, 'PH': 5.718567601795519e-16, 'V': 1.2387651578961043e-18, 'VO(c)': 0.0, 'VO(L)': 0.0, 'MgSiO3(c)': 0.0, 'SiC(c)': 0.0, 'Fe(c)': 0.0, 'Al2O3(c)': 0.0, 'Na2S(c)': 0.0, 'KCl(c)': 0.0, 'Fe(L)': 0.0, 'SiC(L)': 0.0, 'MgSiO3(L)': 0.0, 'H2O(L)': 0.0, 'H2O(c)': 0.0, 'TiO(c)': 0.0, 'TiO(L)': 0.0, 'TiO2(c)': 0.9828721327112375, 'TiO2(L)': 0.0, 'H3PO4(c)': 0.0, 'H3PO4(L)': 0.0, 'Mg2SiO4(c)': 0.0, 'Mg2SiO4(L)': 0.0, 'MgAl2O4(c)': 0.0, 'FeO(c)': 0.0, 'Fe2SiO4(c)': 0.0}
EASY CHEM WARNING: One or more convergence criterianot satisfied! in cond Press, temp 1.0000000000000000 800.00000000000000
What to do is not entirely clear. The best thing likely is to check which species are expected to condense at a given pressure, temperature and composition, and which of these dominate the mass budget, when compared to other species. A useful starting point are papers such as Visser et al. (2010).