Earth
Geothermal gradient
Heat a 40 km column of crust from below at 30 mW/m² while its own radioactivity adds 45 more: the profile curves, steep in the sediments and flatter in the conductive lower crust, and after a few slab time constants of 23 Myr it lands on the closed-form steady state to a hundredth of a degree while the surface heat flow reaches 75 mW/m². Halve the sediments’ conductivity and the whole crust below runs 180 K hotter; double the basal flow and the Moho climbs past 1300 degrees. A uniform half-space after a surface step follows the erf solution. Heat in, made and out balance the stored heat to rounding.
What happens
One-dimensional conduction through a layered crust: eighty finite volumes of half a kilometre with the conductivity and heat production of sediments (2 W/m/K, 1 µW/m³), upper crust (2.5, 2) and lower crust (2.2, 0.5), harmonic face conductivities, the surface held at 10 °C and the base fed 30 mW/m². Explicit stepping at 2 kyr books every face flux and the production, so the stored heat meets what came in, was made and left to rounding; the steady profile is the piecewise-quadratic closed form, and a uniform half-space after a surface step is the erf solution. One second of the experiment is one Myr.
Try this
Watch the profile bend into its steady shape and read the Moho temperature against the closed form; halve the cover’s conductivity, double the basal flow, or start the whole column hot and follow erf.
The model and its limits
Pure conduction with constant properties per layer, no erosion, advection or time-dependent heat production, a fixed basal heat flow. Temperatures are those of this column.
Equations used in the module
Equations ρc·∂T/∂t = ∂/∂z(k·∂T/∂z) + A, q = −k·∂T/∂z k_face = 2k_a·k_b/(k_a + k_b), q(z) = q_b + ∫_z^L A·dz′ T(z) = T_s + ∫_0^z q/k·dz′, q_surface = q_b + Σ A·Δz T(z, t) = T_s + (T_0 − T_s)·erf(z/2√(κt)), τ = 4L²/π²κ
