Skip to content

Commit

Permalink
Merge pull request #298 from FormingWorlds/hn/volatiles
Browse files Browse the repository at this point in the history
Allow setting volatile abundances as both metallicity or concentration
  • Loading branch information
nichollsh authored Nov 28, 2024
2 parents 03c7329 + 6ab57b6 commit fa011c6
Show file tree
Hide file tree
Showing 10 changed files with 201 additions and 73 deletions.
12 changes: 9 additions & 3 deletions input/aragog.toml
Original file line number Diff line number Diff line change
Expand Up @@ -232,11 +232,17 @@ author = "Harrison Nicholls, Tim Lichtenberg"

# Set initial volatile inventory by planetary element abundances
[delivery.elements]
CH_ratio = 1.0 # C/H mass ratio in mantle/atmosphere system
H_oceans = 0.0 # Hydrogen inventory in units of equivalent Earth oceans
# H_oceans = 0.0 # Hydrogen inventory in units of equivalent Earth oceans
H_ppmw = 109.0 # Hydrogen inventory in ppmw relative to mantle mass

CH_ratio = 1.0 # C/H mass ratio in mantle/atmosphere system
# C_ppmw = 0.0 # Carbon inventory in ppmw relative to mantle mass

# NH_ratio = 0.0 # N/H mass ratio in mantle/atmosphere system
N_ppmw = 2.0 # Nitrogen inventory in ppmw relative to mantle mass
S_ppmw = 200.0 # Sulfur inventory in ppmw relative to mantle mass

SH_ratio = 2.0 # S/H mass ratio in mantle/atmosphere system
# S_ppmw = 1.0 # Sulfur inventory in ppmw relative to mantle mass

# Set initial volatile inventory by partial pressures in atmosphere
[delivery.volatiles]
Expand Down
11 changes: 9 additions & 2 deletions input/default.toml
Original file line number Diff line number Diff line change
Expand Up @@ -232,12 +232,19 @@ author = "Harrison Nicholls, Tim Lichtenberg"

# Set initial volatile inventory by planetary element abundances
[delivery.elements]
CH_ratio = 1.0 # C/H mass ratio in mantle/atmosphere system
H_oceans = 6.0 # Hydrogen inventory in units of equivalent Earth oceans
H_ppmw = 0.0 # Hydrogen inventory in ppmw relative to mantle mass
# H_ppmw = 0.0 # Hydrogen inventory in ppmw relative to mantle mass

CH_ratio = 1.0 # C/H mass ratio in mantle/atmosphere system
# C_ppmw = 0.0 # Carbon inventory in ppmw relative to mantle mass

# NH_ratio = 0.0 # N/H mass ratio in mantle/atmosphere system
N_ppmw = 2.0 # Nitrogen inventory in ppmw relative to mantle mass

# SH_ratio = 0.0 # S/H mass ratio in mantle/atmosphere system
S_ppmw = 200.0 # Sulfur inventory in ppmw relative to mantle mass


# Set initial volatile inventory by partial pressures in atmosphere
[delivery.volatiles]
H2O = 30.0 # partial pressure of H2O
Expand Down
14 changes: 10 additions & 4 deletions input/hd63433d.toml
Original file line number Diff line number Diff line change
Expand Up @@ -232,11 +232,17 @@ author = "Harrison Nicholls, Tim Lichtenberg"

# Set initial volatile inventory by planetary element abundances
[delivery.elements]
H_oceans = 8.0 # Hydrogen inventory in units of equivalent Earth oceans
# H_ppmw = 0.0 # Hydrogen inventory in ppmw relative to mantle mass

CH_ratio = 1.0 # C/H mass ratio in mantle/atmosphere system
H_oceans = 8.0 # Hydrogen inventory in units of equivalent Earth oceans, by mass
H_ppmw = 0.0 # Hydrogen inventory in ppmw relative to mantle mass
N_ppmw = 2.01 # Nitrogen inventory in ppmw relative to mantle mass, by mass
S_ppmw = 235.0 # Sulfur inventory in ppmw relative to mass of melt
# C_ppmw = 0.0 # Carbon inventory in ppmw relative to mantle mass

# NH_ratio = 0.0 # N/H mass ratio in mantle/atmosphere system
N_ppmw = 2.01 # Nitrogen inventory in ppmw relative to mantle mass

# SH_ratio = 0.0 # S/H mass ratio in mantle/atmosphere system
S_ppmw = 235.0 # Sulfur inventory in ppmw relative to mantle mass

# Set initial volatile inventory by partial pressures in atmosphere
[delivery.volatiles]
Expand Down
15 changes: 11 additions & 4 deletions input/k218b.toml
Original file line number Diff line number Diff line change
Expand Up @@ -232,11 +232,18 @@ author = "Harrison Nicholls, Tim Lichtenberg"

# Set initial volatile inventory by planetary element abundances
[delivery.elements]
H_oceans = 60.0 # Hydrogen inventory in units of equivalent Earth oceans
# H_ppmw = 0.0 # Hydrogen inventory in ppmw relative to mantle mass

CH_ratio = 0.1 # C/H mass ratio in mantle/atmosphere system
H_oceans = 60.0 # Hydrogen inventory in units of equivalent Earth oceans, by mass
H_ppmw = 0.0 # Hydrogen inventory in ppmw relative to mantle mass
N_ppmw = 2.01 # Nitrogen inventory in ppmw relative to mantle mass, by mass
S_ppmw = 235.0 # Sulfur inventory in ppmw relative to mass of melt
# C_ppmw = 0.0 # Carbon inventory in ppmw relative to mantle mass

# NH_ratio = 0.0 # N/H mass ratio in mantle/atmosphere system
N_ppmw = 2.01 # Nitrogen inventory in ppmw relative to mantle mass

# SH_ratio = 0.0 # S/H mass ratio in mantle/atmosphere system
S_ppmw = 235.0 # Sulfur inventory in ppmw relative to mantle mass


# Set initial volatile inventory by partial pressures in atmosphere
[delivery.volatiles]
Expand Down
32 changes: 19 additions & 13 deletions input/l9859d.toml
Original file line number Diff line number Diff line change
Expand Up @@ -35,7 +35,7 @@ author = "Harrison Nicholls, Tim Lichtenberg"
minimum = 3e2 # yr, minimum time-step
maximum = 1e7 # yr, maximum time-step
initial = 1e3 # yr, inital step size
starspec = 1e8 # yr, interval to re-calculate the stellar spectrum
starspec = 1e9 # yr, interval to re-calculate the stellar spectrum
starinst = 1e2 # yr, interval to re-calculate the instellation
method = "adaptive" # proportional | adaptive | maximum

Expand All @@ -44,7 +44,7 @@ author = "Harrison Nicholls, Tim Lichtenberg"

[params.dt.adaptive]
atol = 0.02 # Step size atol
rtol = 0.08 # Step size rtol
rtol = 0.09 # Step size rtol

# Termination criteria
# Set enabled=true/false in each section to enable/disable that termination criterion
Expand All @@ -53,7 +53,7 @@ author = "Harrison Nicholls, Tim Lichtenberg"
# required number of iterations
[params.stop.iters]
enabled = true
minimum = 5
minimum = 8
maximum = 9000

# required time constraints
Expand Down Expand Up @@ -133,14 +133,14 @@ author = "Harrison Nicholls, Tim Lichtenberg"
module = "agni" # Which atmosphere module to use

[atmos_clim.agni]
p_top = 1.0e-4 # bar, top of atmosphere grid pressure
p_top = 1.0e-5 # bar, top of atmosphere grid pressure
spectral_group = "Honeyside" # which gas opacities to include
spectral_bands = "256" # how many spectral bands?
num_levels = 42 # Number of atmospheric grid levels
chemistry = "eq" # "none" | "eq"
num_levels = 41 # Number of atmospheric grid levels
chemistry = "none" # "none" | "eq"
surf_material = "greybody" # surface material file for scattering
solution_atol = 1e-2 # solver absolute tolerance
solution_rtol = 5e-2 # solver relative tolerance
solution_rtol = 6e-2 # solver relative tolerance
overlap_method = "ee" # gas overlap method

[atmos_clim.janus]
Expand Down Expand Up @@ -181,9 +181,9 @@ author = "Harrison Nicholls, Tim Lichtenberg"
num_levels = 160 # Number of SPIDER grid levels
mixing_length = 2 # Mixing length parameterization
tolerance = 1.0e-9 # solver tolerance
tsurf_atol = 20.0 # tsurf_poststep_change
tsurf_atol = 25.0 # tsurf_poststep_change
tsurf_rtol = 0.01 # tsurf_poststep_change_frac
ini_entropy = 3400.0 # Surface entropy conditions [J K-1 kg-1]
ini_entropy = 3500.0 # Surface entropy conditions [J K-1 kg-1]
ini_dsdr = -4.698e-6 # Interior entropy gradient [J K-1 kg-1 m-1]

[interior.aragog]
Expand Down Expand Up @@ -231,11 +231,17 @@ author = "Harrison Nicholls, Tim Lichtenberg"

# Set initial volatile inventory by planetary element abundances
[delivery.elements]
CH_ratio = 1.0 # C/H mass ratio in mantle/atmosphere system
H_oceans = 0.0 # Hydrogen inventory in units of equivalent Earth oceans, by mass
# H_oceans = 0.0 # Hydrogen inventory in units of equivalent Earth oceans
H_ppmw = 109.0 # Hydrogen inventory in ppmw relative to mantle mass
N_ppmw = 2.01 # Nitrogen inventory in ppmw relative to mantle mass, by mass
S_ppmw = 235.0 # Sulfur inventory in ppmw relative to mass of melt

CH_ratio = 1.0 # C/H mass ratio in mantle/atmosphere system
# C_ppmw = 0.0 # Carbon inventory in ppmw relative to mantle mass

NH_ratio = 0.018 # N/H mass ratio in mantle/atmosphere system
# N_ppmw = 0.0 # Nitrogen inventory in ppmw relative to mantle mass

SH_ratio = 2.16 # S/H mass ratio in mantle/atmosphere system
# S_ppmw = 0.0 # Sulfur inventory in ppmw relative to mantle mass

# Set initial volatile inventory by partial pressures in atmosphere
[delivery.volatiles]
Expand Down
14 changes: 10 additions & 4 deletions input/trappist1c.toml
Original file line number Diff line number Diff line change
Expand Up @@ -232,11 +232,17 @@ author = "Harrison Nicholls, Tim Lichtenberg"

# Set initial volatile inventory by planetary element abundances
[delivery.elements]
H_oceans = 8.0 # Hydrogen inventory in units of equivalent Earth oceans
# H_ppmw = 0.0 # Hydrogen inventory in ppmw relative to mantle mass

CH_ratio = 1.0 # C/H mass ratio in mantle/atmosphere system
H_oceans = 8.0 # Hydrogen inventory in units of equivalent Earth oceans, by mass
H_ppmw = 0.0 # Hydrogen inventory in ppmw relative to mantle mass
N_ppmw = 2.01 # Nitrogen inventory in ppmw relative to mantle mass, by mass
S_ppmw = 235.0 # Sulfur inventory in ppmw relative to mass of melt
# C_ppmw = 0.0 # Carbon inventory in ppmw relative to mantle mass

# NH_ratio = 0.0 # N/H mass ratio in mantle/atmosphere system
N_ppmw = 2.01 # Nitrogen inventory in ppmw relative to mantle mass

# SH_ratio = 0.0 # S/H mass ratio in mantle/atmosphere system
S_ppmw = 235.0 # Sulfur inventory in ppmw relative to mantle mass

# Set initial volatile inventory by partial pressures in atmosphere
[delivery.volatiles]
Expand Down
32 changes: 25 additions & 7 deletions src/proteus/config/_delivery.py
Original file line number Diff line number Diff line change
Expand Up @@ -10,24 +10,42 @@
class Elements:
"""Initial volatile inventory by planetary element abundances.
For hydrogen: only H_oceans or H_ppmw should be used at any one time.
For X in {C, N, S}: only XH_ratio or X_ppmw should be used at any one time.
Attributes
----------
CH_ratio: float
Volatile C/H mass ratio in combined mantle+atmosphere system.
H_oceans: float
Absolute hydrogen inventory, units of equivalent Earth oceans.
H_ppmw: float
Relative hydrogen inventory, ppmw relative to mantle mass.
CH_ratio: float
Carbon metallicity. C/H mass ratio in combined mantle+atmosphere system.
C_ppmw: float
Relative carbon inventory, ppmw relative to mantle mass.
NH_ratio: float
Nitrogen metallicity. N/H mass ratio in combined mantle+atmosphere system.
N_ppmw: float
Relative nitrogen inventory, ppmw relative to mantle mass.
SH_ratio: float
Sulfur metallicity. C/H mass ratio in combined mantle+atmosphere system.
S_ppmw: float
Absolute sulfur inventory, ppmw relative to mantle mass.
"""
CH_ratio: float = field(validator=gt(0))
H_oceans: float = field(validator=ge(0))
H_ppmw: float = field(validator=ge(0))
N_ppmw: float = field(validator=ge(0))
S_ppmw: float = field(validator=ge(0))
H_oceans: float = field(default=0.0, validator=ge(0))
H_ppmw: float = field(default=0.0, validator=ge(0))

CH_ratio: float = field(default=0.0, validator=ge(0))
C_ppmw: float = field(default=0.0, validator=ge(0))

NH_ratio: float = field(default=0.0, validator=ge(0))
N_ppmw: float = field(default=0.0, validator=ge(0))

SH_ratio: float = field(default=0.0, validator=ge(0))
S_ppmw: float = field(default=0.0, validator=ge(0))


@define
Expand Down
112 changes: 86 additions & 26 deletions src/proteus/outgas/calliope.py
Original file line number Diff line number Diff line change
Expand Up @@ -19,11 +19,16 @@

log = logging.getLogger("fwl."+__name__)

# Constants
mass_ocean = ocean_moles * molar_mass['H2']

def construct_options(dirs:dict, config:Config, hf_row:dict):
"""
Construct CALLIOPE options dictionary
"""

invalid = False # Invalid options set by config

solvevol_inp = {}

# Planet properties
Expand All @@ -36,30 +41,89 @@ def construct_options(dirs:dict, config:Config, hf_row:dict):
solvevol_inp["T_magma"] = hf_row["T_magma"]
solvevol_inp['fO2_shift_IW'] = config.outgas.fO2_shift_IW

# Sum hydrogen absolute and relative amounts...
# Calculate hydrogen inventory...

# absolute part
# internally, calliope will convert this to mass in kg as:
# H_kg = H_oceans * number_ocean_moles * molar_mass['H2']
H_abs = float(config.delivery.elements.H_oceans)
# absolute part (H_kg = H_oceans * number_ocean_moles * molar_mass['H2'])
H_abs = float(config.delivery.elements.H_oceans) * mass_ocean

# relative part
# H_kg = H_rel * 1e-6 * M_mantle
# then converted to units of earth oceans, and summed with absolute part
# relative part (H_kg = H_rel * 1e-6 * M_mantle)
H_rel = config.delivery.elements.H_ppmw * 1e-6 * hf_row["M_mantle"]
H_rel /= ocean_moles * molar_mass['H2']

# avoid floating point errors here
if H_abs < 1e-10:
H_abs = 0.0
if H_rel < 1e-10:
H_rel < 0.0

# Elemental inventory
solvevol_inp['hydrogen_earth_oceans'] = H_abs + H_rel
solvevol_inp['CH_ratio'] = config.delivery.elements.CH_ratio
solvevol_inp['nitrogen_ppmw'] = config.delivery.elements.N_ppmw
solvevol_inp['sulfur_ppmw'] = config.delivery.elements.S_ppmw
# use whichever was set (one of these will be zero)
if H_abs < 1.0:
if H_rel < 1.0:
log.error("Hydrogen inventory is unspecified")
invalid = True
else:
H_kg = H_rel
elif H_rel < 1.0:
H_kg = H_abs
else:
log.error("Hydrogen inventory must be specified by H_oceans or H_ppmw, not both")
invalid = True
H_kg = -1 # dummy value

# Calculate carbon inventory (we need CH_ratio for calliope)
CH_ratio = float(config.delivery.elements.CH_ratio)
C_ppmw = float(config.delivery.elements.C_ppmw)
if CH_ratio > 1e-10:
# check that C_ppmw isn't also set
if C_ppmw > 1e-10:
log.error("Carbon inventory must be specified by CH_ratio or C_ppmw, not both")
invalid = True
else:
if C_ppmw < 1e-10:
log.error("Carbon inventory is unspecified")
invalid = True
else:
# calculate C/H ratio for calliope from C_kg and H_kg
CH_ratio = config.delivery.elements.C_ppmw * 1e-6 * hf_row["M_mantle"] / H_kg

# Calculate nitrogen inventory (we need N_ppmw for calliope)
NH_ratio = float(config.delivery.elements.NH_ratio)
N_ppmw = float(config.delivery.elements.N_ppmw)
if NH_ratio > 1e-10:
# check that N_ppmw isn't also set
if N_ppmw > 1e-10:
log.error("Nitrogen inventory must be specified by NH_ratio or N_ppmw, not both")
invalid = True
# calculate N_ppmw
N_ppmw = 1e6 * NH_ratio * H_kg / hf_row["M_mantle"]
else:
if N_ppmw < 1e-10:
log.error("Nitrogen inventory is unspecified")
invalid = True


# Calculate sulfur inventory (we need S_ppmw for calliope)
SH_ratio = float(config.delivery.elements.SH_ratio)
S_ppmw = float(config.delivery.elements.S_ppmw)
if SH_ratio > 1e-10:
# check that S_ppmw isn't also set
if S_ppmw> 1e-10:
log.error("Sulfur inventory must be specified by SH_ratio or S_ppmw, not both")
invalid = True
# calculate S_ppmw
S_ppmw = 1e6 * SH_ratio * H_kg / hf_row["M_mantle"]
else:
if S_ppmw < 1e-10:
log.error("Sulfur inventory is unspecified")
invalid = True

# Volatile abundances are over-specified in the config file.
# The code exits here, rather than above, in case there are multiple
# instances of volatiles being over-specified in the file.
if invalid:
log.error(" a) set X by metallicity, e.g. XH_ratio=1.2 and X_ppmw=0")
log.error(" b) set X by concentration, e.g. XH_ratio=0 and X_ppmw=2.01")
UpdateStatusfile(dirs, 20)
exit(1)

# Pass elemental inventory
solvevol_inp['hydrogen_earth_oceans'] = H_kg / mass_ocean
solvevol_inp['CH_ratio'] = CH_ratio
solvevol_inp['nitrogen_ppmw'] = N_ppmw
solvevol_inp['sulfur_ppmw'] = S_ppmw

# Volatile inventory
for s in vol_list:
Expand All @@ -70,7 +134,8 @@ def construct_options(dirs:dict, config:Config, hf_row:dict):

if (s in ("H2O","CO2","N2","S2")) and not included:
UpdateStatusfile(dirs, 20)
raise RuntimeError(f"Missing required volatile {s}")
log.error(f"Missing required volatile {s}")
exit(1)

return solvevol_inp

Expand All @@ -79,11 +144,6 @@ def calc_target_masses(dirs:dict, config:Config, hf_row:dict):
# make solvevol options
solvevol_inp = construct_options(dirs, config, hf_row)

# warn
if (config.delivery.elements.H_ppmw > 1e-10) \
and (config.delivery.elements.H_oceans > 1e-10):
log.warning("Hydrogen inventory set by summing `H_ppmw` and `H_oceans`")

# calculate target mass of atoms (except O, which is derived from fO2)
if config.delivery.initial == 'elements':
solvevol_target = get_target_from_params(solvevol_inp)
Expand Down
Loading

0 comments on commit fa011c6

Please sign in to comment.