A lightweight thermal modelling tool for physics-based continental heat flow interpolation

Both energy applications, such as assessing one of the controlling factors of conductive geothermal plays, and geodynamics modelling, are influenced by the large uncertainties arising from uneven sampling of the direct observable of the Earth's thermal state, surface heat flow. Heterogeneity in structure and composition of the continental lithosphere complicate the temperature field even in stable provinces in thermal equilibrium. The measurements deviate from what simple relationships with geological and geophysical data predict, requiring more sophisticated schemes such as those based on multivariate inversion (e.g. Mather et al. 2018) and geostatistics (e.g. the similarity method employed by Lucazeau, 2019).

Recently, we aimed at assessing the performance of satellite-gravity-constrained modelling of surface heat flow [1], with the aim of employing the unparalleled spatial uniformity of global gravity models in the fill-in of sparsely sampled surface heat flow data. The model we obtained, in a test area in Central Europe, provided additional information on the lithospheric structure and revealed a satisfactory coherence with the geological features in the area and their controlling effect on the conductive heat transport. That test was based on a fit of radioactive heat production to available heat flow data, based on a misfit linearization and substitution strategy, which we have shown to be independently consistent with available heat production relationships (e.g. Hasterok and Webb, 2017). Furthermore, model validation techniques provide additional metrics on the predictability in areas devoid of heat flow measurements.

To reach those objectives, we developed a finite-difference based solver for the heat equation in conductive, stable lithosphere, relying on the assumption of steady state, 3-D heat conduction from the thermal base of the lithosphere to surface. It allows for non-homogeneous heat production and thermal conductivity, and non-flat upper and bottom boundaries. Concurrent joint forward modelling of the gravity field is also possible.
Through compromise between complexity and approximation, it was designed favouring easy and fast forward modelling, such as in assessing parameter sensitivity and performing grid searches or parameter fitting. Geological models and parameters can be defined using an user-friendly plain text layer-wise definition, which is then turned into a volume, on a rectangular mesh.
Computational requirements are lean: a 75 × 75 × 104 node model such as the one employed in [1] can be forward-modelled on an ordinary workstation in 135 seconds. A direct solver is employed to solve the FD system of linear equations: the Matlab built-in Cholesky decomposition for sparse arrays (Davis, 2006).

Albeit initially developed as an ad-hoc tool for a proof of concept, its ease of use and versatility suggest its potential in other applications. We therefore present the solver and the accompanying tool set, both openly available, along with a set of promising examples.

[1] Pastorutti, A., Braitenberg, C. (2019) "A geothermal application for GOCE satellite gravity data: modelling the crustal heat production and lithospheric temperature field in Central Europe." Geophysical Journal International, doi:10.1093/gji/ggz344