[1]:
import logging
import numpy as np
from jaxmod.constants import EARTH_MASS
from atmodeller import EquilibriumModel, Planet, SpeciesNetwork, debug_logger
logger = debug_logger()
logger.setLevel(logging.INFO)
# For more output use DEBUG
# logger.setLevel(logging.DEBUG)
Atmodeller initialized with double precision (float64)
Atmosphere
This notebook is available at notebooks/atmosphere.ipynb and is easiest to obtain by downloading the source code.
Atmodeller allows you to simultaneuosly solve for surface conditions and an atmospheric P-T profile under equilibrium assumptions. In this example, we follow the model presented in Hakim et al. (2025) for TOI-421b.
[2]:
earth_radius = 6371000 # metre
# Mass and radius of TOI-421b
planet_mass = 6.7 * EARTH_MASS # kg
MEB_radius = 1.65 * earth_radius # metre
# atm_radius = 2.64 * earth_radius # metre
# MEB radius = 1.65 Earth radii in metre for 6.7 Earth masses (for atm_radius = 2.64 Earth radii)
# M-R relation from Hakim et al. (2018) Icarus
# Temperature of TOI-421b
surface_temperature = 3000 # K
top_temperature = 920 # K
We set up arrays for pressure (in bar) and temperature (in K) profile. The “trick” is to use np.nan as a placeholder indicating that the surface pressure should be solved for based on the total volatile mass in the atmosphere, that is, by enforcing mechanical pressure balance at the surface. For all other entires where a pressure value is provided, that value is used directly. The temperature array aligns positionally with the pressure array: the first entry gives the surface temperature,
and the last entry corresponds to the temperature at the final pressure level.
[3]:
temperature = np.array([surface_temperature, 1500, 1000, top_temperature])
pressure = np.array([np.nan, 1, 1e-2, 1e-4])
We then create the model using the gaseous and condensed species we wish to include:
[4]:
species = SpeciesNetwork.create(
(
# Gas species
"H2O_g",
"H2_g",
"O2_g",
"OSi_g",
"H4Si_g",
"CO2_g",
"CO_g",
"CH4_g",
"N2_g",
"NH3_g",
"He_g",
# Condensates
"O2Si_cd",
"O2Si_cd",
"O2Si_cd",
"C_cd",
"CSi_cd",
"Si_cd",
"N4Si3_cd",
)
)
model = EquilibriumModel(species)
[19:32:34 - atmodeller.classes - INFO ] - species_network = ('H2O_g: IdealGas, NoSolubility', 'H2_g: IdealGas, NoSolubility', 'O2_g: IdealGas, NoSolubility', 'OSi_g: IdealGas, NoSolubility', 'H4Si_g: IdealGas, NoSolubility', 'CO2_g: IdealGas, NoSolubility', 'CO_g: IdealGas, NoSolubility', 'CH4_g: IdealGas, NoSolubility', 'N2_g: IdealGas, NoSolubility', 'H3N_g: IdealGas, NoSolubility', 'He_g: IdealGas, NoSolubility', 'O2Si_cd: CondensateActivity, NoSolubility', 'O2Si_cd: CondensateActivity, NoSolubility', 'O2Si_cd: CondensateActivity, NoSolubility', 'C_cd: CondensateActivity, NoSolubility', 'CSi_cd: CondensateActivity, NoSolubility', 'Si_cd: CondensateActivity, NoSolubility', 'N4Si3_cd: CondensateActivity, NoSolubility')
[19:32:34 - atmodeller.classes - INFO ] - Thermodynamic data requires temperatures between 200 K and 4000 K
[19:32:34 - atmodeller.classes - INFO ] - reactions = {0: '1.0 H2_g + 1.0 CO2_g = 1.0 H2O_g + 1.0 CO_g',
1: '4.0 H2_g + 1.0 CO2_g = 2.0 H2O_g + 1.0 CH4_g',
2: '1.0 H2O_g + 1.0 H4Si_g = 3.0 H2_g + 1.0 OSi_g',
3: '1.5 H2_g + 0.5 N2_g = 1.0 H3N_g',
4: '2.0 H2O_g = 2.0 H2_g + 1.0 O2_g',
5: '2.0 H2O_g + 1.0 H4Si_g = 4.0 H2_g + 1.0 O2Si_cd',
6: '2.0 H2O_g + 1.0 H4Si_g = 4.0 H2_g + 1.0 O2Si_cd',
7: '2.0 H2O_g + 1.0 H4Si_g = 4.0 H2_g + 1.0 O2Si_cd',
8: '2.0 H2_g + 1.0 CO2_g = 2.0 H2O_g + 1.0 C_cd',
9: '1.0 H4Si_g + 1.0 CO2_g = 2.0 H2O_g + 1.0 CSi_cd',
10: '1.0 H4Si_g = 2.0 H2_g + 1.0 Si_cd',
11: '3.0 H4Si_g + 2.0 N2_g = 6.0 H2_g + 1.0 N4Si3_cd'}
We set the mass constraints based on a previous Atmodeller calculation that accounted for real-gas behaviour and dissolution into a magma ocean. From that calculation, we extract elemental abundances in the gas phase only (i.e. the atmosphere) and use them as abundance constraints. The values below correspond to a model with solar metallicity and a fully molten mantle.
[5]:
mass_constraints = {
"H": 5.73802837470845e22,
"He": 7.13997e21,
"C": 1.18962697880417e21,
"N": 1.99427e17,
"O": 1.94361960955633e22,
"Si": 3.5648e23,
}
Since we are using the gas masses from the previous calculation directly (which already accounted for gas solubility), we set mantle_melt_fraction=0. The specification of the planet’s mass and surface radius is used only when computing the mechanical pressure balance—that is, for determining the first pressure entry. These parameters play no role when the pressure profile is prescribed directly.
[6]:
planet = Planet(
temperature=temperature,
planet_mass=planet_mass,
mantle_melt_fraction=0,
surface_radius=MEB_radius,
pressure=pressure,
)
Finally, we solve the model:
[7]:
model.solve(state=planet, mass_constraints=mass_constraints)
[19:33:00 - atmodeller.classes - INFO ] - Solve (robust) complete: 4 (100.00%) successful model(s)
[19:33:01 - atmodeller.classes - INFO ] - Multistart summary: 3 (75.00%) models(s) required 1 attempt(s)
[19:33:01 - atmodeller.classes - INFO ] - Multistart summary: 1 (25.00%) models(s) required 2 attempt(s)
[19:33:01 - atmodeller.classes - INFO ] - Solver steps (max) = 253
The output will contain four model evaluations, and you can see that the gas pressure satisfies the prescribed constraints, with the first entry enforcing the mechanical pressure balance at the surface. This workflow therefore provides a convenient way to evaluate the atmspheric chemical signature at different temperatures and pressures while remaining consistent with the planetary surface conditions.
[8]:
# Show the system
model.output.to_dataframes()["state"]
[19:33:01 - atmodeller.output_core - INFO ] - Computing to_dataframes output
[19:33:01 - atmodeller.output_core - INFO ] - Computing asdict output
[8]:
| planet_mass | core_mass_fraction | mantle_melt_fraction | surface_radius | temperature | pressure | mantle_mass | mantle_melt_mass | mantle_solid_mass | surface_area | surface_gravity | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 4.001374e+25 | 0.295335 | 0.0 | 10512150.0 | 3000.0 | 15449.75169 | 2.819629e+25 | 0.0 | 2.819629e+25 | 1.388651e+15 | 24.167502 |
| 1 | 4.001374e+25 | 0.295335 | 0.0 | 10512150.0 | 1500.0 | 1.00000 | 2.819629e+25 | 0.0 | 2.819629e+25 | 1.388651e+15 | 24.167502 |
| 2 | 4.001374e+25 | 0.295335 | 0.0 | 10512150.0 | 1000.0 | 0.01000 | 2.819629e+25 | 0.0 | 2.819629e+25 | 1.388651e+15 | 24.167502 |
| 3 | 4.001374e+25 | 0.295335 | 0.0 | 10512150.0 | 920.0 | 0.00010 | 2.819629e+25 | 0.0 | 2.819629e+25 | 1.388651e+15 | 24.167502 |
If you’re interested in what happens “behind the scenes,” this works because the gas volume is not required during the solution itself; instead, it is computed deterministically during post-processing. This allows the total number of moles to remain fixed—as in the example above—while still accommodating the specification of different pressures, since the gas volume can adjust to maintain thermodynamic consistency with the imposed conditions.