Towards optimisation of geothermal heat recovery An example from the West Netherlands Basin

1% of the potentially recoverable heat will be recovered by is of in-efficient doublet deployment on a ‘first-come, first basis with operational parameters that focus on ob- jectives of small decentralised heat grid demands. Instead, similar to the common-practise approach in the hydrocarbon industry, a regional coordinated ‘masterplan’ approach could be used to increase heat recovery. Utilising numerical simulations for flow and heat transfer in the subsurface, we showed that the heat recovery efficiency could be increased by tens of percentages with such coordinated doublet deployment. Based on cal- culations of the Levelized Costs Of Heat for both deployment strategies, we also show that current financial support schemes do not favour heat recovery optimisation. This study emphasises that although Hot Sedimentary Aquifer resources have the potential to cover a significant part of our energy demand, a radical change in financial support schemes and legislation are required to unlock their true potential.


Introduction
Huge amounts of heat are stored in Hot Sedimentary Aquifers (HSA) in regions all over the world (e.g., [1][2][3][4]). These resources are especially suitable for district heating applications and could play a significant role in our future energy mix when integrated into district heating networks or industrial processes (e.g., [5,6]). With an increasing demand for low-carbon heat, ambitious targets are presented to enhance exploitation of these under-developed resources (e.g., [7][8][9]). It is, however, unclear if these targets are realistic because it is unknown how much heat could be recovered from aquifers. The Recovery efficiency, i.e. the ratio of the total recovered energy and the estimated energy in the resource or recovery factor, is a key performance indicator in hydrocarbon production [10] but rarely mentioned in relation to geothermal energy [11][12][13][14][15][16]. It is mainly discussed in HSA feasibility studies as conjectures [2,[17][18][19] or in relation to hydrocarbon and heat co-production (e.g., [20][21][22]). So far, it has not been an objective of geothermal exploitation strategies. The reason for this is that geothermal resources are often exploited to provide heat to local decentralised grids, instead of optimising heat recovery from the entire resource [23][24][25]. Geothermal resources are often exploited by multiple independent operators with potentially conflicting objectives [26,27].
As a result, exploration licences are issued on a 'first-come, first-served' basis as legislators have no desire to be held responsible for future failure of the individual projects' objectives. Alternatively, doublet deployment could be based on a regional masterplan that aims to optimise heat recovery of the entire HSA resource (e.g., [26,28]), similar to the common-practise approach in the hydrocarbon industry [29][30][31].
This study compares the impact of those two contrasting approaches in doublet deployment on the recoverable heat from (henceforth recovery efficiency) low-enthalpy HSA resources. A case-study is presented from the West Netherlands Basin (WNB) in the Dutch province of Zuid-Holland with an overview of key-parameters of the on-going geothermal exploitation. Based on trends of the geothermal expansion rate, predictions were made on future annual heat production, the recovery efficiency and available production license space in the WNB. In the second part of this study, we compare results of flow and heat transfer simulations for the WNB aquifer exploitation with the current 'first-come, first-served' doublet deployment and the 'masterplan' approach to doublet deployment.
The WNB is selected as the study area because of the availability of geothermal operational data, studies on regional heat demand [32][33][34], as well as extensive geothermal potential studies and Heat in Place (HIP) evaluation [17,35]. Moreover, it was the epicentre of geothermal development in the Netherlands, the country with one of the fastest expansions of geothermal development in Europe in the past decade [36,37]. Current geothermal doublets are operated by independent operators and mainly provide heat for the large greenhouse sector. The main aquifer targets are formed by fluvial sandstones of the Lower Cretaceous Nieuwerkerk Formation at depths between 2 and 3 km (e.g., [38][39][40]). The estimated HIP of this formation is 52,000 PJ [33], which could potentially cover the regional heat demand of 120 PJ/yr the province of Zuid-Holland [41] for decades with a low carbon footprint [42]. After the first doublets in the basin confirmed this significant geothermal potential, a next challenge is to make sure this resource is exploited efficiently increasing the role of geothermal heat production in a future low-carbon energy mix. In Section 2 of this paper, the geothermal development in the WNB is evaluated and predictions for future expansion are made. Section 3 describes how numerical production simulations are employed to explore the scope of optimisation of heat recovery from the main geothermal resource in the WNB. The results of the numerical simulations are presented in Section 4 and discussed in Section 5.

Exploitation between 2007 and 2018
From 2007 to 2018, 13 doublet systems have been realised in the WNB (Fig. 1). Challenges for efficient doublet deployment in the future are evident from this figure. Available space for new operators is limited by the extend of the main aquifer target, a mandatory distance of faults, hydrocarbon fields and the irregular shape and layout of licences that have been granted to independent operators. An overview of key parameters of this development is presented in Table 1. These parameters include doublet capacity in MW th , aquifer name, area of influence, combined screen length as a proxy for net-aquifer thickness, injector-producer spacing, production rates and production temperatures and start date. The parameters values were derived from 'End-Of-Well' reports from the publicly available database Nlog.nl [32], and from the database of the Netherlands Enterprise Agency [33]. Most doublets in the WNB exploit aquifers of the Lower Cretaceous Nieuwerkerk Formation. The True Vertical Depth (TVD) ranges from 1850 to 2480 m (Table 1). An increase of doublet production rate could be recognised in newer doublets. Doublets that have been realised before 2015 were designed for production rates of up to 200 m 3 /h, whereas newer doublets operate at rates between 250 and 360 m 3 /h. As a result, the capacity of the doublet increased up to 18 MW th . Injector-producer spacing (L) ranges from 1 to 2.1 km with an average of 1.6 km. L is derived from 'End-Of-Well' reports or interference test reports. We derived L as from the horizontal distance between top of the production and injection liners of the doublets was used to determine L. The values are provided in kilometres and rounded off to 1 decimal. A course indication of injector-producer spacing is indicated here, whereby we neglected the impact of the slight tilting of the aquifer, and the increasing L along the wellbore because the wells are deviated with a 30-50°angles [32]. The tilting angle of the aquifer could be estimated assuming a 3.5 m/ms seismic velocity at 2.5 km depth [43] and recognising ∼250 ms Two Way Travel-time difference of the top of the Nieuwerkerk over 10 km in one of the WNB fault blocks from Willems et al. [40]. Such a TWT difference would results in some tens of meters over a typical injector producer spacing of 1-2 km. Furthermore, exact calculation of the injector producer spacing at aquifer level is complicated by uncertainty in aquifer thickness of several tens of meters [40]. The associated area of influence (A) that captures the expected extent of the cold-water plume at the moment of cold-water breakthrough is approximated by: In the WNB, A ranges between 2 and 7.2 km 2 , with an average of 5.4 km 2 . The associated energy production per square kilometre is 1.1-4.7 MW/km 2 , with an average of 2.4 MW/ km 2 . Note that for recently drilled doublets not all parameters could be listed. This is because data was not publicly available, or doublets are still in a test/workover phase. Values in italic indicate that they are derived from prognoses instead of measurements.

Future expansion of WNB geothermal exploitation
Between 2007 and 2018, the average expansion rate of the installed capacity was 13 MW th /yr ( Fig. 2A). In the last four years, an increased expansion rate could be recognised of 24 MW th /yr. With these expansion rates, the total installed capacity could be either 540 or 1000 MW th in 2050 ( Fig. 2A). Subsequently, the potential annual heat production in the basin is determined by integrating the total installed capacity in MW th over time and assuming an annual downtime of 20%. The associated annual energy production will be approximately 18 for the expansion rate of 13 MW th /yr and 30 PJ for the expansion rate of 24 MW th /yr in 2050 (Fig. 2B). Schoof et al. [9] defined a target of 50 PJ/yr of geothermal heat production in the Netherlands in 2030. This target could only be met if geothermal exploitation both within and outside the WNB would expand at significantly higher rates, which are shown in Fig. 2A. The recovery of the available geothermal heat is expressed in terms of the recovery efficiency (R, Eq. (1)) and in the Heat Demand Coverage (HDC, Eq. (2)), the percentage of the regional heat that is covered by geothermal exploitation. The recovery efficiency is derived from dividing cumulative heat production (E cu ) by 52,000 PJ, the estimated HIP in the main Lower Cretaceous aquifer target by Kramers et al. [35]: HDC is derived by dividing the total combined installed capacity of all doublets (C T ) [MW th ] by the total annual heat demand (H D ) [MW th ] in the province of Zuid-Holland which was 120 PJ/yr in 2016 according to [41]: In 2018, HSA exploitation provided approximately 4% of Zuid-Holland heat demand (Fig. 2C). In 2050, this could increase to approximately 15-25% depending on the expansion rate trend. Even in the more optimistic expansion trend of 24 MW th /yr, some 0.75-1.2% of the HIP would be recovered (Fig. 2D). Fig. 2C and D underline that a different doublet deployment strategy is required to increase the percentage of geothermal heat in a future low-carbon energy mix and to optimise heat recovery.

Available space for production licenses in the WNB
When the installed capacity will continue to grow with 13 MW th /yr, some 50 doublets will have been realised by 2050. Assuming an that the average injector-producer spacing of the doublets remains 1.6 km (Table 1), it was calculated that these 50 doublets will cover some 225 km 2 . This is approximately 20% of the surface area below which Lower Cretaceous aquifers are expected to have sufficient potential according to Pluymaekers et al. [44]. When the expansion rate of geothermal exploitation is faster, for example 24 MW th /yr like in the past four years, some 80 doublets can be realised by 2050, occupying some 35% of the surface area below which the aquifer has sufficient potential. Note that in both expansion scenarios, not even 1% of the HIP will be recovered in 2050 (Fig. 2D) and pressure interference will become increasingly imminent [45]. Therefore, Fig. 3 shows that a change in doublet deployment and permitting strategy is required to increase the number of doublets and heat recovery efficiency in this basin.

Study area, aquifer model, heat demand
In a second phase of the study, numerical production simulations are employed to compare the recovery efficiency and HDC for different   (Fig. 4). For the sake of simplicity, the aquifer is modelled as a homogeneous and horizontal sandstone layer with constant thickness. For the development of an optimised exploitation strategy, reservoir heterogeneity for flow properties would have to be considered [13], but this is out of the scope of this study. The model consists of a horizontal 100 m thick homogeneous sandstone layer, which is confined between 300 m impermeable over-and under-burden layers providing thermal recharge (Fig. 4). The thickness is derived from averaging aquifer thickness of currently active doublets in the WNB. Two of the four boundaries of the model are formed by the northwest to southeast trending faults, derived from the WNB structural analysis of Duin et al. [46]. The aquifer properties are assumed isotropic and values are derived from Willems et al. [45] ( Table 2). Element size ranged from 0.3 m to 40 m in the aquifer layer and from 40 m to 300 m in the over-and under-burden layers. The total heat demand in the surface area covered by the aquifer model is approximately 80 MWth (Fig. 4). The aquifer model has a surface area of 24.5 km 2 and the average surface heat demand of 3.8 MW th /km 2 (based on heat demand values for horticulture and domestic heating for both municipalities derived from the Dutch National Energy Atlas [34]). The heat demand in our study area is higher than the current geothermal heat production per square kilometre in the province of Zuid-Holland of 2.4 MW th /km 2 .

Doublet deployment scenarios
Heat exploitation from the introduced aquifer model is considered by utilising three doublet deployment scenarios. In the first scenario, 3 doublets are deployed. This scenario represents the currently applied 'first come, first served' deployment strategy. Two doublets are exploiting the Lower Cretaceous aquifer in the Pijnacker fault block. A third exploration license is granted where a new doublet is planned in 2019 [47]. The injector-producer spacing in this scenario ranged from 1200 m for doublet 3, to almost 1700 m for doublet 1 (Fig. 5). Deployment scenarios 2 and 3 represent optimised deployment scenarios according to a 'masterplan' strategy. In scenario 2, nine doublets are deployed with an injector-producer spacing of 800 m. Following Willems et al. [45], the optimal distance (dx) between individual doublets was 1040 m. Doublets were placed in a so called 'checkerboard' configuration whereby adjacent wells have opposite functions. In scenario 3, sixteen doublets are deployed in checkerboard configuration, but here the distance between individual doublets is equal to the injector-producer spacing, which is 800 m. For simplification, the production rate of all doublets is equal to 200 m 3 /h and production commenced at the same moment. These simplifications are made because the aim is to show that the recovery efficiency can be improved with optimised doublet deployment. We do not try to develop an optimised deployment strategy. To recognise interference in the doublet deployment scenarios 1, 2 and 3, two additional simulations with a single doublet and injector-producer spacing of 800 m and 1650 m, respectively, are conducted.

Thermal flow modelling
Numerical production simulations are used to calculate the net energy production and life time of the doublet in the different deployment scenarios. In the numerical modelling procedure, the energy balance is solved for a rigid medium fully saturated with water, with thermal equilibrium between the fluid and solid phases: is the temperature, ρ [kg/m 3 ] and C i [J/(kgK)] are the temperature independent mass density and specific heat capacity, respectively. Subscripts reference to the rock -(r) and brine (w) phase. λ [W/(m/K)] is the thermal conductivity tensor q = −Kµ∇•h is the Darcy velocity vector, with hydraulic conductivity = K k g r w , where k r [m 2 ] is the sandstone or shale permeability, g gravitation acceleration and µ the temperature dependent viscosity. The thermal conductivity tensor is calculated through: λ = (λ eq + (α T )|q|)I + ρ f C f (α L − α T ) qq/|q|. In this equation, |q| is the magnitude of the Darcy velocity vector, I is the identity matrix, α L and α T are the thermal dispersion coefficients in the longitudinal and transversal direction, respectively. The approach of Diersch [48] is followed using the empirical viscosity temperature dependency of Mercer and Pinder [49]. The equivalent heat conductivity, density and the volumetric heat capacity are assumed to be independent of temperature for simplicity and described by: λ eq = (1 − φ)λ r + φλ w and ρC = (1 − φ)ρ r C r + φρ w C w in which φ is the porosity. The pressure field is obtained through solving the continuity equation: is external sinks and sources. The detailed modelling procedure follows the approach explained in Saeid et al. [50,51]. The values of all constant parameters are listed in Table 2. The production simulations yield a production temperature development over time and the required injection and production pressure for the associated production rate and set of parameters. A constant temperature boundary condition of 75°C was assigned to the top and base of the aquifer model. Another temperature boundary condition is assigned to the nodes of the injection well trajectories of 35°C. The two long edges of the model representing the faults flanking the aquifer model are noflow boundaries. The short edges of the model are open-flow boundaries. Fig. 4 presents an overview of these boundary conditions.

Net energy production
The production simulations yield a production temperature development over time and the required injection and production pressure for the associated production rate. The difference between injection and production pressures (ΔP) is used to estimate pump energy losses: , [52]), where Q is the production rate and ε the pump efficiency of 60%. The produced power (E prod ) is estimated by: In Eq. (3), T inj is the injection temperature of 35°C and T prod is the production temperature that changes over time after thermal breakthrough. The net energy (E net ) production is determined by the sum of the produced energy and the pump energy losses. Integration of E net over the production time gives the cumulative energy production. In this calculation, the net energy production is reduced by 20% to account for maintenance down-time. To obtain the recovery efficiency, the cumulative produced heat is divided by the HIP of the aquifer model, following Eq. (1). This HIP is estimated by where V is the bulk rock volume of the active aquifer and C sst is the heat capacity of sandstone (sst) ( Table 2).

Levelized Cost of Heat
The Levelized Cost of Heat (LCOH) is determined for each doublet deployment scenario, using Eq. (4): In this equation, CAPEX t and OPEX t are the respective total Capital and Operational expenses in year t, which are derived from n the parameters listed in Table 3, r is the discount rate and E net,t is the cumulative generated amount of heat in year t, which is derived from the production simulations. A fifteen-year period is chosen because this relates to the maximum duration of the Dutch feed-inn tariff subsidies  [44]. [53]. ESP pump replacement costs have been considered for each injection well in the sixth and eleventh year. Using the LCOH as performance indicator, the financial advantage of optimisation of heat exploitation is quantified for each doublet deployment scenario and compared to the LCOH of conventional gas heating [54].
In our CAPEX calculations we assume that if a single operator exploits the entire aquifer model, this would require an investment for a surface heat distribution grid for selling the produced heat. We neglect these costs for operators with only a single doublet because we assume that in this scenario the geothermal doublet provides heat to an existing decentralised grid. We assume a simplified heat grid that connects all doublets in the three scenarios with associated heat grid costs (C grid ) estimated by: C grid = (N − 1)·dx·C netw . In this equation N is the number of doublets, dx [m] is the doublet distance in each scenario, and C netw is the network costs of 1000 €/m following Daniilidis et al. [56] (Fig. 6). Our heat grids only connect the doublet wellheads (i.e. the surface termination of a wellbore). We neglect additional costs for connecting the grid to heat customers because that would require additional heat consumer analyses, which is out of the scope of this study.

Production temperature development and doublet life time
Results of the numerical production simulations show that with deployment scenario 1, it could take more than some 140 years before the production temperature will drop below uneconomic temperatures for space heating (70°C, [1]). In one of the doublets, the production temperature drops faster than that of the single doublet reference simulation due to its proximity to the north western model boundary. Thermal breakthrough occurs after approximately 90 years in scenario 1. After thermal breakthrough, the net energy production reduces. In scenario 1, it takes between 140 and 200 years for the capacity to drop by 1 MW th . In deployment scenarios 2 and 3, the smaller injector-producer spacing of 800 m and closer doublet distance reduce the thermal breakthrough time significantly to approximately 25 years. Nevertheless, the production temperature does not drop below the minimal required temperature for 45°C for greenhouse heating (T min greenhouse, [1]) in 200 years in scenario 2 and 100 years in scenario 3. The minimal required temperature for space-heating of 70°C is reached in ∼40 years in scenario 3 and in 50 years in scenario 2. The results in Fig. 7 show that not only well spacing (L) but also the doublet distance (dx) affects the life time of the doublets (see Fig. 6 for comparison of L and dx). In scenario 2 and 3, the doublets have the same injector producer spacing (L) of 800 m, but the distance between each doublet (dx) is different. In scenario 3, 16 doublets are deployed with an average doublet distance of 800 m. In scenario 2, 9 doublets are deployed, and the average distance between each doublet is 1040 m. Because of the smaller doublet distance, the production temperature and the associated capacity in MW th reduces faster after thermal breakthrough in  deployment scenario 3. In scenario 2, the capacity of the doublets is reduced with 50% in 100 years. In scenario 3, the capacity of the doublets reduces to 50% of their original value after some 60 years. This shows that, geothermal production licences that target the same aquifer cannot be regarded as stand-alone features because interference is inevitable when multiple operators exploit the same aquifer. The spread of the curves in scenario 2 and 3 is a result of slight variations in well spacing and doublet spacing of up to ± 50 m and the proximity of some doublets to no-flow boundaries.
Higher production rates will increase doublets capacity, but also result in earlier thermal breakthrough and faster reduction of the production temperature thereafter. Fig. 8A-C compare the production temperature development for simulations with production rates of 200 and 400 m 3 /h. Fig. 8D-F compare the associated development of the capacity of each doublet. As a result of the more rapid reduction of production temperature, the capacity of most of the doublets with 400 m 3 /h drops below that of doublets with 200 m 3 /h after 50 years in deployment scenario 2 (Fig. 8E), and after some 40 years in deployment scenario 3 (Fig. 8F). The results in Fig. 8 highlight that an optimised exploitation strategy requires a trade-off between doublet capacity and well density as well as the minimum required production temperature.

Heat recovery efficiency -Pijnacker aquifer model
With the current 'first come, first served' deployment strategy (scenario 1) and doublets operating at 200 m 3 /h, geothermal exploitation could cover almost 40% of the heat demand of the surface area overlying the aquifer model (Fig. 9A). This could be increased to some 60% when the production rate of the three doublets is 400 m 3 /h. An even higher percentage of the regional heat demand could be covered by geothermal exploitation when more doublets are deployed like in scenario 2 and 3. This suggests that HSA exploitation could play a more significant role in the regional energy mix if doublets would be deployed with higher density according to a 'masterplan' approach. Fig. 9B shows that recovery efficiency could be enhanced in a similar  way. with deployment scenario 1 and 200 m 3 /h production rate in each doublet, ∼15% of the HIP could be produced in some 50 years. In simulations with scenarios 2 and 3, 40-65% of the HIP was produced in the same period of time. Because of the large and irregular well spacing in the current 'first come, first served' strategy, most of the heat in the acquirer remains untapped in scenario 1 at the moment of thermal breakthrough (Fig. 9C). The efficiency of heat extraction can be enhanced by denser doublet deployment like in scenarios 2 and 3 ( Fig. 9D and E).

Financial impact of optimised doublet deployment
The Levelized Cost Of Heat (LCOH), calculated over a 15 year period, is compared for different doublet deployment scenarios in Fig. 10. Our calculations result in a LCOH for exploitation with a single doublet of approximately 3 €ct/kWh. This is in between the assumed LCOH with 'today's technology' by Beckers et al. [57] and the Dutch gas heating price according to ECN [54]. The LCOH does not reduce for operators that exploit three doublets at the same time, like in deployment scenario 1. Only a very small reduction in LCOH could be recognised for deployment scenarios 2 and 3. Production rate increase has a much more significant impact on LCOH. The LCOH calculation for the Fig. 7. Development of production temperature in doublet deployment scenario 1, 2 and with a constant production rate of 200 m 3 /h. The injector-producer spacing in scenario 1 is derived from Table 1 and in scenario 2 and 3 this spacing is 800 m for each doublet. Each line represents the production temperature development of a single doublet. Fig. 8. Impact of production flow rate on production temperature development (A-C) and capacity (D-F), for each doublet in the different deployment scenarios (columns). Dotted lines represent simulation results with 400 m 3 /h production rates, solid lines represent production simulations with 200 m 3 /h production rates. simulations with 400 m 3 /h production rates drops to approximately 2.0 €ct/kWh for deployment scenario 1 and single doublet deployment, and 1.9 €ct/kWh for deployment scenario 2 and 3. Note that we included additional 25% CAPEX costs that would be required to upscale well diameter and surface facilities to accommodate for the increased flow rates. These calculations suggest that for current operators in the WNB no financial incentive exists that promotes upscaling exploitation. There is no significant financial advantage of the 'masterplan' doublet deployment that could potentially increase the heat recovery efficiency as was shown in Fig. 9A. For individual operators it is financially more attractive to optimise heat production of a single doublet, then to consider optimisation of the heat recovery efficiency.

Heat recovery efficiency
The pioneering first geothermal projects in the WNB proved the technical potential of HSA exploitation. The next challenge for this new industry is to optimise heat exploitation strategies to utilise the full potential of these resources. This study presents predictions of the current and future role of geothermal heat production in the region. Our results show that with the current expansion rate of the installed geothermal capacity the recovered heat will only be approximately 1% or less of the estimated HIP in 2050. This is much lower than the estimated recovery efficiency values of 20-30%, which are typically used in HSA potential studies [2,19,35,55]. Our simulations do not aim to provide an optimal exploitation strategy, but we aim to show that there is a significant scope for optimisation of the recovery efficiency of heat from HSA resources. With a regional approach to doublet deployment, significantly more heat could be extracted from HSA and geothermal heat exploitation could play a larger role in the low-carbon energy mix of the province of Zuid-Holland and similar regions world-wide than it does today. Our estimation of the recovery efficiency in the WNB could be refined by improving the HIP estimate. In calculation of the HIP in  year production period assuming 20% annual downtime (bluediamond markers) and 40% annual downtime (red, circle markers). For reference, the LCOH value of heat for direct-use applications by Beckers et al. [57] and the heat price of traditional gas heating are indicated. the WNB, Kramers et al. [35] considers neither the uncertainty of the specific heat capacity, porosity and variation of the required injection and production temperature nor the significant impact of thermal recharge on the amount of recoverable heat [58]. Especially because of the significant impact of thermal recharge, we expect that the HIP is actually much larger and therefore our recovery efficiency calculations are still conservative. Other improvements of our recovery efficiency predictions could be made by incorporating fluctuating heat demand and associated variations in production rate in estimations of the total produced heat Fig. 2B). We assume an annual down time of 20%, for example for maintenance, in calculating energy production from the capacity in MW th . This assumption is made because actual annual heat production records are not available. Heat recovery might have been slower in the first 10 years because of extra downtime for unforeseen maintenance, but because of increasing exploitation experience this might be lower in the future. Our annual heat production calculations are based on values from Table 1, which are derived from 'End of Well Reports', and the websites of the geothermal branch-organisation [59] and governmentally managed publicly available databases. Therefore, these table entries represent snapshots of performance indicators and actual values may vary over time. Also, the regional heat demand derived of the Province of Zuid-Holland [41] may be subject to development in the future, which will affect the percentage of heat that geothermal exploitation could cover (Fig. 2C).

Expansion of geothermal exploitation in the WNB
In our analysis, we assumed linear growth trends for WNB geothermal exploitation. Several factors such as increasing demand for low-carbon heat or fluctuating process of conventional gas heating could have an impact on the growth trend of geothermal exploitation in the WNB. Finally, future expansion might be decelerated by the potential impact of interference and availability of space for new production licences. On the other hand, improved understanding of technical challenges may reduce costs for future operators, enhancing financial competitiveness and therefore also the expansion rate. For example, a major current challenge is the casing steel selection. Corrosion of casing-and Electrical Submersible Pump (ESP) systems had significant financial implications in the past years [60]. Improved and consistent hydro-chemical analyses of formation water in future studies are required to improve steel selection, reduce corrosion risk and accelerate exploitation further. Because the above-mentioned simplifications may have negative and positive impact on the recovery efficiency estimate, we expect that the order of magnitude of our estimations is still realistic. This study therefore emphasizes that a radical change is required in doublet deployment rates to make HSA exploitation a more significant contributor to the regional energy mix in Zuid-Holland.

Well placement optimisation
The numerical production simulations in this study are used to explore the potential scope for heat recovery optimisation from HSA resources from a reservoir engineering perspective. HSA exploitation is simulated utilising the current 'first come, first served' doublet deployment and an optimised deployment approach. The results suggest that the heat recovery efficiency could be increased by several tens of percentages when doublet deployment is coordinated according to a masterplan and when smaller injector-producer spacing distances are considered for doublets. The current WNB spacing standard of 1.6 km is derived from exploitation in the Paris Basin, where no thermal breakthrough has been observed, not even after 40 years of exploitation [23]. In our 'masterplan' deployment we use approximately half of the current WNB spacing. This is chosen to emphasize that the current standard might lead to unnecessary long breakthrough times and reduces the possible number of doublets that can exploit the aquifer ( [45,61], Fig. 7). Optimised well spacing and doublet deployment is basin dependent and it should be derived based on detailed geological modelling. This is because the sedimentary facies architecture will reduce thermal breakthrough time when the same production rates are used but also reduce the rate of the production temperature decline after thermal breakthrough [13,61]. Optimal doublet deployment further depends on net aquifer thickness as well as the tilting of the aquifer, lateral thickness variation and aquifer background flow [45,62]. Also, logistical and financial factors must be considered related to the design of surface heat distribution networks, location of consumers as well as minimisation of drilling costs, which could affect wellhead locations at the surface. For simplification, neglected the time it would take to realise doublets. Increased insight of the geological heterogeneities will be obtained with an increasing number of wells. The masterplan should therefore be flexible enough to adopt with the future gained progressive insights on aquifer heterogeneities. In our simulations, doublets commence exploitation at the same time. In reality, it would take several months to realise doublets. This would reduce the cumulative heat production and recovery efficiency and affect the results of Fig. 9. Finally, the use of doublets may not be the most optimal method to extract heat from geothermal aquifers. Alternatively 'triplet' systems could be considered similar to the examples in the Paris Basin and the Soultz-sous-forets EGS system [23,63,64] or even more advanced injector-producer layout could be imagined. Several software tools have been developed that could aid in finding optimised well placement strategies [e.g. 65,66]. The impact of reservoir heterogeneity [13,16,67] on well placement optimisation should be also considered.

Implementation of masterplan deployment
Application of the 'masterplan' deployment approach has so far been hindered by the economic and legislative hurdles. The feed-in tariff schemes and geological risk insurance schemes made a significant contribution to initiate geothermal development in for example France, Germany and the Netherlands [23,53,68]. However, our LCOH calculations (Fig. 10) indicate that little financial incentive exists to promote heat recovery optimisation or the use of the 'masterplan' deployment approach. LCOH does not reduce for exploitation with deployment scenarios 2 and 3. If neighbouring operators would be able to significantly reduce their CAPEX by sharing surface facilities and investments for data acquisition and risk insurance, the LCOH could potentially be reduced with large scale exploitation. Also, larger-scale exploitation with multiple doublets would become attractive if the heat demand or market is far larger than the capacity of a single well. When neighbouring operators are connected to the same heat network, successful doublets could compensate unexpected poor performance of other doublets, which could be a result of local geological heterogeneities. One possible way to promote heat recovery optimisation could be through the creation of regional surface heat distribution networks. This would of course require a large upfront investment but Blom et al. [69] showed that such networks could be profitable as stand-alone infrastructure if users pay a fee for heat traffic. This type of government support was also provided by the Dutch government in the 60s to support gas exploitation after the discovery of the Groningen gas field. A national gas grid was installed connecting every household to the gas grid and ensuring long term domestic gas consumption in the Netherlands [70]. Not only financial support mechanisms but also permitting mechanisms need to be adapted to facilitate and promote heat recovery optimisation. Our production simulations show that interference is inevitable. In Article 42 of the Dutch mining law it is currently stated that neighbouring geothermal operators must come to an agreement in dealing with potential interference [71]. With increasing doublet density, and uncertainty on subsurface flow-paths of reinjected water (e.g., [62]), this will become increasingly challenging. Especially because no accurate methods for monitoring cold-water plume development or the area of influence of doublets are considered.
Interference is also a challenge in the exploitation of other types of geothermal resources. Differently from HSA exploitation, negative interference effects can more readily be solved in Aquifer Thermal Energy Storage (ATES) development [26,28] or high enthalpy geothermal electricity production [14,27], by varying operational parameters or drilling extra (injection) wells to improve sustainable exploitation. The low financial competitiveness of HSA heat production makes it more difficult to justify such measures. Upfront investment for the creation of adequate financial-and legislative incentives as well as coordination of doublet deployment plans are therefore required to enhance HSA exploitation. Our study emphasises that although geological aspects influence considerably operational boundary conditions, the main challenges for efficient HSA exploitation are related to 'man-made' hurdles such as economic competitiveness and conflicting interests of operators.

Conclusion
Results of this study show that exploitation of HSA with a 'firstcome, first-served' approach by individual operators that only develop a single doublet system leads to sub-optimal use of geothermal resources. Geothermal exploitation could cover a significantly larger part of our heat demand, if well placement strategies follow a regionally coordinated 'masterplan' approach. Such an approach could focus on optimising the recovery efficiency. Instead, the commonly applied 'firstcome, first-served' approach aims to cover heat demand of individual operators. So far, geothermal exploitation has been supported by various governments with pilot projects, feed-in tariff schemes and geological insurance schemes. New, tailor-made financial and legislative mechanisms are required to promote more efficient use of geothermal resources. The West Netherlands Basin is used as case study, but in general these results apply to Geothermal resource exploitation worldwide. We conclude that in the WNB: • In 2018, geothermal heat exploitation provided approximately 4% of the total heat demand of the province of Zuid-Holland.
• The installed capacity grew with 24 MW th per year, in the past 4 years. With this expansion rate, geothermal heat exploitation could cover up to 20% of the heat demand in Zuid-Holland in 2050.
• Less than 0.1% of the estimated recoverable heat from Lower Cretaceous aquifers has been recovered so far. With the current expansion rates, this recovery efficiency might not have exceeded 1% by 2050.
The numerical production simulations show that for geothermal resource exploitation in general: • Heat recovery efficiency of up to 30% can be reached, but this requires denser doublet deployment with smaller injector-producer spacing than the current WNB standard of some 1.6 km and coordinated doublet deployment.
• Heat recovery optimisation does not reduce LCOH under current financial support schemes.
• Higher doublet density than the current density in the WNB could still result in sufficient doublet life time but will result in higher interference. Tailor-made legislation is therefore required to facilitate exploitation of a single resource by multiple operators.