diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 903767b3..41d0e9a2 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -69,7 +69,7 @@ jobs:
key: ${{ runner.os }}-${{ steps.os.outputs.release }}-petsc-63b72503
fail-on-cache-miss: true
- name: Install Python requirements
- run: pip install -r py/requirements.txt pytest pytest-timeout ruff gcovr 'proteus-mpl>=1.4.0'
+ run: pip install -r py/requirements.txt pytest pytest-timeout ruff gcovr 'proteus-mpl>=1.5.0'
- name: Build with coverage instrumentation
run: |
make -j CFLAGS_EXTRA="--coverage"
@@ -159,7 +159,7 @@ jobs:
steps:
- uses: actions/checkout@v6
- name: Install Python requirements
- run: pip install -r py/requirements.txt pytest pytest-timeout 'proteus-mpl>=1.4.0'
+ run: pip install -r py/requirements.txt pytest pytest-timeout 'proteus-mpl>=1.5.0'
- name: Run installer
run: ./tools/get_spider.sh
- name: Check the installer produced a runnable binary
diff --git a/.github/workflows/nightly.yml b/.github/workflows/nightly.yml
index 0b5dae32..fae21f76 100644
--- a/.github/workflows/nightly.yml
+++ b/.github/workflows/nightly.yml
@@ -39,7 +39,7 @@ jobs:
./configure --with-fc=0 --with-cxx=0 --with-cc=gcc --download-sundials2 --download-mpich --download-f2cblaslapack PETSC_ARCH=arch-test
make all
- name: Install Python requirements
- run: pip install -r py/requirements.txt pytest pytest-timeout gcovr 'proteus-mpl>=1.4.0'
+ run: pip install -r py/requirements.txt pytest pytest-timeout gcovr 'proteus-mpl>=1.5.0'
- name: Build with coverage instrumentation
run: |
make -j CFLAGS_EXTRA="--coverage"
diff --git a/docs/Validation/atmosphere.md b/docs/Validation/atmosphere.md
index 5202f947..3f3963ae 100644
--- a/docs/Validation/atmosphere.md
+++ b/docs/Validation/atmosphere.md
@@ -1,19 +1,38 @@
# atmosphere.c: Surface energy balance
-**Source under test**: `atmosphere.c`
+## What this page checks
-## Cases
+This page checks the surface radiative energy balance in `atmosphere.c`, which models heat exchange.
+The atmosphere module calculates the thermal flux, which is exchanged between the planetary surface and space.
+Validation means that we test the grey-body radiation routine against the analytical Stefan-Boltzmann law, which provides an exact reference.
+If the surface radiation is wrong, the cooling rate of the magma ocean is incorrect.
+The section below details the validation check, its tolerance, and its test command.
-### Grey-body flux matches Stefan-Boltzmann limit
+## Grey-body flux matches Stefan-Boltzmann limit
-**Reference-pinned test**: `tests/test_atmosphere.py::test_grey_body_flux_matches_stefan_boltzmann_limit`
+**Why this check exists.** Numerical implementations of radiative flux can suffer from incorrect temperature exponents, or from inverted signs. A bug in the surface flux calculation leads to unphysical planetary cooling rates, without raising a runtime error.
-**Reference**: Analytical limit: the grey-body flux F = emissivity sigma (T_surf^4 - T_eqm^4) with emissivity 1 and T_eqm = 273 K (Stefan-Boltzmann).
+**What it tests.** The surface flux $F$ (W/m$^2$) of a grey body depends on the emissivity $\epsilon$ (0 to 1), the Stefan-Boltzmann constant $\sigma$ (W m$^{-2}$ K$^{-4}$), the surface temperature $T_\mathrm{surf}$ and the equilibrium temperature $T_\mathrm{eqm}$ (both in K). The test evaluates $F = \epsilon\sigma(T_\mathrm{surf}^4 - T_\mathrm{eqm}^4)$ with $\epsilon$ = 1 and $T_\mathrm{eqm}$ = 273 K, the values of the blackbody50 run (the standard 50-node test run, `tests/opts/blackbody50.opts`). It compares the reported surface flux against the Stefan-Boltzmann law, evaluated at the reported surface temperature. The test guards against exponent errors, where a T^3 exponent slip changes the flux by orders of magnitude at magma-ocean temperatures. The flux must remain positive (outward), and bounded to the 1e2 to 1e9 W/m2 regime.
-**Tolerance**: rel 1e-6 against the reported surface temperature, using the code constant sigma = 5.670367e-8 W m-2 K-4.
+**How we know the right answer.** The grey-body Stefan-Boltzmann law provides an exact analytical reference for thermal radiation, which governs ideal emitters. An exponent slip to third power changes the flux by orders of magnitude, which easily violates the tolerance of rel 1e-6.
-**Command**: `make`, then `pytest -o addopts= tests/test_atmosphere.py::test_grey_body_flux_matches_stefan_boltzmann_limit`
+**Why it matters.** The atmosphere module computes the surface energy balance, which drives planetary mantle cooling. PROTEUS reads the surface flux `Fatm` as `F_int`, and surface temperature `temperature_surface` as `T_magma`, to couple the interior to the atmosphere.
-**Discrimination guards**: A T^3 exponent slip changes the flux by orders of magnitude at magma-ocean temperatures; the flux is positive (outward) and bounded to the 1e2 to 1e9 W/m2 regime; monotone decline of flux and surface temperature over the first outputs.
+**Known limit.** The test uses its own copy of the value in `constants.c`, sigma = 5.670367e-8 W m-2 K-4. A change of the constant in the code larger than the tolerance fails the test, but a value that is wrong in both places does not. PROTEUS uses this grey-body branch only with `surface_bc_mode = 'grey_body'`; its default is a prescribed flux (the SPIDER option `SURFACE_BC 4`).
-The atmosphere module computes the surface energy balance that drives the cooling. The grey-body branch is pinned against the Stefan-Boltzmann law evaluated from the same output state, and the secular cooling trend is asserted across output times.
+{ width="80%" }
+{ width="80%" }
+
+**How to read the plot.** The x-axis is the surface temperature in K and the y-axis the surface flux in W/m$^2$ on a log scale. The line is the grey-body law $\epsilon\sigma(T^4 - T_\mathrm{eqm}^4)$, with emissivity $\epsilon$ = 1 and equilibrium temperature $T_\mathrm{eqm}$ = 273 K. Here $\sigma$ is the Stefan-Boltzmann constant, and the dashed line is the $T^3$ form, which the test rejects. The circle is the asserted surface flux Fatm of the blackbody50 run at 200 yr, at its reported surface temperature. Its values are shown to 3 significant digits, because the run output moves by about 1e-6 between builds of the PETSc library. A pass means that the circle lies on the line; the test asserts relative 1e-6 on the unrounded values.
+
+??? info "Test details"
+ - Reference: Analytical limit: the grey-body flux F = emissivity sigma (T_surf^4 - T_eqm^4) with emissivity 1 and T_eqm = 273 K (Stefan-Boltzmann).
+ - Tolerance: rel 1e-6 on the flux, computed from the reported surface temperature using the test's copy of the code value sigma = 5.670367e-8 W m-2 K-4.
+ - Test id: `tests/test_atmosphere.py::test_grey_body_flux_matches_stefan_boltzmann_limit`
+ - Command: `make`, then `pytest -o addopts= tests/test_atmosphere.py::test_grey_body_flux_matches_stefan_boltzmann_limit` (add `--validation-plots` to write the plot pair)
+ - Companion test: `tests/test_atmosphere.py::test_surface_cools_and_flux_decays_across_macro_steps` confirms monotone decline of flux and surface temperature over outputs.
+
+## Related code
+
+- `atmosphere.c`: `set_atmosphere_emissivity_and_flux`.
+- `tests/test_atmosphere.py`: `test_grey_body_flux_matches_stefan_boltzmann_limit`, `test_surface_cools_and_flux_decays_across_macro_steps`.
diff --git a/docs/Validation/bc.md b/docs/Validation/bc.md
index ca35e588..7c75b3b5 100644
--- a/docs/Validation/bc.md
+++ b/docs/Validation/bc.md
@@ -1,31 +1,62 @@
# bc.c: Boundary conditions
-**Source under test**: `bc.c`
+## What this page checks
-## Cases
+This page checks the boundary condition (BC) implementations in `bc.c`, which control surface and bottom boundaries.
+The boundary condition routines specify how energy and entropy enter, or leave, the planetary mantle.
+Validation means that we test boundary condition options against analytical identities, and against prescribed inputs.
+If boundary condition dispatch fails, energy conservation is violated at the surface, or at the core-mantle boundary (CMB). The SPIDER options `SURFACE_BC` and `CORE_BC` select the boundary condition at the surface and at the CMB.
+The sections below detail two validation cases, their tolerances, and their test commands.
-### Prescribed-flux BC echoes value and inverts emissivity
+## Constant-entropy BC pins the initial surface
-**Reference-pinned test**: `tests/test_bc.py::test_prescribed_flux_bc_echoes_value_and_inverts_emissivity`
+**Why this check exists.** Fixed boundary values can be overwritten by initial rebalancing steps, or by stencil operations. If the surface boundary assignment fails, the mantle profile starts from an uncalibrated thermodynamic state.
-**Reference**: Analytical identity: with SURFACE_BC 4 the surface flux equals the prescribed value, and the reported emissivity is its grey-body inversion.
+**What it tests.** The test evaluates the constant-entropy boundary option, which writes `ic_surface_entropy` (2550 J/kg/K) onto the surface basic node (the outermost cell boundary). This value is preserved through the steady-state energy solve into the initial output. The test discriminates this setting against the default configuration surface value (2599.7 J/kg/K from the adiabat, the constant-entropy profile, at 2600), which sits about 50 J/kg/K away. Sibling test `test_prescribed_core_flux_slows_cmb_cooling` checks CORE_BC dispatch, CMB entropy difference, and doubled CMB heat flux.
-**Tolerance**: rel 1e-9 on the flux echo; rel 1e-6 on the emissivity inversion.
+**How we know the right answer.** Preserving a specified boundary value is an exact identity. Because the value is assigned directly, it matches to rel 1e-10.
-**Command**: `make`, then `pytest -o addopts= tests/test_bc.py::test_prescribed_flux_bc_echoes_value_and_inverts_emissivity`
+**Why it matters.** Setting fixed surface or core entropies allows controlled numerical experiments and comparisons with simplified interior profiles.
-**Discrimination guards**: The grey-body flux at the same surface temperature differs from the prescribed value by more than 10 percent (dispatch guard), and the emissivity differs from 1 by more than 0.5; the flux is positive and between 1e2 and 1e9 W/m^2. The sibling test `test_mantle_cools_from_the_top_between_surface_and_cmb` (same file) asserts the temperature ordering across the mantle (CMB hotter than surface) and the near-zero surface pressure on the default configuration.
+**Known limit.** For constant-entropy boundary conditions, the surface value is assigned directly rather than determined by a physical flux balance.
-The boundary-condition dispatch is pinned through the prescribed-flux pathway PROTEUS uses in coupled runs: the value handed in through the option must reappear identically in the output, with the emissivity consistently inverted from it.
+{ width="80%" }
+{ width="80%" }
-### Constant-entropy BC pins the initial surface
+**How to read the plot.** The x-axis is the radius and the y-axis the entropy at the basic nodes (the cell boundaries of the mesh) at $t$ = 0. The solid line is the asserted profile of the run with a constant-entropy surface boundary. The dashed line is the default run (context); its surface entropy, 2599.7 J kg$^{-1}$ K$^{-1}$, must differ from 2550 by more than 40. The circle is the asserted surface node, and the hollow square the value 2550 J kg$^{-1}$ K$^{-1}$ of the option `ic_surface_entropy`. A pass means that the circle sits in the square (relative 1e-10) and the deep end of the solid line lies between 2400 and 2700.
-**Reference-pinned test**: `tests/test_bc.py::test_constant_entropy_bc_pins_the_initial_surface`
+??? info "Test details"
+ - Reference: Analytical identity: the constant-entropy boundary writes ic_surface_entropy (2550 J/kg/K) onto the surface basic node, preserved through the steady-state energy solve into the initial output.
+ - Tolerance: rel 1e-10 (assigned, not solved); the core-side entropy between 2400 and 2700 J/kg/K.
+ - Test id: `tests/test_bc.py::test_constant_entropy_bc_pins_the_initial_surface`
+ - Command: `make`, then `pytest -o addopts= tests/test_bc.py::test_constant_entropy_bc_pins_the_initial_surface` (add `--validation-plots` to write the plot pair)
+ - Companion test: `tests/test_bc.py::test_prescribed_core_flux_slows_cmb_cooling` verifies that core heat flux slows cooling at the core-mantle boundary.
-**Reference**: Analytical identity: the constant-entropy boundary writes ic_surface_entropy (2550 J/kg/K) onto the surface basic node, preserved through the steady-state energy solve into the initial output.
+## Prescribed-flux BC echoes value and inverts emissivity
-**Tolerance**: rel 1e-10 (assigned, not solved).
+**Why this check exists.** With a prescribed surface flux, the reported flux must equal the prescribed value and the emissivity must be its grey-body inversion. A wrong boundary-condition (BC) branch or a wrong inversion gives a wrong reported flux or emissivity without an error message.
-**Command**: `make`, then `pytest -o addopts= tests/test_bc.py::test_constant_entropy_bc_pins_the_initial_surface`
+**What it tests.** The test evaluates boundary condition SURFACE_BC 4, where the surface flux equals the prescribed value. The reported emissivity must equal the grey-body inversion value. The test guards against wrong dispatch, confirming that the grey-body flux at the surface temperature differs from the prescribed value by more than 10 percent, and the emissivity differs from 1 by more than 0.5. The flux must be positive, and lie between 1e2 and 1e9 W/m^2. Sibling test `test_mantle_cools_from_the_top_between_surface_and_cmb` asserts temperature ordering across the mantle, and near-zero surface pressure.
-**Discrimination guards**: The default configuration's surface value (2599.7 J/kg/K from the 2600 adiabat) sits about 50 J/kg/K away; positivity and table-range bounds on the core-side value after the steady rebalance. The sibling test `test_prescribed_core_flux_slows_cmb_cooling` (same file) compares against the core-cooling baseline and discriminates the CORE_BC dispatch by the sign of the CMB entropy difference and the doubled CMB heat flux.
+**How we know the right answer.** Under SURFACE_BC 4, reproducing the prescribed flux is an analytical identity. The test checks the flux echo to rel 1e-9, and the inverted emissivity to rel 1e-6.
+
+**Why it matters.** PROTEUS passes the atmospheric radiative flux into SPIDER using option `-SURFACE_BC 4` and `-surface_bc_value`, so SPIDER must echo this flux exactly.
+
+**Known limit.** The test reads only the reported `Fatm` and emissivity; it does not check that the interior total flux at the surface (`Jtot` at node 0, set in `bc.c`) uses the prescribed value. `SURFACE_BC` is read in `atmosphere.c`, `bc.c` and `energy.c`.
+
+{ width="80%" }
+{ width="80%" }
+
+**How to read the plot.** The x-axis is the surface temperature and the y-axis the surface flux on a log scale. The dotted line is the prescribed flux, 1e5 W/m$^2$, and the circle the asserted surface flux Fatm (relative 1e-9) at the reported surface temperature. The solid line is the grey-body law with emissivity 1; the test requires it to be more than 10 % away at that temperature. The dashed line is the grey-body law with the inverted emissivity, about 0.12, asserted to relative 1e-6. A pass means that the circle sits on the dotted and the dashed line, far below the solid one.
+
+??? info "Test details"
+ - Reference: Analytical identity: with SURFACE_BC 4 the surface flux equals the prescribed value, and the reported emissivity is its grey-body inversion.
+ - Tolerance: rel 1e-9 on the flux echo; rel 1e-6 on the emissivity inversion.
+ - Test id: `tests/test_bc.py::test_prescribed_flux_bc_echoes_value_and_inverts_emissivity`
+ - Command: `make`, then `pytest -o addopts= tests/test_bc.py::test_prescribed_flux_bc_echoes_value_and_inverts_emissivity` (add `--validation-plots` to write the plot pair)
+ - Companion test: `tests/test_bc.py::test_mantle_cools_from_the_top_between_surface_and_cmb` verifies that the core-mantle boundary remains hotter than the surface.
+
+## Related code
+
+- `bc.c`: `set_surface_flux_from_atmosphere`, `set_surface_entropy_constant`, `set_cmb_entropy_constant`.
+- `tests/test_bc.py`: `test_prescribed_flux_bc_echoes_value_and_inverts_emissivity`, `test_constant_entropy_bc_pins_the_initial_surface`.
diff --git a/docs/Validation/energy.md b/docs/Validation/energy.md
index f13da5a6..51f1bcbf 100644
--- a/docs/Validation/energy.md
+++ b/docs/Validation/energy.md
@@ -1,43 +1,86 @@
# energy.c: Energy fluxes
-**Source under test**: `energy.c`
+## What this page checks
-## Cases
+This page checks energy flux assembly and internal heating routines in `energy.c`, which control mantle energy transport.
+The energy module calculates 4 heat fluxes, alongside radiogenic and tidal heating. Jcond is conduction, and Jconv is convection, driven by the entropy gradient. Jmix is mixing of melt and solid, active only where both phases exist. Jgrav is the separation of melt from solid under gravity.
+Validation means that we test flux assembly against analytical superposition limits, decay laws, and the balance of the steady-state initial condition (IC).
+If flux calculation or source terms are wrong, mantle cooling rates and core heat extraction are wrong.
+The sections below detail three validation cases, their tolerances, and their test commands.
-### Disabled transport mechanisms carry zero flux
+## Disabled transport mechanisms carry zero flux
-**Reference-pinned test**: `tests/test_energy.py::test_disabled_transport_mechanisms_carry_zero_flux`
+**Why this check exists.** Runtime configuration flags can fail to disconnect disabled flux routines, or leak residual values into total energy balances. A failure in flux superposition produces incorrect thermal gradients without raising an error.
-**Reference**: Analytical limit: a disabled transport mechanism carries zero flux, and with only conduction enabled the total flux equals the conductive flux at every node except the surface node, which carries the grey-body boundary flux.
+**What it tests.** The test evaluates the flux superposition architecture when specific transport terms are switched off. A disabled transport mechanism must hold zero flux at every basic node (cell boundary), within abs 1e-20 W/m^2, which is at least 17 decades below the observed active fluxes. With conduction off, convection gives 9e8 to 2e9 W/m^2 (observed; the test asserts more than 1e3). In a conduction-only run, conduction gives 1.9e-3 to 1.3e-2 W/m^2 (observed). Total flux equals conductive flux to rel 1e-12 at every interior node. The surface node differs by more than 1e3 W/m^2.
-**Tolerance**: abs 1e-20 W/m^2 on the disabled fluxes, at least 17 decades below the active fluxes (convection 9e8 to 2e9 W/m^2 with conduction off; conduction 1.9e-3 to 1.3e-2 W/m^2 in the conduction-only run); rel 1e-12 on the conduction-only identity, and the surface node differs by more than 1e3 W/m^2.
+**How we know the right answer.** Superposition of disjoint flux mechanisms is an analytical identity: an inactive transport term must vanish to rounding precision, and in the conduction-only run the total flux must equal the conductive flux. The magnitude of the convective flux is not compared with a reference.
-**Command**: `make`, then `pytest -o addopts= tests/test_energy.py::test_disabled_transport_mechanisms_carry_zero_flux`
+**Why it matters.** SPIDER models mantle cooling by summing the transport terms, so a leaking term changes every flux of the run.
-**Discrimination guards**: Convection stays active when conduction is off (wrong-path guard); at the surface node the grey-body boundary flux differs from the conductive flux, so the identity holds at the interior nodes only; a scale bound on the full-physics flux. The sibling test in the same file asserts that the total flux is the sum of the components at the interior nodes.
+**Known limit.** For the conduction-only test, the identity holds at interior nodes only, because the surface node transmits the grey-body boundary flux.
-The flux assembly is validated through its superposition structure: each transport term can be switched off through its runtime option, and the total must reproduce the remaining terms to output precision. This pins the dispatch wiring between the options, the individual flux routines, and the total used by the time integration.
+{ width="80%" }
+{ width="80%" }
-### Radiogenic heating pins two isotope decay
+**How to read the plot.** The circles are the largest absolute flux over the basic nodes. The axis is symmetric log, linear between $-10^{-3}$ and $10^{-3}$. Jcond, Jconv, Jmix and Jgrav are the conductive, convective, mixing and separation fluxes, and Jtot is their sum. Tick labels name each case. Values are shown to 2 significant digits, because the run output moves by up to a few 1e-6 between builds of the PETSc library. A pass means that the first 2 circles sit at 0 (1e-20 W/m$^2$) and the last 2 above the dashed line at $10^3$ W/m$^2$.
-**Reference-pinned test**: `tests/test_energy.py::test_radiogenic_heating_pins_two_isotope_decay`
+??? info "Test details"
+ - Reference: Analytical limit: a disabled transport mechanism carries zero flux, and with only conduction enabled the total flux equals the conductive flux at every node except the surface node, which carries the grey-body boundary flux.
+ - Tolerance: abs 1e-20 W/m^2 on the disabled fluxes, at least 17 decades below the observed active fluxes (observed: convection 9e8 to 2e9 W/m^2 with conduction off, conduction 1.9e-3 to 1.3e-2 W/m^2 in the conduction-only run); rel 1e-12 on the conduction-only identity, and the surface node differs by more than 1e3 W/m^2.
+ - Test id: `tests/test_energy.py::test_disabled_transport_mechanisms_carry_zero_flux`
+ - Command: `make`, then `pytest -o addopts= tests/test_energy.py::test_disabled_transport_mechanisms_carry_zero_flux` (add `--validation-plots` to write the plot pair)
+ - Companion test: `tests/test_energy.py::test_total_flux_assembles_from_transport_components` verifies that total flux equals the sum of transport components at interior basic nodes.
-**Reference**: Analytical limit: radiogenic heating is the sum over isotopes of concentration times abundance times heat production, decaying as 2^(-(t - t0)/half_life). A two-isotope configuration with 100-year and 200-year half-lives pins the sum at t = 0 and the weighted decay factor 0.2506 after 200 years.
+## Radiogenic heating pins two isotope decay
-**Tolerance**: rel 1e-6 on the initial sum, rel 1e-4 on the decay factor.
+**Why this check exists.** Radioactive decay models can apply incorrect decay constants, confuse isotope abundances, or miscalculate exponential time evolution. A bug in heating terms leads to wrong internal heat production over geological timescales.
-**Command**: `make`, then `pytest -o addopts= tests/test_energy.py::test_radiogenic_heating_pins_two_isotope_decay`
+**What it tests.** The test evaluates radiogenic heating as the sum of concentration times abundance times heat production over isotopes, decaying as 2^(-(t - t0)/half_life), with t0 the start time of each isotope (0 here). A two-isotope test configuration with 100-year and 200-year half-lives pins the sum at t = 0, and checks the weighted decay factor 0.2506 after 200 years. A dropped isotope shifts the initial sum by 2.4e-3 relative. The two-isotope factor is discriminated from the single-isotope values: 0.25 for al26 alone (2 half-lives in 200 years) and 0.5 for k40 alone (1 half-life). It also differs from the no-decay value 1.0.
-**Discrimination guards**: A dropped isotope shifts the sum by 2.4e-3 relative; the two-isotope decay factor is resolved from the single-isotope values 0.25 and 0.5 and from the no-decay value 1.0.
+**How we know the right answer.** Radioactive decay follows the analytical exponential decay law. The initial sum is checked to rel 1e-6, and the combined decay factor matches to rel 1e-4.
-### Steady-state IC balances interior and surface
+**Why it matters.** Radiogenic decay supplies essential internal heating that delays mantle freezing. PROTEUS reads integrated radiogenic heating `F_radio` (in `ReadSPIDER` of `spider.py`) into whole-planet energy budgets.
-**Reference-pinned test**: `tests/test_energy.py::test_steady_state_ic_balances_interior_and_surface`
+**Known limit.** This test pins decay mechanics for synthetic isotopes, but does not benchmark geological isotope concentrations.
-**Reference**: Analytical limit: the steady-state initial condition equalises the energy flow through every basic node and closes the surface radiation balance sigma (T_surf^4 - teqm^4) with the Stefan-Boltzmann constant of constants.c.
+{ width="80%" }
+{ width="80%" }
-**Tolerance**: Relative energy-flow spread below 1e-8 (observed 2e-11); grey-body pin at rel 1e-5.
+**How to read the plot.** The x-axis is the time in yr and the y-axis the specific radiogenic heating in W/kg on a log scale. The wide pale line is the closed-form sum of the 2 isotopes. The dashed and dotted lines are the isotopes `al26` and `k40` alone, with test half-lives of 100 and 200 yr (context). The dashed line runs inside the wide line, because `al26` gives more than 99 % of the heating. The circles are the asserted Hradio at 0 and 200 yr: the value at 0 to relative 1e-6, and their ratio, about 0.2506, to 1e-4. A pass means that both circles lie on the wide line.
-**Command**: `make`, then `pytest -o addopts= tests/test_energy.py::test_steady_state_ic_balances_interior_and_surface`
+??? info "Test details"
+ - Reference: Analytical limit: radiogenic heating is the sum over isotopes of concentration times abundance times heat production, decaying as 2^(-(t - t0)/half_life). A two-isotope configuration with 100-year and 200-year half-lives pins the sum at t = 0 and the weighted decay factor 0.2506 after 200 years.
+ - Tolerance: rel 1e-6 on the initial sum, rel 1e-4 on the decay factor.
+ - Test id: `tests/test_energy.py::test_radiogenic_heating_pins_two_isotope_decay`
+ - Command: `make`, then `pytest -o addopts= tests/test_energy.py::test_radiogenic_heating_pins_two_isotope_decay` (add `--validation-plots` to write the plot pair)
+ - Companion test: `tests/test_energy.py::test_prescribed_tidal_heating_is_uniform_and_steady` verifies uniform spatial distribution of tidal heating.
-**Discrimination guards**: A T^3 exponent slip misses the flux by roughly a factor of the surface temperature (about 2000 here); sign and scale bounds on the flux.
+## Steady-state IC balances interior and surface
+
+**Why this check exists.** Steady-state solvers can terminate at false equilibria, or fail to balance heat flow across layers. An uncalibrated initial profile introduces unphysical thermal transients at the start of a simulation.
+
+**What it tests.** The test evaluates the steady-state initial condition of the blackbody50 run, which has no internal heating. It equalises the energy flow Etot (W) through every basic node, and closes the surface radiation balance $\sigma(T_\mathrm{surf}^4 - T_\mathrm{eqm}^4)$. Here $\sigma$ is the Stefan-Boltzmann constant of constants.c and $T_\mathrm{eqm}$ = 273 K the equilibrium temperature. A $T^3$ exponent slip changes the surface flux by a factor equal to the surface temperature in kelvin, about 2000 here (about 1960 K).
+
+**How we know the right answer.** In this configuration without internal heating, conservation of energy in steady state requires equal energy flow through every mantle shell, and balanced radiation at the surface. Energy-flow relative spread remains below 1e-8 (observed 2e-11), and the grey-body balance matches to rel 1e-5.
+
+**Why it matters.** The steady-state initial condition starts a run from a profile in energy balance.
+
+**Known limit.** The test uses its own copy of the Stefan-Boltzmann constant in `constants.c`, so a value that is wrong in both places does not fail it; a change in the code larger than the tolerance does.
+
+{ width="80%" }
+{ width="80%" }
+
+**How to read the plot.** Left: the asserted energy flow Etot at every basic node of the steady-state initial condition, against radius. A pass means a flat line, with a relative spread below 1e-8. Right: the solid line is the grey-body law $\sigma(T^4 - 273^4)$. Here $\sigma$ is the Stefan-Boltzmann constant and 273 K the equilibrium temperature. The dashed line is the rejected $T^3$ form. The circle is the asserted surface flux Fatm at the reported surface temperature, which must lie on the solid line (relative 1e-5).
+
+??? info "Test details"
+ - Reference: Analytical limit: the steady-state initial condition equalises the energy flow through every basic node and closes the surface radiation balance sigma (T_surf^4 - teqm^4) with the Stefan-Boltzmann constant of constants.c.
+ - Tolerance: Relative energy-flow spread below 1e-8 (observed 2e-11); grey-body pin at rel 1e-5.
+ - Test id: `tests/test_energy.py::test_steady_state_ic_balances_interior_and_surface`
+ - Command: `make`, then `pytest -o addopts= tests/test_energy.py::test_steady_state_ic_balances_interior_and_surface` (add `--validation-plots` to write the plot pair)
+ - Companion test: `tests/test_energy.py::test_surface_node_flux_is_the_grey_body_atmospheric_flux` verifies that surface node flux equals the atmospheric grey-body value.
+
+## Related code
+
+- `energy.c`: `set_Etot`, `append_Jcond`, `append_Jconv`, `append_Hradio`, `solve_for_steady_state_energy_interior`.
+- `tests/test_energy.py`: `test_disabled_transport_mechanisms_carry_zero_flux`, `test_radiogenic_heating_pins_two_isotope_decay`, `test_steady_state_ic_balances_interior_and_surface`.
diff --git a/docs/Validation/eos.md b/docs/Validation/eos.md
index 0e75067e..10d64f3d 100644
--- a/docs/Validation/eos.md
+++ b/docs/Validation/eos.md
@@ -1,31 +1,62 @@
# eos.c: EOS interface: setup, viscosity, phase boundaries
-**Source under test**: `eos.c`
+## What this page checks
-## Cases
+This page checks the equation of state interface in `eos.c`, which coordinates material property evaluations.
+`eos.c` manages common option parsing for equations of state, loads phase boundaries, and computes composition-dependent viscosity.
+Validation means that we compare boundary interpolations and viscosity shifts against independent implementations and closed-form formulas.
+If phase boundaries or viscosity calculations are wrong, the melt fraction and convection velocity in the mantle are incorrect.
+The sections below detail two validation cases, their tolerances, and their test commands.
-### Phase boundary matches independent table read
+## Phase boundary matches independent table read
-**Reference-pinned test**: `tests/test_eos.py::test_phase_boundary_matches_independent_table_read`
+**Why this check exists.** Table reading routines can swap solidus and liquidus arrays, misplace table indices, or apply incorrect unit conversions. A wrong phase boundary misidentifies the melting state of the mantle without an error message.
-**Reference**: Cross-implementation cross-check of the loaded melt phase boundary against an independent read of liquidus_A11_H13.dat (Andrault et al. 2011 liquidus with the Hirschmann 2013 solidus companion).
+**What it tests.** The phase boundaries are the liquidus, the entropy above which the mantle is fully molten, and the solidus, below which it is fully solid. The test evaluates the loaded liquidus against an independent read of `liquidus_A11_H13.dat` (Andrault et al. 2011 liquidus with the Hirschmann 2013 solidus companion). It checks interpolation at the P = 0 table edge, a mid-mantle node, and an off-node deep point. The solidus entropy at 10 GPa (about 1524 J/kg/K) sits far from the pinned liquidus entropy (about 2128 J/kg/K), guarding against swapped boundaries. Sibling test `test_viscosity_and_conductivity_dispatch_per_phase` pins per-phase transport dispatch with a 19-decade melt to solid contrast.
-**Tolerance**: rel 1e-12 at the P = 0 table edge, a mid-mantle node, and an off-node deep point.
+**How we know the right answer.** The independent reader evaluates the raw table file using Python numpy interpolation. The C interpolation reproduces these reference values to rel 1e-12.
-**Command**: `make tests_c`, then `pytest -o addopts= tests/test_eos.py::test_phase_boundary_matches_independent_table_read`
+**Why it matters.** The phase boundary determines where melting occurs throughout the mantle. SPIDER evaluates this boundary to compute the local and global melt fraction.
-**Discrimination guards**: The solidus value at 10 GPa (about 1524 J/kg/K) sits far from the pinned liquidus value (about 2128 J/kg/K), so a swapped-boundary regression fails; sign and scale guards bound the entropies; a companion test pins the per-phase viscosity and conductivity dispatch with a 19-decade melt/solid contrast.
+**Known limit.** The independent verification reads the same source table file, checking interpolation rather than physical accuracy of the underlying experimental dataset.
-eos.c owns the option parsing shared by all EOS implementations and loads the phase boundaries. The boundary interpolation is compared against numpy.interp on the same file, and the per-phase transport constants configured in the options file are recovered exactly from the evaluation.
+{ width="80%" }
+{ width="80%" }
-### Compositional viscosity matches Spaargaren et al. (2020)
+**How to read the plot.** The x-axis is the pressure in GPa and the y-axis the entropy of the liquidus in J/kg/K. The line is the liquidus table `liquidus_A11_H13.dat` read with numpy (every 10th node, context). The circles are the asserted values of eos.c at 3 probes. These are the table edge at 0 GPa, a table node at 10 GPa, and an off-node point at 65.5 GPa. The cross is the solidus at 10 GPa, 1524.43 J kg$^{-1}$ K$^{-1}$, which must differ from the code value by more than 100. A pass means that the circles lie on the line (relative 1e-12).
-**Reference-pinned test**: `tests/test_eos.py::test_compositional_viscosity_matches_spaargaren`
+??? info "Test details"
+ - Reference: Cross-implementation cross-check of the loaded melt phase boundary against an independent read of liquidus_A11_H13.dat (Andrault et al. 2011 liquidus with the Hirschmann 2013 solidus companion).
+ - Tolerance: rel 1e-12 at the P = 0 table edge, a mid-mantle node, and an off-node deep point.
+ - Test id: `tests/test_eos.py::test_phase_boundary_matches_independent_table_read`
+ - Command: `make tests_c`, then `pytest -o addopts= tests/test_eos.py::test_phase_boundary_matches_independent_table_read` (add `--validation-plots` to write the plot pair)
+ - Companion test: `tests/test_eos.py::test_viscosity_and_conductivity_dispatch_per_phase` verifies that melt and solid phases receive distinct transport properties.
-**Reference**: The piecewise Mg/Si log10 viscosity prefactor of Spaargaren et al. (2020). The test re-encodes the published coefficients independently of eos.c and additionally pins the two branch-pair shifts as hand-derived literals (0.2074 and -1.35803), so a transcription error in either copy of the formula surfaces.
+## Compositional viscosity matches Spaargaren et al. (2020)
-**Tolerance**: abs 1e-9 on the evaluated shifts (closed-form option arithmetic); abs 1e-4 on the literal transcription pins.
+**Why this check exists.** Piecewise empirical formulas can suffer from transcription typos in branch coefficients, incorrect branch selection, or missing signs. An incorrect viscosity prefactor alters convective vigour by orders of magnitude.
-**Command**: `make tests_c`, then `pytest -o addopts= tests/test_eos.py::test_compositional_viscosity_matches_spaargaren`
+**What it tests.** The test evaluates the piecewise Mg/Si log10 viscosity prefactor of Spaargaren et al. (2020). It re-encodes the coefficients independently of `eos.c`, and pins two branch-pair shifts as hand-derived literals (0.2074 and -1.35803). The comment in `eos.c` says the expressions were worked out in an MSc thesis (2018) and are explained in Spaargaren et al. (2020). The test exercises two configurations that shift in opposite directions across more than 1.5 decades, covering all four piecewise branches. Sibling test `test_activation_terms_shape_the_viscosity_profile` checks the pressure and temperature dependence: no shift at P = 0, a positive shift at depth, its damping by the pressure scale, and the sign of the temperature term.
-**Discrimination guards**: The two configurations shift in opposite directions across more than 1.5 decades, covering all four piecewise branches; a collapsed piecewise or a disabled compositional term fails, and the bare-constant edge case pins the term switched off. The activation-term test (same file) pins the zero-pressure limit exactly and the pressure-scale damping direction.
+**How we know the right answer.** The prefactor is a piecewise closed form in Mg/Si. The shifts match option arithmetic to abs 1e-9, and the hand-derived literal transcription pins match to abs 1e-4.
+
+**Why it matters.** Compositional controls on mantle viscosity determine how mantle convective vigor responds to planetary bulk composition.
+
+**Known limit.** The test re-encodes the coefficients as in `eos.c`. The 2 hand-derived literals catch an error that is in both copies at the 4 compositions they probe; a coefficient that affects neither literal is not caught.
+
+{ width="80%" }
+{ width="80%" }
+
+**How to read the plot.** The x-axis is the Mg/Si ratio and the y-axis the log10 viscosity prefactor. The line is the piecewise prefactor of Spaargaren et al. (2020), as the test re-encodes it (context). The hollow squares are its values at the 2 reference compositions, 1.6 and 0.9. The circles are each reference value plus the asserted shift of eos.c (melt log10 viscosity minus the configured 2.0), at 1.4 and 1.2. A pass means that the circles lie on the line (absolute 1e-9), the one at 1.4 above its square and the one at 1.2 below.
+
+??? info "Test details"
+ - Reference: The piecewise Mg/Si log10 viscosity prefactor of Spaargaren et al. (2020). The test re-encodes the coefficients independently of eos.c and additionally pins the two branch-pair shifts as hand-derived literals (0.2074 and -1.35803), so a transcription error in either copy of the formula surfaces.
+ - Tolerance: abs 1e-9 on the evaluated shifts (closed-form option arithmetic); abs 1e-4 on the literal transcription pins.
+ - Test id: `tests/test_eos.py::test_compositional_viscosity_matches_spaargaren`
+ - Command: `make tests_c`, then `pytest -o addopts= tests/test_eos.py::test_compositional_viscosity_matches_spaargaren` (add `--validation-plots` to write the plot pair)
+ - Companion test: `tests/test_eos.py::test_activation_terms_shape_the_viscosity_profile` checks the pressure and temperature dependence of the viscosity (no shift at P = 0, a positive shift at depth, damping by the pressure scale, and the sign of the temperature term).
+
+## Related code
+
+- `eos.c`: `EOSSetUpFromOptions`, `EOSGetPhaseBoundary`, `EOSEvalSetViscosity`.
+- `tests/test_eos.py`: `test_phase_boundary_matches_independent_table_read`, `test_compositional_viscosity_matches_spaargaren`.
diff --git a/docs/Validation/eos_adamswilliamson.md b/docs/Validation/eos_adamswilliamson.md
index 5b59fbdf..e17b8dde 100644
--- a/docs/Validation/eos_adamswilliamson.md
+++ b/docs/Validation/eos_adamswilliamson.md
@@ -1,19 +1,38 @@
# eos_adamswilliamson.c: Adams-Williamson analytic density
-**Source under test**: `eos_adamswilliamson.c`
+## What this page checks
-## Cases
+This page checks the Adams-Williamson model in `eos_adamswilliamson.c`, which builds planetary structure.
+The Adams-Williamson model describes how mantle density increases with depth.
+Validation means that we check computed profiles against the exact analytical formula.
+If the density profile is incorrect, mantle mass, pressure gradients, and shell volumes are wrong.
+The text below details the check, its tolerance, and its test command.
-### Density matches the closed-form profile
+## Density matches the closed-form profile
-**Reference-pinned test**: `tests/test_eos_adamswilliamson.py::test_density_matches_the_closed_form_profile`
+**Why this check exists.** Analytical codes can suffer from sign errors in compressibility, or from wrong fit inputs. A bug in density calculations distorts the pressure field and mesh geometry without an error.
-**Reference**: Analytical limit: rho(P) = rhos - P beta / g in closed form, with the blackbody50 PREM lower-mantle fit parameters.
+**What it tests.** The model has 3 parameters: the surface density $\rho_s$ (kg/m$^3$), a compressibility parameter $\beta$ (1/m) and the gravity $g$ (m/s$^2$, negative by the SPIDER sign convention). The test evaluates the density $\rho(P) = \rho_s - P\beta/g$ at pressure $P$, with the values of the blackbody50 run (the standard 50-node test run, `tests/opts/blackbody50.opts`), fitted to the lower mantle of the Preliminary Reference Earth Model (PREM). It checks the density at the surface (P = 0), at 10 GPa, and at about 135 GPa, the core-mantle boundary pressure. A sign-flipped $\beta$ gives 3967.8 instead of 4190.1 kg/m$^3$ at 10 GPa; the test asserts a gap above 100 kg/m$^3$. Sibling test `test_si_density_invariant_under_nondimensionalisation_choice` confirms that physical SI density outputs stay unchanged under the nondimensionalisation choice.
-**Tolerance**: rel 1e-12 at the surface (P = 0), 10 GPa, and the approximate core-mantle boundary pressure.
+**How we know the right answer.** The Adams-Williamson formula provides an exact closed-form density profile. Numerical values in C agree with the analytical profile to rel 1e-12, reducing at the surface to the exact configured surface density.
-**Command**: `make tests_c`, then `pytest -o addopts= tests/test_eos_adamswilliamson.py::test_density_matches_the_closed_form_profile`
+**Why it matters.** The Adams-Williamson EOS supplies the static structural profile used by SPIDER native mesh code (`parameters.c` and `mesh.c`).
-**Discrimination guards**: A sign-flipped compressibility shifts the 10 GPa density by more than 200 kg/m3; monotonicity under compression and the silicate density range bound the profile; SI invariance under the nondimensionalisation choice is asserted separately.
+**Known limit.** The Adams-Williamson model assumes an idealized profile, rather than a multi-phase EOS.
-The Adams-Williamson EOS provides the static structural profile used by the native mesh. Its density is a closed-form expression, so the evaluation is pinned to rounding accuracy, with the surface limit reducing exactly to the configured surface density.
+{ width="80%" }
+{ width="80%" }
+
+**How to read the plot.** The x-axis is the pressure in GPa and the y-axis the density in kg/m$^3$. The line is the closed form $\rho_s - P\beta/g$, with surface density $\rho_s$, compressibility parameter $\beta$ and gravity $g$. The dashed line is the form with a sign-flipped $\beta$, which the test rejects (context). The circles are the asserted densities of eos_adamswilliamson.c at the surface, at 10 GPa and at 135 GPa. A pass means that the circles lie on the line (relative 1e-12).
+
+??? info "Test details"
+ - Reference: Analytical limit: rho(P) = rhos - P beta / g in closed form, with the blackbody50 PREM lower-mantle fit parameters.
+ - Tolerance: rel 1e-12 at the surface (P = 0), 10 GPa, and the approximate core-mantle boundary pressure.
+ - Test id: `tests/test_eos_adamswilliamson.py::test_density_matches_the_closed_form_profile`
+ - Command: `make tests_c`, then `pytest -o addopts= tests/test_eos_adamswilliamson.py::test_density_matches_the_closed_form_profile` (add `--validation-plots` to write the plot pair)
+ - Companion test: `tests/test_eos_adamswilliamson.py::test_si_density_invariant_under_nondimensionalisation_choice` asserts SI invariance under the nondimensionalisation choice.
+
+## Related code
+
+- `eos_adamswilliamson.c`: `EOSCreate_AdamsWilliamson`, `EOSEval_AdamsWilliamson`.
+- `tests/test_eos_adamswilliamson.py`: `test_density_matches_the_closed_form_profile`, `test_si_density_invariant_under_nondimensionalisation_choice`.
diff --git a/docs/Validation/eos_composite.md b/docs/Validation/eos_composite.md
index ddf91b9d..f7431e57 100644
--- a/docs/Validation/eos_composite.md
+++ b/docs/Validation/eos_composite.md
@@ -1,19 +1,38 @@
# eos_composite.c: Two-phase composite blending
-**Source under test**: `eos_composite.c`
+## What this page checks
-## Cases
+This page checks the two-phase composite equation of state in `eos_composite.c`, which blends melt and solid properties.
+The composite EOS bridges pure melt and solid lookups across the melting interval.
+Validation means that we test the fusion entropy and melt fraction interpolation against exact analytical definitions.
+If the composite blending is wrong, mantle phase fractions, density jumps, and latent heat release are wrong.
+The section below details the validation check, its tolerance, and its test command.
-### Fusion equals liquidus minus solidus and phi is linear in entropy
+## Fusion equals liquidus minus solidus and phi is linear in entropy
-**Reference-pinned test**: `tests/test_eos_composite.py::test_fusion_equals_liquidus_minus_solidus_and_phi_is_linear_in_entropy`
+**Why this check exists.** Composite blending routines can use wrong denominators, miscalculate melting intervals, or flip boundary bounds. A bug in melt fraction interpolation produces wrong phase mixtures throughout the mantle.
-**Reference**: Analytical limit: by construction of the two-phase region, fusion = S_liq - S_sol and phi = (S - S_sol)/fusion between the phase boundaries.
+**What it tests.** At a given pressure, the mantle is fully solid below the solidus entropy $S_\mathrm{sol}$ and fully molten above the liquidus entropy $S_\mathrm{liq}$ (both in J/kg/K). Between them, the melt fraction $\phi$ runs from 0 to 1, and the entropy of fusion is the width of this interval. The test checks the analytical limit in this two-phase region: fusion = $S_\mathrm{liq} - S_\mathrm{sol}$, and $\phi = (S - S_\mathrm{sol})$/fusion at the specific entropy $S$. It tests the melt fraction at one quarter of the melting interval at 10 GPa. A wrong-denominator regression (using liquidus entropy instead of fusion) shifts phi from 0.25 to about 0.07, which is far beyond tolerance. The fusion entropy lies between 100 and 1500 J/kg/K (603.6 J/kg/K at 10 GPa for the dK09 tables). Sibling tests pin the clamp of the melt fraction between 0 and 1 and the pure-phase limits against single-phase lookups.
-**Tolerance**: rel 1e-10 on the fusion identity; rel 1e-6 on the melt fraction at one quarter of the melting interval at 10 GPa.
+**How we know the right answer.** By mathematical construction of the linear entropy interpolation, the melt fraction and fusion entropy satisfy exact definitions. The fusion identity matches to rel 1e-10, and the fractional melt position matches to rel 1e-6.
-**Command**: `make tests_c`, then `pytest -o addopts= tests/test_eos_composite.py::test_fusion_equals_liquidus_minus_solidus_and_phi_is_linear_in_entropy`
+**Why it matters.** The composite EOS determines the local melt fraction and effective thermodynamic properties of partially molten mantle shells.
-**Discrimination guards**: A wrong-denominator regression (liquidus entropy instead of fusion) shifts phi from 0.25 to about 0.07, far beyond tolerance; the fusion entropy lies between 100 and 1500 J/kg/K (603.6 J/kg/K at 10 GPa for the dK09 tables); companion tests pin the [0, 1] clamp and the pure-phase limits against the single-phase lookups.
+**Known limit.** S_liq, S_sol and the fusion entropy come from the code's own outputs, so the test checks internal consistency: the melt fraction is compared with (S - S_sol)/fusion from those values, not with a pinned 0.25.
-The composite EOS blends the melt and solid lookups across the melting interval. The melt fraction definition, the fusion entropy, and the recovery of the pure end-members away from the interval are each pinned; the truncation contract is made visible through the untruncated diagnostic value.
+{ width="80%" }
+{ width="80%" }
+
+**How to read the plot.** The x-axis is the entropy $S$ at 10 GPa and the y-axis the melt fraction. The dotted lines are the solidus and liquidus entropies $S_\mathrm{sol}$ and $S_\mathrm{liq}$ that eos_composite.c reports. The solid line is $(S - S_\mathrm{sol}) / (S_\mathrm{liq} - S_\mathrm{sol})$ between them (context). The dashed line is the form with the liquidus entropy as the denominator, which the test rejects. The circle is the asserted melt fraction at 1675.34 J kg$^{-1}$ K$^{-1}$, about a quarter of the way through the interval. A pass means that the circle lies on the solid line (relative 1e-6), here at about 0.25.
+
+??? info "Test details"
+ - Reference: Analytical limit: by construction of the two-phase region, fusion = S_liq - S_sol and phi = (S - S_sol)/fusion between the phase boundaries.
+ - Tolerance: rel 1e-10 on the fusion identity; rel 1e-6 on the melt fraction at one quarter of the melting interval at 10 GPa.
+ - Test id: `tests/test_eos_composite.py::test_fusion_equals_liquidus_minus_solidus_and_phi_is_linear_in_entropy`
+ - Command: `make tests_c`, then `pytest -o addopts= tests/test_eos_composite.py::test_fusion_equals_liquidus_minus_solidus_and_phi_is_linear_in_entropy` (add `--validation-plots` to write the plot pair)
+ - Companion test: `tests/test_eos_composite.py::test_phase_fraction_bounded_and_truncation_contract_visible` asserts that the melt fraction is clamped between 0 and 1, and the visibility of untruncated values.
+
+## Related code
+
+- `eos_composite.c`: `EOSCreate_Composite`, `EOSEval_Composite`.
+- `tests/test_eos_composite.py`: `test_fusion_equals_liquidus_minus_solidus_and_phi_is_linear_in_entropy`, `test_phase_fraction_bounded_and_truncation_contract_visible`, `test_pure_phase_limits_recover_the_single_phase_lookups`.
diff --git a/docs/Validation/eos_lookup.md b/docs/Validation/eos_lookup.md
index 2be8ff6e..2bed734f 100644
--- a/docs/Validation/eos_lookup.md
+++ b/docs/Validation/eos_lookup.md
@@ -1,19 +1,38 @@
# eos_lookup.c: Tabulated P-S equation of state
-**Source under test**: `eos_lookup.c`
+## What this page checks
-## Cases
+This page checks the tabulated thermodynamic lookup in `eos_lookup.c`, which reads mineral tables.
+`eos_lookup.c` bilinearly reads thermodynamic properties from pressure-entropy grids.
+Validation means that we check C table reads against an independent Python code.
+If table reads fail, the density and the other tabulated properties of the mantle are wrong; the test compares the density only.
+The section below details the check, its tolerance, and its test command.
-### Melt density matches independent bilinear evaluation
+## Melt density matches independent bilinear evaluation
-**Reference-pinned test**: `tests/test_eos_lookup.py::test_melt_density_matches_independent_bilinear_evaluation`
+**Why this check exists.** Table grid reads can suffer from indexing offsets or wrong corner weights. A bug in table reading shifts the density without an error message; the test checks the density column only, so a wrong scale on another column is not caught.
-**Reference**: Cross-implementation cross-check: an independent numpy bilinear evaluation of the shipped melt density table (1TPa-dK09-elec-free set: de Koker and Stixrude 2009 properties as assembled for Bower et al. 2018).
+**What it tests.** The test evaluates the shipped melt density table (1TPa-dK09-elec-free set: de Koker and Stixrude 2009 properties as assembled for Bower et al. 2018). It compares C outputs with an independent Python bilinear check at three off-node (P, S) probe points spanning 7.3 to 89 GPa. The test guards against reading the wrong table: solid-table density at the same probes differs by more than 1 percent. Positive values and the 2000 to 8000 kg/m3 silicate range bound the density scale. Sibling test `test_si_outputs_invariant_under_nondimensionalisation_choice` confirms consistent SI outputs under scaling choices.
-**Tolerance**: rel 1e-9 at three off-node (P, S) probe points spanning 7.3 to 89 GPa.
+**How we know the right answer.** The table is independently parsed with Python numpy routines, sharing no code with the C side. Matching values pin file parsing, column scaling, index search, and bilinear weighting to rel 1e-9.
-**Command**: `make tests_c`, then `pytest -o addopts= tests/test_eos_lookup.py::test_melt_density_matches_independent_bilinear_evaluation`
+**Why it matters.** Lookup tables supply thermodynamic properties to mantle transport and structure routines.
-**Discrimination guards**: The solid-table density at the same probes differs by more than 1 percent (wrong-table guard); positivity and the 2000 to 8000 kg/m3 silicate range bound the scale; a separate test asserts SI outputs are invariant under the choice of nondimensionalisation.
+**Known limit.** The test compares table reads against an independent Python read of the same table file, but does not benchmark the underlying ab initio mineral physics.
-The lookup EOS is validated against a re-implementation of the same mathematical operation that shares no code with the C side: the table file is parsed and bilinearly interpolated with numpy primitives. Agreement pins file parsing, per-column scaling, index search, and weighting in one comparison.
+{ width="80%" }
+{ width="80%" }
+
+**How to read the plot.** The x-axis is the melt density from an independent bilinear read of `density_melt.dat` in numpy. The y-axis is the density from eos_lookup.c, at the 3 off-node probes from 7.3 to 89 GPa. The circles are the asserted C values, and the line is 1:1. The crosses are the solid-table densities at the same probes, which must differ by more than 1 %. A pass means that the circles lie on the line (relative 1e-9).
+
+??? info "Test details"
+ - Reference: Cross-implementation cross-check: an independent numpy bilinear evaluation of the shipped melt density table (1TPa-dK09-elec-free set: de Koker and Stixrude 2009 properties as assembled for Bower et al. 2018).
+ - Tolerance: rel 1e-9 at three off-node (P, S) probe points spanning 7.3 to 89 GPa.
+ - Test id: `tests/test_eos_lookup.py::test_melt_density_matches_independent_bilinear_evaluation`
+ - Command: `make tests_c`, then `pytest -o addopts= tests/test_eos_lookup.py::test_melt_density_matches_independent_bilinear_evaluation` (add `--validation-plots` to write the plot pair)
+ - Companion test: `tests/test_eos_lookup.py::test_thermodynamic_quantities_positive_and_rho_monotonic_in_pressure` asserts positive thermodynamic quantities and monotonic density with pressure.
+
+## Related code
+
+- `eos_lookup.c`: `EOSCreate_Lookup`, `EOSEval_Lookup`, `SetInterp2dValue`.
+- `tests/test_eos_lookup.py`: `test_melt_density_matches_independent_bilinear_evaluation`, `test_thermodynamic_quantities_positive_and_rho_monotonic_in_pressure`.
diff --git a/docs/Validation/ic.md b/docs/Validation/ic.md
index 163987d0..19dc1f06 100644
--- a/docs/Validation/ic.md
+++ b/docs/Validation/ic.md
@@ -1,43 +1,86 @@
# ic.c: Initial condition
-**Source under test**: `ic.c`
+## What this page checks
-## Cases
+This page checks initial condition setup in `ic.c`, which builds starting mantle and volatile states.
+`ic.c` establishes the initial adiabat profile, sets volatile inventories from mass fractions, and converts ocean-mole specifications.
+Validation means that we compare initial profiles and chemical inventories against analytical limits and mass-balance identities.
+If initial conditions are set incorrectly, simulations start from wrong thermal profiles, or from incorrect volatile inventories.
+The sections below detail three validation cases, their tolerances, and their test commands.
-### Prescribed adiabat entropy sets the initial profile
+## Prescribed adiabat entropy sets the initial profile
-**Reference-pinned test**: `tests/test_ic.py::test_prescribed_adiabat_entropy_sets_the_initial_profile`
+**Why this check exists.** Initial condition (IC) routines can misinterpret input entropy parameters or fail to write values to the staggered nodes (the cell centres of the mesh). A bug in IC setup shifts the starting temperature structure, and the later cooling is wrong too.
-**Reference**: Analytical limit: the initial condition prescribes the adiabat entropy (2600 J/kg/K at the top staggered node) with the configured small gradient ramp.
+**What it tests.** An adiabat is a profile of constant specific entropy. The test evaluates the prescribed adiabat entropy (2600 J/kg/K at the top staggered node, the top cell centre) with the configured gradient ramp `ic_dsdr`. The profile mean lies between the adiabat entropy and that value plus |ic_dsdr| times mantle depth. A run with ic_adiabat_entropy = 2400 shifts the mean by 200 J/kg/K to rel 1e-9, guarding against ignored option inputs. Entropy is positive, with its mean between 1e3 and 1e4 J/kg/K. Top-node entropy matches at t = 0 to rel 1e-3 (observed exact). Sibling test `test_initial_adiabat_is_nearly_flat_until_cooling_imprints_structure` checks that subsequent cooling imprints thermal structure.
-**Tolerance**: rel 1e-3 at the top node (observed exact); the profile mean lies between the adiabat entropy and that value plus |ic_dsdr| times the mantle depth; the 2400 J/kg/K run shifts the mean by 200 J/kg/K to rel 1e-9.
+**How we know the right answer.** Prescribing an adiabat provides an analytical limit, which defines an idealized reference. The top-node value, and the mean profile shift, match exact option arithmetic.
-**Command**: `make`, then `pytest -o addopts= tests/test_ic.py::test_prescribed_adiabat_entropy_sets_the_initial_profile`
+**Why it matters.** The initial adiabat sets the starting thermal state of the magma ocean. SPIDER integrates cooling forward from this baseline.
-**Discrimination guards**: A run with ic_adiabat_entropy = 2400 shifts the profile by exactly 200 J/kg/K (wrong-value discrimination); the entropy is positive and its mean between 1e3 and 1e4 J/kg/K. The sibling test `test_initial_adiabat_is_nearly_flat_until_cooling_imprints_structure` (same file) checks that cooling at later outputs breaks the initial flatness, so the IC is not a fixed point of the output pipeline.
+**Known limit.** An error in the entropy gradient of the initial profile (`ic_dsdr`) is caught only if it moves the profile mean outside 2600 to 2628.7 J/kg/K. The check does not cover non-adiabatic initial profiles.
-The initial condition writes the prescribed adiabat onto the mesh. The top-node entropy, the integrated gradient ramp, and the response to a different prescribed value are pinned from the t = 0 output.
+{ width="80%" }
+{ width="80%" }
-### Abundance IC realises the requested inventory
+**How to read the plot.** The x-axis is the radius and the y-axis the entropy at the staggered nodes (the cell centres of the mesh) at $t$ = 0. The solid line is the asserted profile for `ic_adiabat_entropy` = 2600. The dashed line is the rerun with 2400, whose mean must be 200 J kg$^{-1}$ K$^{-1}$ lower (relative 1e-9). The circle is the asserted top staggered node, 2600 to relative 1e-3. The dotted line is the mean of the 2600 profile, and the grey band its asserted range, 2600 to 2628.7. A pass means that the circle sits at 2600, the dotted line inside the band, and the dashed line 200 below the solid one.
-**Reference-pinned test**: `tests/test_ic.py::test_abundance_ic_realises_the_requested_inventory`
+??? info "Test details"
+ - Reference: Analytical limit: the initial condition prescribes the adiabat entropy (2600 J/kg/K at the top staggered node) with the configured small gradient ramp.
+ - Tolerance: rel 1e-3 at the top node (observed exact); the profile mean lies between the adiabat entropy and that value plus |ic_dsdr| times the mantle depth; the 2400 J/kg/K run shifts the mean by 200 J/kg/K to rel 1e-9.
+ - Test id: `tests/test_ic.py::test_prescribed_adiabat_entropy_sets_the_initial_profile`
+ - Command: `make`, then `pytest -o addopts= tests/test_ic.py::test_prescribed_adiabat_entropy_sets_the_initial_profile` (add `--validation-plots` to write the plot pair)
+ - Companion test: `tests/test_ic.py::test_initial_adiabat_is_nearly_flat_until_cooling_imprints_structure` verifies that cooling breaks initial profile uniformity.
-**Reference**: Mass-balance identity through the initial partial-pressure solve: initial_kg equals the requested ppm of the mantle mass per volatile, and the realised reservoirs carry the requested hydrogen and carbon mole totals (conserved by the water and carbon dioxide reactions).
+## Abundance IC realises the requested inventory
-**Tolerance**: rel 1e-9 on the stored inventories and the H and C totals (observed below 1e-15).
+**Why this check exists.** Volatile partition routines can lose element mass during initial solver iterations, or misapply ppm conversions. A bug in volatile setup gives wrong initial volatile budgets in the atmosphere and the interior.
-**Command**: `make`, then `pytest -o addopts= tests/test_ic.py::test_abundance_ic_realises_the_requested_inventory`
+**What it tests.** The test evaluates mass-balance closure through the initial partial-pressure solve. The quantity initial_kg equals requested ppm of mantle mass per volatile. Realised reservoirs hold requested hydrogen and carbon mole totals. Volatile oxygen drifts by about 2e-3 relative because reactions exchange O with the melt fO2 buffer, and realised H2O reservoirs sit about 3.6e-3 relative off bare requests. Discrimination thresholds sit at 1e-4, far below observed signals, while inert networks match to 1e-9.
-**Discrimination guards**: The volatile oxygen total drifts by about 2e-3 relative against the request because the reactions exchange O with the melt's fO2 buffer, and the realised H2O reservoirs sit about 3.6e-3 relative off the bare request; the test thresholds sit at 1e-4, at least twenty-fold below the observed signals, while an inert reaction network would match the request to the 1e-9 the H and C totals meet.
+**How we know the right answer.** Closed-system mass conservation requires total moles of conserved elements to match requested inputs. Stored inventories and H and C totals match to rel 1e-9 (observed below 1e-15).
-### Ocean-moles IC converts moles to mass
+**Why it matters.** SPIDER's own volatile partitioning starts from this allocation. PROTEUS leaves the volatile part of SPIDER unused and computes outgassing with its own modules.
-**Reference-pinned test**: `tests/test_ic.py::test_ocean_moles_ic_converts_moles_to_mass`
+**Known limit.** Volatile oxygen mass is not strictly conserved by volatile species alone because the melt iron buffer exchanges oxygen with dissolved gases.
-**Reference**: The OCEAN_MOLES constant of constants.c (7.68894973907177e22 mol per Earth ocean): initial_kg equals moles times OCEAN_MOLES times the volatile molar mass.
+{ width="80%" }
+{ width="80%" }
-**Tolerance**: rel 1e-9 (closed-form product of stored constants).
+**How to read the plot.** Left: the grey bars are the requested masses of the 4 volatiles, which the test asserts to equal the stored `initial_kg` (relative 1e-9). The circles are the realised reservoirs after the initial partial-pressure solve, on a log scale. In this run (not asserted), the reactions move H$_2$ and CO$_2$ below their bars and H$_2$O and CO above. Right: the realised over the requested moles of H, C and O, with the dashed line at 1. A pass means that H and C sit on the line (relative 1e-9) and O departs from it by more than 1e-4. Here O lies 2.2e-3 above, because the reactions exchange oxygen with the oxygen-fugacity (fO$_2$) buffer of the melt.
-**Command**: `make`, then `pytest -o addopts= tests/test_ic.py::test_ocean_moles_ic_converts_moles_to_mass`
+??? info "Test details"
+ - Reference: Mass-balance identity through the initial partial-pressure solve: initial_kg equals the requested ppm of the mantle mass per volatile, and the realised reservoirs carry the requested hydrogen and carbon mole totals (conserved by the water and carbon dioxide reactions).
+ - Tolerance: rel 1e-9 on the stored inventories and the H and C totals (observed below 1e-15).
+ - Test id: `tests/test_ic.py::test_abundance_ic_realises_the_requested_inventory`
+ - Command: `make`, then `pytest -o addopts= tests/test_ic.py::test_abundance_ic_realises_the_requested_inventory` (add `--validation-plots` to write the plot pair)
+ - Companion test: `tests/test_ic.py::test_abundance_and_ocean_mole_solves_agree` verifies consistency between abundance and ocean-mole initialization paths.
-**Discrimination guards**: Reading the CO2 inventory with the CO molar mass shifts the expected mass by 36 percent; positivity and scale bounds on the CO2 inventory, the smallest in ocean equivalents (0.0087 oceans; H2 is smaller in mass).
+## Ocean-moles IC converts moles to mass
+
+**Why this check exists.** Unit conversions from ocean equivalents can apply wrong molar masses, or miscalculate moles-to-mass constants. A conversion error scales volatile mass by large factors, leading to invalid atmospheric pressures.
+
+**What it tests.** The test evaluates the OCEAN_MOLES constant of constants.c (7.68894973907177e22 mol per Earth ocean). It verifies that initial_kg equals moles times OCEAN_MOLES times volatile molar mass. Reading CO2 inventory with CO molar mass shifts expected mass by 36 percent, which guards against swapped species indices. The test checks positivity and scale bounds on CO2 inventory, the smallest in ocean equivalents (0.0087 oceans; H2 is smaller in mass).
+
+**How we know the right answer.** Moles to mass conversion is a closed-form product of defined constants. The evaluated mass matches to rel 1e-9.
+
+**Why it matters.** Setting volatile budgets in Earth ocean equivalents allows direct comparison with terrestrial water and carbon inventories in planetary models.
+
+**Known limit.** The test uses its own copy of the OCEAN_MOLES constant in `constants.c`, so it checks the conversion, not the definition of the terrestrial ocean mass; a value that is wrong in both places passes.
+
+{ width="80%" }
+{ width="80%" }
+
+**How to read the plot.** The bars are the closed-form masses of the 4 volatiles on a log scale: moles times `OCEAN_MOLES` (moles per Earth ocean) times molar mass. The circles are the stored `initial_kg` of the run, which the test asserts to match to relative 1e-9. The cross is the CO$_2$ inventory with the CO molar mass, about 36 % lower, which the test rejects. A pass means that each circle sits on the top of its bar and the cross lies below the CO$_2$ bar.
+
+??? info "Test details"
+ - Reference: The OCEAN_MOLES constant of constants.c (7.68894973907177e22 mol per Earth ocean): initial_kg equals moles times OCEAN_MOLES times the volatile molar mass.
+ - Tolerance: rel 1e-9 (closed-form product of stored constants).
+ - Test id: `tests/test_ic.py::test_ocean_moles_ic_converts_moles_to_mass`
+ - Command: `make`, then `pytest -o addopts= tests/test_ic.py::test_ocean_moles_ic_converts_moles_to_mass` (add `--validation-plots` to write the plot pair)
+ - Companion test: `tests/test_ic.py::test_abundance_and_ocean_mole_solves_agree` confirms that the ocean-mole and ppm abundance formulations give the same partial pressures to rel 5e-3.
+
+## Related code
+
+- `ic.c`: `set_ic_interior`, `set_ic_atmosphere`, `set_ic_interior_entropy`, `set_ic_atmosphere_from_initial_total_abundance`, `set_ic_atmosphere_from_ocean_moles`.
+- `tests/test_ic.py`: `test_prescribed_adiabat_entropy_sets_the_initial_profile`, `test_abundance_ic_realises_the_requested_inventory`, `test_ocean_moles_ic_converts_moles_to_mass`.
diff --git a/docs/Validation/index.md b/docs/Validation/index.md
index 0b55a6a3..cb68d804 100644
--- a/docs/Validation/index.md
+++ b/docs/Validation/index.md
@@ -1,6 +1,19 @@
# Validation
-Each physics source of SPIDER has at least one `reference_pinned` test that compares the code with a published benchmark, an analytical limit or an independent implementation. The pages of this tab list these cases, one page per C source file, grouped by the process they describe. Each case gives its reference, its tolerance, its test id and the command that runs it. A full-run regression against the frozen blackbody50 reference complements these cases (`tests/test_regression.py`).
+## What this page checks
+
+This page lists the validation cases of SPIDER: each case checks 1 physics source file against a known answer.
+Validation means that we compare model outputs against analytical limits and independent implementations.
+Each physics source file has at least one pinned test, with a stated tolerance and reference value.
+If these components are wrong, the interior thermal evolution that SPIDER computes, alone or inside [PROTEUS](https://proteus-framework.org/PROTEUS), is wrong. PROTEUS is the coupling framework that runs SPIDER as its interior module.
+The table below lists each validation case, organized by physical process.
+A full-run regression against the frozen blackbody50 reference complements these cases (`tests/test_regression.py`).
+
+Terms used on every page:
+
+- **blackbody50**: the standard 50-node test run of SPIDER, set by `tests/opts/blackbody50.opts`, with a grey-body surface (`SURFACE_BC 1`) and no internal heating.
+- **Basic nodes** and **staggered nodes**: the 2 sets of points of the radial mesh. Basic nodes are the cell boundaries, staggered nodes the cell centres.
+- **IC** and **BC**: initial condition and boundary condition.
Tests in the `unit` tier drive the small C test executables in `tests/c` (`make tests_c`); tests in the `smoke` tier run the `spider` binary (`make`) on short configurations. See [Testing SPIDER](../How-to/test.md) for the tiers and the commands.
diff --git a/docs/Validation/interp.md b/docs/Validation/interp.md
index 21ebbfb4..578f2812 100644
--- a/docs/Validation/interp.md
+++ b/docs/Validation/interp.md
@@ -1,19 +1,38 @@
# interp.c: 1-D and 2-D lookup interpolation
-**Source under test**: `interp.c`
+## What this page checks
-## Cases
+This page checks interpolation routines in `interp.c`, which evaluate tabulated material properties.
+`interp.c` provides one-dimensional and two-dimensional linear table interpolation.
+Validation means that we test interpolator behavior against affine data with exact analytical solutions.
+If table interpolation fails, material properties across the mantle are evaluated incorrectly.
+The section below details the validation check, its tolerance, and its test command.
-### Affine data reproduced exactly including last node
+## Affine data reproduced exactly including last node
-**Reference-pinned test**: `tests/test_interp.py::test_affine_data_reproduced_exactly_including_last_node`
+**Why this check exists.** Table interpolators can suffer from index off-by-one errors, boundary truncation bugs, or incorrect linear weighting. An interpolation bug distorts lookups across all equation of state tables.
-**Reference**: Analytical limit: piecewise-linear interpolation reproduces affine data exactly.
+**What it tests.** The test evaluates piecewise-linear interpolation on synthetic affine data y = 2x + 1 at on-node, off-node, and exact-endpoint queries. Nearest-neighbour values at bracketing nodes (5 and 7 around the pinned 6) sit far outside tolerance. Sign and range guards bound the pinned value. A companion test pins the chord value 6.5 on y = x**2, confirming that higher-order schemes would fail.
-**Tolerance**: rel 1e-14 on y = 2x + 1 at on-node, off-node, and exact-endpoint queries.
+**How we know the right answer.** Piecewise-linear interpolation reproduces affine data exactly by mathematical definition. Numerical queries match analytical values to rel 1e-14 across nodes, interior intervals, and endpoints.
-**Command**: `make tests_c`, then `pytest -o addopts= tests/test_interp.py::test_affine_data_reproduced_exactly_including_last_node`
+**Why it matters.** The interpolation layer feeds the phase boundary in `eos.c`, the tabulated lookups in `eos_lookup.c`, and tables in `bc.c`, `ic.c` and `parameters.c`.
-**Discrimination guards**: Nearest-neighbour values at the bracketing nodes (5 and 7 around the pinned 6) sit far outside the tolerance; sign and range guards bound the pinned value; a companion test pins the chord value 6.5 on y = x**2 so a higher-order scheme would fail.
+**Known limit.** The case covers 1-D interpolation with unit scalings; the 2-D bilinear case (`test_bilinear_plane_exact_everywhere_including_grid_corner`) and per-column scaling are not part of it.
-The interpolation layer feeds every EOS table evaluation, so its contract (linear between nodes, clamped outside the table range, exact at the nodes including the last one) is pinned on synthetic tables whose values are known in closed form. Out-of-range clamping rather than extrapolation is asserted on both sides of the table.
+{ width="80%" }
+{ width="80%" }
+
+**How to read the plot.** The hollow squares are the 11 table nodes of $y = 2x + 1$, and the line is this affine function. The circles are the asserted values of interp.c at the 4 queries $x$ = 0, 2.5, 7 and 10 (the last node). The crosses at $x$ = 2.5 are the 2 nearest node values, 5 and 7. The test requires the code value 6 to be more than 0.5 from each. A pass means that the circles lie on the line (relative 1e-14).
+
+??? info "Test details"
+ - Reference: Analytical limit: piecewise-linear interpolation reproduces affine data exactly.
+ - Tolerance: rel 1e-14 on y = 2x + 1 at on-node, off-node, and exact-endpoint queries.
+ - Test id: `tests/test_interp.py::test_affine_data_reproduced_exactly_including_last_node`
+ - Command: `make tests_c`, then `pytest -o addopts= tests/test_interp.py::test_affine_data_reproduced_exactly_including_last_node` (add `--validation-plots` to write the plot pair)
+ - Companion test: `tests/test_interp.py::test_out_of_range_queries_clamp_to_table_edges` asserts out-of-range clamping at table edges.
+
+## Related code
+
+- `interp.c`: `Interp1dCreateAndSet`, `Interp2dCreateAndSet`, `SetInterp1dValue`, `SetInterp2dValue`.
+- `tests/test_interp.py`: `test_affine_data_reproduced_exactly_including_last_node`, `test_out_of_range_queries_clamp_to_table_edges`, `test_curved_data_follows_the_piecewise_linear_contract`.
diff --git a/docs/Validation/matprop.md b/docs/Validation/matprop.md
index d31335b5..a627f3aa 100644
--- a/docs/Validation/matprop.md
+++ b/docs/Validation/matprop.md
@@ -1,19 +1,38 @@
# matprop.c: Blended material properties
-**Source under test**: `matprop.c`
+## What this page checks
-## Cases
+This page checks material property blending in `matprop.c`, which computes transport coefficients throughout the mantle.
+`matprop.c` blends melt and solid properties over the melting interval to determine local viscosity and conductivities.
+Validation means that we check blended viscosity against configured end-member values in pure phases.
+If material property blending is wrong, mantle convective heat transport and viscosity structures are wrong.
+The text below details the check, its tolerance, and its test command.
-### Viscosity recovers configured end members
+## Viscosity recovers configured end members
-**Reference-pinned test**: `tests/test_matprop.py::test_viscosity_recovers_configured_end_members`
+**Why this check exists.** Rheological blending functions can swap the melt and solid limits. A bug in property blending causes the convection solver to use wrong viscosities.
-**Reference**: Analytical limit: away from the rheological transition the blended viscosity recovers the configured end-members (melt_log10visc = 2, solid_log10visc = 21).
+**What it tests.** The test evaluates the analytical limit where blended viscosity recovers configured end-members away from the rheological transition (melt_log10visc = 2, solid_log10visc = 21). It checks two extreme states: the hot initial state of the blackbody50 run (the standard 50-node test run, `tests/opts/blackbody50.opts`) at t = 0 ( where all 50 nodes are molten), and a cold 0-step run started from a 1600 J/kg/K adiabat (where at least 5 nodes are solid). The recovered end-members must differ by more than 10 decades, guarding against swapped phases. Sibling test `test_material_properties_positive_and_conductivity_constant` asserts positive values for density, heat capacity and thermal expansion, and a constant conductivity, alongside configured 4.0 W/m/K conductivity.
-**Tolerance**: abs 0.01 on log10(viscosity) at the strongly molten nodes and abs 0.05 at the strongly solid nodes.
+**How we know the right answer.** Away from the partial-melting interval, blended material properties must approach pure-phase limits asymptotically. Viscosity matches configured values to abs 0.01 on log10(viscosity) at molten nodes and abs 0.05 at solid nodes.
-**Command**: `make`, then `pytest -o addopts= tests/test_matprop.py::test_viscosity_recovers_configured_end_members`
+**Why it matters.** Mantle viscosity governs convective heat transport. PROTEUS reads the viscosity array `visc_b` from SPIDER.
-**Discrimination guards**: The recovered end-members differ by more than 10 decades (swapped-phase guard); all 50 nodes are molten in the hot state and at least 5 solid in the cold one; the viscosity is positive. The sibling test `test_material_properties_positive_and_conductivity_constant` (same file) asserts the positivity of density, heat capacity, conductivity and expansivity at every node and the configured 4.0 W/m/K conductivity.
+**Known limit.** The main test samples only nearly pure nodes, so a wrong weighting inside the melting interval is not caught there. The companion test `test_viscosity_rises_monotonically_across_rheological_transition` checks that the viscosity rises steadily across the transition, by more than 3 decades, but not its shape.
-Material properties are blended across the melting interval from the per-phase EOS evaluations. The viscosity end-members and the viscosity positivity are pinned on two states: the hot initial state of the short blackbody50 run (t = 0, every node molten) and a cold 0-step run started from a 1600 J/kg/K adiabat (solid nodes); the sibling test `test_material_properties_positive_and_conductivity_constant` checks positivity and the constant conductivity on the full run.
+{ width="80%" }
+{ width="80%" }
+
+**How to read the plot.** The x-axis is the radius and the y-axis the log10 viscosity in Pa s at the basic nodes (the cell boundaries of the mesh). The solid line is the hot initial condition, with all 50 nodes molten. It is asserted at the melt end member, the configured melt viscosity of $10^2$ Pa s (dashed, to 0.01 in log10). The dash-dotted line is the cold initial condition at 1600 J kg$^{-1}$ K$^{-1}$. Its circles are the asserted nodes with a melt fraction below $10^{-3}$, at the solid end member $10^{21}$ Pa s (dotted, to 0.05). A pass means that the solid line lies on the dashed line and the circles on the dotted line.
+
+??? info "Test details"
+ - Reference: Analytical limit: away from the rheological transition the blended viscosity recovers the configured end-members (melt_log10visc = 2, solid_log10visc = 21).
+ - Tolerance: abs 0.01 on log10(viscosity) at the strongly molten nodes and abs 0.05 at the strongly solid nodes.
+ - Test id: `tests/test_matprop.py::test_viscosity_recovers_configured_end_members`
+ - Command: `make`, then `pytest -o addopts= tests/test_matprop.py::test_viscosity_recovers_configured_end_members` (add `--validation-plots` to write the plot pair)
+ - Companion test: `tests/test_matprop.py::test_material_properties_positive_and_conductivity_constant` asserts positive values for density, heat capacity and thermal expansion, and a constant conductivity; the main test asserts a positive viscosity.
+
+## Related code
+
+- `matprop.c`: `set_capacitance_staggered`, `set_matprop_basic`.
+- `tests/test_matprop.py`: `test_viscosity_recovers_configured_end_members`, `test_material_properties_positive_and_conductivity_constant`.
diff --git a/docs/Validation/mesh.md b/docs/Validation/mesh.md
index 622d37e1..f03a8b79 100644
--- a/docs/Validation/mesh.md
+++ b/docs/Validation/mesh.md
@@ -1,19 +1,38 @@
# mesh.c: Static mesh and pressure structure
-**Source under test**: `mesh.c`
+## What this page checks
-## Cases
+This page checks static radial mesh and pressure profile construction in `mesh.c`.
+The mesh module sets up the basic nodes (cell boundaries) and staggered nodes (cell centres) of the radial mesh, shell areas, shell volumes, and gravitational pressures.
+Validation means that we test the native pressure profile against the closed-form Adams-Williamson solution.
+If mesh construction fails, shell geometry and pressure fields in the mantle are distorted.
+The section below details the validation check, its tolerance, and its test command.
-### Pressure profile matches Adams-Williamson closed form
+## Pressure profile matches Adams-Williamson closed form
-**Reference-pinned test**: `tests/test_mesh.py::test_pressure_profile_matches_adams_williamson_closed_form`
+**Why this check exists.** Radial grid generation can suffer from sign errors in the hydrostatic integral or a wrong surface radius. An error in mesh setup makes the pressure and volume arrays of every physics module wrong.
-**Reference**: Analytical limit: the Adams-Williamson pressure profile P(r) = -(rhos g / beta)(exp(beta (R - r)) - 1) in closed form.
+**What it tests.** The profile uses the planet radius $R$ = 6371 km, the radius $r$, the surface density $\rho_s$ (kg/m$^3$), the gravity $g$ (m/s$^2$, negative by the SPIDER sign convention) and a compressibility parameter $\beta$ (1/m). The test evaluates the native Adams-Williamson pressure $P(r) = -(\rho_s g/\beta)(\exp(\beta(R - r)) - 1)$ in closed form. It checks the pressures at the basic nodes through the mantle. A sign-flipped beta and the linear small-beta approximation both fall far outside tolerance, confirming proper discrimination. Sibling test `test_mesh_orderings_interleave_and_shell_masses_close` checks geometric node interleaving and shell mass conservation. Sibling test `test_external_linear_density_mesh_replaces_the_aw_profile` confirms that external meshes replace native profiles cleanly.
-**Tolerance**: rel 1e-12 on the basic-node pressures (observed 6e-15; mesh.c evaluates the same closed form, so the agreement is limited by the round trip through the JSON output).
+**How we know the right answer.** Under the Adams-Williamson assumptions, the hydrostatic equation integrates to an exact analytical expression. Computed pressures agree with the closed form to rel 1e-12 (observed 6e-15, limited by precision in JSON serialization).
-**Command**: `make`, then `pytest -o addopts= tests/test_mesh.py::test_pressure_profile_matches_adams_williamson_closed_form`
+**Why it matters.** The mesh module builds the static spatial and pressure framework that every interior routine in SPIDER evaluates. PROTEUS reads radius arrays `radius_b` and pressure arrays `pressure_s` (in `ReadSPIDER` of `spider.py`).
-**Discrimination guards**: Sign and monotonicity of the pressure with depth; the core-mantle boundary pressure is bounded to its physical order of magnitude; a sign-flipped beta and the linear small-beta approximation both fall far outside tolerance; the external-mesh pathway is verified to actually differ from the native profile when supplied.
+**Known limit.** The test evaluates the expected pressure at the code's own radii, so a spacing error is not caught. A wrong core-mantle boundary radius is caught only if it moves the maximum pressure out of 5e10 to 5e11 Pa. The profile assumes an idealised Adams-Williamson compressibility rather than a self-consistent mineral structure.
-The mesh module builds the static pressure and geometry arrays every physics term runs on. The native profile is pinned against the AW closed form and the geometric consistency between basic and staggered nodes (interleaving, shell masses, monotonicity) is asserted, including for an externally supplied non-AW mesh.
+{ width="80%" }
+{ width="80%" }
+
+**How to read the plot.** The x-axis is the radius in km and the y-axis the pressure in GPa. The line is the Adams-Williamson closed form (context). The dashed and dotted lines are the constant-density profile and the form with a sign-flipped compressibility parameter $\beta$. The test rejects both at the core-mantle boundary. The circles are the asserted pressures of mesh.c at the basic nodes (the cell boundaries of the mesh). A pass means that the circles lie on the line (relative 1e-12).
+
+??? info "Test details"
+ - Reference: Analytical limit: the Adams-Williamson pressure profile P(r) = -(rhos g / beta)(exp(beta (R - r)) - 1) in closed form.
+ - Tolerance: rel 1e-12 on the basic-node pressures (observed 6e-15; `mesh.c` evaluates the same closed form through `eos_adamswilliamson.c`, so the agreement is limited by the round trip through the JSON output).
+ - Test id: `tests/test_mesh.py::test_pressure_profile_matches_adams_williamson_closed_form`
+ - Command: `make`, then `pytest -o addopts= tests/test_mesh.py::test_pressure_profile_matches_adams_williamson_closed_form` (add `--validation-plots` to write the plot pair)
+ - Companion test: `tests/test_mesh.py::test_mesh_orderings_interleave_and_shell_masses_close` verifies node interleaving and shell mass conservation.
+
+## Related code
+
+- `mesh.c`: `set_mesh`.
+- `tests/test_mesh.py`: `test_pressure_profile_matches_adams_williamson_closed_form`, `test_mesh_orderings_interleave_and_shell_masses_close`.
diff --git a/docs/Validation/reaction.md b/docs/Validation/reaction.md
index 7a7c72d8..482b6e7d 100644
--- a/docs/Validation/reaction.md
+++ b/docs/Validation/reaction.md
@@ -1,19 +1,38 @@
# reaction.c: Interior chemical reactions
-**Source under test**: `reaction.c`
+## What this page checks
-## Cases
+This page checks volatile chemical reactions in `reaction.c`, which models volatile speciation.
+`reaction.c` calculates chemical equilibrium and volatile transfers between gas pairs.
+Validation means that we check elemental conservation in reacting pairs against exact stoichiometric balances.
+If reaction bookkeeping is incorrect, volatile elements appear or disappear during mantle degassing.
+The text below details the check, its tolerance, and its test command.
-### Elements conserved across each reaction pair
+## Elements conserved across each reaction pair
-**Reference-pinned test**: `tests/test_reaction.py::test_elements_conserved_across_each_reaction_pair`
+**Why this check exists.** Reaction bookkeeping can violate elemental stoichiometry, misapply molar conversions or swap reactant signs. A bug in reaction bookkeeping changes the total carbon or hydrogen mass during a run.
-**Reference**: Analytical conservation: closed-system elemental budgets across the H2O/H2 and CO2/CO reaction pairs of the Bower et al. (2021) three-ocean hydrogen configuration.
+**What it tests.** The test evaluates closed-system elemental budgets for H2O/H2 and CO2/CO reaction pairs of the Bower et al. (2021) three-ocean hydrogen configuration. It checks hydrogen and carbon transfer balances after one 1000-year step. The test guards reaction signs (H2O produced from H2, CO2 from CO in this step), and bounds H2O transfer between 1e16 and 1e20 kg, below 1 percent of H2O inventory. Sibling tests verify reservoir closure to rel 1e-10 (`test_reservoir_closure_for_every_volatile`), and inventory closure to rel 1e-5 (`test_reaction_transfer_balances_the_initial_inventory`).
-**Tolerance**: rel 1e-3 on the hydrogen and carbon transfer balances after one 1000-year step (observed at round-off level, below 1e-15).
+**How we know the right answer.** Stoichiometric conservation in closed chemical systems requires identical element mass before and after reaction. Balances match to rel 1e-3 (observed below 1e-15).
-**Command**: `make`, then `pytest -o addopts= tests/test_reaction.py::test_elements_conserved_across_each_reaction_pair`
+**Why it matters.** In a SPIDER run with volatiles, the reaction module sets the volatile chemistry of the magma ocean and its atmosphere. PROTEUS does not use the reaction part of SPIDER.
-**Discrimination guards**: Sign structure of the transfers (H2O produced from H2, CO2 from CO in this step); the H2O transfer lies between 1e16 and 1e20 kg and below 1 percent of the H2O inventory. Sibling tests in the same file check the per-volatile reservoir closure to rel 1e-10 (`test_reservoir_closure_for_every_volatile`) and the inventory closure to rel 1e-5 (`test_reaction_transfer_balances_the_initial_inventory`).
+**Known limit.** The check covers carbon and hydrogen; volatile oxygen exchanges with the iron buffer of the melt. The test bounds the direction and the size of each transfer, so a wrong equilibrium constant fails only if it reverses a reaction or moves more than 1 % of the water inventory.
-The reaction module exchanges mass between volatile pairs subject to elemental conservation. The reaction_kg bookkeeping is pinned so that the element handed off by one partner reappears exactly in the stoichiometric fraction of the other.
+{ width="80%" }
+{ width="80%" }
+
+**How to read the plot.** The bars are the asserted element masses moved in the first 1000 yr, on a log scale. Values are shown to 2 significant digits, because the run output moves by about 5e-6 between builds of the PETSc library. For hydrogen, the bars are the hydrogen lost by H$_2$ and the hydrogen in the H$_2$O gained. For carbon, they are the carbon lost by CO and the carbon in the CO$_2$ gained. A pass means that the 2 bars of each pair have the same height; the test asserts relative 1e-3 on the unrounded values.
+
+??? info "Test details"
+ - Reference: Analytical conservation: closed-system elemental budgets across the H2O/H2 and CO2/CO reaction pairs of the Bower et al. (2021) three-ocean hydrogen configuration.
+ - Tolerance: rel 1e-3 on the hydrogen and carbon transfer balances after one 1000-year step (observed at round-off level, below 1e-15).
+ - Test id: `tests/test_reaction.py::test_elements_conserved_across_each_reaction_pair`
+ - Command: `make`, then `pytest -o addopts= tests/test_reaction.py::test_elements_conserved_across_each_reaction_pair` (add `--validation-plots` to write the plot pair)
+ - Companion test: `tests/test_reaction.py::test_reservoir_closure_for_every_volatile` checks per-volatile mass closure between melt and gas reservoirs.
+
+## Related code
+
+- `reaction.c`: `get_log10_modified_equilibrium_constant`, `get_reaction_quotient_products`, `get_reaction_quotient_reactants`.
+- `tests/test_reaction.py`: `test_elements_conserved_across_each_reaction_pair`, `test_reservoir_closure_for_every_volatile`, `test_reaction_transfer_balances_the_initial_inventory`.
diff --git a/docs/Validation/rheologicalfront.md b/docs/Validation/rheologicalfront.md
index 02488e55..263a4799 100644
--- a/docs/Validation/rheologicalfront.md
+++ b/docs/Validation/rheologicalfront.md
@@ -1,19 +1,38 @@
# rheologicalfront.c: Rheological front tracking
-**Source under test**: `rheologicalfront.c`
+## What this page checks
-## Cases
+This page checks front tracking in `rheologicalfront.c`, which locates the solidification boundary.
+`rheologicalfront.c` detects the critical melt fraction threshold where mantle rheology shifts from liquid-like to solid-like behavior.
+Validation means that we check the reported front index and depth against the crossing of the melt-fraction profile that the test finds with the same scan, a self-consistency check.
+If front tracking fails, the dynamic boundary between convective magma ocean and solid mantle is misplaced.
+The section below details the validation check, its tolerance, and its test command.
-### Front location is consistent with melt fraction profile
+## Front location is consistent with melt fraction profile
-**Reference-pinned test**: `tests/test_rheologicalfront.py::test_front_location_is_consistent_with_melt_fraction_profile`
+**Why this check exists.** The front depth must follow from the front radius, and the front must lie inside the mantle. A wrong depth formula, an inverted coordinate sign or a front outside the mantle gives a wrong front without an error message.
-**Reference**: Analytical self-consistency: the reported front index matches the phi_critical = 0.4 crossing of the melt-fraction profile, and the depth equals the surface radius minus the front radius.
+**What it tests.** The test evaluates analytical self-consistency: reported front index matches the phi_critical = 0.4 crossing of the melt-fraction profile, and depth equals surface radius minus front radius. It checks that the front depth is positive and bounded by mantle depth, with front pressure located inside the mantle range. Sibling tests pin the un-formed-front limit at t = 0, where the front reports the core-mantle boundary by convention (`test_unformed_front_reports_the_cmb_at_the_fully_molten_start`), and upward migration during bottom-up solidification (`test_front_rises_from_the_cmb_as_the_mantle_solidifies`).
-**Tolerance**: Front index within one node of the recomputed crossing; rel 1e-6 on the depth identity.
+**How we know the right answer.** The test recomputes the front location in Python, with the same rule, from the melt-fraction profile and the geometry arrays. The depth identity matches to rel 1e-6, and the index matches to one node of the recomputed crossing.
-**Command**: `make`, then `pytest -o addopts= tests/test_rheologicalfront.py::test_front_location_is_consistent_with_melt_fraction_profile`
+**Why it matters.** The front block marks where the mantle turns from mostly molten to mostly solid. PROTEUS reads neither its index nor its depth: it takes `Phi_global` from the same output block and `RF_depth` from the dynamic front.
-**Discrimination guards**: Depth positivity and the mantle-depth bound; front pressure within the mantle range. Sibling tests in the same file pin the un-formed-front limit at t = 0, where the front reports the core-mantle boundary by convention (`test_unformed_front_reports_the_cmb_at_the_fully_molten_start`), and the upward migration of the front as the mantle solidifies bottom-up (`test_front_rises_from_the_cmb_as_the_mantle_solidifies`).
+**Known limit.** The test recomputes the crossing with the same rule as `get_crossover_index` and allows one node of difference, so a one-node misplacement passes; the code takes the first node below `phi_critical` without interpolation. PROTEUS takes `RF_depth` from the dynamic front (`rheological_front_dynamic`, the crossing of the convection regime flag at 1.5, between inviscid (2) and viscous (1) convection), which this test does not cover.
-The rheological front marks the transition between the low- and high-viscosity regimes that PROTEUS reads as the dynamic front depth. Its location is recomputed independently from the melt-fraction profile and the geometry arrays of the same output.
+{ width="80%" }
+{ width="80%" }
+
+**How to read the plot.** The y-axis is the radius, shown to 3 significant digits. The squares are the 2 staggered nodes (cell centres) at 1200 yr on each side of the critical melt fraction $\phi_c = 0.4$ (0.403 and 0.3995). The circle is the asserted front radius, the planet radius minus the reported depth. The test asserts the reported depth and pressure against the arrays at the basic nodes (cell boundaries) to relative 1e-6. The melt fractions are not drawn: the first node is only 4.5e-4 below $\phi_c$, so a rounded value is not stable between builds. A pass means that the front node is at most 1 node from the first node below $\phi_c$; here the circle lies between the squares.
+
+??? info "Test details"
+ - Reference: Analytical self-consistency: the reported front index matches the phi_critical = 0.4 crossing of the melt-fraction profile, and the depth equals the surface radius minus the front radius.
+ - Tolerance: Front index within one node of the recomputed crossing; rel 1e-6 on the depth identity.
+ - Test id: `tests/test_rheologicalfront.py::test_front_location_is_consistent_with_melt_fraction_profile`
+ - Command: `make`, then `pytest -o addopts= tests/test_rheologicalfront.py::test_front_location_is_consistent_with_melt_fraction_profile` (add `--validation-plots` to write the plot pair)
+ - Companion test: `tests/test_rheologicalfront.py::test_front_rises_from_the_cmb_as_the_mantle_solidifies` verifies upward front migration as the mantle solidifies.
+
+## Related code
+
+- `rheologicalfront.c`: `get_crossover_index`.
+- `tests/test_rheologicalfront.py`: `test_front_location_is_consistent_with_melt_fraction_profile`, `test_front_rises_from_the_cmb_as_the_mantle_solidifies`, `test_unformed_front_reports_the_cmb_at_the_fully_molten_start`.
diff --git a/docs/Validation/rhs.md b/docs/Validation/rhs.md
index a6b0ae0f..0f8500b7 100644
--- a/docs/Validation/rhs.md
+++ b/docs/Validation/rhs.md
@@ -1,24 +1,38 @@
# rhs.c: Entropy time-derivative assembly
-**Source under test**: `rhs.c`
+## What this page checks
-## Cases
+This page checks the right-hand side time-derivative assembly in `rhs.c`, which governs thermal evolution.
+`rhs.c` compiles flux divergences and internal heating source terms into the rate of change dS/dt of the specific entropy S (J/kg/K) at each node. This right-hand side is the function that the time stepper integrates.
+Validation means that we check reported time derivatives against finite-difference slopes between consecutive output snapshots.
+If right-hand side evaluation is incorrect, the numerical timestepper integrates along an erroneous trajectory.
+The section below details the validation check, its tolerance, and its test command.
-### Reported dS/dt matches finite difference of outputs
+## Reported dS/dt matches finite difference of outputs
-**Reference-pinned test**: `tests/test_rhs.py::test_reported_dsdt_matches_finite_difference_of_outputs`
+**Why this check exists.** The time derivative that SPIDER reports must be the one the time stepper integrates. A unit or scaling slip in the reported value, or a gross sign flip, makes the reported rate disagree with the trajectory without an error message.
-**Reference**: Analytical self-consistency: the reported dS/dt matches the finite-difference slope of consecutive outputs over 10-year macro steps.
+**What it tests.** A macro step is the interval between 2 outputs; here it is 10 years. The test evaluates the reported dS/dt of the blackbody50 run (the standard 50-node test run, `tests/opts/blackbody50.opts`) against the finite-difference slope of the outputs at 10 and 20 yr. At each node, it asserts cooling (dS/dt < 0). It also bounds the largest |dS/dt| between 1e-12 and 1e-4 J/kg/K/s. The run gives about 2e-8 J/kg/K/s; the same rate, per year instead of per second, is about 0.7 and fails the upper bound. Sibling test `test_cooling_sign_structure_and_finiteness` checks that radiogenic and tidal source arrays remain zero and that entropy declines smoothly in the upper molten layer.
-**Tolerance**: rel 2e-2 at every staggered node (observed at most 0.65 percent, midpoint rule over 10 years).
+**How we know the right answer.** In the limit of small timesteps, the discrete difference between output snapshots approaches the instantaneous time derivative. The test asserts that the reported derivative matches the finite-difference rate to rel 2e-2 at the staggered nodes (cell centres). The observed difference is 0.43 to 0.65 percent, the truncation error of the 10-year slope.
-**Command**: `make`, then `pytest -o addopts= tests/test_rhs.py::test_reported_dsdt_matches_finite_difference_of_outputs` (add `--validation-plots` to write the plot pair)
+**Why it matters.** The right-hand side assembly provides the core time derivative that PETSc integrates to advance mantle thermal evolution in SPIDER.
-{ width="80%" }
-{ width="80%" }
+**Known limit.** The trajectory comes from the same right-hand side that reports dS/dt, so a wrong right-hand side gives a consistent wrong trajectory and passes; dS/dt is not compared with an independent physical value. The run has no heating, so a dropped heating term cannot show. Finite differencing over 10-year intervals adds a truncation error of up to 0.65 percent.
-Top: the cooling rate -dS/dt that rhs.c reports, as the mean of the 10 and 20 yr outputs (line), and the slope of the entropy between these outputs (dashed), at every staggered node of the blackbody50 run. Bottom: their relative difference per node (dots) and the tolerance of the assertion (dashed); the largest difference is 0.65 percent.
+{ width="80%" }
+{ width="80%" }
-**Discrimination guards**: Cooling sign at every node; the largest |dS/dt| lies between 1e-12 and 1e-4 J/kg/K/s, so a per-year value fails (a nondimensional value, about 3e-5 here, passes this bound; the finite-difference comparison catches it). The sibling test in the same file checks the zero radiogenic and tidal source arrays and the near-uniform decline in the upper molten layer.
+**How to read the plot.** The y-axis is the largest cooling rate |dS/dt| of the run in J/kg/K/s, on a log scale. The circle is this asserted value, the mean of the 10 and 20 yr outputs, shown to 1 significant digit. The grey band is the asserted window from 1e-12 to 1e-4 J/kg/K/s. The node profiles are not drawn, because they move by about 7e-4 (relative) between the macOS and Linux builds. The main comparison, with the slope of the specific entropy S between the outputs to relative 2e-2 at every node, is not drawn either. A pass means that the circle lies in the grey band.
-rhs.c assembles the entropy tendency from the flux divergence and sources; it is validated as the actual derivative of the trajectory the timestepper follows, by comparing the reported tendency against the slope between closely spaced outputs.
+??? info "Test details"
+ - Reference: Analytical self-consistency: the reported dS/dt matches the finite-difference slope of consecutive outputs over 10-year macro steps.
+ - Tolerance: rel 2e-2 at every staggered node (observed at most 0.65 percent, midpoint rule over 10 years).
+ - Test id: `tests/test_rhs.py::test_reported_dsdt_matches_finite_difference_of_outputs`
+ - Command: `make`, then `pytest -o addopts= tests/test_rhs.py::test_reported_dsdt_matches_finite_difference_of_outputs` (add `--validation-plots` to write the plot pair)
+ - Companion test: `tests/test_rhs.py::test_cooling_sign_structure_and_finiteness` asserts cooling signs and finiteness of time derivatives.
+
+## Related code
+
+- `rhs.c`: `RHSFunction`.
+- `tests/test_rhs.py`: `test_reported_dsdt_matches_finite_difference_of_outputs`, `test_cooling_sign_structure_and_finiteness`.
diff --git a/docs/Validation/twophase.md b/docs/Validation/twophase.md
index 2bb76ea7..efc6a130 100644
--- a/docs/Validation/twophase.md
+++ b/docs/Validation/twophase.md
@@ -1,19 +1,38 @@
# twophase.c: Global melt and solid mass bookkeeping
-**Source under test**: `twophase.c`
+## What this page checks
-## Cases
+This page checks two-phase mass bookkeeping in `twophase.c`, which tracks global magma ocean freezing.
+`twophase.c` integrates local melt fractions over shell masses to compute total liquid and solid mantle masses.
+Validation means that we check global melt fractions against analytical mass-weighted sums and mass closure.
+If two-phase bookkeeping is wrong, global melt fractions and mantle phase inventories are incorrect.
+The text below details the check, its tolerance, and its test command.
-### Global melt fraction is the mass-weighted profile mean
+## Global melt fraction is the mass-weighted profile mean
-**Reference-pinned test**: `tests/test_twophase.py::test_global_melt_fraction_is_the_mass_weighted_profile_mean`
+**Why this check exists.** Mass sum routines can omit geometric shell weighting, confuse node centers with edges, or violate mass balance. A bug in mass bookkeeping produces inaccurate global melt inventories during cooling.
-**Reference**: Analytical identity: the global melt fraction equals Mliq/Mmantle and the mass-weighted mean of the staggered melt-fraction profile.
+**What it tests.** The test evaluates the analytical identity: the global melt fraction equals the liquid mass Mliq over the mantle mass Mmantle, and the mass-weighted mean of the melt fraction at the staggered nodes (cell centres). It checks a part-solidified 1200-year state. An unweighted profile mean differs by more than 1e-3, confirming that the test guards spatial weighting. Global melt fraction must lie between 0 and 1 and, at 1200 yr, between 0.2 and 0.9. Sibling tests verify mantle mass closure Mliq + Msol = Mmantle at every output to rel 1e-12 (`test_mantle_mass_closure_at_every_output_time`), and monotone solidification under cooling (`test_solidification_proceeds_monotonically_under_cooling`).
-**Tolerance**: rel 1e-10 on both identities at the part-solidified 1200-year state.
+**How we know the right answer.** By mass balance, total melt mass divided by mantle mass must identically equal the mass-weighted mean of the radial melt fraction profile. Both identities match to rel 1e-10.
-**Command**: `make`, then `pytest -o addopts= tests/test_twophase.py::test_global_melt_fraction_is_the_mass_weighted_profile_mean`
+**Why it matters.** Global melt fraction and mantle phase masses are primary outputs read by PROTEUS (`M_mantle_liquid`, `M_mantle_solid`, and `Phi_global` in `ReadSPIDER` of `spider.py`).
-**Discrimination guards**: An unweighted profile mean differs by more than 1e-3 (weighting guard); the global melt fraction lies between 0 and 1 and, at 1200 yr, between 0.2 and 0.9. Sibling tests in the same file check the mantle mass closure Mliq + Msol = Mmantle at every output to rel 1e-12 (`test_mantle_mass_closure_at_every_output_time`) and the monotone solidification under cooling with the mantle mass constant (`test_solidification_proceeds_monotonically_under_cooling`).
+**Known limit.** The identity compares 2 code paths that share the code's own `mass_s` and `phi_s`, so a wrong shell mass or melt fraction passes.
-The two-phase bookkeeping integrates the melt distribution into the global masses PROTEUS consumes. Closure, monotone solidification, and the equivalence of the three melt-fraction definitions are pinned across the full frozen-reference run.
+{ width="80%" }
+{ width="80%" }
+
+**How to read the plot.** The circles are the global melt fraction at 1200 yr in 4 forms. The first 3 are the reported `phi_global`, the liquid over mantle mass Mliq / Mmantle, and the mass-weighted mean of the staggered-node melt fraction. The test requires these 3 to agree (relative 1e-10), and the unweighted mean to be more than 1e-3 away. The grey band is the reported value $\pm$ 1e-3. Values are shown to 2 significant digits, because the run output moves by about 2.5e-6 between builds of the PETSc library. A pass means that the first 3 circles lie in the band and the last outside it; the test checks the unrounded values.
+
+??? info "Test details"
+ - Reference: Analytical identity: the global melt fraction equals Mliq/Mmantle and the mass-weighted mean of the staggered melt-fraction profile.
+ - Tolerance: rel 1e-10 on both identities at the part-solidified 1200-year state.
+ - Test id: `tests/test_twophase.py::test_global_melt_fraction_is_the_mass_weighted_profile_mean`
+ - Command: `make`, then `pytest -o addopts= tests/test_twophase.py::test_global_melt_fraction_is_the_mass_weighted_profile_mean` (add `--validation-plots` to write the plot pair)
+ - Companion test: `tests/test_twophase.py::test_mantle_mass_closure_at_every_output_time` asserts mass conservation closure Mliq + Msol = Mmantle across outputs.
+
+## Related code
+
+- `twophase.c`: `set_Mliq`, `set_Msol`, `set_dMliqdt`, `set_rheological_front`.
+- `tests/test_twophase.py`: `test_global_melt_fraction_is_the_mass_weighted_profile_mean`, `test_mantle_mass_closure_at_every_output_time`, `test_solidification_proceeds_monotonically_under_cooling`.
diff --git a/docs/assets/validation/atmosphere_grey_body_flux.svg b/docs/assets/validation/atmosphere_grey_body_flux.svg
new file mode 100644
index 00000000..533907ab
--- /dev/null
+++ b/docs/assets/validation/atmosphere_grey_body_flux.svg
@@ -0,0 +1,1878 @@
+
+
+
diff --git a/docs/assets/validation/atmosphere_grey_body_flux_dark.svg b/docs/assets/validation/atmosphere_grey_body_flux_dark.svg
new file mode 100644
index 00000000..f6de678d
--- /dev/null
+++ b/docs/assets/validation/atmosphere_grey_body_flux_dark.svg
@@ -0,0 +1,1878 @@
+
+
+
diff --git a/docs/assets/validation/bc_constant_entropy_surface.svg b/docs/assets/validation/bc_constant_entropy_surface.svg
new file mode 100644
index 00000000..99087237
--- /dev/null
+++ b/docs/assets/validation/bc_constant_entropy_surface.svg
@@ -0,0 +1,1628 @@
+
+
+
diff --git a/docs/assets/validation/bc_constant_entropy_surface_dark.svg b/docs/assets/validation/bc_constant_entropy_surface_dark.svg
new file mode 100644
index 00000000..07d5b015
--- /dev/null
+++ b/docs/assets/validation/bc_constant_entropy_surface_dark.svg
@@ -0,0 +1,1628 @@
+
+
+
diff --git a/docs/assets/validation/bc_prescribed_flux.svg b/docs/assets/validation/bc_prescribed_flux.svg
new file mode 100644
index 00000000..f9699798
--- /dev/null
+++ b/docs/assets/validation/bc_prescribed_flux.svg
@@ -0,0 +1,1666 @@
+
+
+
diff --git a/docs/assets/validation/bc_prescribed_flux_dark.svg b/docs/assets/validation/bc_prescribed_flux_dark.svg
new file mode 100644
index 00000000..b3a27d93
--- /dev/null
+++ b/docs/assets/validation/bc_prescribed_flux_dark.svg
@@ -0,0 +1,1666 @@
+
+
+
diff --git a/docs/assets/validation/energy_disabled_transport.svg b/docs/assets/validation/energy_disabled_transport.svg
new file mode 100644
index 00000000..a7cdcd1f
--- /dev/null
+++ b/docs/assets/validation/energy_disabled_transport.svg
@@ -0,0 +1,1456 @@
+
+
+
diff --git a/docs/assets/validation/energy_disabled_transport_dark.svg b/docs/assets/validation/energy_disabled_transport_dark.svg
new file mode 100644
index 00000000..33bc34ce
--- /dev/null
+++ b/docs/assets/validation/energy_disabled_transport_dark.svg
@@ -0,0 +1,1456 @@
+
+
+
diff --git a/docs/assets/validation/energy_radiogenic_decay.svg b/docs/assets/validation/energy_radiogenic_decay.svg
new file mode 100644
index 00000000..39034c84
--- /dev/null
+++ b/docs/assets/validation/energy_radiogenic_decay.svg
@@ -0,0 +1,1771 @@
+
+
+
diff --git a/docs/assets/validation/energy_radiogenic_decay_dark.svg b/docs/assets/validation/energy_radiogenic_decay_dark.svg
new file mode 100644
index 00000000..4fa1f6da
--- /dev/null
+++ b/docs/assets/validation/energy_radiogenic_decay_dark.svg
@@ -0,0 +1,1771 @@
+
+
+
diff --git a/docs/assets/validation/energy_steady_state_ic.svg b/docs/assets/validation/energy_steady_state_ic.svg
new file mode 100644
index 00000000..96dbd91b
--- /dev/null
+++ b/docs/assets/validation/energy_steady_state_ic.svg
@@ -0,0 +1,1898 @@
+
+
+
diff --git a/docs/assets/validation/energy_steady_state_ic_dark.svg b/docs/assets/validation/energy_steady_state_ic_dark.svg
new file mode 100644
index 00000000..03a91c8d
--- /dev/null
+++ b/docs/assets/validation/energy_steady_state_ic_dark.svg
@@ -0,0 +1,1898 @@
+
+
+
diff --git a/docs/assets/validation/eos_adamswilliamson_density.svg b/docs/assets/validation/eos_adamswilliamson_density.svg
new file mode 100644
index 00000000..2a886824
--- /dev/null
+++ b/docs/assets/validation/eos_adamswilliamson_density.svg
@@ -0,0 +1,1578 @@
+
+
+
diff --git a/docs/assets/validation/eos_adamswilliamson_density_dark.svg b/docs/assets/validation/eos_adamswilliamson_density_dark.svg
new file mode 100644
index 00000000..2941f191
--- /dev/null
+++ b/docs/assets/validation/eos_adamswilliamson_density_dark.svg
@@ -0,0 +1,1578 @@
+
+
+
diff --git a/docs/assets/validation/eos_composite_melt_fraction.svg b/docs/assets/validation/eos_composite_melt_fraction.svg
new file mode 100644
index 00000000..96edc95b
--- /dev/null
+++ b/docs/assets/validation/eos_composite_melt_fraction.svg
@@ -0,0 +1,1759 @@
+
+
+
diff --git a/docs/assets/validation/eos_composite_melt_fraction_dark.svg b/docs/assets/validation/eos_composite_melt_fraction_dark.svg
new file mode 100644
index 00000000..d0cc956e
--- /dev/null
+++ b/docs/assets/validation/eos_composite_melt_fraction_dark.svg
@@ -0,0 +1,1759 @@
+
+
+
diff --git a/docs/assets/validation/eos_lookup_melt_density.svg b/docs/assets/validation/eos_lookup_melt_density.svg
new file mode 100644
index 00000000..ae3e13d0
--- /dev/null
+++ b/docs/assets/validation/eos_lookup_melt_density.svg
@@ -0,0 +1,1284 @@
+
+
+
diff --git a/docs/assets/validation/eos_lookup_melt_density_dark.svg b/docs/assets/validation/eos_lookup_melt_density_dark.svg
new file mode 100644
index 00000000..ec381c17
--- /dev/null
+++ b/docs/assets/validation/eos_lookup_melt_density_dark.svg
@@ -0,0 +1,1284 @@
+
+
+
diff --git a/docs/assets/validation/eos_phase_boundary.svg b/docs/assets/validation/eos_phase_boundary.svg
new file mode 100644
index 00000000..d3c47af2
--- /dev/null
+++ b/docs/assets/validation/eos_phase_boundary.svg
@@ -0,0 +1,1595 @@
+
+
+
diff --git a/docs/assets/validation/eos_phase_boundary_dark.svg b/docs/assets/validation/eos_phase_boundary_dark.svg
new file mode 100644
index 00000000..4b53baf4
--- /dev/null
+++ b/docs/assets/validation/eos_phase_boundary_dark.svg
@@ -0,0 +1,1595 @@
+
+
+
diff --git a/docs/assets/validation/eos_spaargaren_viscosity.svg b/docs/assets/validation/eos_spaargaren_viscosity.svg
new file mode 100644
index 00000000..4caebd9e
--- /dev/null
+++ b/docs/assets/validation/eos_spaargaren_viscosity.svg
@@ -0,0 +1,1284 @@
+
+
+
diff --git a/docs/assets/validation/eos_spaargaren_viscosity_dark.svg b/docs/assets/validation/eos_spaargaren_viscosity_dark.svg
new file mode 100644
index 00000000..1fccb564
--- /dev/null
+++ b/docs/assets/validation/eos_spaargaren_viscosity_dark.svg
@@ -0,0 +1,1284 @@
+
+
+
diff --git a/docs/assets/validation/ic_abundance_inventory.svg b/docs/assets/validation/ic_abundance_inventory.svg
new file mode 100644
index 00000000..5d4e6585
--- /dev/null
+++ b/docs/assets/validation/ic_abundance_inventory.svg
@@ -0,0 +1,1488 @@
+
+
+
diff --git a/docs/assets/validation/ic_abundance_inventory_dark.svg b/docs/assets/validation/ic_abundance_inventory_dark.svg
new file mode 100644
index 00000000..fd97d075
--- /dev/null
+++ b/docs/assets/validation/ic_abundance_inventory_dark.svg
@@ -0,0 +1,1488 @@
+
+
+
diff --git a/docs/assets/validation/ic_adiabat_profile.svg b/docs/assets/validation/ic_adiabat_profile.svg
new file mode 100644
index 00000000..23d93efc
--- /dev/null
+++ b/docs/assets/validation/ic_adiabat_profile.svg
@@ -0,0 +1,1536 @@
+
+
+
diff --git a/docs/assets/validation/ic_adiabat_profile_dark.svg b/docs/assets/validation/ic_adiabat_profile_dark.svg
new file mode 100644
index 00000000..c53f89db
--- /dev/null
+++ b/docs/assets/validation/ic_adiabat_profile_dark.svg
@@ -0,0 +1,1536 @@
+
+
+
diff --git a/docs/assets/validation/ic_ocean_moles.svg b/docs/assets/validation/ic_ocean_moles.svg
new file mode 100644
index 00000000..b9fb0f10
--- /dev/null
+++ b/docs/assets/validation/ic_ocean_moles.svg
@@ -0,0 +1,1409 @@
+
+
+
diff --git a/docs/assets/validation/ic_ocean_moles_dark.svg b/docs/assets/validation/ic_ocean_moles_dark.svg
new file mode 100644
index 00000000..28a99cbf
--- /dev/null
+++ b/docs/assets/validation/ic_ocean_moles_dark.svg
@@ -0,0 +1,1409 @@
+
+
+
diff --git a/docs/assets/validation/interp_affine.svg b/docs/assets/validation/interp_affine.svg
new file mode 100644
index 00000000..b08969f0
--- /dev/null
+++ b/docs/assets/validation/interp_affine.svg
@@ -0,0 +1,1262 @@
+
+
+
diff --git a/docs/assets/validation/interp_affine_dark.svg b/docs/assets/validation/interp_affine_dark.svg
new file mode 100644
index 00000000..a460591a
--- /dev/null
+++ b/docs/assets/validation/interp_affine_dark.svg
@@ -0,0 +1,1262 @@
+
+
+
diff --git a/docs/assets/validation/matprop_viscosity_end_members.svg b/docs/assets/validation/matprop_viscosity_end_members.svg
new file mode 100644
index 00000000..4d516056
--- /dev/null
+++ b/docs/assets/validation/matprop_viscosity_end_members.svg
@@ -0,0 +1,1602 @@
+
+
+
diff --git a/docs/assets/validation/matprop_viscosity_end_members_dark.svg b/docs/assets/validation/matprop_viscosity_end_members_dark.svg
new file mode 100644
index 00000000..ca4866f9
--- /dev/null
+++ b/docs/assets/validation/matprop_viscosity_end_members_dark.svg
@@ -0,0 +1,1602 @@
+
+
+
diff --git a/docs/assets/validation/mesh_pressure_profile.svg b/docs/assets/validation/mesh_pressure_profile.svg
new file mode 100644
index 00000000..cf9bfe7e
--- /dev/null
+++ b/docs/assets/validation/mesh_pressure_profile.svg
@@ -0,0 +1,1632 @@
+
+
+
diff --git a/docs/assets/validation/mesh_pressure_profile_dark.svg b/docs/assets/validation/mesh_pressure_profile_dark.svg
new file mode 100644
index 00000000..0997202e
--- /dev/null
+++ b/docs/assets/validation/mesh_pressure_profile_dark.svg
@@ -0,0 +1,1632 @@
+
+
+
diff --git a/docs/assets/validation/reaction_element_balance.svg b/docs/assets/validation/reaction_element_balance.svg
new file mode 100644
index 00000000..19123406
--- /dev/null
+++ b/docs/assets/validation/reaction_element_balance.svg
@@ -0,0 +1,1337 @@
+
+
+
diff --git a/docs/assets/validation/reaction_element_balance_dark.svg b/docs/assets/validation/reaction_element_balance_dark.svg
new file mode 100644
index 00000000..ac02f579
--- /dev/null
+++ b/docs/assets/validation/reaction_element_balance_dark.svg
@@ -0,0 +1,1337 @@
+
+
+
diff --git a/docs/assets/validation/rheologicalfront_location.svg b/docs/assets/validation/rheologicalfront_location.svg
new file mode 100644
index 00000000..9fcc254a
--- /dev/null
+++ b/docs/assets/validation/rheologicalfront_location.svg
@@ -0,0 +1,1279 @@
+
+
+
diff --git a/docs/assets/validation/rheologicalfront_location_dark.svg b/docs/assets/validation/rheologicalfront_location_dark.svg
new file mode 100644
index 00000000..af440346
--- /dev/null
+++ b/docs/assets/validation/rheologicalfront_location_dark.svg
@@ -0,0 +1,1279 @@
+
+
+
diff --git a/docs/assets/validation/rhs_dsdt_finite_difference.svg b/docs/assets/validation/rhs_dsdt_finite_difference.svg
index 7cad4ec9..60d8731b 100644
--- a/docs/assets/validation/rhs_dsdt_finite_difference.svg
+++ b/docs/assets/validation/rhs_dsdt_finite_difference.svg
@@ -1,12 +1,12 @@
-