[1]:

import logging

import numpy as np

from atmodeller import (
    ChemicalSpecies,
    EquilibriumModel,
    Planet,
    SpeciesNetwork,
    debug_logger,
    earth_oceans_to_hydrogen_mass,
)
from atmodeller.eos import get_eos_models
from atmodeller.solubility import get_solubility_models
from atmodeller.thermodata import IronWustiteBuffer

logger = debug_logger()
logger.setLevel(logging.INFO)

# For more output use DEBUG
# logger.setLevel(logging.DEBUG)

Atmodeller initialized with double precision (float64)

Basic Usage

This notebook is available at notebooks/basic_usage.ipynb and is easiest to obtain by downloading the source code.

Species and thermodynamic data

The species available in Atmodeller can be found in the thermodata subpackage, where the prefix of the dictionary key denotes the chemical formula in Hill notation and the suffix describes the states of aggregation in accordance with the JANAF convention.

[2]:

# Get all available species
available_species = SpeciesNetwork.available_species()
logger.info("Available species = %s", available_species)

# To create a gas species, for example CO2, where the state of aggregation defaults to 'g':
CO2_g = ChemicalSpecies.create_gas("CO2")

# The unique name that Atmodeller assigns combines the Hill notation and the state of aggregation
logger.info("Species name = %s", CO2_g.data.name)

# Compute the Gibbs energy relative to RT at 2000 K
temperature = 2000.0
gibbs = CO2_g.data.get_gibbs_over_RT(temperature)
logger.info("Gibbs/RT = %s", gibbs)

# Compute the composition
composition = CO2_g.data.composition
logger.info("Composition = %s", composition)

# Access more thermodynamic data
heat_capacity = CO2_g.data.thermo.cp(temperature)
logger.info("Heat capacity = %s", heat_capacity)

# Etc., other methods are available to compute other quantities

[19:33:17 - atmodeller                     - INFO     ] - Available species = ('Al_g', 'AlO_g', 'AlO2_g', 'Al2_g', 'Al2O_g', 'Al2O2_g', 'Al2O3_g', 'Ar_g', 'C_g', 'CN_g', 'CO_g', 'COS_g', 'CO2_g', 'C2H_g', 'C2N_g', 'C3_g', 'C4N2_g', 'CHN_g', 'CH2_g', 'CH3_g', 'CH4_g', 'C2H2_g', 'C2H3_g', 'C2H3N_g', 'C2H4_g', 'C2H5_g', 'C2H6_g', 'Ca_g', 'CaO_g', 'Ca2_g', 'Cl2_g', 'ClH_g', 'Cr_g', 'CrO_g', 'CrO2_g', 'CrO3_g', 'Fe_g', 'FeO_g', 'H_g', 'HO_g', 'HS_g', 'H2_g', 'H2O_g', 'H2O4S_g', 'H2S_g', 'H3N_g', 'H4Si_g', 'He_g', 'K_g', 'KO_g', 'K2_g', 'K2O_g', 'K2O2_g', 'Kr_g', 'Mg_g', 'MgO_g', 'N_g', 'NH_g', 'NH2_g', 'NO_g', 'N2_g', 'Na_g', 'NaO_g', 'Na2_g', 'Na2O_g', 'Na2O2_g', 'Ne_g', 'O_g', 'OS_g', 'OSi_g', 'OTi_g', 'O2_g', 'O2S_g', 'O2Si_g', 'O2Ti_g', 'O3_g', 'O3S_g', 'S2_g', 'Si_g', 'Ti_g', 'Xe_g', 'Al2O3_cd', 'C_cd', 'CSi_cd', 'CaO_cd', 'ClH4N_cd', 'Fe_cd', 'FeO_cd', 'Fe3O4_cd', 'H2O_cd', 'H2O4S_cd', 'MgO_cd', 'MgO3Si_cd', 'Mg2O4Si_cd', 'N4Si3_cd', 'O2Si_cd', 'S_cd', 'Si_cd')
[19:33:17 - atmodeller                     - INFO     ] - Species name = CO2_g
[19:33:18 - atmodeller                     - INFO     ] - Gibbs/RT = [-55.36417507969985]
[19:33:18 - atmodeller                     - INFO     ] - Composition = ImmutableMap({'C': (1, 12.01074, 0.27291212929920894), 'O': (2, 31.99881, 0.7270878707007911)})
[19:33:18 - atmodeller                     - INFO     ] - Heat capacity = [60.334253148740245]

Solubility

Solubility laws are available in the solubility subpackage.

[3]:

solubility_models = get_solubility_models()
logger.info("Solubility models = %s", solubility_models.keys())

CO2_basalt = solubility_models["CO2_basalt_dixon95"]
# Compute the concentration at fCO2=0.5 bar, 1300 K, and 1 bar
# Note that fugacity is the first argument and others are keyword only
concentration = CO2_basalt.concentration(0.5, temperature=1300, pressure=1)
logger.info("Concentration (ppmw) = %s", concentration)

[19:33:18 - atmodeller                     - INFO     ] - Solubility models = dict_keys(['Ar_basalt_jambon86', 'CH4_basalt_ardia13', 'CO2_basalt_dixon95', 'CO_basalt_armstrong15', 'CO_basalt_yoshioka19', 'CO_rhyolite_yoshioka19', 'Cl2_ano_dio_for_thomas21', 'Cl2_basalt_thomas21', 'H2O_ano_dio_newcombe17', 'H2O_basalt_dixon95', 'H2O_basalt_mitchell17', 'H2O_lunar_glass_newcombe17', 'H2O_peridotite_sossi23', 'H2_andesite_hirschmann12', 'H2_basalt_hirschmann12', 'H2_chachan18', 'H2_kite19', 'H2_silicic_melts_gaillard03', 'He_basalt_jambon86', 'Kr_basalt_jambon86', 'N2_basalt_bernadou21', 'N2_basalt_dasgupta22', 'N2_basalt_libourel03', 'Ne_basalt_jambon86', 'S2_andesite_boulliung23', 'S2_basalt_boulliung23', 'S2_sulfate_andesite_boulliung23', 'S2_sulfate_basalt_boulliung23', 'S2_sulfate_trachybasalt_boulliung23', 'S2_sulfide_andesite_boulliung23', 'S2_sulfide_basalt_boulliung23', 'S2_sulfide_trachybasalt_boulliung23', 'S2_trachybasalt_boulliung23', 'Xe_basalt_jambon86', 'NO_SOLUBILITY'])
[19:33:18 - atmodeller                     - INFO     ] - Concentration (ppmw) = 0.22841535272000957

[4]:

N2_basalt = solubility_models["N2_basalt_libourel03"]
# Compute the concentration at fCO2=0.5 bar, 1300 K, and 1 bar
# Note that fugacity is the first argument and others are keyword only
concentration = N2_basalt.concentration(0.20, temperature=1698.15, pressure=1, fO2=10**-16.2)
logger.info("Concentration (ppmw) = %s", concentration)

[19:33:18 - atmodeller                     - INFO     ] - Concentration (ppmw) = 377.1406984833259

[5]:

N2_basalt_dasgupta = solubility_models["N2_basalt_dasgupta22"]
# Compute the concentration at fCO2=0.5 bar, 1300 K, and 1 bar
# Note that fugacity is the first argument and others are keyword only
concentration = N2_basalt_dasgupta.concentration(
    1550, temperature=1773.15, pressure=1708.7, fO2=1.8e-13
)
logger.info("Concentration (ppmw) = %s", concentration)

[19:33:19 - atmodeller                     - INFO     ] - Concentration (ppmw) = 1002.1074478063852

Real gas EOS

Real gas equations of state are available in the eos subpackage.

[6]:

# Get all available EOS models
eos_models = get_eos_models()
logger.info("EOS models = %s", eos_models.keys())

# Get a CH4 model
CH4_eos_model = eos_models["CH4_beattie_holley58"]
# Compute the fugacity at 800 K and 100 bar
fugacity = CH4_eos_model.fugacity(800, 100)
logger.info("Fugacity = %s bar", fugacity)
# Compute the compressibility factor at the same conditions
compressibility = CH4_eos_model.compressibility_factor(800, 100)
logger.info("Compressibility factor = %s", compressibility)
# Etc., other methods are available to compute other quantities

[19:33:19 - atmodeller                     - INFO     ] - EOS models = dict_keys(['H2_chabrier21', 'H2_3000K_chabrier21', 'H2_4000K_chabrier21', 'H2_He_Y0275_chabrier21', 'H2_He_Y0292_chabrier21', 'H2_He_Y0297_chabrier21', 'He_chabrier21', 'CH4_beattie_holley58', 'CO2_beattie_holley58', 'H2_beattie_holley58', 'He_beattie_holley58', 'N2_beattie_holley58', 'NH3_beattie_holley58', 'O2_beattie_holley58', 'CH4_cork_cs_holland91', 'CO_cork_cs_holland91', 'CO2_cork_holland91', 'CO2_cork_holland98', 'CO2_cork_cs_holland91', 'H2_cork_cs_holland91', 'H2O_cork_holland91', 'H2O_cork_holland98', 'H2S_cork_cs_holland11', 'N2_cork_cs_holland91', 'S2_cork_cs_holland11', 'OSi_rk49_connolly16', 'H4Si_rk49_reid87', 'CHN_rk49_reid87', 'H3N_rk49_reid87', 'Ar_cs_saxena87', 'CH4_cs_shi92', 'CO_cs_shi92', 'CO2_cs_shi92', 'COS_cs_shi92', 'H2_shi92', 'H2S_shi92', 'N2_cs_saxena87', 'O2_cs_shi92', 'S2_cs_shi92', 'SO2_shi92', 'H2_vdw_lide05', 'He_vdw_lide05', 'N2_vdw_lide05', 'H4Si_vdw_lide05', 'H2O_vdw_lide05', 'CH4_vdw_lide05', 'H3N_vdw_lide05', 'CHN_vdw_lide05', 'H4Si_wang18', 'CH4_zhang09', 'H2O_zhang09', 'CO2_zhang09', 'H2_zhang09', 'CO_zhang09', 'O2_zhang09', 'C2H6_zhang09'])
[19:33:19 - atmodeller                     - INFO     ] - Fugacity = 103.31104964174332 bar
[19:33:20 - atmodeller                     - INFO     ] - Compressibility factor = 1.0336552056372945

We can also use broadcasting to perform multiple evaluations at once, for example to compute a grid of fugacities:

[7]:

# Define the temperature (K) and pressure (bar) grid
temperature = np.array([1000, 1600])
pressure = np.array([1, 10, 100])

temperature_broadcasted = temperature[:, None]
pressure_broadcasted = pressure[None, :]

# Get a CH4 model
CH4_eos_model = eos_models["CH4_cork_cs_holland91"]
# Compute the fugacity
fugacity = CH4_eos_model.fugacity(temperature_broadcasted, pressure_broadcasted)
logger.info("Fugacity = %s bar", fugacity)
# Compute the compressibility factor at the same conditions
compressibility = CH4_eos_model.compressibility_factor(
    temperature_broadcasted, pressure_broadcasted
)
logger.info("Compressibility factor = %s", compressibility)
# Etc., other methods are available to compute other quantities

[19:33:20 - atmodeller                     - INFO     ] - Fugacity = [[  1.                 10.036513019492777 104.04907255121003 ]
 [  1.                 10.027930271415434 103.07696646185612 ]] bar
[19:33:21 - atmodeller                     - INFO     ] - Compressibility factor = [[1.000406600334726 1.004038918021219 1.039882083041007]
 [1.000310800338152 1.003092414833786 1.030301664157559]]

Model with mass constraints

A common scenario is to calculate how volatiles partition between a magma ocean and an atmosphere when the total elemental abundances are constrained. Planet() defaults to a molten Earth, but the planetary parameters can be changed using input arguments.

[8]:

solubility_models = get_solubility_models()

H2_g = ChemicalSpecies.create_gas("H2")
H2O_g = ChemicalSpecies.create_gas("H2O", solubility=solubility_models["H2O_peridotite_sossi23"])
O2_g = ChemicalSpecies.create_gas("O2")

# This is one way of defining the species collection, and this approach is preferred if you want to
# specify species with solubility laws, as defined for H2O_g above.
species = SpeciesNetwork((H2_g, H2O_g, O2_g))

# Planet has input arguments that you can change. See the class documentation.
planet = Planet()
model = EquilibriumModel(species)

oceans = 1
h_kg = earth_oceans_to_hydrogen_mass(oceans)
o_kg = 6.25774e20
mass_constraints = {"H": h_kg, "O": o_kg}

model.solve(state=planet, mass_constraints=mass_constraints)
output = model.output

# Quick look at the solution
solution = output.quick_look()
logger.info("solution = %s", solution)

# Get complete solution as a dictionary
# solution_asdict = output.asdict()
# logger.info(solution_asdict)

# Get the complete solution as dataframes
# solution_dataframes = output.to_dataframes()

# Write the complete solution to Excel
# output.to_excel("example_mass_constraints")

[19:33:21 - atmodeller.classes             - INFO     ] - species_network = ('H2_g: IdealGas, NoSolubility', 'H2O_g: IdealGas, SolubilityPowerLaw', 'O2_g: IdealGas, NoSolubility')
[19:33:21 - atmodeller.classes             - INFO     ] - Thermodynamic data requires temperatures between 200 K and 6000 K
[19:33:21 - atmodeller.classes             - INFO     ] - reactions = {0: '2.0 H2O_g = 2.0 H2_g + 1.0 O2_g'}
[19:33:28 - atmodeller.classes             - INFO     ] - Solve (robust) complete: 1 (100.00%) successful model(s)
[19:33:28 - atmodeller.classes             - INFO     ] - Multistart summary: 1 (100.00%) models(s) required 1 attempt(s)
[19:33:28 - atmodeller.classes             - INFO     ] - Solver steps (max) = 35
[19:33:30 - atmodeller                     - INFO     ] - solution = {'H2_g': array(15.168753777766486), 'H2_g_activity': array(15.168753777766488), 'H2O_g': array(0.066407208350573), 'H2O_g_activity': array(0.066407208350573), 'O2_g': array(1.59110467117032e-12), 'O2_g_activity': array(1.591104671170319e-12)}

Model with fO2 constraint

Another common scenario is to calculate how volatiles partition between a magma ocean and an atmosphere when fO2 is fixed relative to a buffer and the total elemental abundances are constrained.

[9]:

H2O_g = ChemicalSpecies.create_gas("H2O", solubility=solubility_models["H2O_peridotite_sossi23"])
H2_g = ChemicalSpecies.create_gas("H2")
O2_g = ChemicalSpecies.create_gas("O2")

species = SpeciesNetwork((H2O_g, H2_g, O2_g))

planet = Planet()
model = EquilibriumModel(species)

oceans = 1
h_kg = earth_oceans_to_hydrogen_mass(oceans)
mass_constraints = {"H": h_kg}

# Use the Iron Wustite buffer.  The "-1" argument is the log10 shift relative to the buffer.
fugacity_constraints = {"O2_g": IronWustiteBuffer(-1)}

model.solve(
    state=planet, fugacity_constraints=fugacity_constraints, mass_constraints=mass_constraints
)
output = model.output

# Quick look at the solution
solution = output.quick_look()
logger.info("solution = %s", solution)

[19:33:30 - atmodeller.classes             - INFO     ] - species_network = ('H2O_g: IdealGas, SolubilityPowerLaw', 'H2_g: IdealGas, NoSolubility', 'O2_g: IdealGas, NoSolubility')
[19:33:30 - atmodeller.classes             - INFO     ] - Thermodynamic data requires temperatures between 200 K and 6000 K
[19:33:30 - atmodeller.classes             - INFO     ] - reactions = {0: '2.0 H2O_g = 2.0 H2_g + 1.0 O2_g'}
[19:33:37 - atmodeller.classes             - INFO     ] - Solve (robust) complete: 1 (100.00%) successful model(s)
[19:33:37 - atmodeller.classes             - INFO     ] - Multistart summary: 1 (100.00%) models(s) required 1 attempt(s)
[19:33:37 - atmodeller.classes             - INFO     ] - Solver steps (max) = 10
[19:33:37 - atmodeller                     - INFO     ] - solution = {'H2O_g': array(0.252820417697318), 'H2O_g_activity': array(0.252820417697317), 'H2_g': array(0.774830277331649), 'H2_g_activity': array(0.774830277331647), 'O2_g': array(8.838513516896055e-09), 'O2_g_activity': array(8.838513516896029e-09)}

Defining a network of species

In the above example, a species network was created by first creating the species and then aggregating the species into a SpeciesNetwork. This approach allows for complete generality since each species can also have its own solubility law assigned. However, if you want to create a simple gas network that assumes ideality and no solubility, you can also use the following formulation, where the state of aggregation must be specified after the formula (and separated by an underscore) for all species:

[10]:

species_no_solubility = SpeciesNetwork.create(("H2_g", "H2O_g", "O2_g"))
logger.info("species_no_solubility = %s", species_no_solubility)

[19:33:37 - atmodeller                     - INFO     ] - species_no_solubility = ('H2_g: IdealGas, NoSolubility', 'H2O_g: IdealGas, NoSolubility', 'O2_g: IdealGas, NoSolubility')

Batch calculation

For a batch calculation you can provide arrays to the planet or constraints. All arrays must have the same size because for a batch calculation the array values are aligned by position. Single values will automatically be broadcasted to the maximum array size.

[11]:

solubility_models = get_solubility_models()

H2_g = ChemicalSpecies.create_gas("H2")
H2O_g = ChemicalSpecies.create_gas("H2O", solubility=solubility_models["H2O_peridotite_sossi23"])
O2_g = ChemicalSpecies.create_gas("O2")

species = SpeciesNetwork((H2_g, H2O_g, O2_g))

# Batch temperature and radius, where the entries correspond by position. You could also choose
# to leave one or both as scalars.
surface_temperature = np.array([2000, 2000, 1500, 1500])
earth_radius = 6371000  # m
surface_radius = earth_radius * np.array([1.5, 3, 1.5, 3])

planet = Planet(temperature=surface_temperature, surface_radius=surface_radius)
model = EquilibriumModel(species)

oceans = 1
h_kg = earth_oceans_to_hydrogen_mass(oceans)
o_kg = 6.25774e20
scale_factor = 5
mass_constraints = {
    # We can also batch constraints, as long as we also have a total of 4 entries
    "H": np.array([h_kg, h_kg, h_kg * scale_factor, h_kg * scale_factor]),
    "O": np.array([o_kg, o_kg * scale_factor, o_kg, o_kg * scale_factor]),
}

# Initial solution guess number of moles
initial_log_number_moles = 50

model.solve(
    state=planet,
    initial_log_number_moles=initial_log_number_moles,
    mass_constraints=mass_constraints,
)
output = model.output

# Quick look at the solution
solution = output.quick_look()
logger.info("Quick look = %s", solution)

# Get complete solution as a dictionary
# solution_asdict = output.asdict()
# logger.info(solution_asdict)

# Write the complete solution to Excel
# output.to_excel("example_batch")

[19:33:37 - atmodeller.classes             - INFO     ] - species_network = ('H2_g: IdealGas, NoSolubility', 'H2O_g: IdealGas, SolubilityPowerLaw', 'O2_g: IdealGas, NoSolubility')
[19:33:37 - atmodeller.classes             - INFO     ] - Thermodynamic data requires temperatures between 200 K and 6000 K
[19:33:37 - atmodeller.classes             - INFO     ] - reactions = {0: '2.0 H2O_g = 2.0 H2_g + 1.0 O2_g'}
[19:33:44 - atmodeller.classes             - INFO     ] - Solve (robust) complete: 4 (100.00%) successful model(s)
[19:33:45 - atmodeller.classes             - INFO     ] - Multistart summary: 4 (100.00%) models(s) required 1 attempt(s)
[19:33:45 - atmodeller.classes             - INFO     ] - Solver steps (max) = 35
[19:33:46 - atmodeller                     - INFO     ] - Quick look = {'H2_g': array([3.332489632554158e+00, 3.261151703408906e-05,
       2.697322698724040e+01, 2.614487702886700e+00]), 'H2_g_activity': array([3.332489632554153e+00, 3.261151703408915e-05,
       2.697322698724046e+01, 2.614487702886697e+00]), 'H2O_g': array([0.064492106255406, 0.237741955198083, 0.06418764819667 ,
       0.816541825998788]), 'H2O_g_activity': array([0.064492106255406, 0.237741955198084, 0.06418764819667 ,
       0.816541825998787]), 'O2_g': array([3.109163888478601e-11, 4.412022955811442e+00,
       2.020962748149759e-17, 3.481006703548366e-13]), 'O2_g_activity': array([3.109163888478602e-11, 4.412022955811457e+00,
       2.020962748149771e-17, 3.481006703548367e-13])}

Monte Carlo

Exploring atmospheric compositions in a Monte Carlo model can be achieved with a batch calculation over a range of parameters. Note that in this case the same initial solution is used for all cases.

[12]:

solubility_models = get_solubility_models()

H2_g = ChemicalSpecies.create_gas("H2")
H2O_g = ChemicalSpecies.create_gas("H2O", solubility=solubility_models["H2O_peridotite_sossi23"])
O2_g = ChemicalSpecies.create_gas("O2")

species = SpeciesNetwork((H2_g, H2O_g, O2_g))
planet = Planet()
model = EquilibriumModel(species)

number_of_realisations = 1000
log10_number_oceans = np.random.uniform(0, 3, number_of_realisations)
number_oceans = 10**log10_number_oceans
fO2_min = -3
fO2_max = 3
fO2_log10_shifts = np.random.uniform(fO2_min, fO2_max, number_of_realisations)

oceans = 1
h_kg = earth_oceans_to_hydrogen_mass(number_oceans)
mass_constraints = {"H": h_kg}
fugacity_constraints = {O2_g.name: IronWustiteBuffer(fO2_log10_shifts)}

# Initial solution guess number of moles
initial_log_number_moles = 50 * np.ones(len(species))

model.solve(
    state=planet,
    initial_log_number_moles=initial_log_number_moles,
    mass_constraints=mass_constraints,
    fugacity_constraints=fugacity_constraints,
)
output = model.output

# Quick look at the solution
# solution = output.quick_look()
# logger.info("solution = %s", solution)

# Get complete solution as a dictionary
# solution_asdict = output.asdict()
# logger.info(solution_asdict)

# Write the complete solution to Excel
# output.to_excel("example_monte_carlo")

[19:33:46 - atmodeller.classes             - INFO     ] - species_network = ('H2_g: IdealGas, NoSolubility', 'H2O_g: IdealGas, SolubilityPowerLaw', 'O2_g: IdealGas, NoSolubility')
[19:33:46 - atmodeller.classes             - INFO     ] - Thermodynamic data requires temperatures between 200 K and 6000 K
[19:33:46 - atmodeller.classes             - INFO     ] - reactions = {0: '2.0 H2O_g = 2.0 H2_g + 1.0 O2_g'}
[19:33:55 - atmodeller.classes             - INFO     ] - Solve (robust) complete: 1000 (100.00%) successful model(s)
[19:33:55 - atmodeller.classes             - INFO     ] - Multistart summary: 1000 (100.00%) models(s) required 1 attempt(s)
[19:33:56 - atmodeller.classes             - INFO     ] - Solver steps (max) = 30