Eszter Békési
email: e.bekesi@uu.nl

Maartje Struijk

Damien Bonté

Hans Veldkamp

Jon Limberger

Peter A Fokker

Mark Vrijlandt

Jan-Diederik Van Wees

An updated geothermal model of the Dutch subsurface based on inversion of temperature data

Keywords: 3D thermal field, regional-scale conductive model, pseudo-convection, data assimilation, the Netherlands

 3D subsurface thermal model of the onshore Netherlands  Conductive heat transport  Pseudo-convective approach to approximate hydrothermal convection  Data assimilation using temperature measurements as observations

Introduction

Due to the increasing global energy demand, the shortage of conventional energy sources, and environmental aspects, the exploitation of renewable energy sources are gaining importance quickly. Geothermal energy is a good alternative to non-renewable energy sources, but geothermal energy systems need to satisfy several constraints to become economic. The key parameters of geothermal systems are the reservoir temperature and the fluid flow rate. The latter is strongly dependent on the reservoir permeability that should be sufficiently high to support the flow rate and on the reservoir thickness. Alternatively, favorable reservoir conditions should be present to allow for enhancement of the permeability (e.g. Enhanced Geothermal Systems [START_REF] Breede | A systematic review of enhanced (or engineered) geothermal systems: past, present and future[END_REF])). To assess potential areas for geothermal exploration, large-scale physics-based models integrating geophysics, geology, and geochemistry are required (e.g. [START_REF] Cloetingh | Lithosphere tectonics and thermo-mechanical properties: an integrated modelling approach for Enhanced Geothermal Systems exploration in Europe[END_REF]. However, for the development of a geothermal project, detailed site-specific studies are indispensable in order to minimize the pre-drilling risks associated with uncertainties in the subsurface.

In this paper we present a 3D thermal model of the onshore Netherlands. We describe the thermal state of the subsurface in order to identify thermal anomalies within the basins and basement, highlighting potential areas for geothermal development. Geothermal resources in the Netherlands can be classified as conduction dominated intracratonic basin plays [START_REF] Moeck | A new 'geothermal play type'catalog: Streamlining exploration decision making[END_REF]. Geothermal projects focus on direct heat uses, targeting formations restricted to 2-3 km depth. At those depths, the average geothermal gradient of ~31°C/km prohibits electricity production, since sufficiently high temperatures are only reached at larger depth, below 4-5 km. For clastic reservoirs in the Netherlands, permeability at 4-5 km depth is not high enough for geothermal exploitation. The only exception, could be the Dinantian carbonate platforms with their natural fracture networks. Typically, deep geothermal projects in low to medium enthalpy reservoirs are mostly exploiting granite or sandstone reservoirs. However, in the Netherlands the focus is on the fractured carbonate platforms. The 4-km deep Luttelgeest carbonate platform drilled by the LTG-01 well, shows evidence for the presence of hydrothermal convection [START_REF] Bonté | Subsurface temperature of the onshore Netherlands: new temperature dataset and modelling[END_REF][START_REF] Lipsey | Numerical modelling of thermal convection in the Luttelgeest carbonate platform, the Netherlands[END_REF][START_REF] Van Oversteeg | Fracture permeability assessment in deeply buried carbonates and implications for enhanced geothermal systems: inferences from a detailed well study at Luttelgeest-01, The Netherlands[END_REF]. Higher than average temperatures and the presence of fluids provide favorable conditions for deep geothermal exploitation. Therefore, one of the main aspects of the present study is to identify possible locations where convection in these carbonate platforms might occur and that currently lack temperature measurements.

Investigation of the subsurface temperature distribution in the Netherlands has been carried out for several decades. The first temperature maps of the Dutch subsurface are reported in the Atlas of subsurface temperatures in the European Community [START_REF] Haenel | Atlas of subsurface temperatures in the European Community[END_REF], although the temperature dataset used to construct these maps is not available. In the update of the geothermal atlas [START_REF] Haenel | Atlas of geothermal resources in the European Community[END_REF] a temperature dataset of the Netherlands including measurements from 388 wells is reported. A larger dataset containing 334 bottom hole temperatures (BHT) and 53 drill stem tests (DST) is available in the third edition of the geothermal atlas [START_REF] Hurtig | Geothermal atlas of Europe[END_REF]. The majority of measurements are obtained from the uppermost 3 km of the subsurface, except two values at 5 km depth. Temperature maps in the atlas are only available at European scale and are constructed based on the extrapolation of measured values from shallower depth. Temperature measurements from 464 wells were incorporated into the latest geothermal atlas [START_REF] Hurter | Atlas of geothermal resources in Europe[END_REF]. The maps for the Netherlands in this latter atlas were constructed by [START_REF] Rijkers | Atlas of Geothermal resources in the European Community, the Netherlands[END_REF] who updated the maps published in older versions of the atlases for the Lower Triassic and the Lower Cretaceous formations that have geothermal potential. The temperature maps were generated by kriging, showing significant improvements. Large-scale temperature models for Europe including the Netherlands, calibrated with a compilation of temperature models based on measurements were also constructed by [START_REF] Limberger | Assessing the prospective resource base for enhanced geothermal systems in Europe[END_REF]2018). [START_REF] Verweij | Fluid flow systems analysis on geological time scale in onshore and offshore Netherlands, with special reference to the Broad Fourteens Basin[END_REF] further characterized the temperature distribution in the onshore Netherlands using not only temperature measurements but also a physics-based approach, taking into account calculated estimates of the thermal conductivity and heatflow of the main lithostratigraphical units. Since then the onshore temperature dataset of the Netherlands has been continuously updated (e.g [START_REF] Bonté | Subsurface temperature of the onshore Netherlands: new temperature dataset and modelling[END_REF][START_REF] Boxem | Steady state 1D temperature modeling of the onshore Dutch subsurface[END_REF]. The first 3D temperature model of the onshore Netherlands has been constructed by [START_REF] Bonté | Subsurface temperature of the onshore Netherlands: new temperature dataset and modelling[END_REF]. Here we present an update of their work using the most recent temperature data, an updated geological model and a significantly improved inverse modeling approach.

We modeled the thermal field of the deep subsurface of the Netherlands in 7 steps. We first constructed a physics-based lithospheric-scale thermal model, hereafter referred to as prior model. The resulting thermal model shows a misfit with the observed temperatures, highlighting areas where the steady-state conductive thermal field is perturbed by transient effects and/or convective heat transport (i.e. paleo-temperature fluctuations, groundwater flow). This misfit is reduced with sequential data-assimilation and more detailed modeling procedures, by updating the thermal properties of the layers using ES-MDA (Emerick and Reynolds, 2013a). Ensemble methods such as the Ensemble-Kalman Filter (EnKF) and the Ensemble Smoother (ES) are applicable for large problems having many parameters. Emerick and Reynolds (2013b) showed that in case of non-linear problems, ES-MDA performs better and is computationally more efficient than EnKF. ES-MDA is most commonly applied for history matching in reservoir modeling. For instance, [START_REF] Fokker | Application of an ensemble smoother with multiple data assimilation to the Bergermeer gas field, using PS-InSAR[END_REF] constrained the model parameters of a compacting gas field in the Netherlands with ES-MDA using satellite data. Inversion using ES-MDA to calibrate temperature models was previously applied by [START_REF] Békési | Subsurface temperature model of the Hungarian part of the Pannonian Basin[END_REF] and [START_REF] Limberger | Refining the thermal structure of the European lithosphere by inversion of subsurface temperature data[END_REF]. We followed a similar methodology by applying the ES-MDA to constrain the thermal conductivity of the sedimentary units and heat generation in the upper crust. This work flow yielded our final thermal model hereafter referred to as posterior model, describing the thermal state of the uppermost 10 km of the Dutch subsurface.

Geology and geothermal conditions

The geological structure of the sedimentary cover of the Netherlands was described in detail by [START_REF] Kombrink | New insights into the geological structure of the Netherlands; results of a detailed mapping project[END_REF] (Figure 1). A high resolution 2.5D model (DGM-deep v4.0) representing the main sedimentary units to the base of the Carboniferous, reaching a depth up to 3-4 km in most of the country, is available from the Netherlands Oil and Gas Portal website (www.nlog.nl). A cross section is presented in Figure 2. Detailed descriptions of the units are available from the www.dinoloket.nl website [START_REF] Van Adrichem Boogaert | Stratigraphic nomenclature of the Netherlands, revision and update by RGD and NOGEPA[END_REF]. Units of Devonian, Silurian and Ordovician and older age have seldomly been drilled onshore. The descriptions of those units in the Dutch stratigraphic nomenclator are therefore incomplete and limited due to the lack of data available. . The pre-Variscan basement is unknown in the Netherlands itself, but is considered to consist of Avalonian crust [START_REF] Pharaoh | Tectonic evolution, Petroleum geological atlas of the Southern Permian Basin area[END_REF], which forms the core of the London Brabant Massif south of the Netherlands. [START_REF] Wong | Geology of the Netherlands[END_REF] provide a concise description of the Devonian and younger rocks in the Netherlands. Below, only a short description is given, largely based on [START_REF] Wong | Geology of the Netherlands[END_REF], [START_REF] Kombrink | New insights into the geological structure of the Netherlands; results of a detailed mapping project[END_REF], [START_REF] Bonté | Subsurface temperature of the onshore Netherlands: new temperature dataset and modelling[END_REF] and [START_REF] Smit | Early Carboniferous extension in East Avalonia: 350 My record of lithospheric memory[END_REF].

Figure 1. Early Carboniferous -Late Jurassic structural elements of the Netherlands. The colour coding reflects the remaining sedimentary succession. Areas that have been relatively stable in green and orange, inverted basins in blue (after [START_REF] Kombrink | New insights into the geological structure of the Netherlands; results of a detailed mapping project[END_REF].

In the Devonian and Carboniferous, lithospheric stretching and subsidence associated with the Variscan orogeny enabled the deposition in the Devonian of sandstones and shales in the southern Netherlands [START_REF] Kombrink | The Carboniferous of the Netherlands and surrounding areas; a basin analysis[END_REF] and reefal limestones in the north [START_REF] Van Hulten | Geological factors controlling Early Carboniferous carbonate platform development in the Netherlands[END_REF]. During the Early Carboniferous carbonate ramps in the south, and isolated platforms in the north were formed [START_REF] Reijmer | Seismic stratigraphy of Dinantian carbonates in the southern Netherlands and northern Belgium[END_REF], in response to widespread SW-NE directed extension in an extensional collapse setting [START_REF] Smit | Early Carboniferous extension in East Avalonia: 350 My record of lithospheric memory[END_REF], and forming the structural grain of major fault structures in Late Paleozoic and Mesozoic basin history. In the Late Carboniferous increased siliciclastic input from the Variscan thrust front buried the carbonate platforms. Thick shales filled the deep basin north of the orogen, followed by the deposition of fluvio-deltaic sandstones, shales and coals. The collapse of the Variscan orogen led to the formation of the Southern Permian Basin [START_REF] Doornenbal | Petroleum geological atlas of the Southern Permian Basin area[END_REF][START_REF] Van Wees | On the origin of the southern Permian Basin, Central Europe[END_REF][START_REF] Ziegler | Geological atlas of western and central Europe[END_REF]. Permian deposits, in the Netherlands, are mainly aeolian and fluvial sandstones (Rotliegend Group) overlain by evaporites (Zechstein Group). The evaporites are thick in the north, leading to the occurrence of salt diapirs, walls and domes [START_REF] Ten Veen | Thin-and thick-skinned salt tectonics in the Netherlands; a quantitative approach[END_REF]. Further south, the salt is absent and the Zechstein deposits are developed in a basin margin facies [START_REF] Geluk | Stratigraphy and tectonics of Permo-Triassic basins in the Netherlands and surrounding areas[END_REF]. The breaking up of Pangea in the Triassic caused subsidence and deposition of shales, followed by fluvial and aeolian sandstones in fault-bounded depocentres (Lower Germanic Trias, De Jager (2007). During the deposition of sediments of the Upper Germanic Trias, subsidence gradually increased, and evaporites (Röt), fluvio-lacustrine and marine marls and limestones of the Muschelkalk were deposited (Geluk et al., 2007). Increased sediment input led to the deposition of lagoonal and evaporite sediments of the Keuper Formation. During the Jurassic, smaller fault bounded basins and highs developed, in the Early Jurassic marine shales were deposited (De Jager, 2007). Rifting in the North Sea leading to doming caused erosion and non-deposition [START_REF] Wong | Geology of the Netherlands[END_REF]. In the Late Jurassic and Early Cretaceous, fluvial to shallow marine siliciclastic sediments were deposited in small faultbound basins (Schieland and Rijnland Groups, De Jager (2007), [START_REF] Willems | Doublet deployment strategies for geothermal Hot Sedimentary Aquifer exploitation: Application to the Lower Cretaceous Nieuwerkerk Formation in the West Netherlands Basin[END_REF]). By the end of the Cretaceous, marine shales and carbonates are deposited. In the Late Cretaceous, the northward movement of Africa towards Europe caused the Alpine orogeny, leading to inversion tectonics with uplift and erosion. Increased sediment input led to the deposition of siliciclastic sediments of the Cenozoic North Sea Group. The Quaternary is dominated by fluvio-deltaic, glacial, and coastal depositional siliciclastic sediments. Figure 1 illustrates the tectonic history. It shows areas that, between the Early Carboniferous and Late Jurassic, have remained relatively undisturbed highs in orange, platform areas in green, and inverted basins in blue. The colour coding used in Figure 1 reflects the sedimentary succession that currently exists, and therefore the tectonic history. For instance, in the Friesland Platform area, shown in dark green in the North, Cretaceous sediments are directly overlying Permian Zechstein evaporites, Rotliegendes sandstones and sandstones and shales of the Carboniferous. This is also illustrated in the cross section of Figure 2, which shows the sedimentary units. Note that faults are not shown in Figure 2 because they are not available in the online version of the DGM-deep model, and fault structures were also not incorporated to the thermal model. Also note that the units shown are based on seismic interpretations, and that they do not necessarily correlate fully to the units defined in the ICS timescale. Currently, 25 doublets have been drilled, of which 16 are producing, 4 are suspended, and 5 are not yet producing (source: www.nlog.nl). The energy is used for direct heating purposes, especially greenhouses and in a few cases city heating. This requires temperatures between about 60 and 100 °C. Assuming a geothermal gradient of about 31.3 °C/km and an average surface temperature of 10 °C [START_REF] Bonté | Subsurface temperature of the onshore Netherlands: new temperature dataset and modelling[END_REF], this temperature range exists approximately between 1.5 and 3.0 kilometers depth. With the exception of two doublets that targeted a fault in Carboniferous limestones, all doublets target sandstone reservoirs having sufficient primary porosity and permeability to produce hot water at such rates that stimulation of the reservoir is not required. Figure 2, a generalized cross section through the Netherlands, shows which reservoirs are found at these depths. The main targets for geothermal exploration are the sandstones from the Rotliegend (mainly in the Northern part of the country), Lower Germanic Trias and Upper Jurassic/Lower Cretaceous (mainly in the West Netherlands Basin). To a lesser extent, fractured Devonian quartzites and Carboniferous Limestones are targeted, and Paleogene sandstones.

Subsurface temperature data

Subsurface temperature data from the onshore Netherlands are available on the Dutch Oil and Gas portal (http://www.nlog.nl) and on the Dutch Geothermal Platform (https://geothermie.nl). These temperature measurements are based on Bottom Hole Temperature (BHT) data, Drill-Stem Tests (DST) and Repeat Formation Tests (RFT). Although a large number of BHT data are available for the Dutch onshore, correction is required before they can be used for modeling.

BHTs are recorded as maximum temperatures from well logs and are assumed to be the highest at the bottom of the drilled well. These temperatures may differ from the actual formation temperatures due to the drilling operations. Corrections of BHT data are most commonly based on analytical methods (e.g. [START_REF] Goutorbe | Comparison of several BHT correction methods: a case study on an Australian data set[END_REF] or statistical methods such as the American Association of Petroleum Geologists (AAPG) correction. Numerous modelling methods are also available (e.g. [START_REF] Luijendijk | Thermal state of the Roer Valley Graben, part of the European cenozoic rift system[END_REF][START_REF] Rühaak | 3-D interpolation of subsurface temperature data with measurement error using kriging[END_REF] but remain difficult to apply for large datasets.

The temperature dataset we used in this study is based on the original database from [START_REF] Bonté | Subsurface temperature of the onshore Netherlands: new temperature dataset and modelling[END_REF], extended with new temperature measurements (Figure 3). For correction of the new BHTs we followed the same workflow as [START_REF] Bonté | Subsurface temperature of the onshore Netherlands: new temperature dataset and modelling[END_REF]. The new dataset contains 438 BHTs corrected with the analytical method sampled from 199 wells and 987 statistically corrected BHTs from 401 wells, yielding a total of 1425 BHT temperature values. The total number of raw BHTs is significantly larger than the amount of corrected data, as analytical corrections require multiple measurements (Figure 3c). In addition, the new dataset contains 65 DST and RFT measurements from 36 wells. We also included 17 production temperature measurements from geothermal wells collected from the Dutch Geothermal Platform (https://geothermie.nl).

The new database consisting of 1507 temperature measurements from 505 wells ( Figure 3a,b). The measurements yield an average geothermal gradient of 31 °C/km. However, there are some deep anomalies with high temperatures such as the LTG-01, TJM-02-S2, and WSK-01 wells (see Figure 9 for locations). The temperature field is forward modeled solving the heat equation in 3D, assuming steady-state conditions and conductive heat transfer only:

( ) [1]
where λ is the thermal conductivity, T is the temperature, A is the radiogenic heat generation, and ( )

is the nabla operator. We solve the equation for a low and high resolution hexahedral grid with geometry specifications presented in section 5. Solutions were obtained by a finite-difference approximation using the Preconditioned Conjugate Gradient method (PCG) [START_REF] Limberger | Refining the thermal structure of the European lithosphere by inversion of subsurface temperature data[END_REF]. The boundary conditions at the top were identical for all models: a stationary surface temperature of 8 °C. This value deviates from the average observed yearly surface temperature, which is 10 degrees. We first tested our models with a fixed top boundary condition of 10 °C. Such models systematically overestimated the temperature measurements at depth shallower than 2500 m. To decrease these misfits, we tested our workflow with a lower value of 8 °C, that provided a better fit with the data: the root mean square (RMS) errors at 1000 m and 2000 m (±500 m) depth decreased from 0.8 °C and 0.84 °C to 0.71 °C (Table 3). Therefore, we choose a lower value to our final model to account for the paleo-surface temperatures effected by the recent ice age (e.g. [START_REF] Donders | Strong climate coupling of terrestrial and marine environments in the Miocene of northwest Europe[END_REF][START_REF] Verweij | Reconstruction of basal heat flow, surface temperature, source rock maturity, and hydrocarbon generation in salt-dominated Dutch Basins[END_REF].

The boundary conditions at the bottom were different for the low-and high-resolution models. We chose a fixed temperature at the bottom of the lithosphere with a corresponding value of 1200 °C for the low-resolution models. The heat flow at 10 km depth, obtained from the low-resolution model, was used as boundary condition for the high-resolution models. The vertical edges of all models were assumed to be insulating with a fixed heat flow of zero.

The motivation to use a steady-state model, excluding transient thermal effects, is computational performance for the ES-MDA inversion which requires 100s of model runs. We argue that the use of a steady state model is justified since in most of the Netherlands active tectonic processes are absent or marked by very low sedimentation rates, up to a maximum of 0.1 mm/y in the southeast of the Netherlands, and have resulted in less than 1.5 km of sediments in the last 20 million years. Previously studies, incorporating transient effects, demonstrate that such low sedimentation rates do not deviate considerably from the steady state assumptions (e.g. [START_REF] Van Wees | Probabilistic tectonic heat flow modeling for basin maturation: Assessment method and applications[END_REF].

Thermal convection in the Dinantian carbonate platforms

The previous 3D temperature model from [START_REF] Bonté | Subsurface temperature of the onshore Netherlands: new temperature dataset and modelling[END_REF] could not reproduce the thermal anomaly measured at a depth below 4 km in the LTG-01 well (see Figure 4 for location). [START_REF] Bonté | Subsurface temperature of the onshore Netherlands: new temperature dataset and modelling[END_REF] suggested that higher than expected temperatures and the abnormally low thermal gradient observed within the Dinantian carbonates might be explained by the occurrence of hydrothermal convection. Another possible explanation of the thermal profile observed in the LTG-01 well is the large thermal conductivity contrast between the overlying Westphalian and Namurian layers and the Dinantian carbonates, based on petrophysical analysis in the framework of the SCAN project [START_REF] Veldkamp | Temperature modelling of the Dutch subsurface at the depth of the Dinantian[END_REF]. [START_REF] Van Oversteeg | Fracture permeability assessment in deeply buried carbonates and implications for enhanced geothermal systems: inferences from a detailed well study at Luttelgeest-01, The Netherlands[END_REF] and [START_REF] Lipsey | Numerical modelling of thermal convection in the Luttelgeest carbonate platform, the Netherlands[END_REF] investigated the potential for thermal convection at the LTG-01 well and confirmed that convection is likely to occur in the Dinantian carbonate platform To come to this conclusion, they assessed the minimum permeability (k min ) required for convection based on Rayleigh number analysis following [START_REF] Horton | Convection currents in a porous medium[END_REF] and [START_REF] Lapwood | Convection of a fluid in a porous medium[END_REF]. [START_REF] Lipsey | Numerical modelling of thermal convection in the Luttelgeest carbonate platform, the Netherlands[END_REF] also supported the occurrence of fluid convection with permeability assessment on core samples and numerical modeling.

The Rayleigh number is a dimensionless value that indicates the likelihood of natural or free convection. The equation for the Rayleigh number in a porous medium is written as:

[2]

where is the permeability, g is the gravitational acceleration, is the temperature difference between the top and the bottom of the layer, H is the thickness of the layer and , are the density, the specific heat capacity, the thermal expansion coefficient of the pore fluid, respectively. and λ denote the viscosity of the fluid and the bulk thermal conductivity of the of the medium. If Ra exceeds a certain threshold value, referred to as the critical Rayleigh number (

), convection can take place in the medium. For a horizontal, homogenous isotropic porous medium bounded by fixed temperature conditions, . Equation ( 1) can be rewritten to determine k min using :

[3]

The minimum permeability to initiate hydrothermal convection, k min , for the Dinantian carbonates was calculated to be 3x10 -14 m 2 [START_REF] Van Oversteeg | Fracture permeability assessment in deeply buried carbonates and implications for enhanced geothermal systems: inferences from a detailed well study at Luttelgeest-01, The Netherlands[END_REF] and1.9x10 -14 m 2 [START_REF] Lipsey | Numerical modelling of thermal convection in the Luttelgeest carbonate platform, the Netherlands[END_REF]. The difference between the two values mostly originates from the different values used for the thickness of the convective layer in the Dinantian carbonates, being 600 m and 800 m, respectively. Since the proposed permeability for the interval of the Dinantian by [START_REF] Van Oversteeg | Fracture permeability assessment in deeply buried carbonates and implications for enhanced geothermal systems: inferences from a detailed well study at Luttelgeest-01, The Netherlands[END_REF] was estimated to be 6x10 -14 m 2 and [START_REF] Lipsey | Numerical modelling of thermal convection in the Luttelgeest carbonate platform, the Netherlands[END_REF] reported the possibility of even higher values (up to 10 -13 m 2 ), convection is likely to explain the temperature profile in the LTG-01 well. The extent of the carbonate platforms within the group based on seismic data and well logs, adopted from [START_REF] Kalkman | Ultra-Diepe Geothermie -Casus Almere en Barendrecht[END_REF], is outlined. We assumed that convection may occur in the deeper platforms, marked by the thick black outline.

Based on these findings, we assumed potential hydrothermal convection in the carbonate platforms all over the Netherlands. However, these carbonate platforms are found at significantly different depths in the Netherlands: in the south the platforms are located much shallower than in the north. (Figure 4). is dependent on the pore fluid properties, which show a significant variation with depth. [START_REF] Van Oversteeg | Fracture permeability assessment in deeply buried carbonates and implications for enhanced geothermal systems: inferences from a detailed well study at Luttelgeest-01, The Netherlands[END_REF] analyzed the effect of pressure and temperature on k min and showed that k min strongly decreases with depth. We constructed similar k min -depth curves for different layer thicknesses (Figure 5) to assess if convection could occur in the carbonate platforms in the south of the Netherlands. We calculated fluid properties using a model dependent on temperature, pressure, and salinity (after Van Wees et al. ( 2012)), where the temperature, pressure and salinity dependence of density and viscosity is incorporated after [START_REF] Batzle | Seismic properties of pore fluids[END_REF], and the temperature and salinity dependence of specific heat capacity is based on [START_REF] Grunnberg | Properties of seawater concentrations[END_REF]. The temperature was calculated with a constant geothermal gradient of 31 °C/km, and the pressure was assumed to be hydrostatic. The minimum Rayleigh permeability at larger depth is relatively small: between 3.5-5.5 km k min ranges from 1 to 5 x 10 -14 m 2 (Figure 5a). At shallower depth, however, the values are significantly larger (Figure 5b): at 2 km depth and shallower, k min is larger than ~10 -13 m 2 , which is equal to the highest permeability of the Dinantian [START_REF] Lipsey | Numerical modelling of thermal convection in the Luttelgeest carbonate platform, the Netherlands[END_REF]. We thus concluded that thermal convection is not likely to occur in the carbonate platforms in the south. Our calculations are based on several assumptions: we adopted a horizontal, homogenous isotropic porous medium bounded by fixed temperature conditions. Other factors such as the effect of the extent and geometry of the platform, heterogeneity in permeability can play an important role in the presence/absence of convection [START_REF] Lipsey | Numerical modelling of thermal convection in the Luttelgeest carbonate platform, the Netherlands[END_REF]. These factors need to be taken into account for site-specific studies. Still, our calculations can be used as an approximation to assess the potential for hydrothermal convection in the Dinantian carbonate platforms all over the Netherlands.

Model geometry and thermal properties 5.1. Model geometry

The thermal models have been built in the Dutch coordinate system (EPSG:28992) and their outline has been defined by the political boundaries of the Netherlands. We performed the modelling in two resolutions: the low resolution models extend from the surface up to the depth of the Lithosphere-Asthenosphere Boundary (LAB), while the high resolution models only include the uppermost 10 km of the subsurface. For the low resolution models, temperatures were calculated through a regular 3D grid with a horizontal resolution of 2.5x2.5 km, a vertical resolution of 200 m for the uppermost 10 km, and 3 km down to the LAB. The high resolution models are defined by 1x1 km and 200 m grid spacing in horizontal and vertical directions, respectively.

The model was built as a layered structure that includes the lithospheric mantle, lower crust, upper crust and 14 sedimentary units, summarized in Table 1. For the thickness of the lithosphere we used a constant value of 110 km. We also tested models for the LAB depth after [START_REF] Tesauro | A new thermal and rheological model of the European lithosphere[END_REF] and [START_REF] Artemieva | Lithosphere structure in Europe from thermal isostasy[END_REF], but the thermal effect of the different LAB models in the shallow part of the model (0-10 km) was negligible. The crustal thickness was adopted from the EuCRUST-07 model [START_REF] Tesauro | EuCRUST-07: A new reference model for the European crust[END_REF] with an average thickness of 32 km in the Netherlands. For the subdivision of the sedimentary layers we followed the DGM-deep v4.0 onshore model (adopted from nlog.nl, Figure 2). The sedimentary model is available to the base of the Limburg Group, since the depth of the base of the Carboniferous Limestone Group (Dinantian) is not yet well constrained in all parts of the Netherlands. We estimated the thickness of the Dinantian to be 700 m, which is a rough average of the thicknesses found in the LTG-01 and UHM-02 wells that drilled the entire Dinantian section. It is important to note that the thickness of the Dinantian might be lower in some areas. For example, in the southern part of the Netherlands, where the Dinantian is shallower, the estimated thickness is mostly below 500 m [START_REF] Reijmer | Seismic stratigraphy of Dinantian carbonates in the southern Netherlands and northern Belgium[END_REF]. We selected a constant thickness of 700 m based on the assumption that older sediments with similar thermal properties are likely to be present below the Dinantian.

Following the example of the thermal anomaly found in the Dinantian carbonates in the Luttelgeest-01 (LTG-01) well, we assumed that hydrothermal convection may occur in the carbonate platforms. To distinguish between the basin facies and platform areas of the Dinantian, we treated the carbonate platforms as a separate layer by adopting their geometry after Kalkman (2016) (Figure 6), based on 2D and 3D seismic data and 15 well logs.

The outline of the platforms is marked by larger uncertainty in the northern part of the country, where the platforms are deeply buried (Figure 4). Also, in the South there are multiple wells drilling the Dinantian carbonates, providing constraints for the seismic data. We excluded the platforms in the South from the new layer, assuming that hydrothermal convection is most likely restricted to larger depths.

In the southern part of the Netherlands we included a new layer to represent the deeply buried Palaeozoic sediments with low porosity below the Dinantian (Figure 10, layer 12). We constructed this layer based on crosssections from [START_REF] Duin | Subsurface structure of the Netherlands-results of recent onshore and offshore mapping[END_REF], although its geometry is poorly constrained due to the lack of available data. [START_REF] Kalkman | Ultra-Diepe Geothermie -Casus Almere en Barendrecht[END_REF]. We use scaling of the thermal conductivity to approximate the thermal effects of convection. Outlined in black are the areas below the carbonate platforms for which we scaled the thermal conductivity values.

Thermal properties 5.2.1. Thermal conductivity

We defined the thermal properties of each sedimentary unit based on their lithologies following the methodology of [START_REF] Hantschel | Fundamentals of basin and petroleum systems modeling[END_REF] (Table 1). The bulk matrix thermal conductivity (λ m ) of the lithological components were corrected for pressure conditions and the in-situ temperature after [START_REF] Sekiguchi | A method for determining terrestrial heat flow in oil basinal areas[END_REF]. For shales and carbonates, we took into consideration the change in anisotropy with increasing compaction [START_REF] Hantschel | Fundamentals of basin and petroleum systems modeling[END_REF]. The horizontal matrix conductivity was calculated from λ m using the anisotropy factors of the corresponding lithology. We obtained the bulk thermal conductivity (λ bulk ) for each lithology by taking the geometric mean of k m and the temperature-dependent thermal conductivity of the pore fluid (λ w ) as follows:

( ) [4]
where ϕ is the porosity. If a sedimentary unit consisted of different lithologies, the bulk thermal conductivity within the layer was calculated by taking the harmonic mean of the bulk conductivities.

Because our forward thermal model only solves the steady-state heat equation for conduction, we used a pseudoconvective approach to reproduce the deep thermal anomaly found in the LTG-01 well: we scaled the prior thermal conductivity of the layer corresponding to the carbonate platforms. This method, together with a similar modelling technique, was applied by [START_REF] Békési | Subsurface temperature model of the Hungarian part of the Pannonian Basin[END_REF] for the Pannonian basin. After several model runs to test different thermal conductivity values of the platform layer, we could not reproduce the deep thermal anomaly. Our models managed to fit the convective profile only by the combination of increasing the thermal conductivity of the platforms and decreasing the conductivity of the layer overlain by the Dinantian (Limburg Group). We tried to choose the thermal conductivity of the layer corresponding to the Limburg Group close to its original value, but low enough to fit the temperature measurements in the Dinantian. We tested numerous combinations before attaining acceptable values. We fixed the thermal conductivity of the Limburg Group to 1.5 W m -1 K -1 , and the Dinantian carbonate platforms to 5.0 W m -1 K -1 (Table 1). It is important to note that scaling the thermal conductivities of limestones to unrealistic values of 5 W m -1 K -1 was required to approximate the effect of convection with a purely conductive model. The low thermal conductivity of 1.5 W m -1 K -1 assigned to the Limburg Group is also rather unrealistic, considering its compaction and low porosity. However, Silesian sediment with high coal content might explain such low values. We used these values as prior thermal conductivities before applying the data assimilation.

The thermal conductivity of the upper and lower crust were estimated with the temperature-and pressuredependent relation of [START_REF] Chapman | Thermal gradients in the continental crust[END_REF]. With increasing temperature down to the lithospheric mantle, the contribution of the radiative component of the thermal conductivity increases compared to the contribution of the lattice component. Therefore, the thermal conductivity of the lithospheric mantle was calculated as the sum of the temperature-dependent radiative contribution obtained after [START_REF] Schatz | Thermal conductivity of earth materials at high temperatures[END_REF] and the temperature-and pressure-dependent lattice component following the formula of [START_REF] Xu | Thermal diffusivity and conductivity of olivine, wadsleyite and ringwoodite to 20 GPa and 1373 K[END_REF]. More details on the calculations of the thermal conductivity values of the different units are described by [START_REF] Limberger | Refining the thermal structure of the European lithosphere by inversion of subsurface temperature data[END_REF].

Radiogenic heat production

For the radiogenic heat production of each sedimentary unit, a constant value was selected based on generic values of typical lithologies from [START_REF] Hantschel | Fundamentals of basin and petroleum systems modeling[END_REF]. The heat production of the layers consisting of more than one lithology was calculated as the geometric mean of the corresponding lithologies. We chose 1 μW m -3 as a prior value for heat generation in the upper crust, which was then updated through data assimilation to account for potential under/overestimation of the initial radiogenic heat production. For heat generation in the lower crust and lithospheric mantle, we assigned 0.5 and 0.02 μW m -3 , respectively, also following [START_REF] Hantschel | Fundamentals of basin and petroleum systems modeling[END_REF]. To construct the thermal model of the deep subsurface of the Netherlands, we first adopted the prior thermal properties of the layers listed in Table 1. Initial calculations were made by solving the heat equation in multi-1D (step 1, Figure 7) in the low resolution grid, in order to obtain default values for thermal conductivity in agreement with first order estimates for temperature and pressure. Subsequently, based on these default properties we calculated a prior 3D thermal model (step 2 Figure 7).

ID Layer name

The misfit between modeled and observed temperatures was subsequently reduced by updating the thermal properties of the layers using ES-MDA in step 3-4 (Figure 6). The results are subsequently used in a higher resolution grid and the final posterior thermal model is obtained in step 7, as shown in Figure 7 and described in the following sections.

Inversion

We used inversion of subsurface temperature data to infer the thermal field by varying a selection of the thermal properties of the layers: the radiogenic heat generation of the upper crust and the thermal conductivity of the sedimentary units.

Ensemble Smoother with Multiple Data Assimilation (ES-MDA)

We solved the inverse problem using the Ensemble Smoother (ES), which estimates the parameters in a single step by a global update incorporating all data available (Emerick and Reynolds, 2013a). Ensemble-based methods, such as the Ensemble-Kalman Filter (EnKF) and the ES, are suitable for systems with large numbers of parameters. Emerick and Reynolds (2013b) investigated different ensemble methods and found that in case of non-linear forward models, the Ensemble Smoother with multiple data assimilation (ES-MDA) performs better than the EnKF. Additionally, the ES-MDA is computationally less demanding. Therefore, we applied the ES-MDA for handling the non-linearity between the observations and model parameters.

For the ES, the assimilated model parameter vectors m a are written as (Emerick and Reynolds, 2013a):

( ) ( ) [5] for j = 1, 2, […], N e ,
where N e denotes the number of ensembles, is the cross-covariance matrix between the prior vector of model parameters, m f , and the vector of predicted data, d f . is the auto-covariance matrix of predicted data, and is the covariance matrix of the measurement errors, where N d denotes the number of measurements assimilated.

( ) is an ensemble of vectors with the measurements vector , perturbed normally according to the covariance matrix . An ensemble of solutions is produced by the ensemble update, consistent with the prior statistics, and the mean of the ensemble is taken as the best estimate. In case of the ES-MDA, an ES is applied multiple times, with the output ensemble used as input for the next update. The number of data assimilation steps or iterations, N a must be selected a-priori. The data covariances used for the update steps are increased by a multiplication factor, α i for i=1,2…, N a , and α i must be selected as ∑ . This covariance adjustment to increase ensemble variance is required to reduce filter divergence and to improve performance of the state vector.

Data uncertainty

We used the temperature measurements as observations for the inversion. We assigned uncertainties to the temperature data based on measurement quality and assuming a Gaussian error distribution. We made this assumption after assessing the misfits of the raw and corrected BHT gradients with the average geothermal gradient, that can be approximated with normal distribution (Figure 3c). It is important to note that measurement errors are not necessarily symmetric and they are depth-dependent (e.g. [START_REF] Agemar | Subsurface temperature distribution in Germany[END_REF]. Here we select identical uncertainties for the same measurement types with the same correction method for simplicity.

Uncertainties may be overestimated due to the fact that we choose maximum values assigned to each category. We marked the DSTs and RFTs with an uncertainty of ±8 °C (e.g. [START_REF] Bonté | Subsurface temperature of the onshore Netherlands: new temperature dataset and modelling[END_REF]. For the BHTs, we chose different uncertainties based on the correction method applied to the measurement. We used ±10 °C for BHTs corrected with the analytical method and ±15 °C for values obtained by the statistical method [START_REF] Goutorbe | Comparison of several BHT correction methods: a case study on an Australian data set[END_REF]. The production temperatures of the geothermal wells were treated with an uncertainty of 5 °C. [START_REF] Saeid | An efficient computational model for deep low-enthalpy geothermal systems[END_REF] suggested that the maximum drop in in production temperatures occur within the first years of operation. Additionally, temperatures of the extracted fluids decrease before reaching the wellhead. Since ~half of the production temperatures were obtained within the first year of production (geothermie.nl), we concluded that the uncertainty of ±5 degrees is a good approximation.

During the inversion, only one observation within a grid cell is used. When multiple measurements were present, inversion was limited to the measurement with the lowest uncertainty. As a result, the number of observations in the high-resolution model ( 1284) was larger than for the low-resolution models (1260, Figure 7). Most of the temperature measurements (~85 %) were conducted at depths below 3 km, ensuring a good calibration for the upper part of the model. At larger depths, the amount of controlling points significantly decreases (Figure 3, Table 3.).

Description of the data assimilation procedure

Through the data assimilation, the prior thermal properties including the heat production in the upper crust and thermal conductivity of the sediments were updated to achieve a better fit with the temperature observations. The prior uncertainty in the thermal properties was taken into account by scaling the parameters to a triangular distribution. The spatial variability of the parameters was determined through a spherical variogram. The radius of the variogram is expressed in cells, which corresponds to a certain distance in kilometers.

We obtained the posterior thermal model in 7 steps (Figure 7), where steps 3, 4 and 7 included the parameter updates summarized in Table 2. In each step the ES-MDA was run with 600 ensembles and 4 iterations. We chose 4 iterations following the example of Emerick and Reynolds (2013b). In step 3, we scaled the upper crustal heat generation between 0.6 and 1.6 of the prior value, with a variogram range of 14 cells (corresponding to ~42 km). Then, the thermal conductivity of the sediments was scaled between 0.7 and 1.3, except for the Zechstein Group (layer 8) and the Carboniferous Limestone Group (layers 12 and 13), where a larger variation was allowed. We chose scaling factors of 0.7 and 1.4 for the Zechstein Group (layer 8 ) to account for potential higher thermal conductivity of salts. We scaled the thermal conductivity of the Carboniferous Limestone Group (layers 12 and 13) using a larger range, to allow for more variance approximating the thermal effect of hydrothermal convection (especially in the platform areas). We assumed a smaller-scale lateral variation in the thermal conductivity of the sediments (layer 1-11, 13): in step 4, the range of variations was set to 8 cells (~24 km). In the last modelling step (step 7), further refinement of the thermal conductivity of the sediments was achieved by scaling the values between 0.8 and 1.2 with a variogram range of 30 cells (~30 km).

Results

We show the results of the 7-step modeling workflow, yielding the posterior 3D temperature model of the onshore Netherlands in Figure 8. Temperature slices of 1-6 km depth are presented in 1 km intervals (Figure 8). Modeled temperatures at 1 and 2 km depth show a similar pattern: the hottest area is located in the north, corresponding to the youngest sedimentary units of Neogene to Upper Cretaceous age (Upper North Sea Group, Lower and Middle Sea Group and Chalk Group). At larger depths, the most significant positive anomaly is observed near the location of the LTG-01 well. At 3 km depth, the highest temperatures are modeled near the Belgian border in the south, where Dinantian carbonate platforms are located at relatively shallow depth. Higher temperatures at 4 km depth correlate with the extent of the Dinantian carbonate platforms. This correlation is also observed in the 5 and 6 km depth slices, although the anomalies are more pronounced towards the north and near the LTG-01 well. The Zeeland High and Limburg High appear to be relatively cold at larger depth. To compare our models with the temperature measurements, we constructed 1D profiles at several well locations (Figure 9). We show both the prior and the posterior model results to demonstrate the model improvements after calibration. As prior model we use the uncalibrated model from step 2 in Figure 4.

Prior and posterior model temperatures are almost identical at the locations of the KTG-01, LUT-06, HLH-GT-01,and WWK-01 wells, showing a good fit with the measured temperatures (Figure 9). The thermal conductivity profiles show little variation, suggesting calibration was not necessary to improve the fit. For the remaining profiles, the thermal conductivity of the layers was varied resulting in a better fit of the posterior model with the measurements. For instance, the prior model overestimates temperatures at the location of the MKP-14 well.

Decreasing the thermal conductivity of the Rijnland Group (layer 4) and increasing the thermal conductivity of the Altena Group (layer 6), resulted in a reduced average geothermal gradient, in good agreement with the measurements. In the LTG-01 and WSK-01 wells, the prior model underestimates the temperatures especially below 3 km depth. To fit the high temperatures in the LTG-01 well, an even lower thermal conductivity was required for the Caumer Subgroup (layer 10) and Limburg group (layer 11) to increase the geothermal gradient above the Dinantian carbonate platform (layer 12). Even though the misfit of our prior model with the WSK-01 well measurement is low, a small reduction (~0.2 W m -1 K -1 ) of the conductivity for layer 11 was still necessary to reproduce the measured values. 1. Note that the depth ranges of the prior and posterior thermal conductivities are shifted in depth for some layers as a result of depth and composition differently averaged due to the different horizontal resolution. For locations of the wells see the map in the lower right.

Posterior model temperatures are presented along a section crossing the onshore part of the Netherlands from SW to NE (Figure 10). Temperatures up to 1 km depth show slight variations corresponding to changes in the lithology of the layers. At larger depths, more fluctuations are observed within the sediments. For instance, in the NE, where the thickness of the Zechstein Group (layer 8) is the largest, temperatures are higher above and lower below the salt layer. The geothermal gradient in the top 3 km corresponds to ~30 °C/km. The insulating effect of the sediments is observed at larger depth: temperatures are lower within the basement and higher in the basins in the SW at 3-7 km depth. The heat chimney effect of the Dinantian carbonate platforms (layer 12) approximated by increasing thermal conductivity, results in generally higher temperatures above, and lower temperatures below the layer. The hottest area along the section below 3 km depth corresponds to the location of the LTG-01 well, where temperature measurements suggest the presence of hydrothermal convection. The geothermal gradient varies significantly with larger depth: it is generally lower in the basement and in the highly conductive carbonate platforms, and higher in the sediments with a lower conductivity, especially in the layers overlaying the platforms.

Figure 10. NW-SE cross-section through the onshore Netherlands showing the geometry of the sedimentary units and the depth of the top of the basement (white) used for the modeling. Superimposed are the isotherms from the posterior temperature model. The trace of the cross-section is shown in Figure 9. Layer numbering corresponds to Table 1. The extent of the main Jurassic and Early Cretaceous basins, highs, and platforms along the section is adopted from [START_REF] Kombrink | New insights into the geological structure of the Netherlands; results of a detailed mapping project[END_REF]. Fault structures are not shown in the section as they were not incorporated to the thermal model.

After ES-MDA calibration the posterior model shows an improved fit with the measurements compared to the prior model: temperature observations are closer to the P50 results of the posterior model (Figure 11). On the other hand, some of the data points fall outside the bandwidth of the variation of the ensembles. It suggests that the ensemble variance generated with the data assimilation procedure was not entirely sufficient to reproduce the observed variation within the data, especially at shallow depth.

The mean, median, and RMS misfits normalized to data error of the prior and posterior model are summarized in Table 3. After ES-MDA calibration, overall RMS, including all the measurements, decreases from a value of 0.95 for the prior model to 0.70 for the posterior model. Analyses of the misfits at different depth intervals suggest that the prior model overestimates temperatures at shallower depth (1-2 km), but underestimates temperatures at larger depth. The posterior model shows the same trend, but with lower misfits, especially below 3 km depth. The total mean and median of the prior model are smaller than the corresponding values of the posterior model. This is due to the fact that the positive and negative misfits at different depth cancel out, and it is not representative in terms of the goodness of the fit. Depth 3. Mean, median, and RMS misfits in °C of the prior and posterior models normalized to data error. The misfits are reported at 1-6 km depth, including measurements within a ±500 m interval. The overall misfits are calculated using all the measurements.

Discussion

We constructed 3D physics-based temperature models for the onshore Netherlands calibrated with 1507 temperature measurements. We performed the modeling in seven steps, starting from lithospheric-scale models with two different resolution steps in vertical direction. The advantage of such large-scale temperature models is that they are able to account for the thermo-mechanical state of the whole lithosphere, which has major influence on temperatures at depths relevant for geothermal exploration [START_REF] Cloetingh | Lithosphere tectonics and thermo-mechanical properties: an integrated modelling approach for Enhanced Geothermal Systems exploration in Europe[END_REF]. The boundary conditions of the high-resolution models, extending from the surface to 10 km depth, were derived from these lowresolution lithospheric-scale models. We incorporated temperature measurements for subsequent model steps through an inversion procedure using ES-MDA. Misfits between modeled and observed temperatures decreased during the data assimilation by varying the thermal properties of the sedimentary units and the upper crust. To account for uncertainties of the thermal conductivity in different sedimentary units and the radiogenic heat generation in the upper crust, we scaled the values to triangular distributions. Scaling parameters were determined according to the uncertainties we assigned to each layer. The spatial variation in the properties was introduced through variograms. We selected the size of the variograms based on the wavelengths of perturbations we expected for the thermal properties.

One of the main aspects of our study was to account for the thermal anomaly found in the LTG-01 well. [START_REF] Bonté | Subsurface temperature of the onshore Netherlands: new temperature dataset and modelling[END_REF] suggested that the high temperatures that were measured could be explained by the presence of a magmatic intrusion with high radiogenic heat production formed during the Variscan orogeny [START_REF] Ziegler | Geological atlas of western and central Europe[END_REF]. As an alternative, they proposed the occurrence of hydrothermal convection in the Dinantian carbonate platform situated at depths below 4 km. Based on fracture permeability assessments and numerical modeling, [START_REF] Van Oversteeg | Fracture permeability assessment in deeply buried carbonates and implications for enhanced geothermal systems: inferences from a detailed well study at Luttelgeest-01, The Netherlands[END_REF] and [START_REF] Lipsey | Numerical modelling of thermal convection in the Luttelgeest carbonate platform, the Netherlands[END_REF] concluded that convection is likely to occur in the carbonate platforms at the location of the LTG-01 well. The thermal state of the deep subsurface is important for geothermal exploration that targets the deeply buried Devonian-Carboniferous carbonate formations in the Netherlands. These reservoirs are potentially suitable for industrial heating applications and electricity production. To this end, our goal was to identify possible locations where convection in these carbonate platforms might occur and that currently lack temperature measurements.

We assessed the minimum permeability required for convection in different depth intervals based on Rayleighnumber analysis. Given the permeabilities derived from petrophysical data on the Dinantian carbonates, we found that the occurrence of convection is unlikely in carbonate platforms at depths shallower that 2 km. Locally, karstification or faults and fractures might still lead to sufficient secondary permeability allowing for hydrothermal convection. However, this is currently not supported by observations. We assumed that convection might only occur in Dinantian carbonate platforms below 2 km depth, having sufficient permeabilities due to karstification. We assumed that the carbonate platforms in the Dutch subsurface have similar properties. However, the tectonic setting, sea level and burial/diagenesis history vary for the areas, resulting in different platform geometries and reservoir properties [START_REF] Lipsey | Numerical modelling of thermal convection in the Luttelgeest carbonate platform, the Netherlands[END_REF].

Since our forward model is purely conductive, we approximated the thermal effect of hydrothermal convection by varying the thermal conductivity of the platforms and the overlying layers. We applied this method based on the fact that long term thermal effects can be approximated by a pseudo-convective approach, where the convective layer is marked by a higher than natural thermal conductivity (e.g. [START_REF] Beglinger | Tectonic subsidence history and source-rock maturation in the Campos Basin, Brazil[END_REF][START_REF] Luijendijk | Thermal state of the Roer Valley Graben, part of the European cenozoic rift system[END_REF][START_REF] Schmeling | A scaling law for approximating porous hydrothermal convection by an equivalent thermal conductivity: theory and application to the cooling oceanic lithosphere[END_REF]. A similar approach was applied by [START_REF] Békési | Subsurface temperature model of the Hungarian part of the Pannonian Basin[END_REF] for the Pannonian Basin. They showed that regional fluid flow systems have major influence on the temperature field, resulting in large model uncertainties. We established an a-priori conductivity profile for the carbonate platforms to fit the temperature profile in the LTG-01 well. We allowed prior conductivities to vary in certain bounds in order to obtain a satisfying fit with temperature measurements at shallower depth. Our approach results in relatively high modeled temperatures at larger depth. These show a good correlation with the locations of the carbonate platforms (Figure 8, 4-6 km depth) and could highlight potential target areas for deep geothermal exploration. The application of a conductive model to convective areas has several limitations. For instance, convective cells cannot be properly modeled with conduction, as their extent and geometry cannot be taken into account. The pseudo-convective approach results in larger modeled temperatures throughout the whole convective layer. However, the structure of hydrothermal convection is much more complex (e.g. [START_REF] Guillou-Frottier | Structure of hydrothermal convection in the Upper Rhine Graben as inferred from corrected temperature data and basin-scale numerical models[END_REF]. The temperature distribution inside convection cells is not homogenous: lower temperatures are associated with downwellings, and warmer areas correspond to upwellings. Additionally, hydrothermal convection is more likely to occur locally, and convective cells may extend beyond the carbonate platforms [START_REF] Lipsey | Numerical modelling of thermal convection in the Luttelgeest carbonate platform, the Netherlands[END_REF]. Therefore, our model overestimates temperatures at the depth of the carbonate platforms without the presence of fluid convection and at downwelling zones. Also, the thermal effect of convection beyond the platform areas cannot be captured. Our model is only able to indicate potential locations where higher temperatures might be present. Successful development of the Dinantian carbonates requires site-specific studies taking into account the local geometry of the carbonate platforms and incorporating convective heat transfer.

It is important to note that the LTG-01 well is the only indication for potential convection in the platforms. Another explanation of the thermal anomaly might be the large thermal conductivity contrast between the Numerian and Westpalian sediments and the Dinantian carbonate platforms [START_REF] Veldkamp | Temperature modelling of the Dutch subsurface at the depth of the Dinantian[END_REF]. Our model with a modified thermal conductivity values aims for approximating fluid convection, but the modified thermal conductivities could also account for the case of a strong thermal conductivity contrast. Without the observations from this well, our modeling approach would have most likely been different, with less focus and constraints on the carbonate formations. It demonstrates the need for deep exploration wells that can not only provide more information on deep geological structures, layers, and reservoir properties, but also on the deep thermal structure. The restricted amount of information available at larger depths, leads to a significant increase in uncertainty with depth for our thermal models.

The previous temperature model of the onshore Netherlands was constructed by [START_REF] Bonté | Subsurface temperature of the onshore Netherlands: new temperature dataset and modelling[END_REF]. We incorporated their work in our study and we added several new aspects to the modeling. To improve our model calibration results, we updated the temperature database with more recent temperature measurements, including production temperatures from geothermal wells. We constructed our sedimentary layers based on the sedimentary model of the onshore Netherlands (DGM-deep v4.0). Furthermore, we added new layers corresponding to the Dinantian (Carboniferous) and older Palaeozoic sediments. By taking into account the thermal effect of convection in the Dinantian carbonate platforms, we reproduced the thermal anomaly observed in the LTG-01 well. However, our model only solves the heat equation in steady state, without taking into account transient effects of vertical tectonic motions and paleo surface temperature fluctuations. Vertical motions originate from sedimentation/erosion and large-scale deformations such as lithospheric stretching. [START_REF] Bonté | Subsurface temperature of the onshore Netherlands: new temperature dataset and modelling[END_REF] incorporated vertical motions for the latest stage of basin evolution from 20 Ma until recent, originating from the accumulation of sediments. Since large-scale lithospheric deformations with significant temperature perturbations have no effect on the present-day temperatures in the Netherlands (major tectonism associated with the stretching of the lithosphere and inversion took place over 100 My and ~65 My ago (e.g. [START_REF] Van Wees | Probabilistic tectonic heat flow modeling for basin maturation: Assessment method and applications[END_REF], we concluded that a steady-state temperature model would suffice for the Netherlands.

Our model show a good correspondence with the German temperature maps [START_REF] Agemar | Subsurface temperature distribution in Germany[END_REF] along the eastern border of the Netherlands between 2 and 4 km depth. Large temperature variations (up to ±10 °C) can be observed in 1 to 3 km depth (Figure 8) due to the presence of the highly conductive evaporites within the Zechstein Group (layer 8) mostly in the northern part of the country, where salt diapers are present [START_REF] Geluk | Stratigraphy and tectonics of Permo-Triassic basins in the Netherlands and surrounding areas[END_REF]. The top of the salt layer is marked by high temperatures, for instance in the Friesland Platform, Lauwerszee Trough, Groningen High and Lower Saxony Basin, whereas isotherms are depressed below the highly conductive salts (e.g. Figure 10). The same thermal effect is observed in the westernmost part of the North German Basin (Lower Saxony) [START_REF] Agemar | Subsurface temperature distribution in Germany[END_REF] and in the in the central Polish Basin [START_REF] Zielinski | Thermal effects of Zechstein salt and the Early to Middle Jurassic hydrothermal event in the central Polish Basin[END_REF]. Elevated temperatures above salt layers provide potential target areas for geothermal exploitation in the northeastern part of the Netherlands [START_REF] Daniilidis | Salt intrusions providing a new geothermal exploration target for higher energy recovery at shallower depths[END_REF]. [START_REF] Daniilidis | Salt intrusions providing a new geothermal exploration target for higher energy recovery at shallower depths[END_REF] show that the thermal anomaly is not only dependent on the thickness of the salt layer, but the shape of the salt intrusion also plays an important role. Detailed site-specific studies taking into account the local geology including the potential for a suitable aquifer are necessary to exploit the geothermal potential of salt bodies in sedimentary basin settings. Temperature anomalies attributed to the Silesian sediments with low thermal conductivity can also be observed throughout the Netherlands. These sediments are predominantly shales with relatively high coal content, resulting in an insulating thermal effect (e.g. Figure 10). Modeled temperatures above 3 km depth are lower than average observed temperatures in the West Netherlands Basin, Roer Valley Graben, and Central Netherlands Basin, in agreement with the earlier findings of [START_REF] Bonté | Subsurface temperature of the onshore Netherlands: new temperature dataset and modelling[END_REF]. These pronounced negative anomalies can be explained by the large thickness of the Silesian sediments. The Dinantian carbonate platforms marked by high conductivity appear as positive temperature anomalies at larger depth. Additionally, the carbonate platforms in the south are identified as areas with larger than average temperatures, except for the Limburg High. The coldest areas below 3 km depth correspond to the Zeeland High and Limburg High and are likely caused by the absence of a thick sedimentary cover on top of the highly conductive basement rocks.

Different horizontal mesh resolutions may influence the calculation of thermal properties, resulting in slight variations in predicted temperatures (Figure 9). This effect is negligible in areas where the thermal field is dominantly conductive and no steep temperature gradients are present (Figure 9). On the other hand, within the highly conductive layers including the Dinantian carbonate platforms where convection is approximated by a higher than normal thermal conductivity, modeled temperatures are dependent on the horizontal mesh resolution (Figure 9). The posterior model with a finer horizontal discretization can better capture the variations in the input geological model, providing more reliable results and an improved fit with measured temperatures where steep gradients occur. We did not test different vertical resolution of the models. [START_REF] Kaiser | 3D coupled fluid and heat transport simulations of the Northeast German Basin and their sensitivity to the spatial discretization: different sensitivities for different mechanisms of heat transport[END_REF] concluded that horizontal mesh resolution has a more significant effect on modeled temperatures. They suggested that effects on temperatures modeled with conduction are almost independent from mesh resolution, however, discrepancies may occur at steep gradients or high thermal conductivity contrasts. Our models with a different horizontal discretization agree with their findings. Since the thermal field of the onshore Netherlands is dominantly conductive, we can conclude that both horizontal mesh resolutions are able to approximate the thermal field sufficiently. On the other hand, at locations where fluid convection exists, discrepancies in modeled temperatures may arise.

Calibrating the model with ES-MDA, reduces misfits of the prior model, although the variation in the data is not entirely captured by the ensembles (Figure 11). The required variation within the measurements could be reached by allowing more variation within the model parameters, especially in the shallower part of the model. On the other hand, a larger variation within these thermal properties would not be realistic. We only allowed unrealistic thermal conductivity values in the carbonate platforms and in the Limburg group (overlaying layer) in order to approximate the thermal effect of convection in the deeply buried carbonates. For the rest of the layers, we chose the scaling parameters to ensure that the thermal conductivities remain realistic. It was important especially in the shallow part of the model, where large amounts of measurements are available, many of them having errors up to ±15 °C. Allowing a wider scaling range would provide better fit with the temperature data, although the predicted temperatures would vary significantly in a small spatial proximity, resulting in unrealistic anomalies. We assume that the misfits are partly caused by transient effects or convective heat transfer. Fluid flow along fault structures are commonly identified as a cause of misfits of conductive models (e.g. [START_REF] Freymark | The deep thermal field of the Upper Rhine Graben[END_REF]. Both the prior and posterior models show a remarkable positive misfit at shallow depth: temperatures are systematically overestimated at 1 and 2 km depth (Figure 11, Table 3). The most pronounced misfits are observed in the shallowest depth interval (Table 3, 1000 m ± 400 m). We explain these misfits with the transient thermal effect of the paleo-temperature fluctuations (e.g. [START_REF] Donders | Strong climate coupling of terrestrial and marine environments in the Miocene of northwest Europe[END_REF][START_REF] Verweij | Reconstruction of basal heat flow, surface temperature, source rock maturity, and hydrocarbon generation in salt-dominated Dutch Basins[END_REF]. We attempted to account for the thermal effect of recent glaciation by choosing a lower surface temperature as a boundary condition, although we cannot entirely reproduce it with our steady-state model. Ter [START_REF] Ter Voorde | Weichselian and Holocene climate history reflected in temperatures in the upper crust of the Netherlands[END_REF] also concluded that the misfits with the steady-state thermal profile in shallow depth reflect a transient condition inherited from past climate change. Alternatively, groundwater flows can also explain local anomalies at shallow depth. Measurements errors are also partly responsible for misfits, especially at 1-2 km depth, where a large number of data points are available in close spatial proximity (e.g. in the WNB, CNB-NHP border). Systematic small under-or overestimates might be caused by boundary conditions and/or heat generation in thick layers that have a larger lateral extent (e.g. upper crust or a sedimentary layer that is widespread across the Netherlands. Other local misfits might be caused by uncertainties in model geometry, lithology mixtures (e.g. thermal properties of the sedimentary layers), etc.

Conclusions

We established a 3D high-resolution subsurface temperature model of the onshore Netherlands. One of the most important aspects of this study is the validation of the thermal model with 1507 temperature measurements. We calibrated our model with temperature observations through inversion with ES-MDA. We took into consideration both the data and model uncertainties by assuming a Gaussian distribution for measurement errors and a triangular distribution for scaling the thermal properties. Misfits of the prior model are reduced through the data assimilation procedure: the overall RMS, including all the measurements, decreases from a value of 0.95 °C for the prior model to 0.70 °C for the posterior model. It demonstrates the effectiveness of ES-MDA as a tool for calibrating temperature models, supporting high-resolution external constraints. On the other hand, for areas without temperature data available, especially in the deeper parts of the model, predicted temperatures are strongly dependent on conceptual constraints. Therefore, a reliable geological model and reasonable thermal properties are crucial as modelling input.

By taking into account the thermal effect of convection in the Dinantian carbonate platforms, we reproduced the thermal anomaly observed in the LTG-01 well. Our model reveals areas with potential for hydrothermal convection in the deep carbonate platforms. These locations can be suitable for deep geothermal development of both electricity generation and direct heat uses due to the sufficiently high temperatures and inferred high permeabilities that are required for convection. The temperature model has been incorporated into the updated ThermoGIS project and is available online at thermogis.nl

Figure 2 .

 2 Figure 2. Cross section through the Netherlands showing the the main sedimentary units to the base of the Namurian, adopted from the Netherlands Oil and Gas Portal website (www.nlog.nl). Fault structures are not shown in the section as they were not incorporated to the thermal model. The tectonic evolution of the Netherlands described above has led to a situation where the reservoirs that have been targeted by the oil and gas exploration for more than 70 years are the same a currently targeted by geothermal exploration now. If earlier unsuccessful attempts like the Asten geothermal well of 1988 are ignored, the exploration for geothermal energy started in the Netherlands around 2006, when the first doublet was drilled.

Figure 3 .

 3 Figure 3. Temperature dataset including BHTs, DSTs, RFTs and production temperatures from geothermal wells modified after Bonté et al. (2012) (a). BHT_ICS: BHTs corrected with analytical method, BHT_AAPG: BHTs corrected with statistical method. The temperature measurements were collected from wells drilled in the onshore Netherlands (b). In case of multiple measurement types in a single well, the one with the lower uncertainty is shown in (b). The misfit of raw and corrected BHT gradients with the average Dutch geotherm are shown in (c), with the mean and standard deviation of the fitted Gaussian functions. 4. Conceptual Model 4.1. Forward model

Figure 4 .

 4 Figure 4. Depth of the top of the Carboniferous Limestone Group(Dinantian). The extent of the carbonate platforms within the group based on seismic data and well logs, adopted from[START_REF] Kalkman | Ultra-Diepe Geothermie -Casus Almere en Barendrecht[END_REF], is outlined. We assumed that convection may occur in the deeper platforms, marked by the thick black outline.

Figure 5 .

 5 Figure 5. The dependence of minimum Rayleigh permeability on depth. The curves were constructed by taking into consideration the temperature, pressure, and salinity dependency of water properties following Van Wees et al. (2012). Note the different scales for depth and k min .

Figure 6 .

 6 Figure 6. Map showing the locations of the Dinantian carbonate platforms and other rocks of Dinantian age modified after[START_REF] Kalkman | Ultra-Diepe Geothermie -Casus Almere en Barendrecht[END_REF]. We use scaling of the thermal conductivity to approximate the thermal effects of convection. Outlined in black are the areas below the carbonate platforms for which we scaled the thermal conductivity values.

Figure 7 .

 7 Figure 7. Description of the modeling work flow.

Figure 8 .

 8 Figure 8. Temperature maps of the final posterior model (step 7 in Figure7) at 1-6 km depth. The misfits between modeled temperature and measurements (modeledobserved values) within a ±200 m interval are plotted with circles. 80% of the misfits lie within the ±10 °C interval, where the remaining 20% are dominantly attributed to measurements with the largest uncertainty (±15 °C, BHTs corrected with statistical method). Note the different color scale for the modeled temperatures at various depth.

Figure 9 .

 9 Figure 9. Temperature-depth (black) and thermal conductivity-depth (blue) profiles of the prior model (dashed lines) and posterior model (solid lines) at different well locations. The temperature measurements from the wells are marked by green dots with corresponding error bars. Measurements from wells within 3 km distance are plotted with green crosses. The numbering of the layers corresponds toTable 1. Note that the depth ranges of

Figure 11 .

 11 Figure 11. Distribution of the modeled and observed temperatures ranked by the modeled temperatures at the observation points for the prior (a) and posterior (b) models. Temperature measurements are marked by black dots. The red curves indicate the model results (in (b) the P50 values), the green lines in (b) are the P10 and P90 values. The blue interval in (b) indicates the modeled temperatures for all the 600 ensemble members of the last iteration.

Table 1 .

 1 Lithological composition and prior thermal properties of the layers.

	2	Lower Middle North and Sea Group	Paleogene		50 % Shale -50 % Sandstone
	3	Chalk Group		Upper Cretaceous 100 % Chalk
	4					50 % Shale -25
		Rijnland Group	Lower Cretaceous	% Marl -25 %
						Sandstone
	5	Upper Jurassic Groups	Upper Jurassic	50 % Shale -50 % Sandstone
	6	Altena Group	Lower Middle Jurassic and	75 % Shale -25 % Silt
	7	Lower Upper Germanic Groups	and	Triassic		40 % Shale -30 % Limestone -30 % Sandstone
	8					Zechstein layer
		Zechstein		Lopingian (Late	(Shale to the
		Group		Permian)		South -Salt to the
						North)
	9	Upper Rotliegend Group		Late Guadalupian (Middle Permian)	70 % Sandstone -15 % Shale -15 % Silt
	10 Caumer Subgroup		Westphalian (Carboniferous, Middle Silesian)	60 % Shale -40 % Sandstone
	11	Limburg Group	Namurian Lower Silesian) (Carboniferous,	% Sandstone 60 % Shale -40	1.5
	12 Carboniferous					5
		Limestone Group carbonate	-	Dinantian (Lower Carboniferous)	100 % Limestone
		platforms
	13 Carboniferous	Dinantian (Lower	100 % Limestone Pressure-and
		Limestone		Carboniferous)			temperature-
		Group						dependent
	14 Old Paleozoic sediments	Paleozoic older) (Devonian	and	60 % Shale -40 % Sandstone	based lithotypes Hantschel and on after
								Kauerauf
								(2009)
	-1 Upper crust		-		100 % Granite		Pressure-and	-	1
	-2 Lower crust		-		Anorthosite 100	%	temperature-after dependent	-	0.5
								Chapman
								(1986)
	-3 Lithospheric		-		100 % Peridotite	Pressure-and	-	0.02
		mantle						temperature-
								dependent
								after Schatz

Table 2 .

 2 Overview of the data assimilation procedure.

	Modeling	Description	Prior	Scaling	Variogram	Variogram	Number of
	step		distribution	parameters	model	range	ensemble
							runs and
							iterations
	3	update the	Triangular	A: 0.6 and 1.6	Spherical	14 cells	600
		heat		for layer -1		(~42 km)	ensembles,
		production					4 iterations
		(A) in the
		upper crust
	4	update the		λ V : 0.7 -1.3 for		8 cells	600
		vertical		layer 1-7, 9-10,		(~24 km)	ensembles,
		thermal		14			4 iterations
		conductivity		λ V : 0.7 -1.4 for
		(λ V ) in the		layer 8
		sediments		λ V : 0.6 -1.4 for
				layer 12-13
	7	update		λ V : 0.8 -1.2 for		30 cells	600
		(refine) λ V in		layer 1-13		(~30 km)	ensembles,
		the					4 iterations
		sediments

Acknowledgements

The research leading to these results has received funding from the European Union's Seventh Framework Programme under grant agreement no. 608553 Project IMAGE. We also thank the reviewers for their constructive comments that helped us to improve this paper.

Data availability

Temperature data that we use for the calibration of our model are available at the Netherlands Oil and Gas Portal website (www.nlog.nl). The DGM-deep v4.0 geological model representing the main sedimentary units of the Netherlands and used for our modeling input is published on www.nlog.nl. Our final temperature model is available online at thermogis.nl.