On the connectivity anisotropy in fluvial Hot Sedimentary Aquifers and its influence on geothermal doublet performance

This study ﬁnds that the geothermal doublet layout with respect to the paleo ﬂow direction in ﬂuvial sedimentary reservoirs could signiﬁcantly affect pump energy losses. These losses can be reduced by up to 10% if a doublet well pair is oriented parallel to the paleo ﬂow trend compared to perpendicular. The chance that ﬂow paths are formed perpendicular to this trend strongly depends on the net sandstone volume in the reservoir. Detailed ﬂuvial facies architecture realisations which are used in this study, are generated with a process-based approach utilizing geological data from the Lower Cretaceous Nieuwerkerk Formation in the West Netherlands Basin. Finally, this study emphasizes the importance of detailed facies architecture modelling for the assessment of both risks and production strategies in Hot Sedimentary Aquifers.


Introduction
Hot Sedimentary Aquifers (HSA) are commonly exploited by a doublet system, consisting of a hot-water production and a coldwater reinjection well. Downhole well distance typically is 1000 to 2000 m, and both wells target the same aquifer to maintain pressure support in the reservoir. In fluvial reservoir rocks the doublet connectivity is via a network of permeable fluvial channel sandstone bodies embedded in non-permeable floodplain mudstone. Detailed knowledge of the size, shape, spatial distribution and connectivity of the fluvial sandstone bodies (or: fluvial reservoir architecture) is required to assess the risk of pressure communication loss between the wells and the inherent economic risk of the geothermal energy production projects (Fig. 1).
The effect of the fluvial reservoir architecture on the recovery of hydrocarbons has extensively been studied (e.g. Jones et al., 1995;Larue and Friedmann, 2005;Hovadik, 2006, 2008 Pranter et al., 2007;De Jager et al., 2009). To a much more limited extent, this topic is addressed for geothermal energy production (e.g. Hamm and Lopez, 2012;Crooijmans et al., 2016) and for CO 2 sequestration (e.g. Issautier et al., 2014). Larue and Hovadik (2008) identified 'connectivity' as one of the main parameters that control the recovery efficiency of hydrocarbon reservoirs. Connectivity could be defined as the ratio of the volume of the largest sandstone body cluster and the total sandstone body volume (e.g. Larue and Hovadik, 2006). If the connectivity is high, less isolated clusters occur and therefore fewer wells are required to drain the reservoir (e.g. Geel and Donselaar, 2007). Previous work on connectivity in sedimentary reservoirs identified several main factors that control the chance that sandstone bodies connect: (1) the netsandstone volume or net-to-gross (N/G); (2) the sandstone body geometry, and (3) the range in paleo-flow direction, which determines the reservoir trend (King, 1990;Hovadik, 2006, 2008;Geel and Donselaar, 2007;Ainsworth, 2005;Pranter and Sommer, 2011). Connectivity of reservoir bodies is also influenced by post-depositional faulting (e.g. Bailey et al., 2002), by diagenetic processes, and by depositional permeability heterogeneity within in the sandstone bodies (Willis and Tang, 2010;Henares et al., 2014). To date, studies into the risk assessment of connectivity are dominantly focused on the optimization of hydrocarbon http://dx.doi.org/10.1016/j.geothermics.2016.10.002 0375-6505/© 2016 Elsevier Ltd. All rights reserved.  recovery efficiency. A main goal of connectivity analyses was to identify a N/G threshold below which isolated bodies start to occur. In meandering fluvial reservoirs, this N/G threshold is often recognized between 20 and 30% N/G, depending on the sandstone body geometries (e.g. Larue and Hovadik, 2006). Because in geothermal exploitation well pairs are used, a new directional component in connectivity analyses is required. This can assist engineers to design geothermal doublet systems with the largest possible heat exchange surface between two wells in order to minimize pump energy losses and to delay cold water breakthrough. A conceptual fluvial reservoir model illustrates the difference between the hydrocarbon and geothermal exploitation objectives (Fig. 2). The model contains five wells in an L-pattern with a 500 m spacing and an alignment parallel and perpendicular to the paleo flow direction. In terms of drainable volume, these wells are efficiently placed and intersect most of the sandstone bodies in the reservoir. However, if the wells included geothermal doublets, the distance and orientation of the well pair layout would significantly influence the chance that flow paths are formed between well pairs. Please note that the well spacing in the model is a third to one half of the 1.5 km spacing commonly used in HSA doublets (Lopez et al., 2010;Mottaghy et al., 2011). A larger well spacing would increase the risk of connectivity loss. The chance that sandstone bodies form flow paths parallel and perpendicular to the paleo flow direction (i.e., the connectivity anisotropy) has so far not been investigated.
The West Netherlands Basin (WNB) is an example of an area with fluvial HSA exploitation. Six doublets currently produce from the fluvial Nieuwerkerk Formation (DeVault and Jeremiah, 2002;Van Heekeren and Bakema, 2015). In three of them the doublet layout is parallel to the paleo flow trend. In the other three doublets, the layout is oblique or perpendicular. Productivity and injectivity vary considerably in the WNB (Van Wees et al., 2012). The reduction in injectivity could be related to well layout but also to other factors such as scaling or skin formation. Van Wees et al. (2012) pointed out that unfortunately it is not possible to identify a single cause of this variability because of limited available data. The uncertainty in injectivity and productivity, limits the growth of HSA development. In the Netherlands this is reflected by the fact that approximately 100 exploration licences are granted, while only 14 doublets are actually realised in the past 10 years. Such a gap between HSA potential and actual exploitation exists worldwide (Boxem et al., 2011). Other examples of sedimentary basins with large HSA potential but limited exploitation are the Perth Basin, Australia (Pujol et al., 2015), and the Idaho thrust belt (Welhan, 2016). A better understanding of connectivity anisotropy could reduce the risks associated with HSA exploitation and hence support its growth. Therefore, the first goal of this paper is to evaluate connectivity anisotropy and its dependence on N/G. Secondly, the possible effect of this anisotropy on doublet performance is evaluated. The results should contribute to fluvial HSA development strategies that increase the efficiency and decrease the risks of exploitation.
For this purpose, hundreds of detailed facies architecture realisations have been generated. This stochastic approach, in which reservoir heterogeneities are taken into account, is standard in hydrocarbon exploitation (e.g. Keogh et al., 2007). In contrast, geothermal reservoirs are often modelled as homogeneous layers (e.g. Mottaghy et al., 2011). The realisations are based on a geological dataset of the Lower Cretaceous Nieuwerkerk Formation (DeVault and Jeremiah, 2002). Sediments in this interval were deposited by a syn-rift, meandering fluvial system. Extensional faulting in the Late Jurassic created half-graben structures between southeast to northwest trending faults. These structures guided the paleo-flow direction of the fluvial system. Intervals with different N/G trends are recognized (DeVault and Jeremiah, 2002). N/G trends are determined by the sediment aggradation rate (Shanley and McCabe, 1994;Posamentier and Morris, 2000). As a result of low aggradation rates meander loops have more time to develop. While they migrate laterally, floodplain fines are eroded and a high N/G interval is generated with wide, thick, amalgamated sandstone complexes. In contrast, low N/G intervals with more narrow and isolated sandstone bodies are created as a result of fast aggradation rates. This is caused by frequent flooding and deposition of fine sediments on the floodplain. The varying N/G trends in the Nieuwerkerk Formation, create uncertainty about connectivity of the sandstone bodies between the doublet wells. Our facies modelling approach is similar to the one used by Crooijmans et al. (2016). In our study, the facies realisations are generated with a process-based approach (Karssenberg et al., 2001;Cojan et al., 2004;Karssenberg and Bridge, 2008;Lopez et al., 2009;Grappe et al., 2012). This is different from previous connectivity analyses that used a more standard object-based facies modelling approach (e.g. Keogh et al., 2007). As the facies modelling approach influences both the spatial distribution and the shape of the bodies, it could have an effect on a connectivity analysis (Villamizar et al., 2015). In object-based modelling, the spatial distribution of the sandstone bodies in the models is random. In contrast, in a process-based modelling approach the spatial distribution of sandstone bodies is governed by the simulated sedimentary processes. This creates a more realistic and sedimentologically-based spatial distribution of facies bodies. Another characteristic of the process-based modelling approach is that the geometry of the facies bodies and N/G are related (e.g. Bridge, 2006;Cojan et al., 2004). Previous connectivity studies, however, used an object-based approach where N/G and sediment body geometry were varied independently which could affect the results (e.g., Larue and Hovadik, 2006). With our approach we are able to show that the facies architecture is non-negligible and that detailed geological modelling is required to increase efficiency of HSA exploitation.

Data and methods
Hundreds of facies realisations were generated that vary in N/G. The realisations were analysed in three steps. Firstly, the relation between sandstone body clustering and N/G was determined by a connectivity analysis. Secondly, the connectivity anisotropy was analysed by deriving the equivalent permeability in three directions, parallel, perpendicular and vertical to the paleo flow direction. Finally, well pairs were placed parallel and perpendicular to the paleo flow in the realisations and equivalent permeability and pump energy losses were calculated in steady state production simulations.

Geological dataset
The geological modelling in this study is based on a subsurface dataset of the fluvial Nieuwerkerk Formation in the WNB. The dataset was used to derive a realistic range of heterogeneities to constrain the set of facies realisations. The dataset consists of Gamma-Ray (GR) logs and cores from deep wells, their locations are indicated in Fig. 3. The core study provided thickness ranges of facies bodies which were used as input for the facies modelling. 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 thickness of individual sandstone bodies is approximately 4 m (Fig. 4). Based on the bank-full flow depth, the bank-full flow width was estimated at 40 m (Williams, 1986). Crevasse splay thickness in the cores varied between 0.2 to 0.6 m. (Fig. 4). Furthermore, cores provide porosity-permeability relations for the reservoir property modelling. The gamma ray logs are used to derive N/G ranges of the Nieuwerkerk Formation. GR logs in the WNB show that N/G varies from approximately 10 to 90% (Fig. 3).

Process-based facies modelling
Input parameters for the process-based facies modelling software, Flumy, were (1) channel width and depth, (2) overbank flood deposit thickness, (3) avulsion frequency, (4) flood frequency, (5) maximum overbank flood deposit thickness (H th ) and (6) floodplain topography parameter (henceforth: FT-factor). The thickness of floodplain deposit decreases away from the channel (Fig. 5). The distance at which the thickness decreased exponentially is the FT factor. A high FT factor means that the flood deposit is wide and thick which increases the sediment aggradation rate and decreases the N/G of the realisation. Parameters 1 and 2 were derived from the core analysis and analogues respectively. The other parameters cannot be derived directly from subsurface data. Therefore large ranges of values were used to capture the uncertainty. Avulsion frequency was varied between 200 and 1600 years. This parameter could not be derived from the dataset and hence a large range around a 800 to 1000 year (Törnqvist and Bridge, 2002) mean was used. Avulsion frequency mainly affects the sandstone body width. Flood frequency, H th and the FT factor were the primary controls on N/G. To obtain realisations with N/G values between 10 and 90%, overbank flood frequency was varied between 20 and 200 years, H th between 0.2 to 0.6 m and the FT-factor between 300-900m. During every simulated flood, sediments were deposited on the floodplain with a maximum thickness H th near the channel. In the simulations, sedimentary processes distributed and shaped different facies bodies such as channel lags, point-bars, crevasse splays, mud plugs and floodplain fines. The process-based method implemented in Flumy software is explained in more detail in Cojan et al. (2004), Grappe et al. (2012) and Lopez et al. (2009).

Facies realisations
Six types of facies bodies were distributed in the realisations. The set of realisations were divided in two groups. In realisation Group 1, the paleo flow direction was parallel to the long edge and in Group 2, the paleo flow direction was perpendicular to the long edge of the model (Fig. 6). Our facies realisations have dimensions of 1 km × 2 km × 50 m which relate to a typical area of influence of a HSA geothermal doublet (e.g. Lopez et al., 2010). Grid block size varies from 10 cm near the well bore region to 20 × 20 × 2.5 at 100 m away from the wells. By utilizing this resolution, grid blocks are smaller than the assumed geometries of sandstone bodies (Zhang and Montgomery, 1994;Loughnane et al., 2014). Synthetic GR logs were made by extraction of a facies column from realisations. GR values were assigned to different facies in these columns: Channel Lag 15, Point-Bar 20, Mud Plug 120, Crevasse splay 50, Overbank alluvium 140. These values were derived from the core analysis (Fig. 4).

Porosity and permeability modelling
The facies body types that result from the process-based modelling were divided into two classes, reservoir and non-reservoir. The non-reservoir class includes fine grained facies such as crevasse splays, overbank alluvium and mud plugs. Their assumed permeability and porosity are 5 mD and 10%, respectively. Sandy facies bodies such as point-bars and channel lags were all assumed to be reservoir 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 are 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.0633 e 29.507 , where k is the permeability [mD] and is the porosity [−]. Because of this specific porosity distribution, the arithmetic average sandstone permeability is approximately 1000 mD.

Analysis methods
The set of facies realisations were analysed in three ways. First, the clustering of sandstone bodies was evaluated in a connectivity analysis. Secondly, the equivalent permeability in different directions and between the wells was calculated in steady-state finite-element isothermal production simulations utilizing COM-SOL Multiphysics ® . For this purpose, a pressure difference was applied to opposite model boundaries parallel, perpendicular and vertical to the paleo flow direction. The resulting average Darcy flow velocity was calculated and related to equivalent permeability in all three dimensions. Finally, the formation of flow paths between doublet wells was evaluated in a similar way. More than 2000 doublets wells were placed in the facies realisations. Equivalent permeability between doublet wells was compared for parallel and perpendicular doublet layout.

Connectivity analysis
In all realisations, sandstone body clustering was evaluated as a function of N/G. The connectivity (C) was defined such that it is equal to the ratio of the largest sandstone cluster volume and the total sandstone volume (Fig. 7).
In this work, three definitions for connectivity between grid blocks were compared to obtain more information on how sandstone grid blocks are distributed in the realisations. A first option is to consider two blocks as 'connected' only if they have an adjacent face ( Fig. 8A-1). Secondly, two blocks could be connected if they share an edge (Fig. 8A-2). In this way, not only six but eighteen connections can be formed. Finally, two grid blocks can be considered as connected if they share a corner (Fig. 8A-3) which results in 26 possible connections. In summary, connectivity was calculated for three connectivity definitions defined as 6-, 18-and 26-point connectivity.

Equivalent permeability on realisation scale and in well pairs
In steady-state, finite-element production simulations, a 40 bar pressure difference ( P) was consecutively applied in three dimensions between opposite realisation boundaries. To determine the fluid pressure field, the single-phase steady-state continuity equation was solved with constant viscosity ( ) using Eq. (2).     In this balance, k facies is the permeability field which is assigned to the grid blocks in each realisation and is the water viscosity (0,001 Pas). The detailed modelling procedure follows the approach by Saeid et al., 2014Saeid et al., , 2015. The Darcy flow velocity (v) can be found as, Subsequently, by integrating the fluid flux (q) across the model boundaries on which the pressure difference is applied, the equivalent permeability K can be determined parallel, perpendicular and vertical to the paleo flow direction using Eq. (4) (e.g., Matthaï and Nick, 2009;Bisdom et al., 2016).
where A is the cross-sectional area of the flow and L the distance between the boundaries of the realisation on which the pressure difference is applied. The derived equivalent permeability was related to the N/G of the realisation in the analysis of the results. Finally, equivalent permeability was also determined between doublet well pairs. In all realisations, doublet well pairs were placed at two spacing distances: 800 and 1000 m. For each spacing distance, three different locations were used in each realisation. Well pairs were always placed parallel to the long edge of the models on the central axis. In total, 2100 simulations were carried out.
Large numbers of simulations are required to get statistically meaningful results due to the geological uncertainties associated with random well placement. The simulations yielded a required pressure difference between wells for a 100 m 3 /h production rate. This pressure difference was used to determine equivalent permeability between wells using Eq. (2). Subsequently, the associated pump energy (Watt) was estimated by Eq. (5), where Q is the flow rate and ε the assumed pump efficiency of 60%. The N/G was the main parameter in our analyses and varied between 10% to 90% in the realisations. The equivalent permeability analysis between the well pairs was calculated for different well placement.

Facies architecture analysis
In low N/G reservoirs, impermeable floodplain fines separate the sandstone bodies and form extensive flow baffles perpendicular to the paleo flow direction (Fig. 9A). This is evident from the reservoir example of Fig. 9A with a N/G value below the connectivity threshold (Larue and Hovadik, 2006). Because of this low N/G value, many isolated single storey sandstone bodies occur. At    N/G values above the N/G threshold (Fig. 9B), the sandstone bodies amalgamate, increase in width and form one big cluster with more flow paths. However, still more flow baffles perpendicular to the paleo flow direction can be recognized. Also vertical flow baffles are maintained as indicated by the synthetic GR log. These baffles will decrease vertical permeability and permeability perpendicular to the paleo flow direction. The effect of these baffles on connectivity and equivalent permeability is discussed below.

Connectivity
Results of the connectivity analysis show differences between the three definitions of connectivity (Fig. 10). For the 18-point connectivity and 26-point connectivity definitions results are similar. 6-point connectivity results in approximately 10% lower connectivity value in the same realisation. The difference between 6-point and 18-point connectivity indicates that grid blocks are often close but do not share a face. This will influence production simulations because flow can only occur through grid block faces. A N/G threshold is recognized at 30% N/G. The use of 18-and 26-point connectivity results in a slight shift of this threshold to approximately 25% N/G. Below 30% N/G, the connectivity to N/G relation has a large standard deviation (Fig. 11A). Due to this uncertainty in the connectivity, more realisations are required in the low N/G region to determine a stable average (Fig. 11B). When the N/G is below 10%, more than 25 realisations are required. In the second N/G range from 10% to 15%, the required number of realisations decreases to 15 realisations. Finally, if the N/G is more than 15%, 5 realisations are sufficient.

Equivalent permeability between opposite model boundaries
The relation of equivalent permeability to N/G has a different trend compared to the connectivity analysis (Fig. 12). Equivalent permeability increases between 5 and 100% N/G to 1000 mD which is the average sandstone permeability. This indicates a dependence of fluid flow behaviour on N/G, also above the connectivity threshold. Below 70% N/G, equivalent permeability is higher in the direction parallel to the paleo flow direction, despite the isotropic properties in each grid block (Fig. 12). More flow paths are formed parallel to the paleo flow direction compared to perpendicular. This increases the equivalent permeability. The vertical equivalent permeability behaves differently. Below 30% N/G, only very few vertical sandstone grid block connections are formed from top to base in the realisations. Above this N/G value, the vertical equivalent permeability increases but is lower than equivalent permeability in the horizontal directions. The relations of equivalent permeability in the horizontal directions to N/G are in between the harmonic and arithmetic average permeability curves of the realisations. This means that in every realisation connections are formed between the realisation boundaries.
For a more detailed analysis of the anisotropy, the ratio of perpendicular and parallel equivalent permeability (k per /k par ) and vertical and parallel (k vert /k par ) are determined. These ratios are related to N/G and compared in Fig. 13. Between 10% and 20% N/G, the equivalent permeability is approximately 40% lower in the direction perpendicular compared to parallel to the paleo flow. This anisotropy decreases towards 70% N/G. Above this threshold, equivalent permeability is equal in both horizontal directions. Vertical permeability increases in this range and is equal to the permeability in horizontal direction at 100% N/G.
To determine whether sufficient realisations are used, equivalent permeability values are grouped in 5% N/G bins. The minimal number of realisations for a stable average is five facies realisations per N/G bin (Fig. 14). The variance in the equivalent permeability to N/G relation is significantly lower than the connectivity to N/G relation (Fig. 11).

Equivalent permeability between doublet wells
Comparing equivalent permeability calculations between opposite realisation boundaries and between well pairs, three observations can be made. Firstly, equivalent permeability between doublet wells (Fig. 15A), shows a smaller anisotropy compared to equivalent permeability between realisation boundaries (Fig. 12). Secondly, the equivalent permeability as a function of N/G is lower. Thirdly, the scatter of the results is much larger. These three observations result from the random well placement. In our simulations, it is possible that both wells intersect a different number of sandstone grid blocks or no sandstone grid blocks at all. This would result in unexpected low equivalent permeability, even at high N/G  values. When the simulations are used to estimate pump energy losses, the anisotropy is more clearly recognized (Fig. 15B). Because of the inverse relation between the simulated pressure difference and permeability (Eq. (4)), the effect of doublet layout with respect to the paleo flow trend is more clearly expressed in pump energy losses that relate proportional to pressure (Eq. (5)). Nevertheless, our results show that pump losses are approximately 0.1 to 0.2 MW higher with a perpendicular doublet layout if the N/G is lower than 60%. This would be approximately 10% of the total capacity of a typical WNB HSA doublet.

Connectivity
Our connectivity analysis showed a difference between the use of 6-and 18-point connectivity, especially in lower N/G realisations. It is uncertain to which extend this is a resolution effect. One could imagine that two adjacent grid blocks that connect through an edge are in reality part of a single sandstone body. Higher resolution realisations are required to accurately capture connections in lower N/G reservoirs. In lower N/G reservoirs amalgamation is less frequent and therefore the average surface of connections between sandstone bodies decreases (Bridge, 2006). Smaller objects need higher resolution grid blocks (e.g. Loughnane et al., 2014). This is however not evaluated in this study. The choice for one of the three connectivity definitions depends on the purpose of the study. If models are generated for production simulations, 6-point connectivity relates more to the simulations because flow only occurs through grid block faces. If the realisations are generated for volumetric analyses, the comparison of the three definitions gives information on how the reservoir grid blocks are distributed. This could be used to evaluate the facies modelling.
Despite our process-based facies modelling approach, the connectivity to N/G relation is not significantly different from previous work (e.g. Larue and Hovadik, 2006;Pranter and Sommer, 2011). Therefore, our results do not confirm the expectation of Villamizar et al. (2015) that object-based modelling could result in overestimation of connectivity. The discrepancy between the expectation of Villamizar et al. (2015) and our results could be a consequence of the size of the realisations relative to the size of the individual sandstone bodies. In our 1 km wide models, we describe connectivity of sandstone bodies on a channel belt scale (Donselaar and Overeem, 2008). Villamizar et al. (2015) based their suggestion on studies of Hajek et al. (2010) and Flood and Hampson (2015). These studies recognized autogenic sandstone body clustering in outcrops that are approximately one order of magnitude larger than the channel belt width. Therefore, larger realisations should be used to test this expectation. Connectivity on this larger scale applies to risk assessment of interference between adjacent doublets. Our results apply to connectivity within a single doublet. For evaluation of connectivity on a doublet scale, both facies modelling methods are adequate.

Equivalent permeability on realisation scale
Comparison of our connectivity and equivalent permeability analyses indicates that the connectivity analysis alone is insufficient to determine doublet layout strategies. This analysis is not able to differentiate the potential of HSA reservoirs with different N/G above the connectivity threshold. However, Crooijmans et al. (2016) showed that also above this threshold, doublet life time depends on N/G. In contrast, the equivalent permeability has an increasing trend over the complete N/G range. This shows that on average an increasing number of flow paths is formed when the N/G increases. Furthermore, the equivalent permeability analysis provides directional information on connectivity. Connectivity is the main factor that influences hydrocarbon recovery (Larue and Hovadik, 2008), for heat recovery however additional analyses are required to assess the potential.
To apply the results of our study, three main factors must be taken into account. First, the N/G threshold values and equivalent permeability ratios relate to this specific set of reservoir realisations. They might vary for fluvial reservoirs with different sinuosity, range in paleo flow direction or width to thickness ratio of the sandstone bodies. These are the main parameters that affect connectivity (Larue and Hovadik, 2006). The same workflow but different geological parameters could be applied to assess connectivity anisotropy in other HSA basins. Secondly, the results are affected by simplification of the geological modelling in our study. For example, small-scale internal sandstone body heterogeneities which are smaller than the grid block resolution are neglected. Reservoir properties are assumed isotropic in each grid block. In reality, small-scale sedimentary heterogeneities such as shale drapes, accretion surfaces and bedding planes decrease the permeability perpendicular to the paleo flow direction (e.g. Pranter et al., 2007). This could be accounted for by utilizing anisotropic grid block permeability like in Bierkens and Weerts. (1994). Thirdly, sandstone porosity is randomly assigned to sandstone grid blocks. In reality, grain size heterogeneity within sandstone bodies depends on paleo flow speed, and the proximity to the channel axis and river bends. As a result, the permeability distribution is not random across sandstone bodies (Willis and Tang, 2010). These factors could influence the magnitude of the anisotropy and the N/G threshold above which the anisotropy vanishes.
Finally, these results indicate large risks associated with horizontal wells in contrast to the results of Hamm and Lopez (2012). Because of the low vertical equivalent permeability, the chance is small that flow paths are formed between two wells. To increase this chance, well length should be large which in turn will significantly increase well costs. This is most likely not an attractive strategy because current HSA exploitation with deviated wells is already marginally economic. The vertical equivalent permeability in our results would be higher if the thickness of the realisations was increased. However, it will always remain lower than equivalent permeability in horizontal directions because of frequent vertical flow baffles that are preserved, also in higher N/G aquifers (Fig. 10).

Equivalent permeability between two wells
The anisotropy in equivalent permeability between two wells cannot be clearly recognized compared to equivalent permeability between two opposite model boundaries (Fig. 12). This is a result of geologic uncertainty associated with well placement (Fig. 15A). If in our simulations a constraint had been used for the well location stating that both wells should intersect the same amount of sandstone, the anisotropy in connectivity would have become more clear. However, we chose unconstraint well placement to evaluate the order of magnitude of pump energy losses as a result of realistic well placement. These losses in our calculations are conservative. In reality doublets with high pump energy losses would not be taken into production without any measures to improve injectivity or productivity. Examples of such measures could be (1) continued drilling into a deeper and higher N/G interval, (2) creation of a sidetrack or (3) hydraulic stimulation of the well. Moreover, reservoirs with a 10 to 20% N/G are not likely to be exploited at all. No WNB doublets installed so far encounter only reservoir intervals with N/G lower than 30%, especially not with a small total thickness of 50 m like in our realisations.
The present study focused on the risks associated with perpendicular well layout. However, an advantage of a perpendicular layout could be that longer flow paths are formed which may increase the doublet life time (Hamm and Lopez, 2012). This would allow closer well spacing, reducing well path length and hence drilling costs. Next to reservoir architecture the structural setting is also a doublet layout constraint. For example, fault blocks in the WNB dip perpendicular to the paleo flow direction (DeVault and Jeremiah, 2002). Therefore, a consideration in doublet orientation could be to target the deeper and hotter part of the fault block by a production well, and the more shallow part with an injection well to take advantage of the hydrostatic head within the reservoir. The balance between advantages and disadvantages of these constraints should be analysed in further studies with transient production simulations that provide a basis for net energy optimization. Finally, our results underline the importance of detailed geological modelling. Homogeneous models underestimate the risks related to connectivity. A stochastic approach with detailed modelling of reservoir heterogeneities is required to reduce uncertainties and improve efficiency of HSA doublets.

Conclusions
On the basis of our calculations with detailed model realisations we can conclude that: I In the Nieuwerkerk Formation impermeable facies bodies form significant flow baffles perpendicular to the paleo flow direction, if the N/G is lower than 70%. II Lower pump energy losses can be expected when a geothermal doublet is oriented parallel to the paleo flow direction. III Equivalent permeability between doublet wells has a smaller anisotropy compared to equivalent permeability opposite model boundaries. IV The acquisition of geological data and the use of detailed facies architecture realisations are not negligible. Homogeneous realisations could significantly underestimate the geological risks of geothermal doublets. This study provides a workflow for reservoir engineers to determine the optimal doublet layout in HSA.