MultIHeaTS: A Fast and Stable Thermal Solver for Multilayered Planetary Surfaces

Abstract A fully implicit scheme is proposed for solving the heat equation in 1D heterogeneous media, available as a computationally efficient open-source Python code. The algorithm uses finite differences on an irregular grid and is unconditionally stable due to the implicit formulation. The thermal solver is validated against a stiff analytical solution, demonstrating its robustness in handling stiff initial conditions. Its general applicability for heterogeneous cases is demonstrated through its use in a planetary surface scenario with nonlinear boundary conditions induced by blackbody thermal emission. MultIHeaTS's advantageous stability allows for computation times up to 100 times faster than Spencer’s explicit solver, making it ideal for simulating processes on large timescales. This solver is used to compare the thermal signatures of homogeneous and bilayer profiles on Europa. Results show that homogeneous materials cannot reproduce the thermal signature observed in bilayer profiles, emphasizing the need for multilayer solvers. In order to optimize the scientific return of a space mission, we propose a strategy made of three local time observations that is enough to identify bilayer media, for instance, for the next missions to the Jovian system. A second application of the solver is the estimation of the temperature profile of Europa’s near surface (first 10s m) throughout a 1 million yr simulation with varying orbital parameters. The probability distribution of temperature through depth is obtained. Among its various applications, MultIHeaTS serves as the core thermal solver in a multiphysics simulation model detailed in the companion article by C. Mergny & F. Schmidt.