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.


Introduction: 1D Thermal Model for Planetary Surfaces
The one-dimensional heat equation is routinely solved to calculate the temperature in the shallow subsurface on airless bodies (the Moon, asteroids) or planetary bodies with tenuous atmospheres (Mars).Finite-difference schemes for the heat equation with an implicit or semi-implicit time step are unconditionally stable for linear boundary conditions, when the surface temperature or energy input is a prescribed function of time (Crank & Nicolson 1947;Press et al. 1992).However, the Stefan-Boltzmann radiation law is a nonlinear boundary condition because it involves the surface temperature itself, and it is unclear whether these classical results are still valid.This study is dedicated to the description and analysis of a fast numerical solver for this important application in planetary science.
The governing equations considered here are where Q(t) is the incoming solar irradiance and T(z, t) is the temperature.Variables are listed in Table 1.The Stefan-Boltzmann law, òσT 4 , makes this initial value problem nonlinear, as solutions of the homogeneous equation can no longer be superposed.Two variants of the surface energy balance will also be considered in this work.The thermal properties k, ρ, and c might be functions of z and t.Temperature gradients normal to the surface are often orders of magnitude larger than temperature gradients parallel to the surface, so the one-dimensional heat equation suffices for a great many applications.
The content described in Section 2 was developed by the authors in the time period 2001-2003.We have since used this model for numerous studies, some involving millions of model runs where a 1D heat equation is solved for each pixel (e.g., Aharonson & Schorghofer 2006;Schorghofer et al. 2019;Schorghofer & Williams 2020).Since then, other semi-implicit solvers with a Stefan-Boltzmann law boundary condition have been implemented and briefly described in the literature (e.g., Pohl 2014;Young 2017;Magri et al. 2018), but many 1D thermal models currently used for planetary surfaces use simple explicit time steps (e.g., Spencer et al. 1989;Rozitis & Green 2011;Kieffer 2013) and are not competitive from a computational perspective.There has been a lack of analysis and development for this important application; compare, for example, with the large volume of work on numerical methods to solve Kepler's equation (Taff & Brennan 1989) or the dedicated techniques invented for long-term trajectory integrations (symplectic integrators; Wisdom & Holman 1991), two other standard numerical problems in planetary science.The present work is an overdue effort to advance numerical treatment of the heat equation with radiative cooling on the surface.The desired duration of integration can be long (many solar days over many orbits for many locations), and changes in surface temperature can be rapid on an airless body (at sunrise), challenging the speed and stability of the numerical technique.Section 3 provides validations and benchmarks for our method.In Section 4 we analyze the numerical stability of the method and several approaches to eliminate instabilities.Recently, we developed a Volterra integral equation approach, which is described in Section 5.

Semi-implicit Scheme for Inhomogeneous Medium on
Irregular Grid

Flux-conservative Discretization and Semi-implicit Scheme
Consider grid points at depths z 1 ,K, z N in a direction normal to the surface, with z 1 the first point below the surface.The heat flux is H = −k∂T/∂z.A flux-conservative discretization (e.g., Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
LeVeque 2007) on an irregularly spaced grid is given by ¶ Subscript j refers to position z j .The spatial discretization of the heat Equation (1b) then becomes Introducing the coefficients where superscript n refers to the time step.Hence, a a g g a a g g which leads to the system of equations a a g g a a g g for 1 < j < N.This tridiagonal linear system can be solved in O(N) steps.Moreover, the system is diagonally dominant and can therefore be solved robustly.In the spatial discretization (3), the temperature T j is defined on grid point z j , whereas the conductivity is defined in between points.In Equation (4), the volumetric heat capacity (ρc) j is defined on z j .This is inconvenient because the thermal properties k and ρc, which are use-supplied variables, are evaluated at different depths, offset by half a grid point.In our program implementation, the values for the volumetric heat capacity are shifted by half a grid point, so the thermal properties k and ρc are defined on the same points.When the value of ρc at z j is needed, it is calculated from the average of the two neighboring values.
Although the derivation was made with time-constant thermal parameters k and ρc, it remains applicable if these parameters change slowly with time, because the additional time discretization error will be small.An example is the contribution of radiative heat transfer to the conductivity, which is proportional to T 3 .The time step should be chosen small enough so that the thermal parameters do not change dramatically within a single time step.

Upper Boundary Condition: Prescribed Surface Temperature
Take as a boundary condition a prescribed surface temperature T s = T(0, t).The general formulae (5) and (9) with T 0 = T s and z 0 = 0 yield This is a standard Crank-Nicolson solver for an irregularly spaced grid and inhomogeneous thermal properties.

Upper Boundary Condition: Stefan-Boltzmann Radiation Law
An implicit-in-time discretization of the full problem with nonlinear heat transfer results in a nonlinear algebraic system of equations that could be solved using, e.g., Newton's method.The computational expense of solving such a system at each time step would, however, be large, especially if the properties of the medium are changing in time.Instead, to retain the advantages of an implicit scheme, we linearize the nonlinear boundary condition (1a).
Considering that the flux is evaluated in between grid points (see Equation (3)), it is natural to impose the surface energy balance halfway between z 1 and an imaginary point z 0 = − z 1 , with the surface at z = 0, as illustrated in Figure 1.Introduce the auxiliary quantity T 0 , such that the surface temperature T s = (T 0 + T 1 )/2.The first two grid points are chosen as z 1 = Δz/2 and z 2 = z 1 + Δz = 3z 1 .The coefficients for the upper boundary condition will be derived for this choice of z 2 /z 1 .No restrictions are imposed on z 3 and beyond.
On the surface, The temperature is linearized around a reference temperature This results in the following expression for T 0 : The coefficients (5a) and (5b) involve the imaginary grid point and then α 1 = βk 3/2 and γ 1 = βk 1/2 .The surface temperature is computed as As the reference temperature choose the surface temperature at the current time, = T T r s n .Other choices of T r will be described in Sections 4 and 5.

Upper Boundary Condition with Seasonal Frost
Cover (Mars) On Mars, atmospheric CO 2 can condense (desublimate) seasonally.The surface energy balance (1a) with the latent heat of CO 2 sublimation added is where where H s = − k(∂T/∂z)| z=0 .In addition, adjust the surface albedo and infrared emissivity to account for the difference between bare ground and a frost-covered surface.
(Integration with constant time steps may temporarily result in a slightly negative value for m CO 2 .Energy is conserved more accurately in the long term when the value is allowed to remain negative, although it should be output as zero.)

Lower Boundary Condition
At the lower boundary of the domain, a known value is imposed on the heat flux (1c): (If the z-axis points downward, then an upward heat flux is negative, and the value of H geotherm is positive.)Assume that the position of the hypothetical next grid point is set by Figure 1.Illustration of grid spacing, as used for the derivation with a nonlinear boundary condition.Temperatures T j are defined at the centers of cells, whereas the heat flux H j+1/2 is evaluated in between grid points.The surface is placed halfway from the shallowest grid point because the surface energy balance is imposed on the flux.
In system (9) the equation for j = N becomes g g g g r Alternatively, the geothermal flux could be imposed at depth z N+1/2 .In this case Equation (25) is replaced by with the same definition for γ N .Our current implementation uses Equation (25), but both choices work.

Miscellaneous
Literature sources for the values of thermal conductivity and heat capacity appropriate for planetary surfaces are provided in Appendix A. General thoughts on equilibration times for thermal model calculations are presented in Appendix B.

Analytical Solution for Sinusoidal Surface Temperature
For a sinusoidally varying surface temperature with period P, the solution to the heat equation for a homogeneous material in an infinite half-space is known analytically: where d k p = P is the thermal skin depth.Figure 2 shows a comparison of the analytical solution with the numerical solution.

Asymptotic Solution for Sudden Radiative Cooling
A short-term solution for the heat equation with a Stefan-Boltzmann radiation surface boundary condition was derived by Handelsman & Olmstead (1972).Their nondimensional equations are , where T e is an ambient temperature.In this case, their solution is ( ) for small t, a result also derived in Chambré (1959).After redimensionalizing, the surface temperature is found to change as where T 0 = T(0, 0) is the initial surface temperature.Figure 3 shows that the numerical solver reproduces the expected behavior for this discontinuous change in incoming flux.

Convergence Test
Convergence of some solutions with Δt and Δz has been verified.In the absence of an analytically known solution, the error can be defined as the difference between one numerical solution and another numerical solution at twice the resolution.The difference between the two numerical solutions is defined by a vector norm.Figure 4 demonstrates that the error decreases proportionally to Δt 2 , as expected for a semi-implicit scheme (Crank & Nicolson 1947;Press et al. 1992).An explicit or fully implicit time discretization converges only to first order in Δt.

Test of Energy Conservation
For periodic solutions the heat flux H = −k∂T/∂z, timeaveraged over one period, must be the same at all depths and equal to the heat flux imposed at the bottom boundary.Consider the time average of Equation (1b) over one period.After the solution has equilibrated (has become periodic), and as long as the heat capacity does not vary with time, the time average of H must be constant with depth, even if the thermal properties vary with depth.Figure 5 shows a flux conservation test.The heat flux is conserved even across an interface with a sudden change in thermal properties, something the finite-  difference scheme was not designed for but copes with excellently.

Implementation and Performance
Figure 6 shows the surface temperature over one solar day (sol) for an unobstructed flat surface.The temperature varies smoothly, with a sudden increase at sunrise. Figure 6(b) shows the solution using a time step Δt and 1/10 of that time step; the solutions look almost identical.Plotting the difference between the two solutions (Figure 6(c)) starts to reveal deviations.At sunrise and sunset the solver overestimates the surface temperature.Also shown in Figure 6(c) is a solution using a different reference temperature for the linearization, based on an estimate from the Volterra integral equation, which will be described in detail in Section 5.This reduces the errors to about half.
The thermal model was originally implemented in Modern FORTRAN, and that is also the version benchmarked here.Several implementations are available on GitHub (Schörghofer 2024), along with validation examples.Insight into the relative computational cost of various components is gained by counting, to leading order, the number of floating point operations (FLOPs) based on the expressions in the program code.The tridiagonal solver (Press et al. 1992) takes about 7N FLOPs, where N is the number of grid points.To assign the matrix coefficients and calculate the right-hand side of Equation (9), we count about 22N additional FLOPs in our program implementation.These are necessary because the grid spacing is irregular and the thermal properties are allowed to be inhomogeneous.Our implementation of the upper boundary law (Section 2.3) consumes about 23 FLOPs per time step, which is the equivalent of less than one extra grid point.
Run times were measured on a workstation with a 3.6 GHz Intel Xeon CPU with six cores (but the executable runs on a single thread), running Linux with Gnu FORTRAN compiler version 11.4.0,compiler optimization flag -O2.Our benchmark scenario is for Mars over 6686 sols (one Mars year has 668.60 solar days, so this spans 10 Mars years), with 100 steps per sol and N = 80 grid points.This corresponds to approximately 15-minute time steps over 18.8 Earth years.Total execution time is 1.2 s, or 2 μs per time step.This is measured with the Linux time utility, so it includes start-up cost, trigonometry to calculate the incoming sunlight, and a simple model for atmospheric absorption on Mars, but profiling reveals that these subtasks consume only a few percent of the execution time.This example demonstrates that millions of time steps can be performed in seconds.Taking advantage of (trivial) parallelization on a multicore processor, thermophysical models can be run for millions of pixels over millions of time steps within a few days of compute time.
Profiling revealed that 44% of the time was spent in the tridiagonal solver and an additional 52% in the routines that calculate the coefficients at every time step that serve as input to the tridiagonal solver.Another version was written where the coefficients are precomputed, possible if the thermal properties are time independent.Precomputing the coefficients saves on arithmetic but requires more data movement.This improved execution speed by about one-third.
The C version equals the speed of the FORTRAN version, whereas our current implementations in Matlab and Python are more than an order of magnitude slower.

Numerical Stability
Although the Crank-Nicolson method is unconditionally stable for both a temperature boundary condition, T(0, t) = T s (t), and a linear flux boundary condition, T z (0, t) = Q(t), this is not guaranteed for a nonlinear flux boundary condition like the Stefan-Boltzmann law, kT z (0, t) = òσT 4 (0, t) − Q(t).We have observed a transient instability in cases when the input flux changes abruptly within a single time step, such as a sunrise above an elevated horizon.The above linearization of the Stefan-Boltzmann law (Section 2.3) works well as long as the surface temperature changes slowly, but it fails for large discontinuous changes.
Figure 7 shows results for a horizon elevation of 20°.The discontinuous change in incoming flux causes a major overshoot at the first time step after sunrise.The numerical solution  returns to stable behavior after a few time steps, but the instability is clearly unacceptable.Three approaches to cure the instability will be described in this section.
But first a note about the stability of the explicit scheme.Milton & Goss (1973) derived a stability criterion for radiant heat loss, uniform thermal properties, and uniform grid spacing.After translating their result for the nondimensional equations into dimensional form and into our notation, their result is This is not a rigorous criterion, but they found it to be sufficient, and even generous, for all of their numerical examples.This stability requirement is even worse than for an explicit scheme with a linear boundary condition, where   2).Here sunrise is delayed by a 20°h igh horizon, leading to a discontinuity in the incoming solar flux.All three methods cure the instability, but unphysical oscillations of small magnitude are still discernible, for all three methods.Parameters are as in Figure 6, but with a lower thermal inertia of Γ = 100 J m −2 K −1 s −1/2 and an elevated horizon.
Δt (Δz) 2 /(2κ).We are not aware of any equivalent result for an implicit scheme with radiant heat loss.Williams & Curry (1977) experimented with four methods of incorporating radiant heat loss into a fully implicit solver.A linearization, similar to the one above (Section 2.3), can result in instabilities.An iterative approach, where the tridiagonal linear system is solved repeatedly while a modified Newton method is applied to the surface temperature, turned out to be robust, but it is computationally expensive.They also studied a semi-iterative method, where an iteration is performed in between the two stages of the tridiagonal solver, an interesting approach we have not pursued for the semi-implicit method.
Our approaches to cure the instability are summarized in Table 2.

Flux Smoothing (Reduced Time Step)
One approach to enhance stability is "artificial flux smoothing," where the time step is subdivided into multiple substeps, using linear interpolation of the incoming flux from Q n to Q n+1 .It turns a discontinuous change in Q into a continuous change.This approach does not identify a nascent instability or answer how many substeps are required to cure it, but it reduces the linearization error.
Figure 7 shows the outcome of this approach.When the surface temperature changed by more than 20% (which occurred only once over the last sol shown in the figure), the time step was internally reduced by a factor of five.This eliminates the instability almost entirely.Subsequent steps use the full original time step, because changes in temperature are below the empirically set threshold.
The computational cost is considerable, because the tridiagonal system is re-solved, here five times, if and only if an instability is suspected.For almost the same computational cost, one could reduce the time step from the outset and use fluxes Q for these fractional time steps.However, that approach would not eliminate a discontinuous change in Q, and hence an "artificial" smoothing of Q is preferred.
The procedure could be applied recursively, but recursive functions come with a performance penalty, and we apply flux smoothing at only one level, e.g., with a step 5 times smaller than the original time step.(The number of substeps could be adjusted according to the perceived magnitude of the instability.)

Iterative Predictor-corrector for Reference Temperature
In Section 2.3, the thermal emission is linearized around the reference temperature = T T r s n .If + T s n 1 is far from T s n , a significant error is incurred in the evaluation of the emitted energy.This can be alleviated by repeating the calculation with a new reference temperature T r somewhere in between T s n and . An empirical choice is the geometric mean of the previous reference temperature and the new surface temperature.This predictor-corrector step is iteratively applied until T r is within 20% of + T s n 1 . In our experience, the iterations always converge.
Figure 7 includes a solution with this predictor-corrector, and it successfully dampens the instability.The amplitude of the oscillations is larger than for the flux smoothing method described above, but the computational overhead of this method is smaller than for the flux smoothing method, as one or two iterations often suffice.
The methods described in Sections 4.1 and 4.2 both require that a criterion for a suspected instability is monitored and the temperature profile at time step n is stored instead of being overwritten in place by the new profile at time step n + 1.

Volterra Predictor for Reference Temperature
Another approach to enhance stability is based on the Volterra integral equation and will be described next.Its result when used to choose an improved linearization temperature T r is already included in Figure 7.The outcome is similar to the iterative predictor-corrector approach but requires no additional solves of the tridiagonal system.It is computationally by far the fastest approach among the three stabilizers.Moreover, the Volterra predictor also improves accuracy at all time steps (Figure 6(c)).

Volterra Integral Equation Approach
The initial value problem of a partial differential equation can sometimes be formulated as an integral equation.The onedimensional heat equation can be transformed into a Volterra integral equation.We take advantage of this approach to develop a predictor step for the surface temperature.

Solution to Volterra Integral Equation
Cannon (1984) provides a treatise for the one-dimensional heat equation.Theorem 7.3.1 considers the problem for t > 0. If the functions f and h are piecewise continuous and  is continuous, the solution is of the form Moreover, if  is Lipschitz continuous, i.e., its derivative with respect to t does not exceed a fixed value, then u is unique and f(t) = u(0, t); see also Jumarhon et al. (1996).
On the surface (x = 0), This integral equation for the surface temperature u(0, t) provides an alternative approach to the initial value problem for all times, as long as thermal properties are uniform.We will derive an approximate solution to this integral equation to estimate the surface temperature after a small time step.For this purpose, thermal properties need to be approximately constant only over the shallow depth that heat propagates within a single time step.

Approximation of Solution
To approximate Equations (36) and (37), we only retain terms where the exponent can vanish, because otherwise the exponential factors tend to be small.This leaves only the terms with n = 0 in the sums, solutions with second-order time accuracy and only O(N) FLOPs per time step, where N is the number of grid points.We verified the accuracy of this finite-difference solver and have used it extensively in various planetary science applications.
Linearization of the upper boundary condition can cause instabilities when the input irradiance changes rapidly, as can occur when the Sun rises above an elevated horizon.Three approaches to eliminate the instability are analyzed: artificial smoothing of the input flux using fractional time steps, a predictor-corrector that iteratively adjusts the reference temperature used for the linearization, and a predictor temperature based on an approximate solution of the Volterra integral equation that also adjusts the reference temperature.All three approaches dampen the instability dramatically, with the artificial flux smoothing providing the most robust approach.What remains to be discovered is a numerical method that is guaranteed to be stable for large discontinuous changes in incoming flux.
The Volterra integral equation approach can be used to predict the surface temperature at the next time step.We derived a computationally light Volterra predictor, which is beneficial to both accuracy and stability.

Appendix B Spin-up
Many applications of this thermal model involve periodic solutions set by the orbital period around the Sun.The thicker the domain, the longer it will take for the model to equilibrate, but too thin a domain will yield inaccurate results.A domain depth of 4-5 times the skin depth of the longest period suffices for accuracy but can require, say, 20 periods for thermal equilibration.
An accurate initial guess for the mean temperature is most helpful.When the thermal parameters are constant in time, the time-averaged temperature profile can be calculated from the time-averaged surface temperature.Without geothermal heat, the time-averaged temperature is in fact constant with depth.A suitable initialization temperature can be based on the timeaveraged absorbed flux.A finite temperature amplitude will lead to more radiative cooling than calculated based on this time-averaged flux, but the deviation is large only for small thermal inertia.
The required spin-up time can be further shortened by resetting the temperature profile after the first few orbits.Based on the provisional mean surface temperature T 0 0 , the mean temperature profile is obtained with This profile is then used for re-initialization.This can also be applied "on the fly" when the comparison is made with the temperature time-averaged over the previous period but applied to an instantaneous temperature profile.At very low temperatures, radiative heat loss (σT 4 ) becomes small compared to the heat flux within the thermal skin depth, and it takes a long time to reach equilibrium.In this situation, the surface temperature amplitude is small for the same reason, and the absorbed flux can be time-averaged over one period and used for initialization, a method that has already been pointed out above.

Figure 2 .
Figure 2. Comparison of the numerical solution with the analytical solution for Crank-Nicolson solver with a periodic surface boundary condition.The deviations at the bottom are justified because the analytical solution (28) is for an infinite half-space.Nonequidistant grid points were used in this example.

Figure 3 .
Figure 3. Response of numerical solution to a sudden change in incoming flux compared to the analytically obtained expansion for small times, Equation (29).

Figure 4 .
Figure 4. Convergence with time step Δt for the Crank-Nicolson method with a nonlinear boundary condition.Errors are evaluated as ||T Δt (z, t) − T Δt/2 (z, t)||, where the subscript indicates the time step.The 1-norm is the average absolute deviation between the two numerical solutions, and the max-norm is the maximum absolute deviation, -D DT z t T z t max , ,j t j t j 2 | ( ) ( )|.The black lines have slopes 1 and 2, respectively, and arbitrary prefactors.The rate of convergence is second order.

Figure 5 .
Figure 5. Validation of the conservation of the heat flux (flux-conservative discretization).The ice table at 10 cm depth causes dramatic changes in thermal properties.Left panel: minimum and maximum subsurface temperatures over one Mars year.Middle panel: temperature averaged over one Mars year, which changes linearly within each layer because the thermal conductivity is constant within each layer.Right panel: heat flux averaged over one Mars year, which is preserved across changes in thermal properties and is equal to the heat flux imposed at the bottom boundary of 0.028 W m −2 .

Figure 7 .
Figure 7.Comparison of stabilizers (Table2).Here sunrise is delayed by a 20°h igh horizon, leading to a discontinuity in the incoming solar flux.All three methods cure the instability, but unphysical oscillations of small magnitude are still discernible, for all three methods.Parameters are as in Figure6, but with a lower thermal inertia of Γ = 100 J m −2 K −1 s −1/2 and an elevated horizon.

Table 1
Frequently Used Variables CO 2 is the specific latent heat and m CO 2 is the areal mass density of CO 2 ice.Time integration has to switch between the boundary conditions described in Sections 2.2 and 2.3.Use the boundary condition of Section 2.3 when the surface temperature T s is above the CO 2 frost point temperature or if to the CO 2 frost point temperature based on the atmospheric pressure and calculate the energy difference to update m CO 2 .Specifically, Note.The last column lists the number of additional tridiagonal solves carried out if an instability is anticipated.