Semi-implicit Solver for the Heat Equation with Stefan–Boltzmann Law Boundary Condition

The surface energy balance on an atmosphereless body consists of solar irradiance, subsurface heat conduction, and thermal radiation to space by the Stefan–Boltzmann law. Here we extend the semi-implicit Crank–Nicolson method to this specific nonlinear boundary condition and validate its accuracy. A rapid change in incoming solar flux can cause a numerical instability, and several approaches to dampen this instability are analyzed. A predictor based on the Volterra integral equation formulation for the heat equation is also derived and can be used to improve accuracy and stability. The publicly available implementation provides a fast and robust thermophysical model that has been applied to lunar, Martian, and asteroidal surfaces, on occasion to millions of surface facets or parameter combinations.