Skip to content

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 str "self"
core_heatcap float str "self"
melting_dir str or none none Melting curve folder name in FWL_DATA, for the SPIDER structure module.
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: ":". Tabulated: "PALEOS:iron" (default), "Seager2007:iron". Analytic: "Analytic:iron", "Analytic:MgFeSiO3", etc.
mantle_eos str "PALEOS:MgSiO3" EOS for the mantle layer. Format: ":". Tabulated: "PALEOS:MgSiO3" (default), "PALEOS-2phase:MgSiO3", "Seager2007:MgSiO3", "WolfBower2018:MgSiO3". Analytic: "Analytic:MgSiO3", "Analytic:MgFeSiO3", etc.
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 unified EOS. 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, roughly matching the Stixrude+2014 cryoscopic depression for MgSiO3. Must be in [0.7, 1.0]. Only applies to PALEOS unified EOS; ignored for WolfBower2018 and RTPress100TPa (which use explicit melting curve files). This factor is applied consistently across Zalmoxis (density interpolation), SPIDER (phase boundaries), and the VolatileProfile phi-blending. 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 false 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 shipped with Zalmoxis when those tables are present, which keeps the properties well resolved across the melting-curve discontinuity that a single unified table interpolates through. If the two-phase tables are not available, the property surfaces are built from the unified table alone, and the entropy near the melting curve is less reliable. 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 constant solidus-to-liquidus ratio of the Stixrude (2014)3 MgSiO\(_3\) melting parametrization. The melt fraction then follows from the lever rule between this solidus and liquidus.

With mantle_eos = "PALEOS-2phase:MgSiO3", the solid and liquid tables define the phase boundaries directly. mushy_zone_factor is treated as 1.0 so the solidus coincides with the liquidus, and the latent-heat gap is supplied by the entropy difference between the solid and liquid tables rather than by a fixed temperature depression.

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, mantle-mass fraction, table layout, and the resolved mantle-EOS identity, so a run reuses the cache only when every one of those matches and different equations of state never collide. This 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-10 when left unset.
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.5 Centre of rheological transition in terms of melt fraction. 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.

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), '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, clamping the integrator to 1 yr steps (max_step = 1 yr, versus 100 yr otherwise) for the whole run and stalling it, so 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.
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.
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:

  • Boundary backend assumes a fixed surface state coupling.
  • Interior tidal heating requires an orbit 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

  1. 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. ↩

  2. 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. ↩

  3. Stixrude, L., Melting in super-earths, Philosophical Transactions of the Royal Society A, 372, 20130076, 2014. SciX. ↩