The impact of reduction of doublet well spacing on the Net Present Value and the life time of fluvial Hot Sedimentary Aquifer doublets

This paper evaluates the impact of reduction of doublet well spacing, below the current West Netherlands Basin standard of 1000–1500 m, on the Net Present Value (NPV) and the life time of fluvial Hot Sedimentary Aquifer (HSA) doublets. First, a sensitivity analysis is used to show the possible advantage of such reduction on the NPV. The parameter value ranges are derived from West Netherlands Basin HSA doublet examples. The results indicate that a reduction of well spacing from 1400 to 1000 m could already improve NPV by up to 15%. This effect would be larger in more marginally economic HSA doublets compared to the eywords:


Introduction
Large potential resources of heat are stored in sedimentary rocks. In the Netherlands alone, the Dutch geological survey estimated the total recoverable heat from this type of resource to be approximately 55 times larger than the annual heat consumption CBS). Hot Sedimentary Aquifers (HSA) are especially suitable for 'direct use' or heat production, because they are often found in areas with average thermal gradient (Boxem et al., 2011;Pluymaekers et al., 2012). In these areas temperatures for commercial electricity production (e.g., Shengjun et al., 2011) are found at depths where pore space is generally diminished. However, heating accounts for half of the total energy consumption in, for example, the European Union (European Commission, 2016). Therefore, HSA should be considered as important energy resources. Unfortunately, a large gap currently exists between HSA potential and exploitation. This gap is a result of a combination of high initial investment costs and large uncertainties in both doublet life time and capacity. The height of the initial investment is mainly influenced by a combination of high drilling costs and the replacement of existing fossil fuel based heating networks. These two factors decrease the competitiveness of HSA exploitation with other energy sources and thereby limit its growth. In the Netherlands, for example, two new geothermal projects are realised in HSA each year, while more than 100 exploration licences have been granted since 2007. A reduction of the initial investment costs would reduce the risks for developers and hence stimulate the growth of HSA exploitation. In this paper, we evaluate the effect of well spacing reduction on the Net Present Value (NPV) of a HSA doublet. Also, we discuss the effect such a reduction would have on the life time of a doublet. The current well spacing standard in the West Netherlands Basin (WNB) and Paris Basin (Mijnlieff and Van Wees, 2009;Lopez et al., 2010;Mottaghy et al., 2011;Daniilidis et al., 2016), is 1000-1500 m. Partially, the large distance is used to prevent early cold water breakthrough. Overdesign, however, could lead to unnecessarily long thermal breakthrough time. For example, no thermal breakthrough has yet been reported in the past 40 years of exploitation in the Paris Basin. A well spacing reduction could still result in sufficient life time, while improving the financial situation of a doublet in two ways. First, it may reduce the drilling costs by making the overall well length shorter, at least for the current standard doublet layout with two deviated wells from the same surface location. Second, it could reduce the required pump energy due to shorter flow paths between the wells. Another advantage is the decrease of chance on geological flow baffles between the wells, such as sealed sub-seismic faults (e.g.,  Bailey et al., 2002) or poor sandstone body connectivity (Larue and Hovadik, 2006;Pranter and Sommer, 2011). Finally, more doublets could be realised in the same aquifer which increases the amount of produced geothermal heat (Mijnlieff and Van Wees, 2009). The number of studies on optimisation of geothermal systems is limited and often focus on maximizing energy production (e.g., Sauty et al., 1980;Ekneligoda and Min, 2014;Adams, 2015). Both the possible financial advantage of well spacing reduction and its impact on doublet life time is evaluated in this paper. In the first part, the effect of a variation in well spacing on the Net Present Value (NPV) for a typical WNB doublet is determined. This is derived from a sensitivity analysis in which both production and financial parameters are varied. Parameter ranges are derived from a West Netherlands Basin case study. The analysis is based on finite-element production simulations in ho permeable and impermeable facies bodies in the aquifer could significantly affect doublet life time and capacity (e.g., Pranter et al., 2007;Hamm and Lopez, 2012;Poulsen et al., 2015;Crooijmans et al., 2016). Homogeneous aquifer models do not capture this uncertainty. In our study, detailed fluvial facies architecture realisations are generated utilizing a process-based facies modelling approach (Cojan et al., 2004;Grappe et al., 2012;Hamm and Lopez, 2012) based on a WNB geological dataset. Minimum well spacing is analysed in terms of life time and NPV. The results of this study could be used as an incentive for re-evaluation of HSA well spacing standards. Utilizing detailed facies architecture models, more profound estimates of doublet life time and capacity can be made. This should prevent overdesign and thereby improve the competitiveness of HSA exploitation.

Data and aquifer modelling
The aquifer models in this paper were based on a geological dataset of the fluvial Lower Cretaceous Nieuwerkerk Formation in the West Netherlands Basin (DeVault and Jeremiah, 2002;Jeremiah et al., 2010;Donselaar et al., 2015). This dataset was chosen, because most of the approximately 40 exploration licences target this HSA interval in the WNB. Parameter values, which were used in aquifer modelling and production simulations, were derived from this dataset. Furthermore, production rates and reinjection temperatures were derived from WNB doublet examples. By using this dataset a realistic range of heterogeneities was derived to constraint the set of facies realisations. Two types of aquifer models were created: (1) detailed facies architecture realisations and (2) homogeneous models. The first type of realisations were generated with a process-based approach. The following sections describe the geological data and the modelling approach.

Geological dataset
A subsurface dataset of the fluvial Nieuwerkerk Formation in the WNB formed the basis for the geological modelling in this study. The same dataset and modelling approach as described in Crooijmans et al. (2016) and Willems et al. (2017) was used. The dataset comprised of cores and Gamma-ray (GR) logs of geothermal wells In the WNB. The core study provided thickness ranges of facies bodies which were used as input parameters to generate process-based facies realisations. In approximately 75 m of core in MKP-11 and 25 m in Q13-09, five different types of facies bodies were recognized: floodplain fines, crevasse splays, single-storey channel bodies and amalgamated sandstone complexes. The maximum fining upward sequence that was recognized in the cores was 4 m (Fig. 1). Therefore it was assumed that the paleo bank-full flow depth of the fluvial system that formed the Nieuwerkerk Formation was 4 m. Based on flow depth, the paleo bank-full flow width was estimated at 40 m (Williams, 1986) (Fig. 1). Furthermore, cores provided porosity-permeability relations for the aquifer property modelling. The gamma ray logs were used to derive N/G ranges of the Nieuwerkerk Formation. GR logs in the WNB showed that the Nieuwerkerk Formation significantly varies in thickness and N/G. The thickness varies from 50 to almost 200 m and the N/G the approximately 15-70% in different sections of the aquifer (Fig. 1).

Process-based detailed facies architecture realisations
To generate the detailed facies architecture realisations, a similar approach as in Crooijmans et al. (2016) and Willems et al. (2017) was used. Input parameters for the process-based facies modelling ( Fig. 2) were (1) channel width and depth, (2) maximum overbank flood deposit thickness (H th ), (3) avulsion frequency, (4) flood frequency, and (5) floodplain topography parameter (henceforth: FT-parameter) (Fig. 2). In fluvial systems, the thickness of floodplain deposit decreases away from the channel. The distance at which the thickness decreased exponentially is the FT-parameter (Fig. 2). A high FT-parameter means that the flood deposit is wide and thick, which increases the sediment aggradation rate and decreases the N/G of the realisation. The paleo bank-full flow depth was derived from the core analysis and analogues, respectively. As it was not straightforward derive values of the other parameter from cores such as the flood plain deposit thickness (e.g., Bridge, 2006), values ranges were assumed to capture the uncertainty of the parameter values and to obtain realisations that range in N/G. Flood frequency, maximum flood deposit thickness (H th ) and the FT factor were the primary controls on N/G. To obtain realisations with a wide range of N/G values between 15 and 70%, overbank flood frequency was varied between 20 and 120 years, H th between 0.2 and 0.6 m and the FT-factor between 300 and 600 m. Avulsion frequency is varied from 600 to 1600 years (Törnqvist and Bridge, 2002). During every simulated fluvial flood, sediments were deposited on the floodplain with a maximum thickness H th near the channel (Fig. 2). In the simulations, sedimentary processes distribute and shape different facies bodies such as channel lags, point-bars, crevasse splays, mud plugs and floodplain fines. Realisations have dimensions of 1 km × 2 km × 50 m and the paleo flow direction is parallel to the long edge of the realisations. The process-based Flumy software method was explained in more detail in Cojan et al. (2004), Grappe et al. (2012) and Lopez et al. (2009).

Heterogeneous aquifer models
Facies grid blocks in the realisations were divided into two classes, aquifer and non-aquifer. The non-aquifer class included fine grained facies such as crevasse splays, overbank alluvium and mud plugs. These bodies were all assumed to be relatively impermeable. Their assumed permeability was 5 mD and porosity 10%. Sandy facies bodies such as point-bars and channel lags were all assumed to be aquifer grid blocks. Porosity values were assigned to these blocks based on the core plug porosity data. From this data, a beta distribution correlation function was derived. The distribution characteristics including: mean, standard deviation, skew and kurtosis were equal to 0.28, 0.075, 0.35 and 2.3, respectively. Secondly, the permeability of each grid block was determined by a porosity-permeability relation obtained from petrophysical data of well MKP-11 (TNO, 1977): k = 0.0633e 29.5 . In this equation, k is the permeability [mD] and is the porosity [-].

Homogeneous aquifer models
For the NPV sensitivity analysis, homogeneous aquifer models were created based on the same geological dataset. These models had an average porosity of 28% and permeability between 250 and 2000 mD. This permeability range was derived from WNB HSA well tests (Lingen, 2014) and core plug measurements (TNO, 1977). The thickness of the homogeneous models in the sensitivity analysis varied between 50 and 120 m.

Numerical production simulations
In the production simulations, the aquifer was confined between two 50 m thick impermeable over-and underburden layers that provide thermal recharge. The energy balance was solved for a rigid medium fully saturated with water, in which thermal equilibrium was assumed between the fluid and solid phases: In this balance, t (s) is time, T (K) is the temperature, is the mass density (kg/m 3 ), C w (J/kg K) is the specific heat capacity, (W/mK) is the thermal conductivity, I the identity matrix, and q (m/s) is the Darcy velocity vector. The suffix w refers to the pore fluid and s to the solid matrix. The heat capacity, density and conductivity values for aquifer grid blocks are 730 J/kg K, 2650 kg/m 3 and 2.7 W/mK respectively. For the non-aquifer blocks these values are 950 J/kg K, 2600 kg/m 3 and 2.0 W/mK. The thermal conductivity and the volumetric heat capacity are described in terms of a local volume average. Heat conductivity, density and heat capacity are assumed to be independent of temperature for simplicity and described by = (1 − ) s + w and C = (1 − ) s C s + w C w in which is the porosity. This Darcy flow velocity vector can be determined by: q = (k P)/ , where is k (m 2 ) the intrinsic permeability, the temperature and salinity dependent viscosity explained in Crooijmans et al. (2016), and P (Pa) the pressure. The pressure field is obtained through solving the continuity equation: (∂ w )/∂t + ( w q) = w S, where S (s -1 ) is external sinks and sources. The detailed modelling procedure follows the approach explained in Saeid et al. (2014Saeid et al. ( , 2015. The production simulations yield a production temperature development over time and the required injection and production pressure for the defined production rate. The difference between these pressures for each time step i ( P i ) was used to estimate pump energy losses (E pump,i ) as indicated in Eq. (2) (e.g., Willems et al., 2016), where Q is the constant production rate and ε the pump efficiency. The produced energy (E prod,i ) was estimated by Eq.
(3) on each time step i (e.g., Willems et al., 2016) in which w is the water density of 1050 kg/m 3 and T i the difference between injection and production temperature. The net energy production, or the doublet capacity was determined by the sum of the produced energy and the pump energy losses.

Net Present Value model
A NPV model developed by Van Wees et al. (2010) was utilized to relate production simulations to NPV. Input for the NPV calculations were net energy production in Watt and the economic parameters listed in Table 1. In our study, additional separator costs were included, because in many WNB doublets natural gas co-production occurs. The NPV is the depreciated, discounted, netcumulative income after 15 years. This period was chosen, because it is the maximum duration of the Dutch feed-in tariff scheme (SDE+) for geothermal energy (Van Heekeren andBakema, 2013, 2015). 0.25 MD pump work-over costs were taken into account every five years, which is equal to half of the estimated pump costs. 40% down-time was assumed for maintenance throughout the year.
Because doublet wells are drilled from one surface location, doublet well spacing influences the well length (Fig. 3). In this study, the length of a single well (WL) was approximated by the sum of a vertical (D vert ) and a deviated section (D dev ), WL = D vert + D dev . The vertical section was assumed to be 1500 m. The length of the deviated section, was approximated by: In Eq. (4), True Vertical Depth (TVD) is the total well depth which is 2.2 km and L is the doublet well spacing. Reduction of the well spacing from 1000 m to 800 m could therefore result in a reduction of the total well length by 2%. Because of the assumed well cost per kilometer in Table 1, drilling costs reduce accordingly. For example, typically drilling costs of a 2.2 km deep geothermal doublet is ∼7.2 MD (Van Wees et al., 2010;Nielson and Garg, 2016). Utilizing equation 4, the reduction in spacing from 1000 to 800 m could lead to a reduction of investment of up to 0.16 MD . A reduction from the 1500 m standard to 800 m would decrease the total well length by approximately 9% and the associated well costs by approximately 0.66 MD .

NPV sensitivity analysis
In the sensitivity analysis, the effect of a 100 m and 400 m well spacing variation from a 1000 m base case scenario on the NPV was evaluated. This was compared to the relative impact of other parameters. In this analysis, homogeneous aquifer models were utilized. Two categories of parameters were compared, financial and production related parameters. The financial parameters were varied by ±10% of the base case value. Production related parameter ranges were derived from the geological data and currently active WNB doublets. The parameter ranges are listed in Fig. 6. In the analysis, each parameter was varied individually, while keeping the others at the base case value.

Impact of fluvial facies architecture on life time
In the second part of this paper, the effect of aquifer architecture on the advantage of reduction in well spacing was evaluated. In this part, production simulations were carried out utilizing detailed fluvial facies architecture realisations. In the set of realisations, N/G varied between 15 and 70% N/G. Doublet wells were placed at 400, 600, 800 and 1000 m (Fig. 4A) under the condition that both wells intersect the same amount of sandstone grid blocks, like in Crooijmans et al. (2016). Doublet well pairs were placed parallel to the paleo flow direction as this reduces the required pump energy losses (Willems et al., 2017). For comparison, two homogeneous models were generated in which the thickness was adjusted by 15 and 70% respectively ( Fig. 4B and C). In this way, both model types had the same net-aquifer volume range. Two life time scenarios were considered. In the first scenario A, the life time was reached after the production temperature decreased by 10% of its initial  value. In the second scenario B, the life time was reached when the temperature decreased by 1 • C.

Base case scenario
The thermal capacity of our base case doublet is 4.1 MW (Eqs. (2) and (3)). The associated discounted cumulative cash flow and the net annual income are presented in Fig. 5A. For the base case scenario the NPV is approximately 3.7 MD and represents the cumulative discounted cash flow after 15 years. The associated Internal Rate of Return on Investments (IRR) would be approximately 8 years. The thermal breakthrough moment happens after more than 60 years, even if aquifer thickness is reduced to 50 m (Fig. 5B). After this thermal breakthrough the production temperature only decreases by approximately 1-2 • C per decade.

NPV sensitivity analysis
The result of the sensitivity analysis is presented in Fig. 6. A tornado plot shows the change in NPV in percentages ( NPV) which is defined as: (NPV BC − NPV x )/NPV BC × 100%, as a result of the variation of a single parameter. NPV BC is the base case NPV, and NPV x is the NPV as a result of adjustment of a parameter in the sensitivity analysis. The parameter ranges are presented on the right side of the tornado plot. This analysis indicates that variation of the temperature difference between injection and production water ( T) has the most significant impact on a HSA doublet NPV. A five degree Celsius variation could change the NPV by approximately 50%. A 10 m 3 /h variation of the production rate could vary the NPV by up to 40%. Reducing the well spacing from 1400 to 1000 m could improve the NPV up to 15%. Because well length is estimated by equation 4, the well length decreases by 5.5% if the well is reduced from 1400 to 1000 m. An additional reduction from 1000 to 600 changes well length by 4.2%. Therefore the tornado plot shows a slight asymmetrical result. An additional well spacing reduction from 1000 to 600 m increases the NPV by approximately 13%. The 100 m well spacing variation has smaller impact on the NPV of approximately 4%. Permeability, aquifer thickness and pump efficiency have a limited effect on the NPV. Because of the inverse relation between injection pressure and permeability, the negative effect of a 750 mD permeability reduction on NPV is larger than the positive effect of a 1000 mD permeability increase. Variation of the pump energy efficiency by ±10% did not have a significant impact on NPV (less than 1%). Other parameters that have a significant impact on NPV are CAPEX and the height of the feed-in tariff subsidy. 10% variation of discount rate, heat price and OPEX affect the NPV by approximately 10%. In our analysis, well spacing is the only independent parameter, which is determined prior to drilling. Production rate and T depend on geological uncertainties, surface facility efficiency and on required heat consumption. Also, all financial parameters depend on economy or government policy. Therefore, these results also indicate the uncertainty in NPV as a result of uncertainty in geological and financial circumstances. Note that these results are based on homogeneous aquifer models, the effect of facies and property heterogeneities are not included.
The parameters in the sensitivity analysis affect the NPV of a doublet because they influence the discounted cumulative cash flow. This is because the initial investment costs are changed or the net income is affected. Examples of discounted cumulative cash flow in three scenarios (A, B and C) from the sensitivity analysis are presented in Fig. 7. Scenario A is the base case scenario. In scenario B the well spacing is reduced by 400 m from 1000 to 600 m while all other parameters are kept at their base case value (Fig. 6). As a result, the initial investments costs are reduced while the net income is constant. Therefore, the discounted cumulative cash flow curve is shifted upward, increasing the NPV after 15 years (Fig. 7). In scenario C the production rate is increased from 100 to 110 m 3 /h while all other parameters are kept at their base case values. The higher production rate increased the slope of the discounted cumulative cash flow, which also increases NPV (Fig. 7).
The theoretical advantage of a well spacing reduction on NPV in Fig. 6, applies to a typical WNB HSA doublet. To evaluate the effect of well spacing variation on NPV in doublets in other basins or countries, the sensitivity analysis is repeated for doublets with different drilling costs, production or reinjection temperatures and production rates. These three parameters are chosen because they have the most significant impact on NPV according to our analysis results in Fig. 6. In scenario 1 and 2, drilling costs are varied by 0.3 MD /km from the WNB base case (Table 1). In scenario 3 the T is 5 degrees higher than in our base case scenario. This could be a result of a higher production temperature or higher heat extraction efficiency and affects the net energy production (Eq. (3)). In contrast in scenario 4 the T is 5 degree lower. Finally the production rate is changed by 20 m 3 /h in scenario 5 and 6. In these six scenarios only the base case value of the indicated parameter in scenario 1-6 is varied while all other values are kept at the same value as in Fig. 6. For each scenario, a new NPV value is calculated. Subsequently, the effect of a 100 m well spacing variation from 1000 m on this NPV is calculated. The results are presented in Fig. 8. The results indicate that well spacing variation has greater impact on NPV in areas with higher drilling costs. If the drilling costs are 0.3 MD /km higher than in the WNB, the 10% well spacing reduction changes NPV by approximately 10% which is approximately 5% more compared to the result in Fig. 6. If the doublet capacity is lower because of the 5 degree reduction of T, a 10% well spacing reduction could increase the NPV by up to 15%. Well spacing reduction has an even higher impact on NPV in doublet with lower production rates. In a doublet with an 80 m 3 /h production rate, a 10% well spacing reduction increased the NPV by approximately 25%. In doublets with lower drilling costs, higher T, or higher production rates compared the WNB base case, the impact of well spacing reduction on NPV remains a few percentages. These results imply that reducing the well spacing to improve NPV is more relevant for marginally economic doublets.

Impact of fluvial facies architecture on doublet life time
The results of the production simulations with detailed fluvial facies architecture realisations are compared to simulations in homogeneous aquifer models in Fig. 9. Two life time scenarios and two production rates are compared. Three observations can be made. First, 800 m spacing is sufficient to obtain a 15 years life time in both life time scenarios and a production rate of 100 m 3 /h. If the production rate is 150 m 3 /h, 800 m is only sufficient when T is 7.5 • C. With 600 m spacing the life time exceeds 15 years only with a 100 m 3 /h production rate and with T of 7.5 • C. Secondly, these results also indicate an underestimation of the uncertainty in life time by the homogeneous aquifer models. The range of life time results per well spacing distance is larger, when detailed facies architecture realisations are used. This uncertainty increases for larger well spacing distances. In contrast, it is lower for higher production rates. The uncertainty in life time increases slightly for higher allowed production temperature drop ( T). Thirdly, the results imply that life time estimations are lower in the homogeneous models. The peaks of the distributions of life time of the detailed facies realisations are closer to the life time of the 70% homogeneous aquifer model, especially for larger well spacings.  Fig. 6. T is the difference between injection and production temperatures. For comparison, the effect of 100 m well spacing variation in the base case scenario is presented on top.

Impact of fluvial facies architecture on NPV
The NPV is higher when homogeneous aquifer models are used (Fig. 10A). This could be explained by Fig. 10B. This figure presents the ratio of produced heat (E prod ) and pump energy losses (E pump ) in each doublet scenario. It shows that in the detailed facies realisations, pump energy losses are significantly higher which decreases the net income. Fig. 10B shows that the fluvial facies architecture could significantly influence the net energy production and thereby also the NPV. This is a result of the risk of flow path formation between the wells. Therefore no clear trend is observed between the energy ratio and well spacing in the facies models. In contrast, a slight reduction of the energy ratio is recognized with increasing well spacing for the 70% N/G homogeneous model. The energy ratio is approximately 5% lower for 400 m spacing compared to 1000 m spacing. For the 15% homogeneous model this reduction is not rec- Nevertheless, the effect of well spacing on pump energy losses will not have a significant impact on NPV in this well space range. Note that the ratios in Fig. 10B apply to a single moment in time as pump energy values are time dependent (Fig. 10C and D). This is because of the transient pressure stabilisation in the aquifer and because of the temperature dependent viscosity of the brine in our simulations. Pump energy varies more significantly over time in lower N/G realisations (Fig. 10C) compared to high N/G realisations (Fig. 10D). Heat production remained constant until thermal breakthrough (Fig. 10C). The variation in pump energy losses and heat production were, however, small and varied by less than 0.02 MW in 20 years.
A consequence of reducing the well spacing is a reduction of the maximal possible amount of recoverable heat. Fig. 11 shows NPV and the maximum recoverable energy in homogeneous aquifer models with a N/G of 15, 45 and 70% assuming life time scenario A ( T is 7.5 • C), for two production rates. This comparison shows that an increasing production rate increases NPV but decreases the recoverable energy. In contrast larger well spacing increases the recoverable energy as it increases life time. Larger well spacing is therefore favourable in lower N/G aquifers or in doublet with high production rates. These results suggest that a compromise between financial efficiency and recovery efficiency of a doublet could be obtained depending on the expected sandstone volume, required production rate and required life time. Finally, Fig. 11 confirms that reducing well spacing reduction to improve NPV is more relevant for marginally economic doublets with low production rates. In the 45% and 70% N/G models with 100 m 3 /h production rate, an optimum in NPV is most clearly recognized. No peak in the NPV curves is observed when the life time is too low because of high production rate or low sandstone volume.

Well spacing reduction and the advantage on NPV
Our results show a possible improvement of the NPV of approximately 15%, when well spacing was reduced from 1400 to 1000 m. Even a small reduction of 100 m could already improve NPV by up to 4%. The sensitivity analysis indicates that variation of T (T prod − T inj ) and production rate have a more significant impact on NPV. However, these parameters are to a large extent dependent on (1) geological uncertainties and (2) technological constraints. Geological uncertainties include aquifer temperature of several degrees  and aquifer productivity. Also all financial parameters are influenced by external factors such as economy, government policy and regional rates of equipment and technology. Well spacing is one of the parameters in our analysis that could be designed prior to drilling. So far, well spacing is primarily designed to create sufficient life time (e.g., Motthagy et al., 2011;Daniilidis et al., 2016). Therefore our study takes a more comprehensive look at the effect of well spacing. This becomes increasingly important with more widespread deployment of geothermal doublets. The current spacing standard significantly limits the possible number of doublets in the WNB (e.g., Mijnlieff and Van Wees, 2009). Well spacing reduction would therefore not only improve the financial competitiveness of individual doublets but also enhance exploitation on a regional scale.
Furthermore, our results indicate that well spacing reduction has a stronger impact in doublets with a lower productivity or net energy production rate. This is because the NPV results from a combination of initial investment and the height of the net income. These parameters determine the starting point and slope of the discounted cumulative cash flow (Fig. 5). Optimisation of NPV by reducing well spacing is therefore of particular interest for marginally economic HSA doublets.
Our results are based on the assumption that the minimal required life time is 15 years. This applies to HSA doublets in the Netherlands, because this is the maximum duration of the Dutch feed-in tariff scheme. As a result the well spacing for an optimal NPV is the smallest well distance with a minimal production temperature reduction in this period of time. This is of course related to a specific production rate. Most geothermal systems however, are designed to produce for a much longer time such as 30-50 years, which requires a larger well spacing.
The NPV model by Van Wees et al. (2010) assumes an ideal HSA exploitation scenario. It does not accurately reflect current WNB investment costs and NPV. Firstly, this is a result of neglecting the costs which are associated to years of preparation prior to drilling and consultants. In addition, often unexpected maintenance and workover costs add to the initial investments in the first years. HSA exploitation started approximately 10 years ago in the WNB and currently passes through a learning curve. Therefore, our estimated investment costs apply more to future Dutch doublets that are able to take advantage of the exploitation experience of the first decade of HSA exploitation. Furthermore, the NPV is increased in the base case scenario of the sensitivity analysis (Fig. 5A) because of the relatively high porosity and permeability values of 28% and 1000 mD, respectively. In WNB doublets the NPV might be lower because of the impact of geological heterogeneity on net energy production (e.g., Fig. 10A). Nevertheless, well spacing reduction would be especially relevant for more marginally economic doublets as was shown in Fig. 8.
In our simulations, vertical wells are used. Deviated wells could improve the injectivity and productivity increasing the net energy production and improve NPV. Hence, for optimisation of doublet design, the preferred angle of the deviated well should be determined for a specific aquifer as it depends on the geological model, net sandstone volume and aquifer thickness.

Life time
The main reason for the large current well spacing standard of 1000-1500 m is to avoid early cold water breakthrough. However, our production simulations indicate that this spacing might be overcautious. The life time in our production simulations varies between 60 and 100 years in life time scenario A, and 15-35 years in life time scenario B both with 150 m 3 /h production rate and a 1000 m well spacing (Fig. 9). This shows that despite our conservative aquifer thickness of 50 m, still a considerable life time could be expected. Most current WNB doublets actually have a larger well spacing. However, well spacing optimisation requires a better understanding of the uncertainties in life time. Doublet life time is significantly affected by facies architecture (Hamm and Lopez, 2012;Crooijmans et al., 2016). Our production simulations with detailed facies architecture realisations indicate that homogeneous models underestimate doublet life time. This is because impermeable claystone bodies provide thermal recharge to the cold water plume and flow baffles increase flow path length between the wells. Secondly, homogeneous models overestimate the net energy production of the doublet (Fig. 10B) because flow baffles or barriers could decrease injectivity or productivity. Especially in low N/G aquifers, wells might intersect different sandstone bodies that do not form flow paths between the wells which increases the required pump energy losses. Doublet well pairs were oriented parallel to the paleo flow direction and the orientation of the fluvial sandstone bodies. Perpendicular doublet orientation would increase the required pump energy losses (e.g., Willems et al., 2016) and therefore reduce NPV. In addition, perpendicular doublet orientation would have increased life time uncertainty. This is a result of a higher probability that impermeable claystone bodies form flow baffles or barriers. Flow-barriers can reduce life time when they cut-off permeable parts of the aquifer reducing the net reservoir volume while the flow baffles increase life time when they increase flow path length (e.g., Crooijmans et al., 2016).
Our results show that a statistical approach with multiple realisations is required to capture uncertainties in life time and doublet capacity. The large spread in life time ( Fig. 9) indicates that more realisations would be required to accurately capture this uncertainty in life time and capacity as were used in our study. The goal of our study was, however, to examine how facies heterogeneities influence the theoretical advantage of well spacing reduction on NPV. Our results support the possibility to reduce the well spacing, especially because the net sandstone volume in our models is very conservative. For example, no WNB doublet only encounters a 15-20% N/G aquifer of 50 m thickness. In general, WNB aquifers have a larger thickness and higher N/G percentage. Our conservative sandstone volumes results in low doublet capacities (Fig. 10B) and therefore lower NPV compared to the base case scenario. In reality, doublets with very low injectivity or net sandstone volume would not be taken in production without any measures (e.g., Blöcher et al., 2015). Examples of such a measures are to continue drilling into a higher N/G interval or to increase heat exchange surface by hydraulic stimulation. If we would compensate for this, life time and NPV in our results would increase. The detailed facies architecture realisations in this current study are simplified by neglecting small scale permeability heterogeneities and assuming isotropic aquifer properties in grid blocks. Examples of such small scale heterogeneities are shale drapes, accretion surfaces and bedding planes (e.g., Pranter et al., 2007). These features decrease on average the permeability perpendicular to the paleo flow direction. This could be accounted for adjusting the permeability in different directions in each grid block like in Bierkens and Weerts (1994). Furthermore porosity was randomly distributed amongst the sandstone grid blocks. In reality the porosity of channel lags, point-bars and sand plugs varies more systematically across sandstone bodies as a result of sedimentary process (Willis and Tang, 2010), which could influence flow path formation through the sandstone bodies (Larue and Hovadik, 2008). These heterogeneities could further diffuse the cold water plume. Note that in our simulations flow only occurs through the rock matrix and no fractures are present (e.g. Hardebol et al., 2015;Bisdom et al., 2016). Fractures, especially parallel to the flow between the doublets could decrease the thermal breakthrough time of the cold water plume.
Finally, it could be argued that not only the sandstone volumes in our simulations were conservative, also the assumptions on the minimal production temperature could be considered as such. The minimum required production temperature for economic Dutch HSA exploitation for greenhouse heating is assumed to be 45 • C and 65 • C for district heating purposes . These minimal temperatures are lower than the 67.5 • C in our life time scenario A and much lower than the 74 • C of scenario B. Furthermore, due to likely technological improvements on insulation and technological efficiency in the next decades, lower production temperatures will most likely be sufficient in the future. Our simulations show that production temperature reduces by 0.5-1 • C per year, depending on the net sandstone volume and production rate. The speed of production temperature reduction in our simulations is relatively low compared to previous studies (e.g., Mijnlieff and Van Wees, 2009) because thermal recharge from over and under burden is taken into account (Poulsen et al., 2015). In Fig. 12 the production temperature development over time in homogeneous aquifer models and detailed aquifer architecture models are compared as an example. Doublets in these simulations have a 600 m well spacing and 100 m 3 /h production rate. In the lowest N/G aquifer realisation, the production temperature drops below 60 • C after 60 years. In the highest N/G realisation this production temperature is reached after approximately 140 years. Therefore, Fig. 12 underlines conservative nature of our minimal temperature assumptions.

Conclusion
This paper presents a modelling based study on the advantages and impacts of a well spacing reduction in HSA doublets. Firstly, our results indicate that a reduction of 400 m well spacing could lead to a NPV improvement of up to 15%. The impact of well spacing reduction on NPV is more significant in marginally economic doublets. Secondly, production simulations with detailed facies aquifer architecture derived from WNB geological data are used to evaluate the impact of smaller well spacing on life time. Our results suggest that sufficient life time could be obtained, if the well spacing is reduced below the current 1000-1500 m standard. These current standards aim to avoid early cold water breakthrough. However, our results show potential overdesign of these standards. This results in the negative effect of decreasing the financial competitiveness of geothermal exploitation and reducing the possible number of doublets in a region. Finally, the comparison of production simulations with detailed facies architecture realisations and homogeneous models shows that homogeneous models underestimate the uncertainty in life time and doublet capacity. The aquifer heterogeneities have a significant effect on the flow path formation between the wells. Therefore, a statistical approach with multiple realisation is required to accurately assess HSA potential and to optimise exploitation efficiency.