heat1d

One-dimensional thermal model for planetary science applications

heat1d solves the 1-D heat equation in planetary regolith with depth- and temperature-dependent thermophysical properties, following the formulation of Hayne et al. (2017). Implementations are provided in Python, C, and MATLAB.

Features

• Four numerical solvers: Explicit (Forward Euler), Implicit (Backward Euler), Crank-Nicolson, and Fourier-matrix (frequency domain)
• Temperature-dependent properties: thermal conductivity with radiative $T^3$ term, polynomial heat capacity
• Depth-dependent profiles: exponential density/conductivity transition with custom layer support
• Angle-dependent albedo following Keihm (1984)
• Non-uniform spatial grid with geometrically increasing layer spacing
• Graphical user interface (PySide6) with interactive parameter panel, run manager, and 5 plot types
• JPL Horizons/SPICE integration for precise ephemeris-driven illumination
• Eclipse modeling for satellite bodies (e.g., Moon in Earth's shadow)
• PSR crater modeling: bowl-shaped permanently shadowed regions (Ingersoll & Svitek, 1992)
• Sloped surfaces: tilted-surface insolation with self-shadowing (Braun & Mitchell, 1983) and indirect heating from surrounding flat terrain
• YAML configuration system shared between Python and C implementations
• Validation suite against Apollo heat flow data and Diviner radiometer observations
• C implementation for the highest computational performance (same YAML config, all solvers, 24 validation tests)

Quick Start

Installation

The heat1d package on PyPI is an obsolete release that is no longer updated, so don't use pip install heat1d. Install the current version directly from GitHub (requires Python ≥ 3.9). The Python package lives in the python/ subdirectory of this repository:

pip install "heat1d @ git+https://github.com/phayne/heat1d.git#subdirectory=python"

For the graphical user interface (PySide6):

pip install "heat1d[gui] @ git+https://github.com/phayne/heat1d.git#subdirectory=python"

For development, clone the repository and do an editable install:

git clone https://github.com/phayne/heat1d.git
cd heat1d/python
pip install -e ".[dev,gui]"

Command Line

# Default: Moon equator, 1 diurnal cycle
heat1d

# Moon at 45 N with implicit solver, 2 output days
heat1d --lat 45 --solver implicit --ndays 2

# Load a YAML config file, override latitude
heat1d moon.yaml --lat 30

# Use JPL Horizons for real ephemeris
heat1d --use-spice --lat 0 --lon 0 --start-time "2024-06-15 12:00" --ndays 1

# Run the validation suite
heat1d --validate

# Launch the GUI
heat1d-gui

Python API

import numpy as np
from heat1d import Model, Configurator, planets

config = Configurator(solver="crank-nicolson")
model = Model(planet=planets.Moon, lat=np.deg2rad(30), ndays=1, config=config)
model.run()

# Surface temperature vs local time
print(f"Peak T: {model.T[:, 0].max():.1f} K")
print(f"Min  T: {model.T[:, 0].min():.1f} K")

# Access depth grid and full T(time, depth) array
print(f"Grid: {model.N_z} layers to {model.profile.z[-1]:.2f} m depth")

Supported Planetary Bodies

The following bodies have built-in orbital and thermophysical parameters and can be selected by name with --planet (CLI), planet: name: (YAML), or planets.<Name> (Python):

Body	Type	Surface properties	Horizons ID
Moon (default)	Satellite of Earth	Lunar regolith (Hayne et al., 2017)	301
Mercury	Planet	Lunar regolith defaults	199
Europa	Galilean satellite	Porous ice regolith	502
Ganymede	Galilean satellite	Porous ice regolith	503
Small bodies, e.g. Bennu, Patroclus	Asteroids, incl. near-Earth, main-belt, and Jupiter Trojan (Lucy mission targets)	Lunar regolith defaults	e.g. 2101955 (Bennu)

The other built-in small bodies are defined in python/heat1d/planets.py.

heat1d --planet Europa --lat 0

model = Model(planet=planets.Bennu, lat=0, ndays=1)

--lon (or lon in Python, in radians) is always east longitude. Without an ephemeris, longitudes are referenced so that 0° E faces the Sun at perihelion, which is the IAU definition for Mercury. Mercury's 3:2 spin-orbit resonance then puts its hot poles (~700 K) at 0° and 180° E and its warm poles (~570 K) at 90° and 270° E:

heat1d --planet Mercury --lat 0 --lon 90    # warm pole

To choose the orbital position (season) directly, use --true-anomaly (degrees, 0 = perihelion; true_anomaly in Python, in radians), which places local noon at that point in the orbit. It cannot be combined with --lon or --use-spice:

heat1d --planet Bennu --true-anomaly 180    # noon at aphelion

Any property can be overridden in the planet: section of a YAML config file: albedo, emissivity, ks, kd, rhos, rhod, H, cp0, cpCoeff, Qb, S, rAU, day, year, eccentricity, obliquity, and Lp. Changing rAU rescales the solar flux S by the inverse square of distance unless S is also given. This also makes bodies whose entries carry only orbital data (Pluto, Triton, etc.) usable by supplying their thermophysical properties; running one without them stops with an error listing what is missing.

Other airless bodies can be modeled with planet: name: Custom, which starts from lunar regolith properties, combined with --use-spice --body-id <ID> to fetch the body's orbit, rotation period, radius, and albedo from JPL Horizons:

# ceres.yaml
planet:
  name: "Custom"

heat1d ceres.yaml --use-spice --body-id 2000001 --lat 0 --lon 0 \
    --start-time "2025-01-01 00:00" --ndays 1

Example Output

Temperature vs. Depth and Local Time

The thermal wave penetrates only a few centimeters into lunar regolith. Below the thermal skin depth (~7 cm, accounting for temperature-dependent conductivity and heat capacity), diurnal variations are strongly damped.

Multi-Latitude Comparison

Peak daytime temperatures decrease with latitude due to lower solar incidence angles, while nighttime temperatures converge as all surfaces radiate to the same cold sky.

Solvers

heat1d implements four numerical methods for the 1-D heat equation. Three are finite-difference time-stepping schemes; the fourth solves in the frequency domain:

Solver	Method	Stability	Time Accuracy	Typical Steps/Day	Relative Speed
Explicit	Forward Euler	CFL-limited	$O(\Delta t)$	~830	1x
Implicit	Backward Euler + TDMA	Unconditional	$O(\Delta t)$	~24	~35x
Crank-Nicolson	Semi-implicit + TDMA	Unconditional	$O(\Delta t^2)$	~24	~35x
Fourier-matrix	Frequency domain	N/A (periodic)	Spectral	N/A	~1000x

Time-Stepping Solvers

The explicit scheme (Hayne et al. 2017, Eq. A17) is a straightforward forward Euler discretization. It updates each grid node from its neighbors at the current time step, requiring no linear algebra. However, the Courant-Friedrichs-Lewy (CFL) stability condition ($\Delta t \le \frac{\Delta z^2}{2\kappa}$) limits the time step to ~3000 seconds for typical lunar parameters, resulting in ~830 steps per synodic day.

The implicit and Crank-Nicolson schemes evaluate spatial derivatives at the new time level (fully or partially), producing a tridiagonal linear system at each step. This system is solved in O(N) time using the Thomas algorithm (TDMA). Because they are unconditionally stable, both schemes can take time steps ~35x larger than the explicit scheme. The Crank-Nicolson scheme averages explicit and implicit contributions, achieving second-order accuracy in time at the same computational cost as the first-order implicit scheme.

Both the implicit and Crank-Nicolson solvers support adaptive time-stepping via Richardson extrapolation (step doubling), automatically adjusting the time step to maintain a user-specified temperature accuracy.

Fourier-Matrix Solver

The Fourier-matrix solver eliminates time-stepping entirely by solving the periodic steady-state in the frequency domain. It decomposes the diurnal surface flux into Fourier harmonics and propagates each frequency through the subsurface using 2×2 transmission matrices. Nonlinear surface radiation is handled via Newton iteration in the time domain, using a circulant admittance matrix constructed from the frequency-domain impedance.

An outer iteration loop captures the solid-state greenhouse effect (thermal pumping): the nonlinear $T^3$ dependence of thermal conductivity produces a net downward heat flux that elevates subsurface temperatures. The solver computes the exact time-averaged rectification flux and adjusts the equilibrium temperature profile accordingly.

This approach is ~1000× faster than time-stepping because it solves the periodic steady state directly, without equilibration orbits. A complete lunar diurnal cycle is solved in ~100 ms. It is also the default equilibration solver -- even when using time-stepping methods for output, the Fourier solver initializes the temperature profile to the converged periodic state, eliminating multi-orbit spin-up.

The Fourier-matrix approach builds on the harmonic decomposition of periodic lunar thermal models introduced by Linsky (1966), the thermal quadrupole (transfer matrix) formalism of Pipes (1957) and Maillet et al. (2000), and the harmonic balance technique for nonlinear periodic circuits (Kundert & Sangiovanni-Vincentelli, 1986). The solid-state greenhouse (thermal pumping) correction follows the analysis of Linsky (1966), Section V. For detailed equations and derivations, see the Numerical Methods documentation.

PSR Crater Modeling

heat1d can model temperatures on the floors of permanently shadowed regions (PSRs) in bowl-shaped polar craters, following the analytical framework of Ingersoll & Svitek (1992). The crater floor receives no direct sunlight but is heated by scattered sunlight and thermal infrared re-emitted by the illuminated crater walls.

The crater geometry is parameterized by the depth-to-diameter ratio $d/D$, which determines the surface area ratio $f = 4(d/D)^2 / [1 + 4(d/D)^2]$ and the crater half-angle $\beta = \arccos(1 - 2f)$. A PSR exists when the maximum solar elevation angle at the given latitude is less than $\beta$, setting a minimum $d/D$ that increases with distance from the pole. For the Moon ($\epsilon \approx 1.54$°), even shallow craters ($d/D \approx 0.01$) can host PSRs at the poles, while at 70° latitude only the deepest simple craters ($d/D \gtrsim 0.1$) qualify.

The model accounts for cavity radiation trapping (enhanced effective emissivity) and computes the absorbed flux from wall scattering. PSR mode works with all time-stepping solvers and adaptive timestepping, but not with the Fourier-matrix solver (which automatically falls back to implicit).

# PSR crater floor at 85°N with d/D = 0.2
heat1d --lat 85 --psr-d-D 0.2 --ndays 1

model = Model(planet=planets.Moon, lat=np.deg2rad(85), ndays=1, psr_d_D=0.2)
model.run()

For the full derivation and formulas, see the Boundary Conditions documentation.

Sloped Surfaces

heat1d can model temperatures on an isolated tilted surface surrounded by flat terrain, given the slope angle and azimuth (degrees clockwise from north: 0=N, 90=E). Direct insolation follows the Braun & Mitchell (1983) geometry with two independent occlusion tests: the direct beam is zero whenever the sun is behind the tilted surface (self-shadowing) or below the local horizon — the tangent plane of the spherical body — so sunrise/sunset timing accounts for occlusion by the body itself, even for steep slopes that would otherwise geometrically "see" the sun at night.

By default, sloped runs also include ground heating — the indirect flux from thermal emission and reflected sunlight of the surrounding flat terrain, seen by the slope at view factor $\sin^2(s/2)$. This is computed from an automatic flat-surface companion run at the same latitude, whose equilibrated diurnal cycle supplies a periodic indirect-flux table applied during both equilibration and output. The feature works with the Python and C backends, all solvers, and the JPL Horizons illumination mode (which uses the queried solar azimuth/elevation).

# 25° south-facing slope at 30°N, with ground heating (default)
heat1d --lat 30 --slope 25 --slope-az 180

# Direct insolation only
heat1d --lat 30 --slope 25 --slope-az 180 --no-ground-heating

model = Model(planet=planets.Moon, lat=np.deg2rad(30), ndays=1,
              slope=np.deg2rad(25), slope_az=np.deg2rad(180))
model.run()

See the Sloped Surfaces documentation for the geometry, view factors, and approximations.

JPL Horizons Ephemeris (SPICE)

With --use-spice, heat1d drives the surface boundary condition with the actual solar geometry at a specific site and time instead of an idealized diurnal cycle. It queries the JPL Horizons system for the Sun's azimuth, elevation, and distance as seen from the given latitude and longitude. This captures the real solar distance, seasonal declination, and (for satellites) eclipses by the parent body. A network connection is required.

# Apollo 15 landing site (26.13 N, 3.63 E), one lunar day starting 2024-06-15 12:00 UTC
heat1d --use-spice --lat 26.13 --lon 3.63 --start-time "2024-06-15 12:00" --ndays 1

# Equator at the sub-Earth point through the total lunar eclipse of 2025-03-14,
# using an explicit stop time and a fine output step to resolve the eclipse
heat1d --use-spice --lat 0 --lon 0 \
    --start-time "2025-03-13 00:00" --stop-time "2025-03-15 00:00" --output-dt 0.05

In the second example the Moon passes through Earth's shadow near local noon, and the surface temperature drops from ~380 K to ~175 K before recovering.

Notes:

• Times are UTC and must be given as "YYYY-MM-DD HH:MM".
• Set the run length with --ndays (in diurnal cycles) or --stop-time.
• The model is first equilibrated with the built-in Fourier solver; the Horizons ephemeris drives the output period only.
• Because the run starts at an arbitrary UTC epoch, the local_time_hr output column begins at that epoch's local time rather than at noon (0 = local noon, 12 = local midnight).
• Eclipses are detected automatically for satellites; disable with --no-eclipses, or set the parent body with --parent-body-id.
• --body-id overrides the Horizons body ID (e.g. 301 for the Moon, 2101955 for Bennu).

See the CLI Reference for all options.

Graphical User Interface

The GUI provides interactive control over all model parameters with real-time visualization:

• Parameter panel with tabs for body selection, simulation settings, material properties, and numerical options
• Run manager for tracking and comparing multiple simulation runs
• Five plot types: diurnal curves, depth profiles, heatmaps, flux, and combined views
• Publication-quality export (PDF, PNG, SVG) and CSV data export
• Parameter sweeps across any model variable (linear or log spacing)
• Custom depth profiles via the layer editor for non-standard regolith structures
• JPL Horizons search for finding any solar system body by name

Install with the [gui] extra (see Installation) and launch with heat1d-gui.

C Implementation

The c/ directory contains a standalone C implementation optimized for batch runs and high performance. It supports all four solvers, reads the same YAML configuration files as the Python version, and can model any airless planetary body (Moon, Mercury, asteroids, icy satellites, etc.) by configuring the physical and orbital parameters.

Building

cd c/
brew install fftw libyaml   # macOS; see c/README.md for Linux packages
make
make test                    # 24 validation checks

Running

Legacy positional arguments (backward compatible):

# Moon equator, TI=55, H=0.06, highland albedo
./heat1d 0 55 0.06 0.12 > temperature.txt

# Apollo 15 site (26 N, mare albedo, Fourier solver)
./heat1d 26 55 0.06 0.06 3 > apollo15.txt

# Implicit solver with adaptive timestepping
./heat1d 0 55 0.06 0.12 2 480 480 1.0 > implicit_adaptive.txt

YAML configuration (shares config files with Python):

# Use the Python example config
./heat1d --config ../python/heat1d/examples/moon_default.yaml --ti 55

# Override latitude and albedo for an Apollo site
./heat1d --config ../python/heat1d/examples/moon_default.yaml \
    --lat 26 --ti 55 --albedo 0.06

# Inspect the full configuration
./heat1d --config moon.yaml --verbose 2>config_dump.txt

Output is temperature data to stdout (one row per time step, one column per depth layer), plus loctime.txt and profile_z_dz_rho_k.txt in the working directory. See c/README.md for complete documentation.

Theory

heat1d solves the 1-D heat equation in porous planetary regolith:

$$\rho c_p \frac{\partial T}{\partial t} = \frac{\partial}{\partial z} \left( K \frac{\partial T}{\partial z} \right)$$

with the following depth- and temperature-dependent material properties from Hayne et al. (2017):

• Density: exponential transition from surface (1100 kg/m³) to depth (1800 kg/m³) with scale height $H$
• Thermal conductivity: phonon (contact) conductivity plus a radiative $T^3$ component: $K = K_c \left[1 + \chi (T/350)^3 \right]$
• Heat capacity: 4th-order polynomial in temperature following Hemingway et al. (1981)
• Surface boundary: radiative equilibrium with angle-dependent albedo $A(\theta) = A_0 + a(\theta/45{^\circ})^3 + b(\theta/90{^\circ})^8$
• Bottom boundary: constant geothermal heat flux (0.018 W/m² for the Moon)

For full derivations, see the Theory documentation and Hayne et al. (2017).

Documentation

Detailed documentation is organized by topic:

Topic	Description
Theory	Heat equation physics and flux formulation
Thermophysical Properties	Density, conductivity, and heat capacity models
Boundary Conditions	Surface energy balance and bottom heat flux
Numerical Methods	Solver equations, stability, and comparison
Spatial Grid	Non-uniform grid construction and skin depth
Initialization	Temperature profile initialization
Equilibration	Convergence to periodic steady state
Validation	Comparison with Apollo and Diviner data
CLI Reference	Command-line options and examples
API Reference	Python class and function documentation
C Implementation	Building, running, and configuring the C code
MATLAB Implementation	MEX-file version for MATLAB

Validation

The model is validated against lunar temperature measurements from the Apollo 15 and 17 heat flow experiments and the Diviner Lunar Radiometer Experiment. All 8 validation checks pass within the published uncertainties across all four solvers.

heat1d --validate

See the Validation documentation for details.

Citation

If you use heat1d in your research, please cite:

Hayne, P. O., Bandfield, J. L., Siegler, M. A., Vasavada, A. R., Ghent, R. R., Williams, J.-P., Greenhagen, B. T., Aharonson, O., Elder, C. M., Lucey, P. G., & Paige, D. A. (2017). Global regolith thermophysical properties of the Moon from the Diviner Lunar Radiometer Experiment. Journal of Geophysical Research: Planets, 122, 2371-2400. doi:10.1002/2017JE005387

License

MIT License. See LICENSE for details.

Authors

• Paul O. Hayne - University of Colorado Boulder (paul.hayne@lasp.colorado.edu)
• K.-Michael Aye - Package maintainer