Chemistry Schemes#
The vertical chemistry routines calculate the volume mixing ratio (VMR) of each species, as defined as the fraction of the local total number density of the atmosphere.
Exo Skryer provides several vertical chemistry calculation functions in the vert_chem module.
Several schemes are included in Exo Skryer, from simple constant profiles to chemical equilibrium to quenching timescale approximation. In many schemes we assume an H2-He dominated atmosphere, forcing the total VMR to satisfy
For table-interpolated chemistry (fastchem_grid_jax), returned species are
the mapped active-opacity set, and VMRs are used directly from the CE table
without additional renormalization in vert_chem.
Background Gas: H2 and He Filling#
After all trace-species VMRs \(X_i\) are determined, the remaining fraction of the atmosphere is filled with H2 and He following the scheme of Welbanks et al. (2021), ensuring \(\sum_i X_i = 1\):
where the He to H2 ratio is fixed at the solar value from Asplund et al. (2021):
Atomic H and Free-Electron Fraction#
For ultra-hot Jupiter atmospheres where H2 dissociation is significant, the filling scheme can be extended to include atomic H. A H/H2 VMR ratio is defined as a retrievable quantity
with the background filler fraction
and the solar He/H ratio from Asplund et al. (2021):
Defining \(N_{\rm H}\) as the dimensionless abundance of hydrogen nuclei in the filler,
the H2, H, and He VMRs are then
Note that in certain circumstances \(X_{\rm H} > X_{\rm H_2}\), which should be taken into account when setting prior bounds for the H/H2 ratio.
This scheme is combined with a free-electron number density fraction \(f_{\rm e^-}\) as a separate retrieved parameter,
where \(n_{\rm e^-}\) [cm-3] is the electron number density and \(n_{\rm tot}\) [cm-3] the total background gas number density. Together, the H/H2 and \(f_{\rm e^-}\) parameters enable a consistent recovery of H- free–free opacity, an important continuum source in ultra-hot Jupiter atmospheres.
params:
- { name: log_10_H_over_H2, dist: uniform, low: -6, high: 0, transform: logit, init: -3 }
- { name: log_10_ne_over_ntot, dist: uniform, low: -6, high: 0, transform: logit, init: -4 }
Constant VMR#
The constant VMR profile assumes a constant value for each species as given by the sampled (or delta) parameter, log_10_f_x.
import numpy as np
import matplotlib.pyplot as plt
from exo_skryer.vert_Tp import Modified_Milne
from exo_skryer.vert_chem import constant_vmr
from exo_skryer.data_constants import bar
nlev = 100
p_bot = np.log10(100.0)
p_top = np.log10(1e-4)
p_lev = np.logspace(p_bot, p_top, nlev) * bar
params_tp = {
"T_int": 1200.0,
"T_ratio": 0.333,
"log_10_g": 4.5,
"log_10_k_ir": -2.0,
"log_10_p_t": 0.0,
"beta": 0.55
}
T_lev, T_lay = Modified_Milne(p_lev, params_tp)
p_lay = (p_lev[1:] - p_lev[:-1]) / np.log(p_lev[1:] / p_lev[:-1])
species_order = ("H2O", "CO", "CO2", "CH4")
chem_kernel = constant_vmr(species_order)
params = {
"log_10_f_H2O": -3.0,
"log_10_f_CO": -4.0,
"log_10_f_CO2": -8.0,
"log_10_f_CH4": -6.0,
}
vmr_lay = chem_kernel(p_lay, T_lay, params, nlev - 1)
fig, ax = plt.subplots(figsize=(10, 5))
for key in ("H2O", "CO", "CO2", "CH4"):
if key in vmr_lay:
ax.semilogy(vmr_lay[key], p_lay / bar, label=key)
ax.set_xlabel("VMR", fontsize=16)
ax.set_ylabel("pressure [bar]", fontsize=16)
ax.set_title("Constant VMR Chemistry", fontsize=14)
ax.legend(fontsize=10)
ax.set_xscale('log')
ax.invert_yaxis()
plt.tight_layout()
(Source code, png, hires.png, pdf, svg)
Example YAML Configuration:
physics:
vert_chem: constant_vmr
params:
- { name: log_10_f_H2O, dist: uniform, low: -9, high: -1, transform: logit, init: -3 }
- { name: log_10_f_CO, dist: uniform, low: -9, high: -1, transform: logit, init: -4 }
- { name: log_10_f_CO2, dist: uniform, low: -9, high: -1, transform: logit, init: -8 }
- { name: log_10_f_CH4, dist: uniform, low: -9, high: -1, transform: logit, init: -6 }
Constant VMR (CLR Parameterization)#
constant_vmr_clr parameterizes composition in centered-log-ratio (CLR)
space and applies a softmax transform so trace-species plus filler are always
non-negative and normalized.
Aliases: constant_clr, clr.
This scheme supports two parameter styles:
native CLR parameters:
clr_<species>standard abundance parameters:
log_10_f_<species>
When log_10_f_* is provided, Exo Skryer converts those values to CLR
coordinates internally before applying the softmax transform.
Example YAML Configuration:
physics:
vert_chem: constant_vmr_clr
params:
- { name: clr_H2O, dist: uniform, low: -12.0, high: 2.0, transform: logit, init: -4.0 }
- { name: clr_CO, dist: uniform, low: -12.0, high: 2.0, transform: logit, init: -5.0 }
- { name: clr_CH4, dist: uniform, low: -12.0, high: 2.0, transform: logit, init: -6.0 }
Alternative (log10 VMR inputs, internally converted to CLR):
physics:
vert_chem: constant_vmr_clr
params:
- { name: log_10_f_H2O, dist: uniform, low: -12.0, high: -1.0, transform: logit, init: -4.0 }
- { name: log_10_f_CO, dist: uniform, low: -12.0, high: -1.0, transform: logit, init: -5.0 }
- { name: log_10_f_CH4, dist: uniform, low: -12.0, high: -1.0, transform: logit, init: -6.0 }
Chemical Equilibrium with rate JAX#
We can also use the semi-analytical chemical equilibrium scheme, Reliable Analytic Thermochemical Equilibrium (rate), from Cubillos et al (2019). This was converted into JAX compatible python from the original python code found on GitHub.
import numpy as np
import matplotlib.pyplot as plt
from pathlib import Path
from exo_skryer.rate_jax import load_nasa9_cache
from exo_skryer.vert_Tp import Modified_Milne
from exo_skryer.vert_chem import CE_rate_jax
from exo_skryer.data_constants import bar
try:
root = Path(__file__).resolve().parent
except NameError:
root = Path.cwd().resolve()
for _ in range(5):
if (root / "NASA9").is_dir():
break
root = root.parent
load_nasa9_cache(str(root / "NASA9"))
nlev = 100
p_bot = np.log10(1000.0)
p_top = np.log10(1e-8)
p_lev = np.logspace(p_bot, p_top, nlev) * bar
params_tp = {
"T_int": 1200.0,
"T_ratio": 0.333,
"log_10_g": 4.5,
"log_10_k_ir": -2.0,
"log_10_p_t": 0.0,
"beta": 0.55
}
T_lev, T_lay = Modified_Milne(p_lev, params_tp)
p_lay = (p_lev[1:] - p_lev[:-1]) / np.log(p_lev[1:] / p_lev[:-1])
params = {"M_to_H": 0.0, "C_to_O": 0.55}
vmr_lay = CE_rate_jax(p_lay, T_lay, params, nlev - 1)
fig, ax = plt.subplots(figsize=(10, 5))
for key in ("H2O", "CO", "CH4", "NH3", "HCN", "CO2"):
if key in vmr_lay:
ax.semilogy(vmr_lay[key], p_lay / bar, label=key)
ax.set_xlabel("VMR", fontsize=16)
ax.set_ylabel("pressure [bar]", fontsize=16)
ax.set_title("CE Chemistry (RateJAX)", fontsize=14)
ax.legend(fontsize=10)
ax.set_xscale('log')
ax.invert_yaxis()
plt.tight_layout()
(Source code, png, hires.png, pdf, svg)
Example YAML Configuration:
physics:
vert_chem: CE_rate_jax
params:
- { name: M_to_H, dist: uniform, low: -1.0, high: 2.0, transform: logit, init: 0.0 }
- { name: C_to_O, dist: uniform, low: 0.1, high: 1.5, transform: logit, init: 0.55 }
Chemical Equilibrium with FastChem 5D Grid Interpolation#
This backend interpolates a precomputed FastChem grid over
(temperature, pressure, M_to_H, C_to_O) using JAX
RegularGridInterpolator.
Interpolation is performed in log-space coordinates
(log10(T), log10(P), M/H[dex], log10(C/O)) for both legacy linear NPZ
grids and log10 NPZ grids. VMR and MMW tables are interpolated in log10 space
and converted back to linear values.
The chemistry key is fastchem_grid_jax (aliases: ce_fastchem_grid,
fastchem_ce_grid). Legacy ce / fastchem_jax aliases are also
supported and route to this backend.
When physics.opac_special enables H- continuum:
bf requires mapped
H-in the CE gridff requires mapped
Hande-in the CE grid
Initialization raises if required species are missing.
The backend also provides interpolated mean molecular weight to the vert_mu
stage via an internal cache key (used automatically by vert_mu: dynamic and
vert_mu: auto).
Example YAML Configuration:
physics:
vert_chem: fastchem_grid_jax
fastchem_grid_jax:
grid_path: ../../FastChem/fastchem_grid_5d_log10.npz
solver:
mode: vmap
bounds:
mode: clip
species_map:
# Optional per-species mapping override
# CO: C1O1
Chemical Equilibrium with Atmodeller#
This backend solves equilibrium chemistry with atmodeller using an explicit
species network.
The chemistry key is atmodeller.
Example YAML Configuration:
physics:
vert_chem: atmodeller
atmodeller:
species_network:
- H2O_g
- CO_g
- CO2_g
- CH4_g
- NH3_g
- H2_g
- H_g
- He_g
solver:
multistart: 5
atol: 1.0e-5
rtol: 1.0e-5
max_steps: 128
params:
- { name: M_to_H, dist: uniform, low: -1.0, high: 2.0, transform: logit, init: 0.0 }
- { name: C_to_O, dist: uniform, low: 0.1, high: 1.5, transform: logit, init: 0.55 }
Chemical Equilibrium with Element Potentials JAX#
This backend solves thermochemical equilibrium from NASA-9 Gibbs free energies using the element-potentials formulation (Lagrange multipliers on element budgets).
The chemistry kernel key is easychem_jax (alias: easychem). Species are configured explicitly in a dedicated
YAML block, and all listed species must exist as <species>.txt files in
data.nasa9.
Solver mode can be configured:
scan: warm-starts layer-to-layer (typically more robust)vmap: cold-start parallel solves (often faster on GPU)
Example YAML Configuration:
physics:
vert_chem: easychem_jax
easychem_jax:
species: [H2O, CO, CO2, CH4, NH3, HCN, H2, He]
elements: [H, He, C, N, O]
e_ref: H
p0_bar: 1.0
solver:
mode: scan
max_steps: 64
tol: 1.0e-11
throw: false
prefer_chord: true
params:
- { name: M_to_H, dist: uniform, low: -1.0, high: 2.0, transform: logit, init: 0.0 }
- { name: C_to_O, dist: uniform, low: 0.1, high: 1.5, transform: logit, init: 0.55 }
Quenching Timescale Approximation#
import numpy as np
import matplotlib.pyplot as plt
from pathlib import Path
from exo_skryer.rate_jax import load_nasa9_cache
from exo_skryer.vert_Tp import Modified_Milne
from exo_skryer.vert_chem import quench_approx
from exo_skryer.data_constants import bar
try:
root = Path(__file__).resolve().parent
except NameError:
root = Path.cwd().resolve()
for _ in range(5):
if (root / "NASA9").is_dir():
break
root = root.parent
load_nasa9_cache(str(root / "NASA9"))
nlev = 100
p_bot = np.log10(1000.0)
p_top = np.log10(1e-8)
p_lev = np.logspace(p_bot, p_top, nlev) * bar
params_tp = {
"T_int": 1200.0,
"T_ratio": 0.333,
"log_10_g": 4.5,
"log_10_k_ir": -2.0,
"log_10_p_t": 0.0,
"beta": 0.55
}
T_lev, T_lay = Modified_Milne(p_lev, params_tp)
p_lay = (p_lev[1:] - p_lev[:-1]) / np.log(p_lev[1:] / p_lev[:-1])
params = {"M_to_H": 0.0, "C_to_O": 0.55, "log_10_Kzz": 8.0, "log_10_g": 4.5}
vmr_lay = quench_approx(p_lay, T_lay, params, nlev - 1)
fig, ax = plt.subplots(figsize=(10, 5))
for key in ("H2O", "CO", "CH4", "NH3", "HCN", "CO2"):
if key in vmr_lay:
ax.semilogy(vmr_lay[key], p_lay / bar, label=key)
ax.set_xlabel("VMR", fontsize=16)
ax.set_ylabel("pressure [bar]", fontsize=16)
ax.set_title("Quench Approx Chemistry", fontsize=14)
ax.legend(fontsize=10)
ax.set_xscale('log')
ax.invert_yaxis()
plt.tight_layout()
Example YAML Configuration:
physics:
vert_chem: quench_approx
params:
- { name: M_to_H, dist: uniform, low: -1.0, high: 2.0, transform: logit, init: 0.0 }
- { name: C_to_O, dist: uniform, low: 0.1, high: 1.5, transform: logit, init: 0.55 }
- { name: log_10_Kzz, dist: uniform, low: 5.0, high: 11.0, transform: logit, init: 8.0 }
- { name: log_10_g, dist: uniform, low: 2.0, high: 4.0, transform: logit, init: 3.0 }