Interior structure and energetics¶
PROTEUS separates the interior into two coupled subsystems: structure (hydrostatic equilibrium, density profile, planet radius) and energetics (thermal evolution, melt fraction, heat flux). Each has its own module selection and parameters.
Submodule documentation: Aragog | SPIDER | Zalmoxis. See also Model description.
Interior structure [interior_struct]¶
The structure module computes the planet's density profile, radius, and core-mantle boundary position by solving the hydrostatic equilibrium equation with an equation of state (EOS).
| Parameter | Type | Default | Description |
|---|---|---|---|
core_frac |
float | 0.325 |
Fraction of the planet's interior radius corresponding to the core. Must be > 0 and < 1. |
core_frac_mode |
str | "mass" |
How core_frac is interpreted. 'radius': fraction of planet radius. 'mass': fraction of total planet mass. Only 'radius' is supported when module = 'spider'. The zalmoxis module always interprets core_frac as a mass fraction and ignores this flag (a warning is emitted at runtime if 'radius' is set with module = 'zalmoxis'). Choices: "radius", "mass". |
module |
str or none | "zalmoxis" |
Module for solving the planet's interior structure. Choices: none, "dummy", "spider", "zalmoxis". |
core_density |
float or str | "self" |
Density of the planet's core [kg m-3]. Set to 'self' for self-consistent calculation by Zalmoxis (requires module = 'zalmoxis'). |
core_heatcap |
float or str | "self" |
Specific heat capacity of the planet's core [J kg-1 K-1]. Set to 'self' for self-consistent calculation by Zalmoxis (requires module = 'zalmoxis'). |
melting_dir |
str or none | none |
Melting curve name in FWL_DATA. Required for the SPIDER structure module and for any run without a PALEOS table set; not read with module = 'zalmoxis' and a single PALEOS mantle EOS. With module = 'dummy' and a PALEOS mantle EOS, Aragog requires it but the solve uses the PALEOS set. |
eos_dir |
str or none | none |
EOS folder name in FWL_DATA, for the SPIDER structure module. |
Zalmoxis [interior_struct.zalmoxis]¶
Zalmoxis solves the full hydrostatic structure with tabulated EOS (PALEOS,
Chabrier, or Seager 2007). It supports 2-layer (core + mantle) and 3-layer
(core + mantle + ice/water) models. The newton_* parameters apply when
outer_solver = "newton"; the update_* and mesh_* parameters control
the structure update triggers (Zalmoxis structure updates are expensive and
decoupled from the main coupling loop); the miscibility_* parameters
belong to the experimental binodal-aware mode.
Equation of state
| Parameter | Type | Default | Description |
|---|---|---|---|
core_eos |
str | "PALEOS:iron" |
EOS for the core layer. Format: " |
mantle_eos |
str | "PALEOS:MgSiO3" |
EOS for the mantle layer. Format: " |
ice_layer_eos |
str or none | none |
EOS for the ice/water layer (3-layer model). 'none' for 2-layer model (core + mantle only). Tabulated: "PALEOS:H2O", "Seager2007:H2O". Analytic: "Analytic:H2O". |
mushy_zone_factor |
float | 0.8 |
Cryoscopic depression factor controlling the width of the mushy zone (partially molten region) in the PALEOS EOS family. Defines the solidus as T_sol = T_liq * mushy_zone_factor. 1.0 = sharp phase boundary (no mushy zone). 0.8 = solidus at 80% of the liquidus temperature, the constant solidus-to-liquidus ratio of Stixrude+2014 for MgSiO3, applied here to the PALEOS liquidus. Must be in [0.7, 1.0]. Applies to the PALEOS EOS family (PALEOS, PALEOS-2phase, PALEOS-API, PALEOS-API-2phase); ignored for WolfBower2018 and RTPress100TPa (which use explicit melting curve files) and for Seager2007/Analytic (no derived solidus). Must be >= 0.7 and <= 1.0. |
mantle_mass_fraction |
float | 0 |
Fraction of the planet's interior mass corresponding to the mantle. Required for 3-layer models (with ice layer) and for T-dependent 2-layer models (WolfBower2018, RTPress100TPa) where it partitions mass between core and mantle layers. Must be >= 0 and < 1. |
dry_mantle |
bool | true |
Structure EOS assumes a dry mantle. Set False for melt-fraction-aware dissolved-volatile mixing in the mantle density (per-shell volatile profile); the dissolved mass then stays inside the interior mass target. |
Grid and solver
| Parameter | Type | Default | Description |
|---|---|---|---|
num_levels |
int | 150 |
Number of Zalmoxis radius layers. |
outer_solver |
str | "newton" |
Outer mass-radius solver: 'newton' (recommended) or 'picard'. Choices: "picard", "newton". |
use_jax |
bool | true |
Use the JAX backend for the structure solver. |
use_anderson |
bool | true |
Anderson Type-II Picard acceleration on the density loop. |
solver_tol_outer |
float | 0.003 |
Relative tolerance for mass convergence (outer loop). Must be > 0. |
solver_tol_inner |
float | 0.0001 |
Relative tolerance for density convergence (inner loop). Must be > 0. |
solver_max_iter_outer |
int | 100 |
Max iterations for mass convergence (outer loop). Must be >= 10. |
solver_max_iter_inner |
int | 100 |
Max iterations for density convergence (inner loop). Must be >= 10. |
Newton solver tuning (used when outer_solver = "newton")
| Parameter | Type | Default | Description |
|---|---|---|---|
newton_max_iter |
int | 30 |
Maximum Newton iterations (outer_solver = 'newton'). Must be >= 5. |
newton_tol |
float | 0.0001 |
Newton convergence tolerance (outer_solver = 'newton'). Must be > 0. |
newton_relative_tolerance |
float | 1e-09 |
Integrator relative tolerance for the Newton path. Must be > 0. |
newton_absolute_tolerance |
float | 1e-10 |
Integrator absolute tolerance for the Newton path. Must be > 0. |
Structure update triggers
| Parameter | Type | Default | Description |
|---|---|---|---|
update_dphi_abs |
float | 0.05 |
Re-solve the structure when the global melt fraction changes by this absolute amount since the last solve. Must be > 0 and < 1. |
update_dtmagma_frac |
float | 0.05 |
Re-solve the structure when T_magma changes by this fraction. Must be > 0 and < 1. |
update_dw_comp_abs |
float | 0.05 |
Re-solve when the relative dissolved-volatile (H2O or H2) mantle mass fraction changes by this amount. Must be > 0 and < 1. |
update_interval |
float | 1000000000.0 |
Maximum time between structure updates [yr]; effectively disabled at the default. Must be >= 0. |
update_min_interval |
float | 0 |
Minimum time between structure updates [yr]; prevents thrashing. Must be >= 0. |
update_stale_ceiling |
float | 25000.0 |
Time since the last successful re-solve after which a trigger refires [yr]; 0 disables. Must be >= 0. |
mesh_max_shift |
float | 0.05 |
Maximum fractional radius shift per structure update. Must be > 0 and < 1. |
mesh_convergence_interval |
float | 10.0 |
Convergence relaxation time after a mesh update [yr]. Must be > 0. |
Initialisation
| Parameter | Type | Default | Description |
|---|---|---|---|
equilibrate_init |
bool | true |
Equilibrate structure and composition before the main loop. |
equilibrate_max_iter |
int | 15 |
Maximum equilibration iterations. Must be >= 1. |
equilibrate_tol |
float | 0.01 |
Equilibration convergence tolerance. Must be > 0. |
P-S entropy lookup tables
| Parameter | Type | Default | Description |
|---|---|---|---|
lookup_nP |
int | 1350 |
Number of pressure points in SPIDER P-S tables generated from PALEOS. Must be >= 100. |
lookup_nS |
int | 280 |
Number of entropy points in SPIDER P-S tables generated from PALEOS. Must be >= 50. |
Miscibility (experimental, not production-ready)
| Parameter | Type | Default | Description |
|---|---|---|---|
global_miscibility |
bool | false |
Enable the H2-silicate binodal-aware radial structure. True is rejected at config load: it requires the binodal handoff on the Zalmoxis side, which is not yet implemented. |
miscibility_max_iter |
int | 10 |
Maximum miscibility iterations. Must be >= 1. |
miscibility_tol |
float | 0.01 |
Miscibility convergence tolerance. Must be > 0. |
How the PALEOS mantle EOS is applied
The MgSiO\(_3\) mantle EOS resolves through two distinct paths depending on the <source> prefix of mantle_eos.
With mantle_eos = "PALEOS:MgSiO3" (the default), the hydrostatic structure solve uses the PALEOS unified MgSiO\(_3\) table for the density profile.
The phase-specific property surfaces used by Aragog (density, heat capacity, thermal expansion, adiabatic gradient) and the pressure-entropy lookup tables are built from the PALEOS two-phase solid and liquid tables, which keeps the properties well resolved across the melting-curve discontinuity that a single unified table interpolates through.
The liquidus is the analytic PALEOS curve (Belonoshko et al. 2005 below 2.55 GPa, Fei et al. 2021 above, in Simon-Glatzel form), and the solidus is derived as \(T_\mathrm{sol}(P) = f\,T_\mathrm{liq}(P)\) with \(f\) the mushy_zone_factor (default 0.8).
The default is the constant solidus-to-liquidus ratio (\(\approx 0.809\)) of the Stixrude (2014)3 MgSiO\(_3\) melting parametrization, applied to the PALEOS liquidus instead of the Stixrude liquidus.
The derived solidus is therefore a constant depression of the PALEOS liquidus, not the Stixrude solidus curve: its absolute value lies well above the Stixrude solidus at low to moderate pressure (about 2900 K against 1700 K at 20 GPa, and about 10% higher at 140 GPa) and below it above roughly 200 GPa.
The melt fraction then follows from the lever rule between this solidus and liquidus.
With mantle_eos = "PALEOS-2phase:MgSiO3", the SPIDER/Aragog entropy tables use the separate solid and liquid PALEOS tables, which supply the latent-heat entropy gap across the melting curve directly rather than through a single interpolated unified table.
Their phase boundaries follow the same construction as the unified case: the liquidus is the analytic PALEOS curve, and the solidus is derived as \(T_\mathrm{sol}(P) = f\,T_\mathrm{liq}(P)\) with \(f\) the mushy_zone_factor.
A run with a PALEOS EOS requires the table of each PALEOS layer, with the MgSiO\(_3\) unified table only when a layer is PALEOS:MgSiO3. A mantle with a PALEOS component also requires the MgSiO\(_3\) two-phase solid and liquid pair: the high-resolution pair for PALEOS-2phase:MgSiO3-highres, and the PALEOS-API cache for the PALEOS-API family.
When one is missing, the run stops before the first solve with an error that names each missing file and the command that fetches them, proteus get interiordata --config-path <config.toml>; it does not fall back to another table set.
With dry_mantle = false it also requires the tables of the dissolved volatiles (PALEOS:H2O and Chabrier:H).
A resumed SPIDER run that keeps its pressure-entropy tables needs only the tables of its layers. With temperature_mode = "liquidus_super" it also needs the two-phase pair that its anchor reads. An Aragog resume always needs the two-phase pair, which it reads when the mesh changes.
A mixture mantle (a + in mantle_eos) takes its melting curves, in the structure and in the energetics, and its SPIDER or Aragog tables from its MgSiO\(_3\) component: with PALEOS:MgSiO3 it gets exactly the curves and tables of a PALEOS:MgSiO3 mantle, and with WolfBower2018:MgSiO3 the eos_dir P-S tables (by default Wolf and Bower 2018) and the melting_dir curves. The other components enter only the structure density.
A PALEOS mantle with no MgSiO\(_3\) component (for example mantle_eos = "PALEOS:H2O", "PALEOS:iron", or a mixture of the two) uses the density of its own material in the structure, while the PALEOS MgSiO\(_3\) melting curves and the two-phase MgSiO\(_3\) tables of its PALEOS family stand in for the energetics.
With SPIDER or Aragog, the energetics must follow a PALEOS, WolfBower2018 or RTPress100TPa EOS; a mantle whose energetics would follow a Seager2007, Chabrier or analytic EOS is rejected when the configuration loads.
temperature_mode = "liquidus_super" solves the initial adiabat on the MgSiO\(_3\) two-phase tables for every mantle, so those tables are required at the start of the run.
The run logs one warning that names the cases that apply at the start of each run segment, so a resumed run logs it again.
mushy_zone_factor and the two-phase structure solve
In a PROTEUS-coupled run, load_zalmoxis_solidus_liquidus_functions builds the mushy_zone_factor * liquidus solidus described above and passes it into the Zalmoxis structure solve. The same curve pair sets the SPIDER/Aragog table boundaries, the adiabatic gradient in the mushy zone, and the two-phase density (PALEOS-2phase, PALEOS-API-2phase), so mushy_zone_factor acts consistently in the tables and in the structure.
A standalone Zalmoxis run uses its own rock_solidus and rock_liquidus keys (Stixrude 2014 by default). With those defaults a two-phase mantle does not depend on mushy_zone_factor; with no unified PALEOS layer elsewhere (for example a Seager2007:iron core), a value below 1.0 is rejected at validation. A unified PALEOS layer such as a PALEOS:iron core does honor the factor in its own density, so validation accepts it, but it does not change the two-phase mantle. Setting rock_liquidus = "PALEOS-liquidus" selects the PALEOS liquidus and derives the solidus as mushy_zone_factor * liquidus, as in the coupled case.
For a 1 \(M_\oplus\) planet with a Seager2007:iron core and a PALEOS-2phase:MgSiO3 mantle in a standalone Zalmoxis run with a linear temperature profile, the default Stixrude (2014) curves give a radius that does not depend on mushy_zone_factor, whereas rock_liquidus = "PALEOS-liquidus" makes the resolved radius respond to mushy_zone_factor, so mushy_zone_factor = 0.8 and 1.0 give different radii. Run this standalone case to quantify the shift for a given planet and table set.
Two-phase table versions
Two versions of the PALEOS two-phase MgSiO\(_3\) tables are in circulation: the set shipped in the Zalmoxis data directory, and the finer-grid set on Zenodo that the reference-data manifest fetches.
They are generated from the same PALEOS EOS but on different grids, so for exact reproducibility record which set a run used (the shipped tables resolve through the Zalmoxis material registry; the fetched tables land under FWL_DATA).
The wet mantle carries dissolved water only
With dry_mantle = false and the production outgassing module (CALLIOPE), the only species blended into the mantle density is H2O. This is the intersection of the two sides of the coupling, not an arbitrary restriction: a species needs both a dissolved mass from the outgassing chemistry and an EOS table in the structure solver, and under CALLIOPE H2O is the only species with both. CALLIOPE also dissolves CO2, CO, CH4, N2, and S2, but these have no structure-side EOS mapping. Their dissolved mass stays inside the interior mass target, and because the silicate blend fraction is defined as the remainder 1 - w_H2O, those kilograms are packed at silicate density: the mass budget stays exact while the compressibility of these species is misrepresented. The radius error scales as the species' dissolved mass fraction times its density contrast against silicate, negligible at typical inventories. If a composition ever pushes an unmapped species to percent-level dissolved mass fractions, the fix is an EOS mapping for it, not a bookkeeping change. H2 has an EOS mapping (Chabrier:H) but no solubility law in CALLIOPE; a dissolved H2 mass can arise from the atmodeller outgassing module (its H2 solubility law), from the H2-silicate binodal partitioning (outgas.h2_binodal), or from the radial miscibility solve (interior_struct.zalmoxis.global_miscibility). The last two are rejected at config load until that physics is production ready, and the atmodeller wet combination has not been through the wet-path conservation verification. Every verified configuration therefore exercises the wet path with water alone, which is exactly the case covered by the dissolved-water conservation verification. Whoever lifts either gate or runs atmodeller wet must extend the wet-path verification to H2 first: the per-shell binodal suppression removes dissolved H2 from the structure without returning it to the atmosphere, so water-grade mass closure does not carry over.
P-S entropy lookup tables
The lookup_nP/lookup_nS tables are derived from PALEOS EOS data at
runtime and cached. They provide the entropy-to-temperature mapping used by
Aragog and SPIDER. By default each run derives its own copy of these tables
under its output data/ directory. Set the PROTEUS_PS_CACHE_DIR
environment variable to an absolute path to instead share one derived copy
across runs: tables are stored in a subdirectory keyed by pressure ceiling,
resolution, mushy zone factor, table layout, the resolved mantle-EOS identity,
and the Zalmoxis table generator, so a run reuses the cache only when every one
of those matches and different equations of state never collide. The generator
identity is the Zalmoxis version plus a digest of two source files,
zalmoxis/eos_export.py and zalmoxis/melting_curves.py; a new version or an
edit to either file starts a new subdirectory. The digest does not cover other
Zalmoxis modules (for example the PALEOS-API cache modules) or the contents of
the EOS data files, which enter the key only through their resolved paths. A
resumed run (proteus start -r) keeps the energetics P-S tables it already
uses, in its own data/ directory or in the shared cache, whatever differs from
the current settings, including a changed mantle EOS or a new table generator,
unless they were built for another planet mass (another P_max, after an impact
or on a walk-back past one); those are replaced by the tables of the current mass.
The first time a process reads them, it logs one warning that names the
difference, or the reason the current key is not checked (for example a
PALEOS-API mantle EOS, which is not resolved on resume, or a mantle EOS that is
no longer PALEOS). The structure solve uses the current melting curves. A
resumed run with no kept tables (no completion marker with both phase-boundary
files) logs a warning and continues on the tables of the current key, which it
builds if they are absent. Sharing tables is the mitigation for the
per-run disk duplication that a grid or batch of same-EOS runs would
otherwise incur, since all such runs then read one shared copy. The cache
directory is not size-limited or auto-pruned; it grows with the number of
distinct pressure, resolution, and EOS combinations run against it, so point
it at storage sized for the campaign.
Interior energetics [interior_energetics]¶
The energetics module evolves the mantle thermal state (temperature,
entropy, melt fraction) and computes the interior heat flux that drives
atmospheric evolution. The section carries the shared transport, heating,
rheology, coupling-limit, and material-property parameters used by the
solver modules; the const_* parameters enable a constant-properties mode
that bypasses the EOS tables with an analytical
\(T(S) = T_\mathrm{ref} \exp((S - S_\mathrm{ref}) / C_p)\) relationship for
controlled parity tests.
| Parameter | Type | Default | Description |
|---|---|---|---|
module |
str | "aragog" |
Module for simulating the magma ocean. Choices: "spider", "aragog", "dummy", "boundary". |
num_levels |
int | 80 |
Number of radial grid levels for the energetics domain. Must be >= 40. |
rtol |
float | -1.0 |
Relative numerical tolerance for the interior ODE solver. SPIDER: -ts_sundials_rtol (used internally via atol_sf scaling). Aragog: scipy solve_ivp rtol. The deprecated aliases num_tolerance and [interior_energetics.spider].tolerance_rel copy into this field. Resolves to 1e-8 for Aragog and 1e-10 for SPIDER when left unset. With interior_struct.module = 'spider' it also sets the relative tolerance of the secant solve for the interior radius. |
num_tolerance |
float | -1.0 |
Deprecated alias for rtol; emits a DeprecationWarning when set. |
atol |
float | 1e-10 |
Absolute numerical tolerance for the interior ODE solver. SPIDER: -ts_sundials_atol (scaled by atol_sf at runtime). Aragog uses [interior_energetics.aragog].atol_temperature_equivalent instead because its state variable is entropy, not temperature, and a direct entropy-scale atol would be unintuitive to tune. Must be > 0. |
flux_guess |
float | -1 |
Initial heat flux guess [W m-2]. When < 0 (default), computed automatically as sigma * T_magma^4. Set to a positive value to prescribe a specific initial flux. Set to 0 for zero initial flux. |
surface_bc_mode |
str | "flux" |
Surface boundary condition for SPIDER/Aragog: 'flux' (prescribed F_atm from the atmosphere module) or 'grey_body' (native grey-body boundary condition computed inside the interior solver). Choices: "flux", "grey_body". |
Transport physics
| Parameter | Type | Default | Description |
|---|---|---|---|
trans_conduction |
bool | true |
Include conductive heat transfer. |
trans_convection |
bool | true |
Include convective heat transfer (mixing length theory). |
trans_grav_sep |
bool | true |
Include gravitational separation (Stokes settling). |
trans_mixing |
bool | true |
Include the chemical mixing flux. |
Heating
| Parameter | Type | Default | Description |
|---|---|---|---|
heat_tidal |
bool | false |
Include tidal heating?. |
heat_radiogenic |
bool | true |
Include radiogenic heat production?. |
radio_tref |
float | 4.567 |
Reference age for setting radioactive decay [Gyr]. Must be >= 0. |
radio_Al |
float | 0.0 |
Concentration (ppmw) of aluminium-26 at reference age t=radio_tref; 1.23 is the canonical early solar system value. Must be >= 0. |
radio_Fe |
float | 0.0 |
Concentration (ppmw) of iron-60 at reference age t=radio_tref. Must be >= 0. |
radio_K |
float | 310.0 |
Concentration (ppmw) of potassium-40 at reference age t=radio_tref. Must be >= 0. |
radio_U |
float | 0.031 |
Concentration (ppmw) of uranium at reference age t=radio_tref. Must be >= 0. |
radio_Th |
float | 0.124 |
Concentration (ppmw) of thorium-232 at reference age t=radio_tref. Must be >= 0. |
Rheology and convection
| Parameter | Type | Default | Description |
|---|---|---|---|
rfront_loc |
float | 0.4 |
Centre of rheological transition in terms of melt fraction. SPIDER receives it as -phi_critical. Aragog centres its viscosity blend, its reported rheological front and the melt-fraction ramp of kappah_floor on it; the boundary module uses it as the transition point of its viscosity models. SPIDER and Aragog use melt mass fraction (0.4 by volume is about 0.375 by mass at a 10 percent melt-solid density contrast); the boundary module takes its melt fraction linearly in temperature between solidus and liquidus and its rheological transition width from phase_transition_width, not rfront_wid. Must be > 0 and < 1. |
rfront_wid |
float | 0.2 |
Width of rheological transition in terms of melt fraction. Must be > 0 and < 1. |
grain_size |
float | 0.001 |
Crystal settling grain size [m]. Must be > 0. |
mixing_length |
str | "nearest" |
Mixing-length scale: 'nearest' (distance to nearest boundary) or 'constant' (a quarter of the mantle depth). Choices: "nearest", "constant". |
kappah_floor |
float | 10.0 |
Eddy diffusivity floor [m2 s-1]; prevents mixing-length transport from freezing out. Must be >= 0. |
Coupling limits
| Parameter | Type | Default | Description |
|---|---|---|---|
tmagma_atol |
float | 20.0 |
Maximum absolute change in T_magma per PROTEUS step [K]. Must be >= 0. |
tmagma_rtol |
float | 0.02 |
Maximum relative change in T_magma per PROTEUS step. Must be >= 0. |
tmagma_tides_step |
float | 10.0 |
Maximum change in T_magma allowed when tides are active [K]. Must be >= 0. |
Ultra-thin boundary layer
| Parameter | Type | Default | Description |
|---|---|---|---|
param_utbl |
bool | false |
Enable the ultra-thin boundary layer parameterisation. |
param_utbl_const |
float | 1e-07 |
Ultra-thin boundary layer scaling constant [K-1]. Must be > 0. |
Hydrostatic EOS (Adams-Williamson)
| Parameter | Type | Default | Description |
|---|---|---|---|
adams_williamson_rhos |
float | 4078.95095544 |
Adams-Williamson surface density [kg/m^3]. Matches SPIDER -adams_williamson_rhos and Aragog _MeshParameters.surface_density. Must be > 0. |
adams_williamson_beta |
float | 1.1115348931000002e-07 |
Adams-Williamson density gradient [1/m]. Matches SPIDER -adams_williamson_beta. Aragog derives its own via bulk modulus, so this value applies to SPIDER only. Must be > 0. |
adiabatic_bulk_modulus |
float | 260000000000.0 |
Adiabatic bulk modulus [Pa] used by Aragog's Adams-Williamson EOS (_MeshParameters.adiabatic_bulk_modulus). SPIDER derives its own. Must be > 0. |
Phase material properties
| Parameter | Type | Default | Description |
|---|---|---|---|
melt_log10visc |
float | 2.0 |
log10 viscosity of molten silicate [Pa s]. Matches SPIDER -melt_log10visc (2.0 = 1e2 Pa s). |
solid_log10visc |
float | 22.0 |
log10 viscosity of solid silicate [Pa s]. Matches SPIDER -solid_log10visc (22.0 = 1e22 Pa s). Shared by SPIDER and Aragog so both apply the same solid-phase rheology; a mis-set value diverges both solvers' solid-phase rheology by the same factor. |
melt_cond |
float | 4.0 |
Thermal conductivity of molten silicate [W/m/K]. Matches SPIDER -melt_cond. Must be > 0. |
solid_cond |
float | 4.0 |
Thermal conductivity of solid silicate [W/m/K]. Matches SPIDER -solid_cond. Must be > 0. |
eddy_diffusivity_thermal |
float | 1.0 |
Multiplier on the internally-computed thermal eddy diffusivity. SPIDER: -eddy_diffusivity_thermal (1.0 default). |
eddy_diffusivity_chemical |
float | 1.0 |
Multiplier on the internally-computed chemical eddy diffusivity. SPIDER: -eddy_diffusivity_chemical (1.0 default). |
latent_heat_of_fusion |
float | 4000000.0 |
Latent heat of fusion of silicate [J/kg]. Aragog uses this as a scalar in _PhaseMixedParameters. SPIDER derives it per-(P,S) from dS * T_fus via the EOS tables; the SPIDER derivation is more physically correct but this scalar is good to ~10% at Earth-mantle conditions. TODO: switch Aragog to SPIDER's derivation once the EntropyEOS exposes a dS_fus(P) method. Must be > 0. |
phase_transition_width |
float | 0.1 |
Width [fraction] of the mushy-zone transition in Aragog's _PhaseMixedParameters. Sets the width of the phase boundary in Aragog's mixed-phase blending (viscosity, thermal conductivity etc.). Distinct from [interior_energetics.spider].matprop_smooth_width which is SPIDER's analogous knob for its own solver. Must be > 0 and < 1. |
Core thermal model
| Parameter | Type | Default | Description |
|---|---|---|---|
core_tfac_avg |
float | 1.147 |
Core T_avg / T_cmb ratio from adiabatic gradient (Bower+2018 Table 2). Used by Aragog's _BoundaryConditionsParameters.tfac_core_avg. SPIDER derives its own internally. Must be > 0. |
Diagnostics
| Parameter | Type | Default | Description |
|---|---|---|---|
write_flux_diagnostics |
bool | false |
When True, Aragog's NetCDF output includes per-component flux decomposition (Jcond_b, Jconv_b, Jgrav_b, Jmix_b) and basic-node state variables (dSdr_b, eddy_diff_b, phi_basic_b, T/cp/rho_basic_b). Adds ~10 fields per snapshot; default False keeps output compact. Useful for diagnosing T_core and CMB-closure behaviour near phi=0. SPIDER path ignores this flag (uses SPIDER's own JSON output which already includes Jcond_b, Jconv_b, Jgrav_b, Jmix_b). |
Constant-properties mode
| Parameter | Type | Default | Description |
|---|---|---|---|
const_properties |
bool | false |
Enable constant-properties mode: bypass the EOS tables and use the analytical T(S) = T_ref * exp((S - S_ref) / Cp) relationship, for controlled parity tests. |
const_rho |
float | 4000.0 |
Constant density [kg/m3]. Must be > 0. |
const_Cp |
float | 1000.0 |
Constant heat capacity [J/kg/K]. Must be > 0. |
const_alpha |
float | 1e-05 |
Constant thermal expansivity [1/K]. Must be > 0. |
const_cond |
float | 4.0 |
Constant thermal conductivity [W/m/K]. Must be > 0. |
const_log10visc |
float | 2.0 |
Constant log10 dynamic viscosity [Pa.s]. |
const_T_ref |
float | 3500.0 |
Reference temperature for T(S) = T_ref * exp((S-S_ref)/Cp) [K]. Must be > 0. |
const_S_ref |
float | 3000.0 |
Reference entropy for T(S) [J/kg/K]. No positivity constraint since entropy reference states can be zero or negative. |
Aragog [interior_energetics.aragog]¶
Aragog is the default interior thermal evolution module. It solves the mantle energy equation with CVODE (SUNDIALS) using JAX-derived analytic Jacobians for robust convergence.
| Parameter | Type | Default | Description |
|---|---|---|---|
mass_coordinates |
bool | true |
Whether to use mass coordinates in the model. Default is True. Uses uniform spacing in mass coordinate space, giving larger cells at the surface where density is lower, matching SPIDER's mesh. |
backend |
str | "jax" |
ODE backend selector. Default 'jax'. - 'jax' : CVODE with JAX-derived RHS and JAX analytic Jacobian (jax.jacrev). Recommended for production. - 'numpy' : CVODE with numpy RHS and CVODE finite-difference Jacobian. Available for development and SPIDER-parity comparisons; less robust than 'jax' at production resolution because the FD-Jacobian noise can trip Aragog's T_core-jump retry guard at tight tolerances. The diffrax direct-JAX integration path is research-only and gated on a code-level flag in proteus.interior_energetics.aragog. Choices: "numpy", "jax". |
atol_temperature_equivalent |
float | 1e-08 |
Effective temperature-scale absolute tolerance [K] for Aragog's ODE integrator. Default 1e-8 matches SPIDER's atol=rtol=1e-8 setting; this tight tolerance avoids a marginal-stability bifurcation at the first dt jump after equilibration. Must be > 0. |
core_bc |
str | "energy_balance" |
Core-mantle boundary condition mode. Default 'energy_balance'. Valid values: - 'quasi_steady': alpha-factor heat-flux partition; gives about -19% T_core offset vs SPIDER. - 'energy_balance': SPIDER bit-parity BC with dSdr_cmb as a new state variable (mirrors SPIDER bc.c:76-131). - 'gradient': gradient-based state with two boundary entropies as state variables. - 'bower2018': experimental, do not use for production. Choices: "quasi_steady", "energy_balance", "gradient", "bower2018". |
phase_smoothing |
str | "tanh" |
Phase-boundary smoothing for Jgrav and Jmix: 'tanh' (SPIDER parity) or 'cubic_hermite'. Choices: "tanh", "cubic_hermite". |
separation_viscosity |
str | "mixture" |
Drag viscosity for the gravitational-separation velocity v_rel = |dRho| g F(phi) / eta. 'melt' (fixed single-phase liquid viscosity, SPIDER parity) keeps separation active below the rheological transition in coupled caps-off runs, which collapses the CMB temperature; 'mixture' (rheological-transition-blended bulk viscosity) ties the drag viscosity to the same solid-fraction rise that stiffens the bulk rheology, so separation locks up at the same melt fraction instead, which is why the default here is 'mixture' while Aragog's own default stays 'melt' for SPIDER parity. The regime boundaries are the porosities where adjacent permeability laws cross (Bower et al. 2018, section 2.1, Eqs. 13a to 13c). Choices: "melt", "mixture". |
solver_method |
str | "cvode" |
ODE solver: 'cvode' (SUNDIALS, SPIDER parity; needs scikits-odes-sundials and stops the run at setup without it), 'radau' (scipy), 'bdf' (scipy). Choices: "cvode", "radau", "bdf". |
scalar_gravity_override |
bool | false |
Scalar-gravity comparison knob. When True, the external mesh file that Zalmoxis writes has its gravity column overwritten with a uniform scalar (the surface value from hf_row['gravity']) before Aragog reads it, so Aragog's per-node gravity path interpolates to that scalar everywhere. False by default; set True only when running a paired scalar-gravity comparison. |
phi_step_cap |
float | 0.0 |
Per-call melt-fraction step cap. When > 0 and any staggered cell is in or near the two-phase window at solve() entry, a CVODE root function (and the equivalent scipy event) returns control at the exact time the larger of the global mass-weighted |ΔΦ| and the maximum single-cell |Δφ| reaches this cap. Off by default: the schema default 0.0 resolves to no cap, because on a benign freezing-front crossing the root function slices the coupled step into many small ones and drives the reported CMB heat flux briefly negative where the uncapped run stays positive, so the cap is a debugging control, not a production setting. Set a positive value to enable it; -1.0 is the explicit off spelling. An explicit 0.0 is rejected at load, since it cannot be told apart from the unset default; any other negative, NaN, or infinity is rejected too. |
temperature_step_cap |
float | 0.0 |
Per-call per-cell temperature step cap [K]. Shares the same root function as phi_step_cap and fires on the maximum single-cell |ΔT| since solve() entry. When enabled it bounds the per-cell temperature change on the solid adiabat just below the solidus, where the melt-fraction cap cannot act because a fully solid cell's melt fraction no longer moves. Off by default (schema default 0.0 resolves to no cap); on a benign freezing-front crossing the caps slice the coupled step into many small ones and drive the reported CMB heat flux briefly negative, so they are a debugging control, not a production setting. Set a positive value to enable it; -1.0 is the explicit off spelling. An explicit 0.0 is rejected at load, since it cannot be told apart from the unset default; any other negative, NaN, or infinity is rejected too. |
entropy_step_cap |
float | 0.0 |
Per-call per-cell entropy step cap [J/kg/K], in the native solver variable; same role as temperature_step_cap without an EOS lookup in the root function. Off by default (schema default 0.0 resolves to no cap); on a benign freezing-front crossing the caps slice the coupled step into many small ones and drive the reported CMB heat flux briefly negative, so they are a debugging control, not a production setting. Set a positive value to enable it; -1.0 is the explicit off spelling. An explicit 0.0 is rejected at load, since it cannot be told apart from the unset default; any other negative, NaN, or infinity is rejected too. |
phase_boundary_entropy_margin |
float | 200.0 |
Phase-boundary proximity band [J/kg/K] within which a staggered cell counts as near a solidus or liquidus crossing, tightening the integrator max_step so CVODE resolves the stiff two-phase RHS across the boundary. This is a solver-accuracy control, not a cosmetic step-size setting: at the default it reproduces the fixed band and the converged trajectory is unchanged, but lowering it below the default can under-resolve a real phase crossing and shift the converged state, because CVODE's local error control can accept an over-large step across the near-discontinuous RHS. Keeping the default reproduces current behaviour, and modestly widening the band does not move a converged result because tighter steps only refine an adaptive integrator; but a value orders of magnitude above the default makes every cell count as near a boundary at all times, so every call runs under the phase-boundary step control: with phase_boundary_cap = 'fixed' that clamps the integrator to 1 yr steps (versus 100 yr otherwise) and stalls the run, and with 'rate' every call runs as CVODE segments with max_step between 1 and 100 yr. Keep the band of order a few hundred J/kg/K. Default 200.0, matching Aragog's own default; a positive value is required (0 or negative is not a valid disabled state for a proximity band). Must be > 0. |
phase_boundary_cap |
str | "rate" |
Phase-boundary step-size policy: 'fixed' tightens max_step to 1 yr whenever any cell is near or inside the two-phase band; 'rate' (default) uses event-driven CVODE segments with max_step scaled by approach time to the nearest phase boundary, clipped to [1, 100] yr. 'rate' needs the CVODE solver and a core_bc other than 'gradient'; otherwise Aragog logs it and uses 1 yr. Its accuracy is verified at rtol 1e-8; Aragog warns once per solver at an rtol above 1e-7. Choices: "fixed", "rate". |
tolerance_struct |
float | 100.0 |
Absolute mass tolerance [kg] for the secant solver in determine_interior_radius. Default 100 kg; pairs with Spider's matching field so both backends drive the same outer-loop convergence criterion. Must be > 0. |
SPIDER [interior_energetics.spider]¶
SPIDER is the C interior module. It requires PETSc and produces results consistent with Bower et al. (2018) 1.
| Parameter | Type | Default | Description |
|---|---|---|---|
solver_type |
str | "bdf" |
SUNDIALS integrator choice. Choices: "adams", "bdf". |
tolerance_rel |
float | -1.0 |
Deprecated alias for Interior.rtol. Set interior_energetics.rtol at the top level instead. |
matprop_smooth_width |
float | 0.01 |
Melt-fraction window width for smoothing material properties across the solidus/liquidus. Passed to SPIDER as -matprop_smooth_width and to Aragog via _PhaseMixedParameters. Must be > 0 and < 1. |
tolerance_struct |
float | 100.0 |
Absolute mass tolerance [kg] for the interior-radius secant solver. Must be > 0. |
log_output |
bool | true |
Write SPIDER solver log output. |
Dummy [interior_energetics.dummy]¶
A parameterised cooling model with prescribed solidus and liquidus.
| Parameter | Type | Default | Description |
|---|---|---|---|
mantle_tliq |
float | 2700.0 |
Mantle liquidus temperature [K]. Must be >= 0. |
mantle_tsol |
float | 1700.0 |
Mantle solidus temperature [K]. Must be >= 0. |
mantle_rho |
float | 4550.0 |
Mantle density [kg m-3] for the dummy density profile and the fallback mantle-mass estimate. When the interior structure provides the masses, the mantle mass is taken as M_int - M_core instead. Must be > 0. |
mantle_cp |
float | 1792.0 |
Mantle specific heat capacity [J kg-1 K-1]. Must be >= 0. |
heat_internal |
float | 0.0 |
Fixed internal heating rate (e.g. radiogenic) [W kg-1]. Tidal heating is handled separately via heat_tidal and is added on top of this value. Must be >= 0. |
Boundary [interior_energetics.boundary]¶
A 0-D box model for the mantle thermal evolution based on Schaefer et al. (2016) 2, with prescribed solidus and liquidus and parameterised convective heat transport.
| Parameter | Type | Default | Description |
|---|---|---|---|
T_solidus |
float | 1420.0 |
Mantle solidus temperature [K]. Must be >= 0. |
T_liquidus |
float | 2020.0 |
Mantle liquidus temperature [K]. Must be > 0. |
critical_rayleigh_number |
float | 1100.0 |
Critical Rayleigh number for onset of convection [-]. Must be > 0. |
nusselt_exponent |
float | 0.33 |
Nusselt-Rayleigh scaling exponent [-]. Must be > 0. |
silicate_heat_capacity |
float | 1200.0 |
Silicate heat capacity [J/kg/K]. Must be > 0. |
core_density |
float | 10738.0 |
Core density [kg/m^3]. Must be > 0. |
core_shear |
float | 0.1 |
Core shear modulus [Pa]. Must be > 0. |
core_bulk |
float | 500000000000.0 |
Core bulk modulus [Pa]. Must be > 0. |
atm_heat_capacity_const |
bool | true |
Always use fallback atmosphere heat capacity?. |
atm_heat_capacity |
float | 17000.0 |
Used as fallback for atmosphere heat capacity when layer-specific value is not available [J/kg/K]. Must be > 0. |
silicate_density |
float | 4103.0 |
Silicate density [kg/m^3]. Default taken from Fei et. al. 2021 (https://ui.adsabs.harvard.edu/abs/2021NatCo..12..876F). Must be > 0. |
thermal_conductivity |
float | 4.2 |
Thermal conductivity [W/m/K]. Must be > 0. |
thermal_diffusivity |
float | 1e-06 |
Thermal diffusivity [m^2/s]. Must be > 0. |
thermal_expansivity |
float | 2e-05 |
Thermal expansivity [1/K]. Must be > 0. |
viscosity_model |
int | 2 |
Viscosity parameterisation model. Choices: 1, 2, 3. |
dynamic_viscosity |
float | 3800000000.0 |
Reference dynamic viscosity [Pa s] for Arrhenius solid mantle model. Must be > 0. |
activation_energy |
float | 350000.0 |
Activation energy [J/mol] for Arrhenius solid mantle model. Must be > 0. |
creep_parameter |
float | 26.0 |
Creep parameter [-] for Arrhenius solid mantle model. Must be > 0. |
viscosity_prefactor |
float | 0.00024 |
Viscosity prefactor [Pa s] for Vogel-Fulcher-Tammann magma ocean model. Must be > 0. |
viscosity_activation_temp |
float | 4600 |
Activation temperature [K] for Vogel-Fulcher-Tammann magma ocean model. Must be > 0. |
logging |
bool | false |
Whether to create diagnostic CSV data files from boundary interior module. |
Constraints¶
Cross-field constraints enforced when the config file loads:
- Require
interior_struct.melting_dirfor Aragog without a PALEOS table set. - Boundary backend assumes a fixed surface state coupling.
- Check the mantle EOS that the dummy structure hands to SPIDER or Aragog.
- Require a PALEOS or temperature-dependent energetics key for SPIDER or Aragog.
- Interior tidal heating requires an tides module to be enabled.
- Aragog requires at least one energy transport term to be enabled.
- Validate Boundary backend's solidus/liquidus ordering.
- Dummy interior requires the liquidus to sit above the solidus.
- SPIDER requires at least one energy transport term to be enabled.
- Validate Zalmoxis EOS format strings and reject unsupported miscibility options.
See also: Interior modules | Structure module | Melting curves
-
Bower, D.J., Sanan, P. & Wolf, A.S., Numerical solution of a non-linear conservation law applicable to the interior dynamics of partially molten planets, Physics of the Earth and Planetary Interiors, 274, 49-62, 2018. SciX. ↩
-
Schaefer, L., Wordsworth, R.D., Berta-Thompson, Z. & Sasselov, D., Predictions of the atmospheric composition of GJ 1132b, The Astrophysical Journal, 829, 63, 2016. SciX. ↩
-
Stixrude, L., Melting in super-earths, Philosophical Transactions of the Royal Society A, 372, 20130076, 2014. SciX. ↩