An updated 3D Temperature model of the Netherlands

Subsurface temperature is a critical parameter when assessing the geothermal energy potential of a region. Regardless of how favorable an aquifer may be in terms of porosity, permeability, or depth, geothermal exploitation is not economically viable if temperatures are insufficient. As part of the ThermoGIS project, we produce nationwide estimates of geothermal energy potential for the entire Netherlands, which requires a high-resolution statistically robust model of subsurface temperature.

This research adopts a combined physics-based and data-driven approach to estimate the three-dimensional temperature field beneath the Netherlands down to a depth of 10 km. The model is discretized on a regular mesh with a horizontal resolution of 1 × 1 km in longitude and latitude, and a variable vertical resolution that averages between 10 and 30 m in the upper 5 km, and 200 m down to 10 km. This represents a five-fold increase in resolution compared to the most recent published temperature model of the Netherlands (Bekesi et al., 2020).

We first construct a three-dimensional lithological model of the Netherlands comprising 101 distinct litho-stratigraphic layers. Based on expert stratigraphic knowledge, lithological compositions are assigned to each layer. These layers are then populated with thermal conductivity and radiogenic heat production values derived from standard reference data (Hantschel and Kauerauf, 2009), yielding an initial prior model.

Using these prior conductivity and radiogenic heat fields, we solve the three-dimensional steady-state heat diffusion equation using centered finite differences. The model parameters are subsequently updated using Ensemble Smoother Multiple Data Assimilation (Emerick and Reynolds, 2013) to match a high-quality dataset of mostly corrected bottom-hole temperature measurements, and DST and geothermal production temperatures.

A key innovation distinguishing this model from previous temperature models of the Netherlands (Bonte et al., 2012; Bekesi et al., 2020) is the use of efficient numerical solvers combined with a more accurate and detailed lithological model, enabling an order-of-magnitude increase in spatial resolution. Our model was written entirely in python and the code will be made open source upon publication.