Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

Economic and fault stability analysis of geothermal field development in direct-use hydrothermal reservoirs

Economic and fault stability analysis of geothermal field development in direct-use hydrothermal... A.Daniilidis@tudelft.nl Faculty of Civil Engineering The installed capacity of geothermal systems for direct use of heat is increasing world- and Geosciences, section wide. As their number and density is increasing, the their interaction with subsurface Reservoir Engineering, Delft University of Technology, faults becomes more important as they could lead to safety risks from induced seismic- Stevinweg 1, 2628CN Delft, ity. Assessment and management of such risks is essential for the further development The Netherlands and extension of geothermal energy for heating. At the same time, the economic out- put of geothermal systems can be marginal and is hence often supported by subsidy schemes. A combined assessment of fault stability and economic output could help operators to balance economic and safety aspects, but this is currently not common practice. In this study we present a methodology to assess field development plans based on fault stability and Net Present Value (NPV ) using reservoir simulations of a fluvial, heterogeneous sandstone representative of the majority of direct-use Dutch geothermal systems. We find that the highest friction coefficient leading to exceed- ance of the Mohr–Coulomb failure criteria in this sandstone is 0.17; such values could be encountered in clay-rich fault gouges. Similar or lower fault permeability compared to the reservoir results in no changes and an increase respectively of both NPV and fault stability with larger Fault-to-Well Distance (FWD). Fault permeability higher than the reservoir permeability results in a minor increase in NPV with smaller FWD. Our results demonstrate that a combined analysis of thermal, hydraulic, mechanical and economic assessment supports a responsible and viable development of geothermal resources at a large scale. The importance of a high spatial density of supporting stress data will be essential for a better understanding and quantification of economic and fault stability effects of geothermal operations. Keywords: Direct-use geothermal, NPV, Economic, Fault stability, Field development Introduction Installed capacity and generated energy of geothermal systems for direct-use of heat is on the rise (Lund and Toth 2020) as they form a promising energy source that can deliver renewable, low-carbon, baseload energy for industrial and domestic usage. Fur- ther acceleration of installed capacity of these systems is envisioned in several countries (e.g. Knapek et  al. (2013), Schoof et  al. (2018), Guglielmetti et  al. (2019), Moscariello et al. (2020)). A higher density of direct-use geothermal operations can lead to potential © The Author(s) 2021. This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/. Zaal et al. Geotherm Energy (2021) 9:12 Page 2 of 26 interference between systems (Willems et  al. 2017b; Daniilidis et  al. 2021a) and more interaction with structural subsurface elements such as faults (Daniilidis et  al. 2020b). The proper assessment, evaluation and mitigation of safety risks is essential for local communities and the environment. One of these risks is that geothermal operations induce seismic events. Concerns related to the these or other disturbances could delay or even inhibit development of geothermal projects (Gaucher et  al. 2015; Buijze et  al. 2019b; Stauffacher et  al. 2015). The required acceleration and expansion of geothermal operations depend on the assessment and management of their potential safety risks (Buijze et al. 2019b; Knoblauch et al. 2019). This is essential for social acceptance of the increased use of geothermal energy. Additionally, the economic output of direct-use geothermal energy is subject to both subsurface and economic uncertainties (Daniilidis et al. 2017; Compernolle et al. 2019) and often requires subsidies to achieve profitability (Daniilidis et  al. 2017; Eyerer et  al. 2020). Increasing the number of geothermal developments requires a stable yield of prospective projects and low financial risks (Schoof et  al. 2018). Financial risks can be reduced by increasing the probability of success and avoiding low-yield projects. It is therefore important to investigate the profitability of a potential reservoir upfront for better investment decisions. For example, maximising flow rates has shown to be eco - nomically beneficial for operators as they have been demonstrated to lead to increased Net Present Value output even though this means an accelerated cooling of the reservoir (Daniilidis et al. 2020b, 2021a, b). Previous studies have found decreased heat production with a low permeability fault in close proximity and parallel to the geothermal doublet (Liang et  al. 2018). Interac- tion of direct-use geothermal operations with pre-faults in conduction-dominated geo- thermal settings has shown significant effects on energy generation, system lifetime and economic output (Daniilidis et al. 2016, 2020b). The proximity to a sealing fault not only influences these three success factors, it also affects the available reservoir volume for heat exchange (Daniilidis et al. 2016, 2020b). High fault permeability combined with the relative fault block positions can also have positive effects on energy and system life - time when connection to deeper parts of the reservoir is achieved through the fault (Daniilidis et al. 2020b). Moreover, depending on well spacing and fault permeability, the presence and proximity of faults can lead to positive interference between two doublets (Daniilidis et al. 2021a). Due to the low magnitude and small number of earthquakes related to geothermal production in sedimentary settings (Evans et  al. 2012; Buijze et  al. 2019b) where con- duction is the dominant heat transfer mechanism (Moeck 2014) seismic hazards from direct-use geothermal operations have not been addressed previously. However, in sedimentary systems induced seismicity related to gas production has been the focus of several studies (Buijze et al. 2019a; Van den Bogert 2015; Candela et al. 2019). There - fore, the focus on fault stability of pre-existing faults is more relevant for these geologi- cal settings. Conversely, coupled thermo-hydraulic-mechanical models have been used to describe fracture deformation under cold water injection and related stress changes (Salimzadeh et  al. 2018) or to assess the potential of Enhanced Geothermal Systems (EGS) (Lepillier et al. 2019). Such models have been further used to identify Coulomb- stress changes attributed to fluid injection as the main trigger for seismic events in an Zaal  et al. Geotherm Energy (2021) 9:12 Page 3 of 26 EGS setting (Troiano et al. 2013). More recently, the impact of re-injection of cold water into supercritical geothermal systems in terms of fault stability has been analysed (Pari- sio et al. 2019). Nonetheless, for conduction-dominated fields the combined influences of specific reservoir conditions and operational settings on the profitability and on the occurrence of fault instability of pre-existing faults, which might result in induced seis- micity, remains under-explored. This work addresses for the first time a single, unified assessment of fault stability and economic performance. We present a novel methodology that combines existing approaches for analysis of these two aspects. We assess the influence of fault presence and proximity to a direct-use geothermal doublet in sedimentary, conduction-domi- nated geothermal fields on economic output and fault stability for a range of fault prop - erties and operational choices. The Delft Advanced Research Terra Simulator (DARTS) (Khait 2019) is used to perform reservoir simulations of the Delft sandstone, a heteroge- neous fluvial sandstone, one of the main reservoirs for geothermal energy production in the Southern part of The Netherlands (Mijnlieff 2020). The simulation results are used to assess the economic output of a typical doublet based on an updated economic model that includes the current renewable-energy subsidy scheme for The Netherlands. The simulations resolve the pore-pressure change, while the fault stability is assessed based on the Mohr–Coulomb failure criterion. We simulate pre-existing faults positioned at a varying distance from the wells of a geothermal doublet, each fault with a different constant permeability. Operational choices considered include the flow rate and the re- injection temperature of the geothermal doublet. Our final analysis compares the gener - ated Net Present Value (NPV) and identifies the range of friction coefficients for which the Mohr–Coulomb stresses are exceeded, leading to potential fault instability. The out - line of the paper is as follows: “Geological setting” section provides a short introduction to the geological setting of the study area, “Delft Advanced Research Terra Simulator (DARTS)”,  “Fault-stability model” and “Economic model” sections describe the simula- tion methodology, fault stability model and economic model, respectively. “Results” section  provides the results of injection pressure, doublet power, spatial distribution of pressure and temperature, NPV, fault stability and combined NPV and fault stability. “Discussion” section places the findings in the broader context of geothermal develop - ment and “Conclusions” section gives the conclusions. Geological setting The target reservoir of this work is the Delft Sandstone Member (DSSM), which is part of the Nieuwerkerk formation in the West Netherlands Basin (WNB). The WNB is a 60 km wide trans-tensional basin, located in the southwestern parts of the Neth- erlands and contains stratigraphic formations of the Late Jurassic/Early Cretaceous (Pluymaekers et al. 2012; Bonté et al. 2012). The DSSM is layered and consists of fine to coarse gravelly sandstone sequences deposited by a north-west flowing meander - ing river and bands of silt- and claystones deposited as floodplains in a lower coastal plain setting (Van Adrichem  Boogaert and Kouwe 1993; Donselaar et  al. 2015). The deposits in the WNB have large thickness variations, which are the result of exten- sional and compressional faulting. Because of the present faults three main anticlinal and synclinal structures can be classified: the Delft High, the Pijnacker High and the Zaal et al. Geotherm Energy (2021) 9:12 Page 4 of 26 Fig. 1 Typical SW-NE geological cross-section through the WN, after van Balen et al. (2000). The location of the hypothetical reservoir considered in this study is tentatively indicated with a red rectangle. For the exact location of the cross-section, see van Balen et al. (2000) Fig. 2 Schematic of the analysis. Input data used are detailed in Table 1 Vrijenban Syncline which lies in between the two highs (Gilding 2010). The highest part of the Vrijenban Syncline near the Delft High is found at 1900 m true vertical depth in which the top of the DSSM is at 2200 m true vertical depth. The syncline is bounded by two main normal faults: the Delft main boundary fault in the southwest and the Pijnacker main boundary fault in the northwest. Additionally, the area con- tains several inverse normal faults and reverse faults. Movements of the two main faults and the syncline blocks along the fault reactivated and inversed the main faults, which also created new reverse normal faults in the area (Gilding 2010). A typi- cal SW-NE cross-section through the WNB, based on the model of van Balen et  al. (2000), is shown in Fig. 1. The location of the hypothetical reservoir considered in this study is indicated with a red rectangle in the Schieland Sandstone. Methods The present analysis is based on three discrete parts: reservoir simulations using the Delft Advanced Research Terra Simulator (DARTS), economic simulations using an economic model and fault-stability analysis with a model based on the Coulomb fail- ure criterion. An overview of the different parts and interconnections is presented in Fig. 2. Zaal  et al. Geotherm Energy (2021) 9:12 Page 5 of 26 Delft Advanced Research Terra Simulator (DARTS) The Delft Advanced Research Terra Simulator (DARTS) is a Python and C++ based dynamic reservoir simulator, developed at Delft University of Technology (Khait 2019; Voskov 2020) and has demonstrated fast and accurate performance for geothermal applica- tions (Wang et al. 2019, 2020). It is based on the framework of Operator-Based Lineariza- tion (OBL) (Voskov 2017), which makes it fit for the simulation of complex multiphase flow and transport systems (Khait 2019; Khait and Voskov 2018b). For the simulation of geothermal reservoirs three fundamental equations are used, which are the conservation of mass, conservation of energy and Darcy’s law for fluid flow in porous media (O’Sullivan et al. 2001). The equations are described in Voskov (2017), Khait and Voskov (2018b). The conservation of mass with n components and n phases is given c p by Eq.  1, in which the three terms denote respectively the change in mass with time, the mass flux due to flow and a sink- and source term: n  n n p p p φ x ρ s − ∇ · x ρ v + x ρ q˜ = 0, cj j j cj j j cj j j ∂t (1) j=1 j=1 j=1 c = 1, . . . , n , in which φ is the porosity, x the mole fraction of component c in phase j, s the phase cj j saturations, ρ the phase molar density, n the number of phases present, n the number j p c of components, q ˜ the phase rate per unit volume and v the phase velocity. j j The equation for the conservation of energy accounts for the stored energy, convection, conduction and sink and source terms: n  n n p p p φ ρ s U + (1 − φ)U − ∇ · h ρ v +∇ · (κ∇T ) + x h ρ q˜ = 0, j j j r j j j cj j j j ∂t j=1 j=1 j=1 (2) in which U is the phase internal energy, U the internal energy of the rock, h the phase j r j enthalpy, κ the thermal conductivity and ∇T the temperature gradient. Thermal convection is controlled by Darcy’s Law: rj v =− K (∇p − g ∇d) , j = 1, . . . , n , j j j p (3) in which K is the permeability tensor, k the relative permeability, µ the phase viscosity, rj j p the vector of pressures in phase j, g is the gravity term and ∇d is the vector of depths. j j All fluid properties are a function of both pressure and temperature according to the IAPWS97 formulation (Romera 2020). Equations  1,  2, 3 are transformed into Eqs.  4 and  5 of flow and transport for general multi-component fluids by applying finite-volume discretization and a backward Euler approximation in time (Voskov 2017; Khait and Voskov 2018b): Zaal et al. Geotherm Energy (2021) 9:12 Page 6 of 26 n   n    n  n p n+1 p n p p l l l l V φ x jρ s − φ x jρ s − �t x ρ Ŵ �ψ + �t ρ x q = 0, c = 1, ... , n , c j j c j j p cj j c cj j j j=1 j=1 l j=1 j=1 (4) and n   n p p n+1 n V φ ρ s U + (1 − φ)U − φ ρ s U + (1 − φ)U j j j r j j j r j=1 j=1 (5) n  n p p l l l l l l − �t h ρ Ŵ �ψ + Ŵ �T + �t h ρ q = 0, j j j j j j c j=1 j=1 in which V is the control volume, q , which is equal to q ˜ V is a source of phase j, �ψ is j j the difference in pressures between the blocks connected via interface l, T is the dif- ference in temperature between blocks connected via interface l, Ŵ , which is equal to l l Ŵ k rj is the phase transmissibility, Ŵ is the constant geometrical part of transmissibility, l l l l Ŵ , equal to Ŵ κ , equal to φ s γ − κ + κ the thermal transmissibility. r r c j=1 j j Equations 4 and 5 are approximated in time using a fully implicit method, which intro- duces non-linearity to the system (Voskov 2017). Linearisation in DARTS is done using the Operator-Based Linearisation approach (OBL) (Voskov 2017). OBL simplifies the complicated non-linear physics by parameterising the different operators in physical space with a limited number of supporting points, providing a piece-wise representation of the operators (Wang et al. 2019; Khait and Voskov 2018a, b). During simulations the operators are evaluated based on multi-linear interpolations, which reduces computa- tion time and improves the simulation of non-linear systems. The successful application of OBL in thermal multi-component and multi-phase flow simulations has been dem - onstrated in recent publications (Voskov 2017; Khait and Voskov 2018a, b; Wang et  al. 2019, 2020, Daniilidis et al. 2020a, Saeid et al. 2020). The well index is calculated accord - ing to the Peaceman well equations (Peaceman 1983), no skin effect is considered and flow inside the wells only considers buoyancy and does not include frictional pressure and heat transfer. Fault‑stability model The fault-stability model is based on the Mohr–Coulomb theory. Fault stability depends on the regional stress field (stress tensor), the orientation of the fault surface, its rock frictional characteristics (e.g. cohesive strength and friction coefficient) and the ratio between shear- and normal stresses acting on the fault (Morris et  al. 1996). A fault is considered to be stable if it does not reach the Mohr–Coulomb failure criterion, as described in Amonton’s law (Jaeger et al. 2009): τ = C + µ(σ − P), f ft (6) in which τ is the shear strength of the rock (MPa), C the cohesion of the rock (MPa), µ the apparent coefficient of internal friction (−) (Byerlee 1978), σ the total stress nor- ft mal to the failure plane (MPa), P the pore pressure (MPa) and σ − P the effective shear ft stress σ . ft Zaal  et al. Geotherm Energy (2021) 9:12 Page 7 of 26 According to Amonton’s law the stability or failure of a fault is determined by the ratio of shear to normal stresses acting on a surface and this also defines the slip tendency (Morris et  al. 1996). Failure is likely to occur when the resolved shear stress (τ ) e quals or exceeds the frictional sliding resistance, also called the failure shear stress (τ ) ( Jae- ger et  al. 2009). Hence the tendency for a pre-existing fault to fail is given by Eq.  7, in which the failure stress is proportional to the normal stress acting on the surface (C = 0, assuming no cohesion develops due to healing): τ ≥ τ = µ(σ ). f (7) ft Equation 6 gives the failure criterion describing the shear strength of the fault. In com- bination with the Mohr circle, the failure criterion determines the fault stability. We assume that failure occurs if the Mohr circle touches or intersects the line that describes the relation between the normal stress and the effective shear stress for a given friction coefficient, that is, where the effective shear stress (Eq.  7) equals or exceeds the failure shear stress of the fault. The maximum pore pressure values are computed based on the corresponding minimum required friction coefficient to reach the failure point and is visualised in Fig.  3a. In our setup, the pore pressure values used are the ones from the cells next to the fault plane. In this analysis we do not distinguish between seismic or aseismic slip, but only the exceedance of the failure criterion at any point of the fault and identification of the corresponding pore pressure change psimilar to for example Jansen et  al. (2019); Hettema (2020); Moeck et  al. (2009); Parisio et  al. (2019)] . In our study, poroelastic effects are not included. In the reservoir simulations (see “Reservoir model, sensitivity study and input parame- ters” section) we assume the presence of a vertical fault in a normal faulting regime. There - fore, failure in this case does not take place on the σ –σ plane, but rather on the σ –σ 1 3 2 3 Fig. 3 Mohr circles for a friction coefficient of 0.5 (black dashed line) showing the initial stress state (black circle) and the effective stress state (orange circle) with the pore pressure increase required for a critically dipping fault along the σ –σ plane (a). Pore pressure increase required for a vertical fault to reach failure 1 3 along the σ –σ plane (b) and effect of thermal stresses for a laterally confined reservoir where the thermal 2 3 stress acts on the horizontal components of the stress tensor with a T of 49.2 K and rock properties according to Table 3 (c). Pore pressure increase required to reach failure along the σ –σ plane of a thermally 2 3 stressed reservoir (d) Zaal et al. Geotherm Energy (2021) 9:12 Page 8 of 26 plane that defines the angle of the vertical fault plane with respect to the horizontal compo - nents of the stress tensor (Fig. 3b) resulting in strike-slip failure Markou and Papanastasiou (2018). Thermal stresses can cause significant changes in the stress tensor, especially around the injector well where the largest temperature differences occur. In this study we consider thermal stresses for a laterally constrained reservoir, therefore only acting on the horizontal components of the stress tensor (Fig. 3c). The thermal stress is evaluated based on (Jaeger et al. 2009): α E�T � = , (8) (1 − ν) in which � is the thermal stress change (MPa), α is the rock linear thermal expansion σ t coefficient, E is the rock Young’s modulus, T is the temperature difference (K) and ν is the rock Poisson’s ratio (−). The combined effect of thermal stresses and pore pressure is visualised in Fig. 3d. Based on the failure criterion, the critical fault strike angles can be evaluated, similar to Soltanzadeh and Hawkes (2009). Economic model The economic model investigates the development of a geothermal project and is imple - mented as outlined in previous research (Daniilidis et  al. 2017, 2020b). Multiple stages are included in the development among which the research phase, exploratory phase, development- and production phases and the abandonment phase. All phases contrib- ute to the total costs of a geothermal project, depending on multiple factors such as the geological uncertainties, government policies and required equipment (Daniilidis et  al. 2017; Schoof et al. 2018). The costs of the development phase include the drilling costs of two wells, which form the doublet and are a function of depth. For each well the drill- ing costs are computed as (TNO 2019): W = 375, 000 + 1150MD + 0.3MD , (9) cost in which MD is the measured depth (m) and the doublet costs are in Euro (€). Experi- ence from the first drilled well is estimated to result in decreased drilling costs of the second well by means of a learning curve reducing the costs by 7.5% (Lukawski et  al. 2014). The Net Present Value is based on the project expenses and revenues according to: NPV = � , (10) (1 + r) where C is the net cashflow generated at time t, r the monthly discount rate (%), t the elapsed number of time periods since the project start (t  =  0 is the time at which the wells are drilled). The net cashflow ( C ) in € is computed according to: C = I + I − CapEx − OpEx − OpEx − AbEx, t heat subsidy annual pumping (11) where I is the income from the produced heat, I is the subsidy income, CapEx heat subsidy are the capital expenditures, OpEx the annual operational expenses, OpEx annual pumping Zaal  et al. Geotherm Energy (2021) 9:12 Page 9 of 26 the pumping expenses and AbEx the abandonment cost. The CapEx consist of the total investment costs of the project and include the exploration, drilling costs (W ), con- cost struction and unforeseen construction costs associated with preparation, drilling and installation of the equipment and drill site. Equipment includes the injector and pro- ducer pumps, the heat exchanger and the screens and filters installed in the well. The OpEx consist of all costs made during production and include an annual component as a percentage of CapEx and the expenditures for running the injector and producer pumps, according to the electricity price. The annual OpEx include the maintenance- and work - over costs, rent of facilities, staff payment and insurance costs and amounts to 4% of the CapEx per anum for the annual costs for inspection, maintenance, workovers and costs for staffing (de Boer et al. 2016) and 1% of the CapEx per anum for heat network mainte - nance to ensure the heat is transported to the district heating system; the district heating itself is considered as pre-existing and its maintenance is not included in this analysis. These percentages remain constant during the project’s lifetime. The income not gener - ated during maintenance and workovers are added to the fixed OpEx in the model to take into account the effects of the downtime (Daniilidis et al. 2017). For the insurance the model incorporates a premium paid by the operator, which is a one-time payment of 7% of the maximum compensation available for a typical geothermal reservoir at 2000 m depth (Rijksdienst voor Ondernemend Nederland 2019a). The Weighted Average Cost of Capital (WACC) is used as the discount rate for the NPV. This is a weighted average of the costs of capital and debt. Lenders and equity holders expect to receive a certain return on the capital or funds they provided. This expected return is given by the WACC and is calculated as : E D WACC = Re + Rd(1 − Tc), (12) V V in which Re is the cost of equity, Rd is the cost of debt, V the total market value of the project’s financing (E + D) and Tc the corporate tax rate, E the market value of equity (0.3) and D the market value of debt (0.7). In the Netherlands, projects that may be sub- sidized with the SDE+ regulation are financed in a 70/30 debt over equity ratio. This all results in a WACC value of 5.4%. Dutch subsidy The project revenues include the income generated from the delivered heat and from the Sustainable Duurzame Energieproductie (SDE+) subsidy. In the Netherlands the heat price is determined as 90% of the Dutch market price of gas (TTF) (Ministerie van Econ- omische Zaken en Klimaat 2019) and this heat price determines the generated income and the subsidy amount. The TTF price is based on the High Heating Value (HHV) of gas. The incorporated gas price in the model is based on statistically historical data (Schoots and Hammingh 2019). The SDE+ subsidy provides financial compensation for up to 6000 full-load production hours over a maximum period of 15 years (Rijksdienst voor Ondernemend Nederland 2019b) for projects with depths between 500 and 4000 m. The subsidy payouts are based on the yearly amount of power required for produc - tion, which is fixed, the market price of thermal energy and the full-load production Zaal et al. Geotherm Energy (2021) 9:12 Page 10 of 26 hours. Based on this, the model takes into account a final subsidy amount of €0.033 per kWh and this is implemented on a yearly basis for 15 years. Energy evaluation The produced energy is computed according to (Daniilidis et al. 2017): E = q�TC ρ �t, produced f f (13) where q is the volumetric flow rate (m /s), T the temperature difference between the produced and injected water (K), C the specific heat capacity of the fluid at constant pressure (J/(kgK)), ρ the density fluid (kg/m ) and t the time interval (h). Heat losses at surface and in the heat exchangers are not taken into account in this energy evaluation. The required pumping energy is determined as the product of the flowrate and the pressure difference between the wells, integrated over time (van’t Spijker and Ungemach 2016): q�P E = �t, pump (14) where q is the volumetric flow rate ( m /s ), P the pressure difference between injector and producer wells (Pa) and η the pump efficiency taken here as 55%. The governing equations are implemented in the model as presented by (Daniilidis et  al. 2017, 2020b) and use production fluid temperature and well pressure from the DARTS simulations. Reservoir model, sensitivity study and input parameters The temperature, pressure, NPV and fault stability are analysed using a 3D heteroge - neous reservoir model for the Delft area (Fig.  4) comprised of approx. 100 k elements. Fig. 4 Map view of the 3D reservoir model of the heterogeneous Delft Sandstone. The injection well is indicated in blue and the production well indicated red, well spacing is 990 m. Each cell is 30 m by 30 m by 2.5 m (x, y, z). The vertical fault is indicated by the purple line (shown here with a permeability of 0.75 mD), has a width of 30 m and extends throughout the domain along the x-axis. The different fault permeability values used in the analysis are listed in Table 3 Zaal  et al. Geotherm Energy (2021) 9:12 Page 11 of 26 The geological setting has been described in “ Geological setting” section and references therein. The reservoir model represents a heterogeneous fluvial system with varying net-to-gross (N/G). It includes different fluvial deposits that have been generated by the process-based model Flumy and are classified as reservoir (points bars, sand plugs and channel lags) and non-reservoir (crevasse splays, overbank alluvium and mud plugs). Using petrophysical well data, the following permeability-porosity relationship is used, after Willems et al. (2017b): 29.507φ K = 0.0633e . (15) The reservoir model, as well as the porosity and permeability distribution have been described extensively in Willems et al. (2017b) and have also been used in previous stud- ies (Wang et  al. 2019; Shetty et  al. 2018; Daniilidis et  al. 2020a). The top and bottom boundary conditions act as infinitely large domains with a fixed temperature and pro - vide conductive thermal recharge to the reservoir (Daniilidis et al. 2020a). We perform a sensitivity study consisting of a number of cases, each with a different choice of reservoir properties and/or operational settings. The input parameters that are the same for each of the cases simulated with the reservoir model are given in Table 1. The well locations remain constant throughout the analysis and the assumed vertical fault is positioned with respect to the wells, changing the Fault-to-Well Distance (FWD); the fault strike is parallel to the plane defined by the wells and therefore both wells have the same distance to the fault. Using a vertical fault is considered as an acceptable approximation in many fault stability studies Jansen et  al. (2019), Van  den Bogert (2015) for reservoir faults in normal faulting regimes where faults are often steeply dipping. Table 1 Input data for the reservoir model Parameter Value Unit nx × ny × nz 60 × 40 × 42 – dx × dy × dz 30 × 30 × 2.5 m Permeability range 1–3400 mD Porosity range 5–30 % Location producer 450,600 m (X,Y) Location injector 1440,600 m (X,Y) Fault permeability See Table 2 mD Fault thickness 30 m −15 2 Overburden permeability 9.9 × 10 (0.001) m (mD) Top reservoir depth 2200 m Model reservoir thickness 105 m Fault-wells distance see Table 2 m Uniform initial pressure 20.1 MPa Uniform initial temperature 345.75 K Flowrate prod/inj 4200 m /day Injection temperature see Table 2 K Simulation time 30 × 365 days Thermal conductivity rock 100 kJ/m/day/K Heat capacity rock 2200 kJ/m /K Zaal et al. Geotherm Energy (2021) 9:12 Page 12 of 26 Table 2 Varying parameters of the reservoir model. The lowest fault permeability of 0.75 mD is referred to as a sealing fault and the highest fault permeability of 750 mD as transmissive fault in “Results” section Parameter Value Unit Flowrate 150, 225, 300, 375 m /h Fault permeability 0.75, 7.5, 75, 750 mD Injection temperature 303.15, 308.15, 313.15 K Fault-to-Well Distance (FWD) 60, 90, 180, 360 m Table 3 Input data for the fault model Parameter Value Unit σ 45 MPa σ 34 MPa σ 30 MPa Young’s modulus 16 GPa Poisson’s ratio 0.25 – Cohesion coefficient 0 MPa Inital pore pressure 20.1 MPa The stress tensor values are based on the work of Mechelse (2017) Table 4 Economic model input Parameter Value Unit Occurrence Source Production time 30 Years – – No. of doublets 1 – – – No. of wells 2 – – – Exploration costs 460,000 € Once van den Bosch et al. (2013) ESP costs 800,000 € Every 5 years van’t Spijker and Ungemach (2016) Power facility 1,500,000 € Once van den Bosch et al. (2013) Drill location 300,000 € Once van den Bosch et al. (2013) Unforeseen costs 5 % of construction costs Once Daniilidis et al. (2017) ESP downtime 15 days Every 5 years Daniilidis et al. (2017) ESP efficiency 55 % – van ’t Spijker and Ungemach (2016) Electricity price 0.071 €/kWh – CBS (2019) Gas price 0.016 €/kWh – - Fixed OpEx 5 % of CapEx Annually de Boer et al. (2016) Drill insurance 7 % of maximum com- Once Rijksdienst voor Ondernemend pensation available Nederland (2019a) $/€conversion 0.90 – – - AbEx single well 1,275,500 € Once Osundare et al. (2018) Market value of equity 0.3 – – Rijksdienst voor Ondernemend Nederland (2019b) Market value of debt 0.7 – – Rijksdienst voor Ondernemend Nederland (2019b) Cost of equity 14.5 % – Rijksdienst voor Ondernemend Nederland (2019b) Cost of debt 2 % – Rijksdienst voor Ondernemend Nederland (2019b) Corporate tax 25 % – Rijksoverheid (2019) Zaal  et al. Geotherm Energy (2021) 9:12 Page 13 of 26 For each of the cases in the sensitivity study, a simulation is performed for a combina- tion of parameters for both reservoir conditions and operation options. These param - eters are provided in Table 2. The input parameters of the fault stability model are given in Table 3. The inputs of the economic model are all based on Dutch financial regulations. The input used in the model is given in Table 4. Results In this section, we present the results of the sensitivity study as described in “Reservoir model, sensitivity study and input parameters” section. This study is performed with the models described in “Delft Advanced Research Terra Simulator (DARTS)”, “Fault-stabil- ity model” and “Economic model” sections following the workflow as presented in Fig. 2. In the following, we will subsequently address the resulting reservoir pressure and tem- perature as well as the NPV and fault-stability analysis for the reservoir and operational conditions under consideration. This section concludes with the combined results of the NPV and the fault stability. Injection pressure Injection pressure increases over time with most significant changes taking place within the first 3 years as pressure dissipates inside the reservoir and stabilises after about 20 years of production (Fig.  5). The injection rate defines the level at which the plateau occurs, with higher rates resulting in increased injection pressure. The distance to the fault influences the final pressure reached; with an FWD of 60 m away from the wells an increase of up to 2 MPa compared to a distance of 390 m for the lowest fault perme- ability. With an FWD of 390 m the fault permeability has no effect on the maximum injection pressure encountered. However, when the fault is closer to the wells the lower fault permeability of 0.75 mD results in higher injection pressure. This effect is less pro - nounced for a fault permeability of 7.5 mD and is minimal for the fault permeability of 75 mD; this indicates that a fault permeability of 75 mD is close to the average reservoir permeability and therefore the distance of the fault is not significant as the fluid passing the fault behaves hydraulically similar to the fluid in the reservoir layers. A fault perme - ability of 750 mD results in a decrease of injection pressure when the fault is positioned closer to the wells, since the fault provides less resistance to flow than the reservoir layers. Doublet power Doublet power shows a slight decrease over time due to production temperature decrease (Fig. 6). The power starts decreasing after 15 years for the low rate of 150 m /h with minor reductions in the order of 0.3 MW and after circa 7 years for the high rate of 375 m /h with reductions of up to 3 MW. Low fault permeability of 0.75 mD leads to earlier power drop with a FWD of 60 m. As the FWD is increasing the power loss over time is decreasing. The discrepancy caused by the FWD is decreasing as the rates decrease and the fault permeability increases. With a fault permeability of 75 mD the discrepancy caused by the FWD is minimised, similar to the injection pressure (Fig. 5). A fault permeability of 750 mD results in an earlier power drop at circa 5 years for the Zaal et al. Geotherm Energy (2021) 9:12 Page 14 of 26 Fig. 5 Injection pressure over time for fault permeability of 0.75 mD, 7.5 mD, 75 mD and 750 mD, FWD of 60 m, 90 m, 180 m and 360 m and injection temperature of 303.15 K FWD of 60 m. With the increase of FWD the onset of the power drop is delayed and is also less pronounced. Even though the power drop decreases as the FWD increases for 3 3 both high (375 m /h) and low (m /h) rates, the mechanism is different depending on the fault permeability (see also Fig. 8). For low fault permeability (0.75 mD) the volume available to to cold water plume is reduced and reaches the producer well earlier. Con- trary to this, for high fault permeability (750 mD) the fault acts as conduit connecting the injector to the producer well. Spatial distribution of temperature and pressure changes The pressure distribution inside the reservoir is similar for a case with a transmis - sive fault and a case with a sealing fault with an FWD of 390 m (Fig.  7a, b). For both cases, the maximum and minimum pressure values are observed at the location of the injector and producer wells respectively. At an FWD of 60 m, a sealing fault limits the available volume over which pressure can dissipate, resulting in pressure build- up around the injector and along the fault plane close to the injector (Fig.  7c). This creates a large pressure differential across the fault plane. Similarly, a lower pressure is observed around the producer and along the fault close to the producer. With the fault positioned 60 m away from the wells, the transmissive fault allows for the dis- sipation of pressure due to its higher permeability relative to the sealing-fault case (Fig. 7d). This implies that the pressure difference across the fault (north block) is less pronounced in this case compared to a sealing fault at the same distance, and so is the pressure difference between the areas of the injector and producer wells. This leads to relatively lower pressure values around the injector and a weaker pressure contrast between the injector and the fault, compared to case b with the a transmissive fault at an FWD of 390 m. The likely explanation for this is that the fault provides a fluid pathway with less hydraulic resistance and the pressure dissipates across the fault. Figure  8 shows the temperature distribution inside the reservoir in the presence of a transmissive and a sealing fault with an FWD of 60 m and 390 m. A sealing fault at an FWD of 390 m restricts the extent of the cold plume on the north side of the fault (Fig.  8a) compared to a transmissive fault (Fig.  8b). In the latter case the increased fault permeability allows the cold plume to permeate the fault, reducing the temper- ature on both sides of the fault. As such, it allows the cold plume to reach further towards the north of the reservoir. At an FWD of 60 m, a sealing fault allows mostly conductive heat transfer towards the north fault block while convective heat transfer is limited (Fig.  8c). This results in a reduced extension of the cold plume to the north and a restriction of the plume Zaal  et al. Geotherm Energy (2021) 9:12 Page 15 of 26 Fig. 6 Doublet power over time for fault permeability of 0.75 mD, 7.5 mD, 75 mD and 750 mD, FWD of 60 m, 90 m, 180 m and 360 m and injection temperature of 303.15 K Fig. 7 Map view at a depth of 2177.5 m showing the pressure distribution in the DSSM reservoir after 30 years of production with an FWD of 390 m (a, b) and 60 m (c, d) for a sealing (a, c) and a transmissive fault (b, d). Data shown use an injection temperature of 303.15 K and a flowrate of 9000 m /day . Injector and producer wells are marked with a blue and red point respectively and the fault position is indicated in black around the injector where the higher pressure gradients occur. In this location, the cold plume initiates and favours conductive transfer. A transmissive fault at an FWD of 60 m (Fig. 8d) diverts a large part of the flow and results in an extension of the cold plume along the fault direction and towards the producer, leading to earlier thermal breakthrough. In this case, the connection between the two fault blocks is strong and resembles that of the FWD of 390 m, albeit slightly restricted due to the diversion through the fault plane. The vertical reservoir profile along the injector for an FWD of 390 m reveals only minor pressure differences between the case with a sealing fault (Fig.  9a) and the case with a transmissive fault (Fig. 9c) as illustrated in Fig. 9b, d respectively. The sealing fault with an FWD of 390 m results in slightly higher pressure closer to the fault plane. At Zaal et al. Geotherm Energy (2021) 9:12 Page 16 of 26 Fig. 8 Map view at a depth of 2177.5 m showing the temperature distribution in the DSSM reservoir after 30 years of production with an FWD of 390 m (a, b) and 60 m (c, d) for a sealing (a, c) and a transmissive fault (b, d). Data shown use an injection temperature of 303.15 K and a flowrate of 9000 m /day. Injector and producer wells are marked with a blue and red point respectively and the fault position is indicated in black. The notation on the bottom right shows the volume percentage that exhibits a �T > 1.5 K after 30 years of simulation compared to the initial reservoir temperature. The percentage is calculated throughout the whole reservoir volume and it is not limited to the horizontal cross-section shown here an FWD of 60 m the sealing fault significantly increases the pressure between the well and the fault but also around the injector away from the fault, towards the south of the domain (Fig.  9f ) across the full length of the well. The higher pressure in the south in this case compared to the case of the transmissive fault extends all the way to the edge of the domain. The case with a permeable fault at an FWD of 60 m exhibits lower pressure values around the injector, even compared to the case where the fault is located further away (Fig.  9h). The relatively higher permeability of the fault in this case offers an easy pathway to the fluid and reduces pressure buildup. NPV The generated NPV after 30 years of production for the cases with the sealing fault (0.75 mD) shows minor variations at the lowest flow rate (150 m /h) for each injection tem- perature depending on the fault distance to the wells in these cases (Fig.  10). We find that for the lowest flow rate of 150 m /h an injection temperature of 308.15 K or lower is required to achieve positive NPV values. With higher rates, the NPV becomes more sensitive to changes in FWD. For the lowest injection temperature of 303.15 K, increas- ing the rate from 150 to 375 m /h (a factor 2.5) increases the NPV from 5 to 22.5 M€ Zaal  et al. Geotherm Energy (2021) 9:12 Page 17 of 26 (a factor 4.5). In the case of the highest rate of 375 m /h the NPV for an FWD of 60 m is equivalent to maintaining an FWD of 390 m and selecting an injection temperature that is higher by 5 K (308.15 K instead of 303.15 K). The NPV differences caused by the fault proximity are smaller for cases with a higher fault permeability of 7.5 mD across all injection temperatures and flow rates investigated. At a fault permeability of 75 mD the fault distance does not seem to affect the generated NPV and differences in flow rate and injection temperature have a relatively stronger effect. Cases with a higher fault perme - ability to 750 mD show small NPV benefits as the fault distance to the wells is larger. The discrepancy between the cases with a sealing (0.75 mD) and a transmissive fault (750 mD) with lower FWD can be attributed to the lower pumping costs due to lower pres- sure (see also Fig.  9). Possibly, the relatively larger reservoir volume available for heat exchange in the high fault-permeability cases also contributes to a higher NPV, since cold water is averted away from the producer well through the fault and as such pen- etrates a larger volume—64% of volume affected compared to 56% for the transmissive and sealing case respectively at an FWD of 60 m (see also Fig. 8). Fault stability The pore pressure increase required to reach failure shows a similar trend for the ini - tial stress state and the thermal stress states of the injection temperatures considered (Fig.  11a). The σ and σ are equally affected by the thermal stress and therefore the 2 3 shapes of the pore pressure curves remain similar, with only the values of the pore pres- sure leading to failure is increasing. The critical angle between the maximum horizon - tal stress S (which is σ ) and the fault strike shows that the angle range of 25° to H 2 max 43° is critical depending on friction coefficient, with higher friction coefficients leading to smaller angles (Fig.  11b). The thermal stresses slightly shift the critical angle but the effect is minor. The maximum encountered pore pressure derived from the reservoir simulations is higher for cases with a small Fault-to-Well Distance and higher flow rates (Fig.  12). Lower injection temperature further increases the maximum pore pressure with up to an additional 1.5 MPa for the highest flow rate and shortest distance between the sealing fault and the wells. For the highest rates, the effect of the increase in injection tempera - ture by 10 K is almost equivalent to the effect of a decrease in flow rate of 75 m /h. For an injection temperature of 303.15 K and assuming a friction coefficient of 0.15 for the sealing fault, the Mohr–Coulomb failure criterion for pore pressure is exceeded with a 3 3 60 m and 90 m FWD with a rate of 375 m /h and with a 60 m FWD for the 300 m /h rate. A higher injection temperature of 313.15 K would require an even lower friction coefficient of 0.12 for the sealing fault to fail at similar distances and rates. It should be noted that experimental results show a pure quartz gouge to have a friction coefficient between 0.66 and 0.75 (Tembe et al. 2010) while friction coefficient values below 0.3 are related to clay rich minerals (Ikari et  al. 2009). For higher fault permeability cases, the failure criterion is only exceeded when the fault is closer to the wells (FWD below 100 m) assuming the highest rates and a friction coefficient of 0.15 and 0.12 for the injection temperature of 303.15 K and 313.15 K respectively. For a fault permeability of 750 mD friction coefficients below 0.15 and 0.12 are required to reach failure for 303 K and 313 K accordingly. Zaal et al. Geotherm Energy (2021) 9:12 Page 18 of 26 Fig. 9 Vertical cross-section of the reservoir at the injection well location showing the permeability (left column) and corresponding pressure distribution (right column) in the DSSM reservoir. A sealing fault at an FWD of 390 m (a) and 60 m (c) from the injector well and a transmissive fault at an FWD of 390 m (b) and 60 m (d) distance from the injector well are shown. The injector well is marked by a blue line. Data shown for an 303.15 K injection temperature and a flowrate of 9000 m /day Combined NPV and fault stability analysis Combining the NPV results and the results of the maximum pore pressure increase for the highest considered injection temperature cases (313.15 k) shows that with a seal- ing fault (0.75 mD) a lower FWD results both in a lower NPV and a higher maximum pore pressure difference (Fig.  13). The trend is similar for all rates but the differences are larger for cases with larger flow rates. The higher rates and shorter distances lead to an exceedance of the failure criterion with a friction coefficient of 0.17 and just below 0.14 for an injection temperature of 303.15 K and 313.15 K respectively. When fault perme- ability is 75 mD, NPV is no longer affected by the distance to fault and only the rates affect the NPV and maximum pore pressure difference. The trend changes with a trans - missive fault (750 mD) and we observe higher NPV values and maximum pore pressure differences for lower FWDs. The NPV changes, however are minor. Lower injection temperature of 303.15 K leads to higher NPV values and also higher maximum pore pressure differences. Discussion Our results support previous findings of increased NPV for higher flow rates and lower temperatures (Daniilidis et  al. 2016, 2020b, 2021a; Willems et  al. 2017a). For the flow rate ranges considered here and assuming current heat and electricity prices, the addi- tional revenues from high flow rates outweigh the costs of increased pumping require - ments. The NPV results also reflect that either option of increasing the rates or lowering the injection temperature can be used to improve the generated economic value when all other parameters remain the same (Daniilidis et al. 2016). In our simulations, keeping the FWD large is beneficial for both NPV and fault stability when fault permeability is higher than most reservoir permeability values (750 mD). Nonetheless, higher rates have a stronger positive effect on the NPV compared to the effects of a larger FWD. A smaller FWD results in a lower NPV and therefore there is no financial incentive for operators Zaal  et al. Geotherm Energy (2021) 9:12 Page 19 of 26 Fig. 10 Comparison of the generated NPV for the range of injection temperature, distance between the fault and the wells, flowrates and fault permeability considered to knowingly position doublets closer to a fault. If faults are of sub-seismic resolution and their existence or location is not known, smaller FWDs might of course still occur. However, if fault permeability is higher than the assumed reservoir permeability, only minor benefits can be achieved by reducing the FWD. Previous studies have also sug - gested that effect of the FWD to system lifetime and NPV could be optimised if the hydraulic behaviour of both the fault and reservoir is well characterised (Daniilidis et al. 2021a). A combined assessment of thermal, hydraulic, mechanical and economic aspects can be foreseen as a future research path for direct-use geothermal developments. With respect to fault stability we identify a friction coefficient of 0.17 and just below 0.14 for an injection temperature of 303.15 K and 313.15 K respectively as the high- est values at which the pore-pressure change from operations results in exceeding the Mohr–Coulomb failure criterion. The exceedance only occurs when the highest flow rate (375 m /h), short FWD (60 m) and lowest fault permeability (0.75 mD) are used. Temperature effects prove important even without the inclusion of poroelasticity as the temperature effect on fluid density and viscosity alone is sufficient to differentiate the interaction with the Mohr–Coulomb failure criterion. Additionally, compared to the ini- tial stress state, the thermally induced stress further reduces the distance to failure. A difference of up to 3 MPa for the same friction coefficient can be observed in the pore pressure increase required to reach failure for the temperature range considered (303.15 K to 313.15 K); this difference however is minimized at lower friction coefficients. The pore pressure increase itself depends strongly on the fault permeability. Additionally, our analysis does not examine whether the increased pore pressure would result in fault slip as the Mohr–Coulomb criterion can also be exceeded aseismically. Therefore, our evaluation of fault stability serves as a worst-case estimate within the methodological choices presented in “Fault-stability model” section and the initial stress values consid- ered. A determination of the initiation point of seismic slip that includes aseismic effects would be needed to assess actual seismic risk. It should be noted that even though it was not the focus of this work, the high injection rates of 375 m /h result in pressure levels in excess of σ (30 MPa), which could lead to shear failure of the reservoir rock close to the injector well. Experimental research has identified friction coefficients below 0.35 in faults that are rich in clay minerals (Ikari et  al. 2009). These faults exhibit only stable slip unless sig - nificant mineralogy changes are incurred (Ikari et  al. 2011). Other experimental stud - ies report friction coefficients below 0.3 for clay content above 70% (Tembe et al. 2010; Takahashi et  al. 2007). Pure quartz gouge has been experimentally estimated to have Zaal et al. Geotherm Energy (2021) 9:12 Page 20 of 26 Fig. 11 Pore pressure required to reach failure along the σ –σ plane for the initial stress state and the 2 3 thermal stress states considered for a range of friction coefficients (a). The critical fault strike angle with respect to S for the different stress states for the range of friction coefficients (b) max a friction coefficient between 0.66 and 0.75 depending on the grain size (Tembe et  al. 2010). While the target reservoir used in this study is a sandstone, the Delft Sandstone is known to exhibit large ranges of Net-to-Gross and hence also contains significant amounts of shale and therefore clay minerals (Donselaar et al. 2015; Donselaar 2016). As a result, the presence of clay minerals inside the fault gouge cannot be ruled out. Moreo- ver, fluid circulation in close proximity to the fault could initiate or intensify migration of fines inside the fault plane and aid the onset of clay smearing (Vrolijk et  al. 2016). Nonetheless, since clay-smearing potential is proportional to fault throw and inversely proportional to fault-zone thickness (Vrolijk et al. 2016), faults exhibiting clay smearing should be easier to characterise and considered in the design of well trajectories. The risk of induced seismicity increases with the injected volume but also with the proximity of the injection point to the fault (Zang et al. 2014; Davis and Frohlich 1993). Similarly, in our study we find that the MC failure criterion is exceeded earlier with higher flow rates, smaller FWD and lower fault permeability. Stress induced by fluid injection can alter the stress distribution around faults (Altmann et  al. 2014). While a low-permeability fault results in sharper pressure (dP) and temperature differences (dT), a high permeability fault could lead to a much larger volume being perturbed by dP and dT and this could further increase the potential for seismicity. Nonetheless, induced seismicity also has a temporal component and cessation of operations does not mean an instant reduction in seismic risk (Gaucher et al. 2015). Fault reactivation related to fluid injection or extraction is primarily a result of changes in pore pressure, but poroelas- tic and chemical effects can also play a role and should be included for a more detailed understanding. However, by considering only vertical faults, the poroelastic effects are minimized (Hettema 2020). Fault proximity is relevant for the expected NPV of geothermal operations, and thus an important aspect for operators to consider, even if the regulation framework does not impose any restrictions on placement of the wells relative to faults or does not require the monitoring and study of mechanical effects. As such, this economical aspect pro - vides an additional argument for the characterisation of faults present in the area of interest. Establishing the friction coefficient for fault surfaces and fault transmissivity is crucial to enable a more robust assessment of seismic risk from geothermal operations. For most geothermal regions, a more detailed spatial analysis of the stress field mag - nitude and orientation, at a resolution higher than presently available, is crucial. This should be complemented with fault-gouge friction coefficient characterisation in the lab. Zaal  et al. Geotherm Energy (2021) 9:12 Page 21 of 26 Recent studies in the Molasse basin in Germany have also pointed towards the need for further stress-field characterisation (Seithel et al. 2019). Also in The Netherlands, obser - vational efforts in support of a stress-field characterisation should be intensified, as the current stress-data density impairs a more detailed assessment. The methodology presented in this paper could also be adapted and used for EGS or fractured hydrothermal systems. In these cases, a different reservoir characterization and model would need to be developed [similar to e.g. Lepillier et  al. (2019)], the con- version of the heat to electricity would need to be modelled and an economic model for the combined electricity- and heat production would need to be developed. Moreover, poroelastic effects would need to be included as their impact is much stronger in the geological settings where geothermal systems are used for electricity production. Even if more data on both the friction coefficient and stress orientation and magni - tude are available, the evaluation of all potential development options and stress condi- tions remains a major undertaking requiring considerable simulation efforts. Additional aspects such as the well positioning within the considered reservoir permeability field could further increase complexity of the required studies. Moreover, fault stability could also be assessed for faults with varying dip values in relation to both the wells and the stress field. The inclusion of mechanical and eventually chemical coupling in the sim - ulation can severely increase the computational expense (Pandey et  al. 2018). Recent improvements in computational efficiency, for example with the use of GPUs in DARTS (Khait et al. 2020), open up possibilities for a more detailed simulation, including these processes but also allowing for uncertainty quantification on multiple scales. This type of simulators will also open the pathway to studies of coupled thermo-hydro-mechanical– chemical processes. Fig. 12 Comparison of the maximum encountered pressure against fault to well distance for the range of flowrates, fault permeability and the highest and lowest injection temperature considered. The horizontal lines show the failure pressures as result of friction coefficients 0.15, 0.17, 0.19 and 0.12, 0.14, 0.16 for the injection temperature of 303.15 K and 313.15 K respectively Zaal et al. Geotherm Energy (2021) 9:12 Page 22 of 26 Fig. 13 Combined results for NPV and maximum encountered pore pressure difference as a function of FWD, flowrate and fault permeability. Both the highest and lowest injection temperature are considered. The horizontal lines show the pore pressure for failure-criterion exceedance assuming friction coefficients of 0.15, 0.17, 0.19 and 0.12, 0.14, 0.16 for the injection temperature of 303.15 K and 313.15 K respectively Conclusions This study presents a methodology to combined the assessment of the NPV and fault- stability analysis for direct-use geothermal energy production that is demonstrated in the Dutch context and is directly applicable to the WNB. THe methodology presented can be applied with minor changes to other sedimentary, direct-use systems and can be adapted to other settings with different geological and economic conditions. Reservoir simulations with a heterogeneous reservoir model representing the fluvial Delft Sand - stone considered the effect of a range of flow rates and injection temperatures. By plac - ing a vertical fault at varying distance from the fixed well positions and with different fault-permeability values the role of these conditions was evaluated. An NPV analysis was carried out using an updated economic model of the Dutch geothermal subsidy context. With simulations we resolved the pressure field inside the reservoir and close to the fault. The resulting pore pressure change was used as input for the Mohr–Coulomb failure criterion to identify the friction coefficient that results in unstable behaviour. Our main findings from these evaluations for a heterogeneous geothermal reservoir include: • both NPV and fault stability improve with a larger Fault-to-Well Distance (FWD) when fault permeability is lower than the average reservoir permeability, • fault permeability higher than the average reservoir permeability results in relatively small NPV benefits when placing the wells closer to the fault, • fault permeability at similar levels as the average reservoir permeability results in no changes in the generated NPV when choosing a different FWD, • a sealing fault limits the reservoir volume for fluid flow but heat conduction across the fault remains significant over the simulated period of 30 years, • compared to a case with a transmissive fault, a case with a sealing fault at close prox- imity to the injection well significantly alters the pressure field both between the injector and the fault but also further away from the injector inside the reservoir domain, Zaal  et al. Geotherm Energy (2021) 9:12 Page 23 of 26 • the fluid density and viscosity changes caused by a re-injection temperature decrease of 10 K (313.15 K to 303.15 K) are sufficient to cause changes in fault stability by shifting the maximum pore pressure encountered by up to 1 MPa compared to a case with higher temperature, and • the highest friction coefficient leading to exceedance of the Mohr–Coulomb failure criterion under the geological assumptions of our study was 0.17 and just below 0.14 for an injection temperature of 303.15 K and 313.15 K respectively. Such low values of friction coefficient can be encountered in clay rich fault gouges. • the presence of a transmissive fault can drastically alter the resulting shape of the cold-water plume by diverting and connecting different channelised parts of the res - ervoir. Future directions of our research include the analysis of the fault impact with variable well positions to better explore the impact of heterogeneous permeability and its inter- action with the fault hydraulic properties. Additionally, the combination of a heteroge- neous reservoir and heterogeneous fault properties will be further explored for its effects on both NPV and fault stability. We expect that the combined assessment of thermal, hydraulic, mechanical and economic aspects will become increasingly relevant as geo- thermal operations expand and the density of direct-use geothermal systems increases. A comprehensive stress-field characterisation with increased spatial density of support - ing data will be essential for a better understanding and quantification of the economic and fault-stability effects of geothermal operations. Acknowledgements The authors would like to acknowledge David Bruhn for providing feedback on the manuscript draft and Anne Pluymak- ers and Phil Vardon for their insights on geomechanics. Additionally, we would like to thank our three reviewers for their detailed, critical and constructive feedback and comments. Authors’ contributions CZ developed the fault stability model and implemented the economic model, performed the simulations in DARTS and contributed to the visualisation and the writing of the manuscript. AD conceptualised the idea for the work, supervised CZ and contributed to the development of the fault stability and economic model, preparing the visualisations, perform- ing the simulations in DARTS and writing the manuscript. FV supervised CZ and provided feedback on the DARTS simu- lations, fault stability and economic analysis, contributed to the visualisations and to writing the manuscript. All authors of this article contributed to reviewing the manuscript. All authors read and approved the final manuscript. Funding This study is partly supported by the DESTRESS project, funded by the European Union’s Horizon 2020 research and innovation program under Grant agreement ID 691728. FV was funded by the Delft Technology Fellowship of Delft University of Technology. Availability of data and materials The dataset generated and/or analysed during the current study is available from the corresponding author upon reasonable request. Declarations Competing interests The authors declare that they have no competing interests. Received: 18 December 2020 Accepted: 8 April 2021 References Altmann J, Müller B, Müller T, Heidbach O, Tingay M, Weißhardt A. Pore pressure stress coupling in 3d and consequences for reservoir stress states and fault reactivation. Geothermics. 2014;52:195–205. https:// doi. org/ 10. 1016/j. geoth ermics. 2014. 01. 004. Zaal et al. Geotherm Energy (2021) 9:12 Page 24 of 26 Bonté D, van Wees J-D, Verweij J. Subsurface temperature of the onshore Netherlands: new temperature dataset and modelling. Neth J Geosci Geol en Mijnbouw. 2012;91(4):491–515. https:// doi. org/ 10. 1017/ s0016 77460 00003 54. Buijze L, Bogert PAJ, Wassing BBT, Orlic B. Nucleation and arrest of dynamic rupture induced by reservoir depletion. J Geophys Res Solid Earth. 2019a;124(4):3620–45. https:// doi. org/ 10. 1029/ 2018j b0169 41. Buijze L, van Bijsterveldt L, Cremer H, Paap B, Veldkamp H, Wassing BB, van Wees J-D, van Yperen GC, ter Heege JH. Review of induced seismicity in geothermal systems worldwide and implications for geothermal systems in the Netherlands. Neth J Geosci. 2019b;98:E13. https:// doi. org/ 10. 1017/ njg. 2019.6. Byerlee J. Friction of rocks. Pure Appl Geophys. 1978;116(4–5):615–26. https:// doi. org/ 10. 1007/ bf008 76528. Candela T, Osinga S, Ampuero J-P, Wassing B, Pluymaekers M, Fokker PA, Wees J-D, Waal HA, Muntendam-Bos AG. Depletion-induced seismicity at the groningen gas field: Coulomb rate-and-state models including differential compaction effect. J Geophys Res Solid Earth. 2019;124(7):7081–104. https:// doi. org/ 10. 1029/ 2018j b0166 70. CBS. Aardgas en elektriciteit, gemiddelde prijzen van energieverbruikers. 2019. https:// opend ata. cbs. nl/ statl ine/#/ CBS/ nl/ datas et/ 81309 NED/ table? froms tatweb. Accessed 02 Nov 2019. Compernolle T, Welkenhuysen K, Petitclerc E, Maes D, Piessens K. The impact of policy measures on profitability and risk in geothermal energy investments. Energy Econ. 2019;84:104524. https:// doi. org/ 10. 1016/j. eneco. 2019. 104524. Daniilidis A, Doddema L, Herber R. Risk assessment of the Groningen geothermal potential: from seismic to reservoir uncertainty using a discrete parameter analysis. Geothermics. 2016;64:271–88. https:// doi. org/ 10. 1016/j. geoth ermics. 2016. 06. 014. Daniilidis A, Alpsoy B, Herber R. Impact of technical and economic uncertainties on the economic performance of a deep geothermal heat system. Renew Energy. 2017;114:805–16. https:// doi. org/ 10. 1016/j. renene. 2017. 07. 090. Daniilidis A, Khait M, Saeid S, Vokov D, Bruhn D. A high performance framework for the optimization of geothermal systems, comparing energy production and economic output. In: Proceedings world geothermal congress. 2020a;491–515. Daniilidis A, Nick HM, Bruhn DF. Interdependencies between physical, design and operational parameters for direct use geothermal heat in faulted hydrothermal reservoirs. Geothermics. 2020b;86:101806. https:// doi. org/ 10. 1016/j. geoth ermics. 2020. 101806. Daniilidis A, Nick HM, Bruhn DF. Interference between geothermal doublets across a fault under subsurface uncertainty; implications for field development and regulation. Geothermics. 2021a;91:102041. https:// doi. org/ 10. 1016/j. geoth ermics. 2021. 102041 Daniilidis A, Saeid S, Doonechaly NG. The fault plane as the main fluid pathway: Geothermal field development options under subsurface and operational uncertainty. Renewable Energy. 2021b;171:927–946. https:// doi. org/ 10. 1016/j. renene. 2021. 02. 148 Davis SD, Frohlich C. Did (or will) fluid injection cause earthquakes?—criteria for a rational assessment. Seismol Res Lett. 1993;64(3–4):207–24. https:// doi. org/ 10. 1785/ gssrl. 64.3- 4. 207. de Boer S, Bourdon E, Butterman G, et al. Evaluatie Algemeen instrumentarium geothermie. 2016. Donselaar M, Groenenberg R, Gilding D. Reservoir geology and geothermal potential of the delft sandstone member in the west Netherlands basin. In: Proceedings world geothermal congress. 2015;2015:1–9. Donselaar M. Reservoir architecture modelling for geothermal energy production—case study of the delft sandstone member, west Netherlands. In: 78th EAGE conference and exhibition 2016. EAGE Publications BV, 2016;1–5. https:// doi. org/ 10. 3997/ 2214- 4609. 20160 0596. Evans KF, Zappone A, Kraft T, Deichmann N, Moia F. A survey of the induced seismic responses to fluid injection in geothermal and CO2 reservoirs in Europe. Geothermics. 2012;41:30–54. https:// doi. org/ 10. 1016/j. geoth ermics. 2011. 08. 002. Eyerer S, Schifflechner C, Hofbauer S, Bauer W, Wieland C, Spliethoff H. Combined heat and power from hydrothermal geothermal resources in Germany: an assessment of the potential. Renew Sustain Energy Rev. 2020;120:109661. https:// doi. org/ 10. 1016/j. geoth ermics. 2011. 08. 002. Gaucher E, Schoenball M, Heidbach O, Zang A, Fokker PA, van Wees J-D, Kohl T. Induced seismicity in geothermal reser- voirs: a review of forecasting approaches. Renew Sustain Energy Rev. 2015;52:1473–90. https:// doi. org/ 10. 1016/j. rser. 2015. 08. 026. Gilding D. Heterogeneity determination of the delft subsurface for heat flow modelling. Master’s thesis, Applied Earth Sciences, Delft University of Technology. 2010. Guglielmetti L, Moscariello A, Meyer M, de Bono CN, Berthoud NA. Geothermie 2020: Exploration and development of geothermal energy in geneva. In: 81st EAGE conference and exhibition 2019 Workshop Programme. European Association of Geoscientists & Engineers. 2019;1–4. https:// doi. org/ 10. 3997/ 2214- 4609. 20190 1962. Hettema M. Analysis of mechanics of fault reactivation in depleting reservoirs. Int J Rock Mech Mining Sci. 2020;129:104290. Ikari MJ, Saffer DM, Marone C. Frictional and hydrologic properties of clay-rich fault gouge. J Geophys Res. 2009. https:// doi. org/ 10. 1029/ 2008j b0060 89. Ikari MJ, Marone C, Saffer DM. On the relation between fault strength and frictional stability. Geology. 2011;39(1):83– 6. https:// doi. org/ 10. 1130/ g31416.1. Jaeger JC, Cook NG, Zimmerman R. Fundamentals of rock mechanics. Hoboken: Wiley; 2009. Jansen JD, Singhal P, Vossepoel FC. Insights from closed-form expressions for injection- and production-induced stresses in displaced faults. J Geophys Res Solid Earth. 2019;124(7):7193–212. https:// doi. org/ 10. 1029/ 2019J B0179 32. Khait M, Voskov D. Adaptive parameterization for solving of thermal/compositional nonlinear flow and transport with buoyancy. SPE J. 2018a;23(02):522–34. https:// doi. org/ 10. 2118/ 182685- pa. Khait M, Voskov D. Operator-based linearization for efficient modeling of geothermal processes. Geothermics. 2018b;74:7–18. https:// doi. org/ 10. 1016/j. geoth ermics. 2018. 01. 012. Khait M. Delft advanced research terra simulator: general purpose reservoir simulator with operator-based lineariza- tion. 2019. http:// resol ver. tudel ft. nl/ uuid: 5f0f9 b80- a7d6- 488d- 9bd2- d68b9 d7b4b 87. Zaal  et al. Geotherm Energy (2021) 9:12 Page 25 of 26 Khait M, Voskov D, Zaydullin R. High performance framework for modelling of complex subsurface flow and trans- port applications. In: Conference proceedings, ECMOR XVII. 2020;2020:1–18. https:// www. earth doc. org/ conte nt/ papers/ 10. 3997/ 2214- 4609. 20203 5188. Knapek E, Dilger G, Hegele H. Deutschlands zukunft mit erdwärme. Tech. rep., Wirtschaftsforum Geothermie, GtV- Bundesverband Geothermie. 2013. Knoblauch TA, Trutnevyte E, Stauffacher M. Siting deep geothermal energy: Acceptance of various risk and benefit scenarios in a Swiss-German cross-national study. Energy Policy. 2019;128:807–16. https:// doi. org/ 10. 1016/j. enpol. 2019. 01. 019. Lepillier B, Daniilidis A, Gholizadeh ND, Bruna P-O, Kummerow J, Bruhn D. A fracture flow permeability and stress dependency simulation applied to multi-reservoirs, multi-production scenarios analysis. Geotherm Energy. 2019;7(1):1–16. https:// doi. org/ 10. 1186/ s40517- 019- 0141-8. Liang X, Xu T, Feng B, Jiang Z. Optimization of heat extraction strategies in fault-controlled hydro-geothermal reser- voirs. Energy. 2018;164:853–70. https:// doi. org/ 10. 1016/j. energy. 2018. 09. 043. Lukawski M, Anderson B, Augustine C, et al. Cost analysis of oil, gas, and geothermal well drilling. J Pet Sci Eng. 2014;118:1–14. Lund J, Toth A. Direct utilization of geothermal energy 2020 worldwide review. In: Proceedings world geothermal congress. 2020;2020:1–39. Markou N, Papanastasiou P. Petroleum geomechanics modelling in the eastern mediterranean basin: analysis and application of fault stress mechanics. Oil Gas Sci Technol Revue d’IFP Energies nouvelles. 2018;73:57. https:// doi. org/ 10. 2516/ ogst/ 20180 34. Mechelse E. The in-situ stress field in the Netherlands: Regional trends, local deviations and an analysis of the stress regimes in the northeast of the Netherlands. Master’s thesis, Applied Earth Sciences, Delft University of Technol- ogy. 2017. Mijnlieff HF. Introduction to the geothermal play and reservoir geology of the Netherlands. Neth J Geosci. 2020. https:// doi. org/ 10. 1017/ njg. 2020.2. Ministerie van Economische Zaken en Klimaat, Staatscourant van het Koninkrijk der Nederlanden, nr 9735. 2019. https:// zoek. offic ieleb ekend makin gen. nl/ stcrt- 2019- 9735. Accessed 05 Nov 2019. Moeck I, Kwiatek G, Zimmermann G. Slip tendency analysis, fault reactivation potential and induced seismicity in a deep geothermal reservoir. J Struct Geol. 2009;31(10):1174–82. https:// doi. org/ 10. 1016/j. jsg. 2009. 06. 012. Moeck IS. Catalog of geothermal play types based on geologic controls. Renew Sustain Energy Rev. 2014;37:867–82. https:// doi. org/ 10. 1016/j. rser. 2014. 05. 032. Morris A, Ferrill DA, Henderson DB. Slip-tendency analysis and fault reactivation. Geology. 1996;24(3):275. https:// doi. org/ 10. 1130/ 0091- 7613% 281996% 29024% 3C0275% 3Asta afr% 3E2.3. co% 3B2 Moscariello A, Guglielmetti L, Omodeo-Salé S, De Haller A, Eruteya OE, Lo H, Clerc N, Makhloufi Y, Do Couto D, De Oliveira GF, Perozzi L, De Oliveira F, Hollmuller P, Quiquerez L, Nawratil De Bono C, Martin F, Meyer M. Heat pro- duction and storage in western switzerland: advances and challenges of intense multidisciplinary geothermal exploration activities, an 8 years progress report. In: Proceedings world geothermal congress. 2020;2020:1–10. O’Sullivan MJ, Pruess K, Lippmann MJ. State of the art of geothermal reservoir simulation. Geothermics. 2001;30(4):395–429. https:// doi. org/ 10. 1016/ S0375- 6505(01) 00005-0. Osundare O, Teodoriu C, Faclone G, et al. Estimation of plugging and abandonment costs based on different EU regulations with application to geothermal wells. In: 43rd workshop on geothermal reservoir engineering. 2018;1–13. Pandey S, Vishal V, Chaudhuri A. Geothermal reservoir modeling in a coupled thermo-hydro-mechanical-chemical approach: a review. Earth-Sci Rev. 2018;185:1157–69. https:// doi. org/ 10. 1016/j. earsc irev. 2018. 09. 004. Parisio F, Vilarrasa V, Wang W, Kolditz O, Nagel T. The risks of long-term re-injection in supercritical geothermal systems. Nat Commun. 2019. https:// doi. org/ 10. 1038/ s41467- 019- 12146-0. Peaceman D. Interpretation of well-block pressures in numerical reservoir simulation with nonsquare grid blocks and anisotropic permeability. Soc Pet Eng J. 1983;23(03):531–43. https:// doi. org/ 10. 2118/ 10528- pa. Pluymaekers M, Kramers L, van Wees J-D, Kronimus A, Nelskamp S, Boxem T, Bonté D. Reservoir characterisation of aqui- fers for direct heat production: Methodology and screening of the potential reservoirs for the Netherlands. Neth J Geosci Geol en Mijnbouw. 2012;91(4):621–36. https:// doi. org/ 10. 1017/ s0016 77460 00004 1x. Rijksdienst voor Ondernemend Nederland. Regeling nationale EZ subsidies - Risico’s dekken voor Aardwarmte. Zevende en achtste openstelling: Handleiding Garantiestelling tegen het risico van misboring. 2019a. Rijksdienst voor Ondernemend Nederland. SDE+ najaar. 2019b. www. rvo. nl/ sites/ defau lt/ files/ 2019/ 10/ Broch ure% 20SDE% 20naj aar% 202019. pdf. Accessed 04 Nov 2019. Rijksoverheid. 2019. Wijzigingen vennootschapsbelasting, belastingplan. 2020. https:// www. rijks overh eid. nl/ onder wer- pen/ belas tingp lan/ belas tingw ijzig ingen- voor- onder nemers/ venno otsch apsbe lasti ng. Accessed 06 Nov 2019. Romera J. iapws.iapws97 module. 2020. https:// iapws. readt hedocs. io/ en/ latest/ iapws. iapws 97. html. Accessed 27 Apr Saeid S, Yang W, Daniilidis A, Khait M, Vokov D, Bruhn D. Lifetime and energy prediction of geothermal systems: Uncer- tainty analysis in highly heterogeneous geothermal reservoirs (Netherlands). In: Proceedings world geothermal congress. 2020;2020:1–8. Salimzadeh S, Paluszny A, Nick HM, Zimmerman RW. A three-dimensional coupled thermo-hydro-mechanical model for deformable fractured geothermal systems. Geothermics. 2018;71:212–24. https:// doi. org/ 10. 1016/j. geoth ermics. 2017. 09. 012. Schoof F, van der Hout M, van Zanten J, et al. Masterplan Aardwarmte in Nederland. Tech. rep., Platform Geothermie, DAGO, Stichting Warmtenetwerk, EBN. 2018. https:// www. ebn. nl/ wp- conte nt/ uploa ds/ 2018/ 05/ 20180 529- Maste rplan- Aardw armte- in- Neder land. pdf. Schoots K, Hammingh P. Klimaat- en Energieverkenning 2019. Tech. rep., Den Haag: Planbureau voor de Leefomgeving. 2019. https:// www. pbl. nl/ sites/ defau lt/ files/ downl oads/ pbl- 2019- klima at- en- energ iever kenni ng- 2019- 3508. pdf. Accessed 04 Oct 2019. Zaal et al. Geotherm Energy (2021) 9:12 Page 26 of 26 Seithel R, Gaucher E, Mueller B, Steiner U, Kohl T. Probability of fault reactivation in the bavarian molasse basin. Geother- mics. 2019;82:81–90. https:// doi. org/ 10. 1016/j. geoth ermics. 2019. 06. 004. Shetty S, Voskov D, Bruhn D. Numerical strategy for uncertainty quantification in low enthalpy geothermal projects. In: 43rd workshop on geothermal reservoir engineering. 2018;1–16. Soltanzadeh H, Hawkes CD. Assessing fault reactivation tendency within and surrounding porous reservoirs during fluid production or injection. Int J Rock Mech Mining Sci. 2009;46(1):1–7. Stauffacher M, Muggli N, Scolobig A, Moser C. sep, Framing deep geothermal energy in mass media: the case of Switzer - land. Technol Forecast Soc Change. 2015;98:60–70. https:// doi. org/ 10. 1016/j. techf ore. 2015. 05. 018. Takahashi M, Mizoguchi K, Kitamura K, Masuda K. Eec ff ts of clay content on the frictional strength and fluid transport property of faults. J Geophys Res. 2007. https:// doi. org/ 10. 1029/ 2006j b0046 78. Tembe S, Lockner DA, Wong T-F. Eec ff t of clay content and mineralogy on frictional sliding behavior of simulated gouges: Binary and ternary mixtures of quartz, illite, and montmorillonite. J Geophys Res. 2010. https:// doi. org/ 10. 1029/ 2009j b0063 83. TNO. Economic model ( Thermogis). 2019. https:// www. therm ogis. nl/ en/ econo mic- model. Accessed 02 Nov 2019. Troiano A, Giuseppe MGD, Troise C, Tramelli A, Natale GD. A coulomb stress model for induced seismicity distribution due to fluid injection and withdrawal in deep boreholes. Geophys J Int. 2013;195(1):504–12. https:// doi. org/ 10. 1093/ gji/ ggt229. Van Adrichem Boogaert H, Kouwe W, Stratigraphic nomenclature of the Netherlands, revision and update by Rijks Geolo- gische Dienst (RGD) and Netherlands oil and gas exploration and production association (NOGEPA). Mededelingen Rijks Geologische Dienst 50. 1993. van Balen R, van Bergen F, de Leeuw C, Pagnier H, Simmelink H, van Wees J, Verweij J. Modelling the hydrocarbon generation and migration in the west Netherlands basin, the Netherlands. Neth J Geosci Geol en Mijnbouw. 2000;79(1):29–44. Van den Bogert P. Impact of various modelling options on the onset of fault slip and fault slip response using 2-dimen- sional finite-element modelling. Restricted, Report No. SR 15. 2015. van den Bosch R, Flipse B, Vorage R. Stappenplan Winning Aardwarmte voor Glastuinbouw. Tech. rep., Kas als energie- bron. 2013. https:// www. kasal sener giebr on. nl/ conte nt/ resea rch/ 14899_ Stapp enplan_ Aardw armte_ 12- 2013. pdf. Accessed 06 Oct 2019. van ’t Spijker H, Ungemach P. Definition of Electrosubmersible Pump (ESP) design and selection workflow. 2016. https:// www. kasal sener giebr on. nl/ conte nt/ user_ upload/ Defin ition_ of_ ESP_ desig n-_ final_ versi on_ 20160 629_- public_ v2. pdf. Accessed 02 Nov 2019. Voskov DV. Operator-based linearization approach for modeling of multiphase multi-component flow in porous media. J Comput Phys. 2017;337:275–88. https:// doi. org/ 10. 1016/j. jcp. 2017. 02. 041. Voskov D. DARTS—delft advanced terra reservoir simulator. 2020. https:// darts. citg. tudel ft. nl. Accessed 09 July 2020. Vrolijk PJ, Urai JL, Kettermann M. Clay smear: review of mechanisms and applications. J Struct Geol. 2016;86:95–152. https:// doi. org/ 10. 1016/j. jsg. 2015. 09. 006. Wang Y, Khait M, Voskov D, et al. Benchmark test and sensitivity analysis for geothermal applications in the Netherlands. In: 44th workshop on geothermal reservoir engineering. 2019;1–11. Wang Y, Voskov D, Khait M, Bruhn D. An efficient numerical simulator for geothermal simulation: a benchmark study. Appl Energy. 2020;264:114693. https:// doi. org/ 10. 1016/j. apene rgy. 2020. 114693. Willems C, Nick H, Goense T, Bruhn D. The impact of reduction of doublet well spacing on the net present value and the life time of fluvial hot sedimentary aquifer doublets. Geothermics. 2017a;68:54–66. https:// doi. org/ 10. 1016/j. geoth ermics. 2017. 02. 008. Willems CJ, Nick HM, Weltje GJ, Bruhn DF. An evaluation of interferences in heat production from low enthalpy geother- mal doublets systems. Energy. 2017b;135:500–12. https:// doi. org/ 10. 1016/j. energy. 2017. 06. 129. Zang A, Oye V, Jousset P, Deichmann N, Gritto R, McGarr A, Majer E, Bruhn D. Analysis of induced seismicity in geothermal reservoirs—an overview. Geothermics. 2014;52:6–21. https:// doi. org/ 10. 1016/j. geoth ermics. 2014. 06. 005. Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geothermal Energy Springer Journals

Economic and fault stability analysis of geothermal field development in direct-use hydrothermal reservoirs

Loading next page...
 
/lp/springer-journals/economic-and-fault-stability-analysis-of-geothermal-field-development-EhifquCgc0
Publisher
Springer Journals
Copyright
Copyright © The Author(s) 2021
eISSN
2195-9706
DOI
10.1186/s40517-021-00193-0
Publisher site
See Article on Publisher Site

Abstract

A.Daniilidis@tudelft.nl Faculty of Civil Engineering The installed capacity of geothermal systems for direct use of heat is increasing world- and Geosciences, section wide. As their number and density is increasing, the their interaction with subsurface Reservoir Engineering, Delft University of Technology, faults becomes more important as they could lead to safety risks from induced seismic- Stevinweg 1, 2628CN Delft, ity. Assessment and management of such risks is essential for the further development The Netherlands and extension of geothermal energy for heating. At the same time, the economic out- put of geothermal systems can be marginal and is hence often supported by subsidy schemes. A combined assessment of fault stability and economic output could help operators to balance economic and safety aspects, but this is currently not common practice. In this study we present a methodology to assess field development plans based on fault stability and Net Present Value (NPV ) using reservoir simulations of a fluvial, heterogeneous sandstone representative of the majority of direct-use Dutch geothermal systems. We find that the highest friction coefficient leading to exceed- ance of the Mohr–Coulomb failure criteria in this sandstone is 0.17; such values could be encountered in clay-rich fault gouges. Similar or lower fault permeability compared to the reservoir results in no changes and an increase respectively of both NPV and fault stability with larger Fault-to-Well Distance (FWD). Fault permeability higher than the reservoir permeability results in a minor increase in NPV with smaller FWD. Our results demonstrate that a combined analysis of thermal, hydraulic, mechanical and economic assessment supports a responsible and viable development of geothermal resources at a large scale. The importance of a high spatial density of supporting stress data will be essential for a better understanding and quantification of economic and fault stability effects of geothermal operations. Keywords: Direct-use geothermal, NPV, Economic, Fault stability, Field development Introduction Installed capacity and generated energy of geothermal systems for direct-use of heat is on the rise (Lund and Toth 2020) as they form a promising energy source that can deliver renewable, low-carbon, baseload energy for industrial and domestic usage. Fur- ther acceleration of installed capacity of these systems is envisioned in several countries (e.g. Knapek et  al. (2013), Schoof et  al. (2018), Guglielmetti et  al. (2019), Moscariello et al. (2020)). A higher density of direct-use geothermal operations can lead to potential © The Author(s) 2021. This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/. Zaal et al. Geotherm Energy (2021) 9:12 Page 2 of 26 interference between systems (Willems et  al. 2017b; Daniilidis et  al. 2021a) and more interaction with structural subsurface elements such as faults (Daniilidis et  al. 2020b). The proper assessment, evaluation and mitigation of safety risks is essential for local communities and the environment. One of these risks is that geothermal operations induce seismic events. Concerns related to the these or other disturbances could delay or even inhibit development of geothermal projects (Gaucher et  al. 2015; Buijze et  al. 2019b; Stauffacher et  al. 2015). The required acceleration and expansion of geothermal operations depend on the assessment and management of their potential safety risks (Buijze et al. 2019b; Knoblauch et al. 2019). This is essential for social acceptance of the increased use of geothermal energy. Additionally, the economic output of direct-use geothermal energy is subject to both subsurface and economic uncertainties (Daniilidis et al. 2017; Compernolle et al. 2019) and often requires subsidies to achieve profitability (Daniilidis et  al. 2017; Eyerer et  al. 2020). Increasing the number of geothermal developments requires a stable yield of prospective projects and low financial risks (Schoof et  al. 2018). Financial risks can be reduced by increasing the probability of success and avoiding low-yield projects. It is therefore important to investigate the profitability of a potential reservoir upfront for better investment decisions. For example, maximising flow rates has shown to be eco - nomically beneficial for operators as they have been demonstrated to lead to increased Net Present Value output even though this means an accelerated cooling of the reservoir (Daniilidis et al. 2020b, 2021a, b). Previous studies have found decreased heat production with a low permeability fault in close proximity and parallel to the geothermal doublet (Liang et  al. 2018). Interac- tion of direct-use geothermal operations with pre-faults in conduction-dominated geo- thermal settings has shown significant effects on energy generation, system lifetime and economic output (Daniilidis et al. 2016, 2020b). The proximity to a sealing fault not only influences these three success factors, it also affects the available reservoir volume for heat exchange (Daniilidis et al. 2016, 2020b). High fault permeability combined with the relative fault block positions can also have positive effects on energy and system life - time when connection to deeper parts of the reservoir is achieved through the fault (Daniilidis et al. 2020b). Moreover, depending on well spacing and fault permeability, the presence and proximity of faults can lead to positive interference between two doublets (Daniilidis et al. 2021a). Due to the low magnitude and small number of earthquakes related to geothermal production in sedimentary settings (Evans et  al. 2012; Buijze et  al. 2019b) where con- duction is the dominant heat transfer mechanism (Moeck 2014) seismic hazards from direct-use geothermal operations have not been addressed previously. However, in sedimentary systems induced seismicity related to gas production has been the focus of several studies (Buijze et al. 2019a; Van den Bogert 2015; Candela et al. 2019). There - fore, the focus on fault stability of pre-existing faults is more relevant for these geologi- cal settings. Conversely, coupled thermo-hydraulic-mechanical models have been used to describe fracture deformation under cold water injection and related stress changes (Salimzadeh et  al. 2018) or to assess the potential of Enhanced Geothermal Systems (EGS) (Lepillier et al. 2019). Such models have been further used to identify Coulomb- stress changes attributed to fluid injection as the main trigger for seismic events in an Zaal  et al. Geotherm Energy (2021) 9:12 Page 3 of 26 EGS setting (Troiano et al. 2013). More recently, the impact of re-injection of cold water into supercritical geothermal systems in terms of fault stability has been analysed (Pari- sio et al. 2019). Nonetheless, for conduction-dominated fields the combined influences of specific reservoir conditions and operational settings on the profitability and on the occurrence of fault instability of pre-existing faults, which might result in induced seis- micity, remains under-explored. This work addresses for the first time a single, unified assessment of fault stability and economic performance. We present a novel methodology that combines existing approaches for analysis of these two aspects. We assess the influence of fault presence and proximity to a direct-use geothermal doublet in sedimentary, conduction-domi- nated geothermal fields on economic output and fault stability for a range of fault prop - erties and operational choices. The Delft Advanced Research Terra Simulator (DARTS) (Khait 2019) is used to perform reservoir simulations of the Delft sandstone, a heteroge- neous fluvial sandstone, one of the main reservoirs for geothermal energy production in the Southern part of The Netherlands (Mijnlieff 2020). The simulation results are used to assess the economic output of a typical doublet based on an updated economic model that includes the current renewable-energy subsidy scheme for The Netherlands. The simulations resolve the pore-pressure change, while the fault stability is assessed based on the Mohr–Coulomb failure criterion. We simulate pre-existing faults positioned at a varying distance from the wells of a geothermal doublet, each fault with a different constant permeability. Operational choices considered include the flow rate and the re- injection temperature of the geothermal doublet. Our final analysis compares the gener - ated Net Present Value (NPV) and identifies the range of friction coefficients for which the Mohr–Coulomb stresses are exceeded, leading to potential fault instability. The out - line of the paper is as follows: “Geological setting” section provides a short introduction to the geological setting of the study area, “Delft Advanced Research Terra Simulator (DARTS)”,  “Fault-stability model” and “Economic model” sections describe the simula- tion methodology, fault stability model and economic model, respectively. “Results” section  provides the results of injection pressure, doublet power, spatial distribution of pressure and temperature, NPV, fault stability and combined NPV and fault stability. “Discussion” section places the findings in the broader context of geothermal develop - ment and “Conclusions” section gives the conclusions. Geological setting The target reservoir of this work is the Delft Sandstone Member (DSSM), which is part of the Nieuwerkerk formation in the West Netherlands Basin (WNB). The WNB is a 60 km wide trans-tensional basin, located in the southwestern parts of the Neth- erlands and contains stratigraphic formations of the Late Jurassic/Early Cretaceous (Pluymaekers et al. 2012; Bonté et al. 2012). The DSSM is layered and consists of fine to coarse gravelly sandstone sequences deposited by a north-west flowing meander - ing river and bands of silt- and claystones deposited as floodplains in a lower coastal plain setting (Van Adrichem  Boogaert and Kouwe 1993; Donselaar et  al. 2015). The deposits in the WNB have large thickness variations, which are the result of exten- sional and compressional faulting. Because of the present faults three main anticlinal and synclinal structures can be classified: the Delft High, the Pijnacker High and the Zaal et al. Geotherm Energy (2021) 9:12 Page 4 of 26 Fig. 1 Typical SW-NE geological cross-section through the WN, after van Balen et al. (2000). The location of the hypothetical reservoir considered in this study is tentatively indicated with a red rectangle. For the exact location of the cross-section, see van Balen et al. (2000) Fig. 2 Schematic of the analysis. Input data used are detailed in Table 1 Vrijenban Syncline which lies in between the two highs (Gilding 2010). The highest part of the Vrijenban Syncline near the Delft High is found at 1900 m true vertical depth in which the top of the DSSM is at 2200 m true vertical depth. The syncline is bounded by two main normal faults: the Delft main boundary fault in the southwest and the Pijnacker main boundary fault in the northwest. Additionally, the area con- tains several inverse normal faults and reverse faults. Movements of the two main faults and the syncline blocks along the fault reactivated and inversed the main faults, which also created new reverse normal faults in the area (Gilding 2010). A typi- cal SW-NE cross-section through the WNB, based on the model of van Balen et  al. (2000), is shown in Fig. 1. The location of the hypothetical reservoir considered in this study is indicated with a red rectangle in the Schieland Sandstone. Methods The present analysis is based on three discrete parts: reservoir simulations using the Delft Advanced Research Terra Simulator (DARTS), economic simulations using an economic model and fault-stability analysis with a model based on the Coulomb fail- ure criterion. An overview of the different parts and interconnections is presented in Fig. 2. Zaal  et al. Geotherm Energy (2021) 9:12 Page 5 of 26 Delft Advanced Research Terra Simulator (DARTS) The Delft Advanced Research Terra Simulator (DARTS) is a Python and C++ based dynamic reservoir simulator, developed at Delft University of Technology (Khait 2019; Voskov 2020) and has demonstrated fast and accurate performance for geothermal applica- tions (Wang et al. 2019, 2020). It is based on the framework of Operator-Based Lineariza- tion (OBL) (Voskov 2017), which makes it fit for the simulation of complex multiphase flow and transport systems (Khait 2019; Khait and Voskov 2018b). For the simulation of geothermal reservoirs three fundamental equations are used, which are the conservation of mass, conservation of energy and Darcy’s law for fluid flow in porous media (O’Sullivan et al. 2001). The equations are described in Voskov (2017), Khait and Voskov (2018b). The conservation of mass with n components and n phases is given c p by Eq.  1, in which the three terms denote respectively the change in mass with time, the mass flux due to flow and a sink- and source term: n  n n p p p φ x ρ s − ∇ · x ρ v + x ρ q˜ = 0, cj j j cj j j cj j j ∂t (1) j=1 j=1 j=1 c = 1, . . . , n , in which φ is the porosity, x the mole fraction of component c in phase j, s the phase cj j saturations, ρ the phase molar density, n the number of phases present, n the number j p c of components, q ˜ the phase rate per unit volume and v the phase velocity. j j The equation for the conservation of energy accounts for the stored energy, convection, conduction and sink and source terms: n  n n p p p φ ρ s U + (1 − φ)U − ∇ · h ρ v +∇ · (κ∇T ) + x h ρ q˜ = 0, j j j r j j j cj j j j ∂t j=1 j=1 j=1 (2) in which U is the phase internal energy, U the internal energy of the rock, h the phase j r j enthalpy, κ the thermal conductivity and ∇T the temperature gradient. Thermal convection is controlled by Darcy’s Law: rj v =− K (∇p − g ∇d) , j = 1, . . . , n , j j j p (3) in which K is the permeability tensor, k the relative permeability, µ the phase viscosity, rj j p the vector of pressures in phase j, g is the gravity term and ∇d is the vector of depths. j j All fluid properties are a function of both pressure and temperature according to the IAPWS97 formulation (Romera 2020). Equations  1,  2, 3 are transformed into Eqs.  4 and  5 of flow and transport for general multi-component fluids by applying finite-volume discretization and a backward Euler approximation in time (Voskov 2017; Khait and Voskov 2018b): Zaal et al. Geotherm Energy (2021) 9:12 Page 6 of 26 n   n    n  n p n+1 p n p p l l l l V φ x jρ s − φ x jρ s − �t x ρ Ŵ �ψ + �t ρ x q = 0, c = 1, ... , n , c j j c j j p cj j c cj j j j=1 j=1 l j=1 j=1 (4) and n   n p p n+1 n V φ ρ s U + (1 − φ)U − φ ρ s U + (1 − φ)U j j j r j j j r j=1 j=1 (5) n  n p p l l l l l l − �t h ρ Ŵ �ψ + Ŵ �T + �t h ρ q = 0, j j j j j j c j=1 j=1 in which V is the control volume, q , which is equal to q ˜ V is a source of phase j, �ψ is j j the difference in pressures between the blocks connected via interface l, T is the dif- ference in temperature between blocks connected via interface l, Ŵ , which is equal to l l Ŵ k rj is the phase transmissibility, Ŵ is the constant geometrical part of transmissibility, l l l l Ŵ , equal to Ŵ κ , equal to φ s γ − κ + κ the thermal transmissibility. r r c j=1 j j Equations 4 and 5 are approximated in time using a fully implicit method, which intro- duces non-linearity to the system (Voskov 2017). Linearisation in DARTS is done using the Operator-Based Linearisation approach (OBL) (Voskov 2017). OBL simplifies the complicated non-linear physics by parameterising the different operators in physical space with a limited number of supporting points, providing a piece-wise representation of the operators (Wang et al. 2019; Khait and Voskov 2018a, b). During simulations the operators are evaluated based on multi-linear interpolations, which reduces computa- tion time and improves the simulation of non-linear systems. The successful application of OBL in thermal multi-component and multi-phase flow simulations has been dem - onstrated in recent publications (Voskov 2017; Khait and Voskov 2018a, b; Wang et  al. 2019, 2020, Daniilidis et al. 2020a, Saeid et al. 2020). The well index is calculated accord - ing to the Peaceman well equations (Peaceman 1983), no skin effect is considered and flow inside the wells only considers buoyancy and does not include frictional pressure and heat transfer. Fault‑stability model The fault-stability model is based on the Mohr–Coulomb theory. Fault stability depends on the regional stress field (stress tensor), the orientation of the fault surface, its rock frictional characteristics (e.g. cohesive strength and friction coefficient) and the ratio between shear- and normal stresses acting on the fault (Morris et  al. 1996). A fault is considered to be stable if it does not reach the Mohr–Coulomb failure criterion, as described in Amonton’s law (Jaeger et al. 2009): τ = C + µ(σ − P), f ft (6) in which τ is the shear strength of the rock (MPa), C the cohesion of the rock (MPa), µ the apparent coefficient of internal friction (−) (Byerlee 1978), σ the total stress nor- ft mal to the failure plane (MPa), P the pore pressure (MPa) and σ − P the effective shear ft stress σ . ft Zaal  et al. Geotherm Energy (2021) 9:12 Page 7 of 26 According to Amonton’s law the stability or failure of a fault is determined by the ratio of shear to normal stresses acting on a surface and this also defines the slip tendency (Morris et  al. 1996). Failure is likely to occur when the resolved shear stress (τ ) e quals or exceeds the frictional sliding resistance, also called the failure shear stress (τ ) ( Jae- ger et  al. 2009). Hence the tendency for a pre-existing fault to fail is given by Eq.  7, in which the failure stress is proportional to the normal stress acting on the surface (C = 0, assuming no cohesion develops due to healing): τ ≥ τ = µ(σ ). f (7) ft Equation 6 gives the failure criterion describing the shear strength of the fault. In com- bination with the Mohr circle, the failure criterion determines the fault stability. We assume that failure occurs if the Mohr circle touches or intersects the line that describes the relation between the normal stress and the effective shear stress for a given friction coefficient, that is, where the effective shear stress (Eq.  7) equals or exceeds the failure shear stress of the fault. The maximum pore pressure values are computed based on the corresponding minimum required friction coefficient to reach the failure point and is visualised in Fig.  3a. In our setup, the pore pressure values used are the ones from the cells next to the fault plane. In this analysis we do not distinguish between seismic or aseismic slip, but only the exceedance of the failure criterion at any point of the fault and identification of the corresponding pore pressure change psimilar to for example Jansen et  al. (2019); Hettema (2020); Moeck et  al. (2009); Parisio et  al. (2019)] . In our study, poroelastic effects are not included. In the reservoir simulations (see “Reservoir model, sensitivity study and input parame- ters” section) we assume the presence of a vertical fault in a normal faulting regime. There - fore, failure in this case does not take place on the σ –σ plane, but rather on the σ –σ 1 3 2 3 Fig. 3 Mohr circles for a friction coefficient of 0.5 (black dashed line) showing the initial stress state (black circle) and the effective stress state (orange circle) with the pore pressure increase required for a critically dipping fault along the σ –σ plane (a). Pore pressure increase required for a vertical fault to reach failure 1 3 along the σ –σ plane (b) and effect of thermal stresses for a laterally confined reservoir where the thermal 2 3 stress acts on the horizontal components of the stress tensor with a T of 49.2 K and rock properties according to Table 3 (c). Pore pressure increase required to reach failure along the σ –σ plane of a thermally 2 3 stressed reservoir (d) Zaal et al. Geotherm Energy (2021) 9:12 Page 8 of 26 plane that defines the angle of the vertical fault plane with respect to the horizontal compo - nents of the stress tensor (Fig. 3b) resulting in strike-slip failure Markou and Papanastasiou (2018). Thermal stresses can cause significant changes in the stress tensor, especially around the injector well where the largest temperature differences occur. In this study we consider thermal stresses for a laterally constrained reservoir, therefore only acting on the horizontal components of the stress tensor (Fig. 3c). The thermal stress is evaluated based on (Jaeger et al. 2009): α E�T � = , (8) (1 − ν) in which � is the thermal stress change (MPa), α is the rock linear thermal expansion σ t coefficient, E is the rock Young’s modulus, T is the temperature difference (K) and ν is the rock Poisson’s ratio (−). The combined effect of thermal stresses and pore pressure is visualised in Fig. 3d. Based on the failure criterion, the critical fault strike angles can be evaluated, similar to Soltanzadeh and Hawkes (2009). Economic model The economic model investigates the development of a geothermal project and is imple - mented as outlined in previous research (Daniilidis et  al. 2017, 2020b). Multiple stages are included in the development among which the research phase, exploratory phase, development- and production phases and the abandonment phase. All phases contrib- ute to the total costs of a geothermal project, depending on multiple factors such as the geological uncertainties, government policies and required equipment (Daniilidis et  al. 2017; Schoof et al. 2018). The costs of the development phase include the drilling costs of two wells, which form the doublet and are a function of depth. For each well the drill- ing costs are computed as (TNO 2019): W = 375, 000 + 1150MD + 0.3MD , (9) cost in which MD is the measured depth (m) and the doublet costs are in Euro (€). Experi- ence from the first drilled well is estimated to result in decreased drilling costs of the second well by means of a learning curve reducing the costs by 7.5% (Lukawski et  al. 2014). The Net Present Value is based on the project expenses and revenues according to: NPV = � , (10) (1 + r) where C is the net cashflow generated at time t, r the monthly discount rate (%), t the elapsed number of time periods since the project start (t  =  0 is the time at which the wells are drilled). The net cashflow ( C ) in € is computed according to: C = I + I − CapEx − OpEx − OpEx − AbEx, t heat subsidy annual pumping (11) where I is the income from the produced heat, I is the subsidy income, CapEx heat subsidy are the capital expenditures, OpEx the annual operational expenses, OpEx annual pumping Zaal  et al. Geotherm Energy (2021) 9:12 Page 9 of 26 the pumping expenses and AbEx the abandonment cost. The CapEx consist of the total investment costs of the project and include the exploration, drilling costs (W ), con- cost struction and unforeseen construction costs associated with preparation, drilling and installation of the equipment and drill site. Equipment includes the injector and pro- ducer pumps, the heat exchanger and the screens and filters installed in the well. The OpEx consist of all costs made during production and include an annual component as a percentage of CapEx and the expenditures for running the injector and producer pumps, according to the electricity price. The annual OpEx include the maintenance- and work - over costs, rent of facilities, staff payment and insurance costs and amounts to 4% of the CapEx per anum for the annual costs for inspection, maintenance, workovers and costs for staffing (de Boer et al. 2016) and 1% of the CapEx per anum for heat network mainte - nance to ensure the heat is transported to the district heating system; the district heating itself is considered as pre-existing and its maintenance is not included in this analysis. These percentages remain constant during the project’s lifetime. The income not gener - ated during maintenance and workovers are added to the fixed OpEx in the model to take into account the effects of the downtime (Daniilidis et al. 2017). For the insurance the model incorporates a premium paid by the operator, which is a one-time payment of 7% of the maximum compensation available for a typical geothermal reservoir at 2000 m depth (Rijksdienst voor Ondernemend Nederland 2019a). The Weighted Average Cost of Capital (WACC) is used as the discount rate for the NPV. This is a weighted average of the costs of capital and debt. Lenders and equity holders expect to receive a certain return on the capital or funds they provided. This expected return is given by the WACC and is calculated as : E D WACC = Re + Rd(1 − Tc), (12) V V in which Re is the cost of equity, Rd is the cost of debt, V the total market value of the project’s financing (E + D) and Tc the corporate tax rate, E the market value of equity (0.3) and D the market value of debt (0.7). In the Netherlands, projects that may be sub- sidized with the SDE+ regulation are financed in a 70/30 debt over equity ratio. This all results in a WACC value of 5.4%. Dutch subsidy The project revenues include the income generated from the delivered heat and from the Sustainable Duurzame Energieproductie (SDE+) subsidy. In the Netherlands the heat price is determined as 90% of the Dutch market price of gas (TTF) (Ministerie van Econ- omische Zaken en Klimaat 2019) and this heat price determines the generated income and the subsidy amount. The TTF price is based on the High Heating Value (HHV) of gas. The incorporated gas price in the model is based on statistically historical data (Schoots and Hammingh 2019). The SDE+ subsidy provides financial compensation for up to 6000 full-load production hours over a maximum period of 15 years (Rijksdienst voor Ondernemend Nederland 2019b) for projects with depths between 500 and 4000 m. The subsidy payouts are based on the yearly amount of power required for produc - tion, which is fixed, the market price of thermal energy and the full-load production Zaal et al. Geotherm Energy (2021) 9:12 Page 10 of 26 hours. Based on this, the model takes into account a final subsidy amount of €0.033 per kWh and this is implemented on a yearly basis for 15 years. Energy evaluation The produced energy is computed according to (Daniilidis et al. 2017): E = q�TC ρ �t, produced f f (13) where q is the volumetric flow rate (m /s), T the temperature difference between the produced and injected water (K), C the specific heat capacity of the fluid at constant pressure (J/(kgK)), ρ the density fluid (kg/m ) and t the time interval (h). Heat losses at surface and in the heat exchangers are not taken into account in this energy evaluation. The required pumping energy is determined as the product of the flowrate and the pressure difference between the wells, integrated over time (van’t Spijker and Ungemach 2016): q�P E = �t, pump (14) where q is the volumetric flow rate ( m /s ), P the pressure difference between injector and producer wells (Pa) and η the pump efficiency taken here as 55%. The governing equations are implemented in the model as presented by (Daniilidis et  al. 2017, 2020b) and use production fluid temperature and well pressure from the DARTS simulations. Reservoir model, sensitivity study and input parameters The temperature, pressure, NPV and fault stability are analysed using a 3D heteroge - neous reservoir model for the Delft area (Fig.  4) comprised of approx. 100 k elements. Fig. 4 Map view of the 3D reservoir model of the heterogeneous Delft Sandstone. The injection well is indicated in blue and the production well indicated red, well spacing is 990 m. Each cell is 30 m by 30 m by 2.5 m (x, y, z). The vertical fault is indicated by the purple line (shown here with a permeability of 0.75 mD), has a width of 30 m and extends throughout the domain along the x-axis. The different fault permeability values used in the analysis are listed in Table 3 Zaal  et al. Geotherm Energy (2021) 9:12 Page 11 of 26 The geological setting has been described in “ Geological setting” section and references therein. The reservoir model represents a heterogeneous fluvial system with varying net-to-gross (N/G). It includes different fluvial deposits that have been generated by the process-based model Flumy and are classified as reservoir (points bars, sand plugs and channel lags) and non-reservoir (crevasse splays, overbank alluvium and mud plugs). Using petrophysical well data, the following permeability-porosity relationship is used, after Willems et al. (2017b): 29.507φ K = 0.0633e . (15) The reservoir model, as well as the porosity and permeability distribution have been described extensively in Willems et al. (2017b) and have also been used in previous stud- ies (Wang et  al. 2019; Shetty et  al. 2018; Daniilidis et  al. 2020a). The top and bottom boundary conditions act as infinitely large domains with a fixed temperature and pro - vide conductive thermal recharge to the reservoir (Daniilidis et al. 2020a). We perform a sensitivity study consisting of a number of cases, each with a different choice of reservoir properties and/or operational settings. The input parameters that are the same for each of the cases simulated with the reservoir model are given in Table 1. The well locations remain constant throughout the analysis and the assumed vertical fault is positioned with respect to the wells, changing the Fault-to-Well Distance (FWD); the fault strike is parallel to the plane defined by the wells and therefore both wells have the same distance to the fault. Using a vertical fault is considered as an acceptable approximation in many fault stability studies Jansen et  al. (2019), Van  den Bogert (2015) for reservoir faults in normal faulting regimes where faults are often steeply dipping. Table 1 Input data for the reservoir model Parameter Value Unit nx × ny × nz 60 × 40 × 42 – dx × dy × dz 30 × 30 × 2.5 m Permeability range 1–3400 mD Porosity range 5–30 % Location producer 450,600 m (X,Y) Location injector 1440,600 m (X,Y) Fault permeability See Table 2 mD Fault thickness 30 m −15 2 Overburden permeability 9.9 × 10 (0.001) m (mD) Top reservoir depth 2200 m Model reservoir thickness 105 m Fault-wells distance see Table 2 m Uniform initial pressure 20.1 MPa Uniform initial temperature 345.75 K Flowrate prod/inj 4200 m /day Injection temperature see Table 2 K Simulation time 30 × 365 days Thermal conductivity rock 100 kJ/m/day/K Heat capacity rock 2200 kJ/m /K Zaal et al. Geotherm Energy (2021) 9:12 Page 12 of 26 Table 2 Varying parameters of the reservoir model. The lowest fault permeability of 0.75 mD is referred to as a sealing fault and the highest fault permeability of 750 mD as transmissive fault in “Results” section Parameter Value Unit Flowrate 150, 225, 300, 375 m /h Fault permeability 0.75, 7.5, 75, 750 mD Injection temperature 303.15, 308.15, 313.15 K Fault-to-Well Distance (FWD) 60, 90, 180, 360 m Table 3 Input data for the fault model Parameter Value Unit σ 45 MPa σ 34 MPa σ 30 MPa Young’s modulus 16 GPa Poisson’s ratio 0.25 – Cohesion coefficient 0 MPa Inital pore pressure 20.1 MPa The stress tensor values are based on the work of Mechelse (2017) Table 4 Economic model input Parameter Value Unit Occurrence Source Production time 30 Years – – No. of doublets 1 – – – No. of wells 2 – – – Exploration costs 460,000 € Once van den Bosch et al. (2013) ESP costs 800,000 € Every 5 years van’t Spijker and Ungemach (2016) Power facility 1,500,000 € Once van den Bosch et al. (2013) Drill location 300,000 € Once van den Bosch et al. (2013) Unforeseen costs 5 % of construction costs Once Daniilidis et al. (2017) ESP downtime 15 days Every 5 years Daniilidis et al. (2017) ESP efficiency 55 % – van ’t Spijker and Ungemach (2016) Electricity price 0.071 €/kWh – CBS (2019) Gas price 0.016 €/kWh – - Fixed OpEx 5 % of CapEx Annually de Boer et al. (2016) Drill insurance 7 % of maximum com- Once Rijksdienst voor Ondernemend pensation available Nederland (2019a) $/€conversion 0.90 – – - AbEx single well 1,275,500 € Once Osundare et al. (2018) Market value of equity 0.3 – – Rijksdienst voor Ondernemend Nederland (2019b) Market value of debt 0.7 – – Rijksdienst voor Ondernemend Nederland (2019b) Cost of equity 14.5 % – Rijksdienst voor Ondernemend Nederland (2019b) Cost of debt 2 % – Rijksdienst voor Ondernemend Nederland (2019b) Corporate tax 25 % – Rijksoverheid (2019) Zaal  et al. Geotherm Energy (2021) 9:12 Page 13 of 26 For each of the cases in the sensitivity study, a simulation is performed for a combina- tion of parameters for both reservoir conditions and operation options. These param - eters are provided in Table 2. The input parameters of the fault stability model are given in Table 3. The inputs of the economic model are all based on Dutch financial regulations. The input used in the model is given in Table 4. Results In this section, we present the results of the sensitivity study as described in “Reservoir model, sensitivity study and input parameters” section. This study is performed with the models described in “Delft Advanced Research Terra Simulator (DARTS)”, “Fault-stabil- ity model” and “Economic model” sections following the workflow as presented in Fig. 2. In the following, we will subsequently address the resulting reservoir pressure and tem- perature as well as the NPV and fault-stability analysis for the reservoir and operational conditions under consideration. This section concludes with the combined results of the NPV and the fault stability. Injection pressure Injection pressure increases over time with most significant changes taking place within the first 3 years as pressure dissipates inside the reservoir and stabilises after about 20 years of production (Fig.  5). The injection rate defines the level at which the plateau occurs, with higher rates resulting in increased injection pressure. The distance to the fault influences the final pressure reached; with an FWD of 60 m away from the wells an increase of up to 2 MPa compared to a distance of 390 m for the lowest fault perme- ability. With an FWD of 390 m the fault permeability has no effect on the maximum injection pressure encountered. However, when the fault is closer to the wells the lower fault permeability of 0.75 mD results in higher injection pressure. This effect is less pro - nounced for a fault permeability of 7.5 mD and is minimal for the fault permeability of 75 mD; this indicates that a fault permeability of 75 mD is close to the average reservoir permeability and therefore the distance of the fault is not significant as the fluid passing the fault behaves hydraulically similar to the fluid in the reservoir layers. A fault perme - ability of 750 mD results in a decrease of injection pressure when the fault is positioned closer to the wells, since the fault provides less resistance to flow than the reservoir layers. Doublet power Doublet power shows a slight decrease over time due to production temperature decrease (Fig. 6). The power starts decreasing after 15 years for the low rate of 150 m /h with minor reductions in the order of 0.3 MW and after circa 7 years for the high rate of 375 m /h with reductions of up to 3 MW. Low fault permeability of 0.75 mD leads to earlier power drop with a FWD of 60 m. As the FWD is increasing the power loss over time is decreasing. The discrepancy caused by the FWD is decreasing as the rates decrease and the fault permeability increases. With a fault permeability of 75 mD the discrepancy caused by the FWD is minimised, similar to the injection pressure (Fig. 5). A fault permeability of 750 mD results in an earlier power drop at circa 5 years for the Zaal et al. Geotherm Energy (2021) 9:12 Page 14 of 26 Fig. 5 Injection pressure over time for fault permeability of 0.75 mD, 7.5 mD, 75 mD and 750 mD, FWD of 60 m, 90 m, 180 m and 360 m and injection temperature of 303.15 K FWD of 60 m. With the increase of FWD the onset of the power drop is delayed and is also less pronounced. Even though the power drop decreases as the FWD increases for 3 3 both high (375 m /h) and low (m /h) rates, the mechanism is different depending on the fault permeability (see also Fig. 8). For low fault permeability (0.75 mD) the volume available to to cold water plume is reduced and reaches the producer well earlier. Con- trary to this, for high fault permeability (750 mD) the fault acts as conduit connecting the injector to the producer well. Spatial distribution of temperature and pressure changes The pressure distribution inside the reservoir is similar for a case with a transmis - sive fault and a case with a sealing fault with an FWD of 390 m (Fig.  7a, b). For both cases, the maximum and minimum pressure values are observed at the location of the injector and producer wells respectively. At an FWD of 60 m, a sealing fault limits the available volume over which pressure can dissipate, resulting in pressure build- up around the injector and along the fault plane close to the injector (Fig.  7c). This creates a large pressure differential across the fault plane. Similarly, a lower pressure is observed around the producer and along the fault close to the producer. With the fault positioned 60 m away from the wells, the transmissive fault allows for the dis- sipation of pressure due to its higher permeability relative to the sealing-fault case (Fig. 7d). This implies that the pressure difference across the fault (north block) is less pronounced in this case compared to a sealing fault at the same distance, and so is the pressure difference between the areas of the injector and producer wells. This leads to relatively lower pressure values around the injector and a weaker pressure contrast between the injector and the fault, compared to case b with the a transmissive fault at an FWD of 390 m. The likely explanation for this is that the fault provides a fluid pathway with less hydraulic resistance and the pressure dissipates across the fault. Figure  8 shows the temperature distribution inside the reservoir in the presence of a transmissive and a sealing fault with an FWD of 60 m and 390 m. A sealing fault at an FWD of 390 m restricts the extent of the cold plume on the north side of the fault (Fig.  8a) compared to a transmissive fault (Fig.  8b). In the latter case the increased fault permeability allows the cold plume to permeate the fault, reducing the temper- ature on both sides of the fault. As such, it allows the cold plume to reach further towards the north of the reservoir. At an FWD of 60 m, a sealing fault allows mostly conductive heat transfer towards the north fault block while convective heat transfer is limited (Fig.  8c). This results in a reduced extension of the cold plume to the north and a restriction of the plume Zaal  et al. Geotherm Energy (2021) 9:12 Page 15 of 26 Fig. 6 Doublet power over time for fault permeability of 0.75 mD, 7.5 mD, 75 mD and 750 mD, FWD of 60 m, 90 m, 180 m and 360 m and injection temperature of 303.15 K Fig. 7 Map view at a depth of 2177.5 m showing the pressure distribution in the DSSM reservoir after 30 years of production with an FWD of 390 m (a, b) and 60 m (c, d) for a sealing (a, c) and a transmissive fault (b, d). Data shown use an injection temperature of 303.15 K and a flowrate of 9000 m /day . Injector and producer wells are marked with a blue and red point respectively and the fault position is indicated in black around the injector where the higher pressure gradients occur. In this location, the cold plume initiates and favours conductive transfer. A transmissive fault at an FWD of 60 m (Fig. 8d) diverts a large part of the flow and results in an extension of the cold plume along the fault direction and towards the producer, leading to earlier thermal breakthrough. In this case, the connection between the two fault blocks is strong and resembles that of the FWD of 390 m, albeit slightly restricted due to the diversion through the fault plane. The vertical reservoir profile along the injector for an FWD of 390 m reveals only minor pressure differences between the case with a sealing fault (Fig.  9a) and the case with a transmissive fault (Fig. 9c) as illustrated in Fig. 9b, d respectively. The sealing fault with an FWD of 390 m results in slightly higher pressure closer to the fault plane. At Zaal et al. Geotherm Energy (2021) 9:12 Page 16 of 26 Fig. 8 Map view at a depth of 2177.5 m showing the temperature distribution in the DSSM reservoir after 30 years of production with an FWD of 390 m (a, b) and 60 m (c, d) for a sealing (a, c) and a transmissive fault (b, d). Data shown use an injection temperature of 303.15 K and a flowrate of 9000 m /day. Injector and producer wells are marked with a blue and red point respectively and the fault position is indicated in black. The notation on the bottom right shows the volume percentage that exhibits a �T > 1.5 K after 30 years of simulation compared to the initial reservoir temperature. The percentage is calculated throughout the whole reservoir volume and it is not limited to the horizontal cross-section shown here an FWD of 60 m the sealing fault significantly increases the pressure between the well and the fault but also around the injector away from the fault, towards the south of the domain (Fig.  9f ) across the full length of the well. The higher pressure in the south in this case compared to the case of the transmissive fault extends all the way to the edge of the domain. The case with a permeable fault at an FWD of 60 m exhibits lower pressure values around the injector, even compared to the case where the fault is located further away (Fig.  9h). The relatively higher permeability of the fault in this case offers an easy pathway to the fluid and reduces pressure buildup. NPV The generated NPV after 30 years of production for the cases with the sealing fault (0.75 mD) shows minor variations at the lowest flow rate (150 m /h) for each injection tem- perature depending on the fault distance to the wells in these cases (Fig.  10). We find that for the lowest flow rate of 150 m /h an injection temperature of 308.15 K or lower is required to achieve positive NPV values. With higher rates, the NPV becomes more sensitive to changes in FWD. For the lowest injection temperature of 303.15 K, increas- ing the rate from 150 to 375 m /h (a factor 2.5) increases the NPV from 5 to 22.5 M€ Zaal  et al. Geotherm Energy (2021) 9:12 Page 17 of 26 (a factor 4.5). In the case of the highest rate of 375 m /h the NPV for an FWD of 60 m is equivalent to maintaining an FWD of 390 m and selecting an injection temperature that is higher by 5 K (308.15 K instead of 303.15 K). The NPV differences caused by the fault proximity are smaller for cases with a higher fault permeability of 7.5 mD across all injection temperatures and flow rates investigated. At a fault permeability of 75 mD the fault distance does not seem to affect the generated NPV and differences in flow rate and injection temperature have a relatively stronger effect. Cases with a higher fault perme - ability to 750 mD show small NPV benefits as the fault distance to the wells is larger. The discrepancy between the cases with a sealing (0.75 mD) and a transmissive fault (750 mD) with lower FWD can be attributed to the lower pumping costs due to lower pres- sure (see also Fig.  9). Possibly, the relatively larger reservoir volume available for heat exchange in the high fault-permeability cases also contributes to a higher NPV, since cold water is averted away from the producer well through the fault and as such pen- etrates a larger volume—64% of volume affected compared to 56% for the transmissive and sealing case respectively at an FWD of 60 m (see also Fig. 8). Fault stability The pore pressure increase required to reach failure shows a similar trend for the ini - tial stress state and the thermal stress states of the injection temperatures considered (Fig.  11a). The σ and σ are equally affected by the thermal stress and therefore the 2 3 shapes of the pore pressure curves remain similar, with only the values of the pore pres- sure leading to failure is increasing. The critical angle between the maximum horizon - tal stress S (which is σ ) and the fault strike shows that the angle range of 25° to H 2 max 43° is critical depending on friction coefficient, with higher friction coefficients leading to smaller angles (Fig.  11b). The thermal stresses slightly shift the critical angle but the effect is minor. The maximum encountered pore pressure derived from the reservoir simulations is higher for cases with a small Fault-to-Well Distance and higher flow rates (Fig.  12). Lower injection temperature further increases the maximum pore pressure with up to an additional 1.5 MPa for the highest flow rate and shortest distance between the sealing fault and the wells. For the highest rates, the effect of the increase in injection tempera - ture by 10 K is almost equivalent to the effect of a decrease in flow rate of 75 m /h. For an injection temperature of 303.15 K and assuming a friction coefficient of 0.15 for the sealing fault, the Mohr–Coulomb failure criterion for pore pressure is exceeded with a 3 3 60 m and 90 m FWD with a rate of 375 m /h and with a 60 m FWD for the 300 m /h rate. A higher injection temperature of 313.15 K would require an even lower friction coefficient of 0.12 for the sealing fault to fail at similar distances and rates. It should be noted that experimental results show a pure quartz gouge to have a friction coefficient between 0.66 and 0.75 (Tembe et al. 2010) while friction coefficient values below 0.3 are related to clay rich minerals (Ikari et  al. 2009). For higher fault permeability cases, the failure criterion is only exceeded when the fault is closer to the wells (FWD below 100 m) assuming the highest rates and a friction coefficient of 0.15 and 0.12 for the injection temperature of 303.15 K and 313.15 K respectively. For a fault permeability of 750 mD friction coefficients below 0.15 and 0.12 are required to reach failure for 303 K and 313 K accordingly. Zaal et al. Geotherm Energy (2021) 9:12 Page 18 of 26 Fig. 9 Vertical cross-section of the reservoir at the injection well location showing the permeability (left column) and corresponding pressure distribution (right column) in the DSSM reservoir. A sealing fault at an FWD of 390 m (a) and 60 m (c) from the injector well and a transmissive fault at an FWD of 390 m (b) and 60 m (d) distance from the injector well are shown. The injector well is marked by a blue line. Data shown for an 303.15 K injection temperature and a flowrate of 9000 m /day Combined NPV and fault stability analysis Combining the NPV results and the results of the maximum pore pressure increase for the highest considered injection temperature cases (313.15 k) shows that with a seal- ing fault (0.75 mD) a lower FWD results both in a lower NPV and a higher maximum pore pressure difference (Fig.  13). The trend is similar for all rates but the differences are larger for cases with larger flow rates. The higher rates and shorter distances lead to an exceedance of the failure criterion with a friction coefficient of 0.17 and just below 0.14 for an injection temperature of 303.15 K and 313.15 K respectively. When fault perme- ability is 75 mD, NPV is no longer affected by the distance to fault and only the rates affect the NPV and maximum pore pressure difference. The trend changes with a trans - missive fault (750 mD) and we observe higher NPV values and maximum pore pressure differences for lower FWDs. The NPV changes, however are minor. Lower injection temperature of 303.15 K leads to higher NPV values and also higher maximum pore pressure differences. Discussion Our results support previous findings of increased NPV for higher flow rates and lower temperatures (Daniilidis et  al. 2016, 2020b, 2021a; Willems et  al. 2017a). For the flow rate ranges considered here and assuming current heat and electricity prices, the addi- tional revenues from high flow rates outweigh the costs of increased pumping require - ments. The NPV results also reflect that either option of increasing the rates or lowering the injection temperature can be used to improve the generated economic value when all other parameters remain the same (Daniilidis et al. 2016). In our simulations, keeping the FWD large is beneficial for both NPV and fault stability when fault permeability is higher than most reservoir permeability values (750 mD). Nonetheless, higher rates have a stronger positive effect on the NPV compared to the effects of a larger FWD. A smaller FWD results in a lower NPV and therefore there is no financial incentive for operators Zaal  et al. Geotherm Energy (2021) 9:12 Page 19 of 26 Fig. 10 Comparison of the generated NPV for the range of injection temperature, distance between the fault and the wells, flowrates and fault permeability considered to knowingly position doublets closer to a fault. If faults are of sub-seismic resolution and their existence or location is not known, smaller FWDs might of course still occur. However, if fault permeability is higher than the assumed reservoir permeability, only minor benefits can be achieved by reducing the FWD. Previous studies have also sug - gested that effect of the FWD to system lifetime and NPV could be optimised if the hydraulic behaviour of both the fault and reservoir is well characterised (Daniilidis et al. 2021a). A combined assessment of thermal, hydraulic, mechanical and economic aspects can be foreseen as a future research path for direct-use geothermal developments. With respect to fault stability we identify a friction coefficient of 0.17 and just below 0.14 for an injection temperature of 303.15 K and 313.15 K respectively as the high- est values at which the pore-pressure change from operations results in exceeding the Mohr–Coulomb failure criterion. The exceedance only occurs when the highest flow rate (375 m /h), short FWD (60 m) and lowest fault permeability (0.75 mD) are used. Temperature effects prove important even without the inclusion of poroelasticity as the temperature effect on fluid density and viscosity alone is sufficient to differentiate the interaction with the Mohr–Coulomb failure criterion. Additionally, compared to the ini- tial stress state, the thermally induced stress further reduces the distance to failure. A difference of up to 3 MPa for the same friction coefficient can be observed in the pore pressure increase required to reach failure for the temperature range considered (303.15 K to 313.15 K); this difference however is minimized at lower friction coefficients. The pore pressure increase itself depends strongly on the fault permeability. Additionally, our analysis does not examine whether the increased pore pressure would result in fault slip as the Mohr–Coulomb criterion can also be exceeded aseismically. Therefore, our evaluation of fault stability serves as a worst-case estimate within the methodological choices presented in “Fault-stability model” section and the initial stress values consid- ered. A determination of the initiation point of seismic slip that includes aseismic effects would be needed to assess actual seismic risk. It should be noted that even though it was not the focus of this work, the high injection rates of 375 m /h result in pressure levels in excess of σ (30 MPa), which could lead to shear failure of the reservoir rock close to the injector well. Experimental research has identified friction coefficients below 0.35 in faults that are rich in clay minerals (Ikari et  al. 2009). These faults exhibit only stable slip unless sig - nificant mineralogy changes are incurred (Ikari et  al. 2011). Other experimental stud - ies report friction coefficients below 0.3 for clay content above 70% (Tembe et al. 2010; Takahashi et  al. 2007). Pure quartz gouge has been experimentally estimated to have Zaal et al. Geotherm Energy (2021) 9:12 Page 20 of 26 Fig. 11 Pore pressure required to reach failure along the σ –σ plane for the initial stress state and the 2 3 thermal stress states considered for a range of friction coefficients (a). The critical fault strike angle with respect to S for the different stress states for the range of friction coefficients (b) max a friction coefficient between 0.66 and 0.75 depending on the grain size (Tembe et  al. 2010). While the target reservoir used in this study is a sandstone, the Delft Sandstone is known to exhibit large ranges of Net-to-Gross and hence also contains significant amounts of shale and therefore clay minerals (Donselaar et al. 2015; Donselaar 2016). As a result, the presence of clay minerals inside the fault gouge cannot be ruled out. Moreo- ver, fluid circulation in close proximity to the fault could initiate or intensify migration of fines inside the fault plane and aid the onset of clay smearing (Vrolijk et  al. 2016). Nonetheless, since clay-smearing potential is proportional to fault throw and inversely proportional to fault-zone thickness (Vrolijk et al. 2016), faults exhibiting clay smearing should be easier to characterise and considered in the design of well trajectories. The risk of induced seismicity increases with the injected volume but also with the proximity of the injection point to the fault (Zang et al. 2014; Davis and Frohlich 1993). Similarly, in our study we find that the MC failure criterion is exceeded earlier with higher flow rates, smaller FWD and lower fault permeability. Stress induced by fluid injection can alter the stress distribution around faults (Altmann et  al. 2014). While a low-permeability fault results in sharper pressure (dP) and temperature differences (dT), a high permeability fault could lead to a much larger volume being perturbed by dP and dT and this could further increase the potential for seismicity. Nonetheless, induced seismicity also has a temporal component and cessation of operations does not mean an instant reduction in seismic risk (Gaucher et al. 2015). Fault reactivation related to fluid injection or extraction is primarily a result of changes in pore pressure, but poroelas- tic and chemical effects can also play a role and should be included for a more detailed understanding. However, by considering only vertical faults, the poroelastic effects are minimized (Hettema 2020). Fault proximity is relevant for the expected NPV of geothermal operations, and thus an important aspect for operators to consider, even if the regulation framework does not impose any restrictions on placement of the wells relative to faults or does not require the monitoring and study of mechanical effects. As such, this economical aspect pro - vides an additional argument for the characterisation of faults present in the area of interest. Establishing the friction coefficient for fault surfaces and fault transmissivity is crucial to enable a more robust assessment of seismic risk from geothermal operations. For most geothermal regions, a more detailed spatial analysis of the stress field mag - nitude and orientation, at a resolution higher than presently available, is crucial. This should be complemented with fault-gouge friction coefficient characterisation in the lab. Zaal  et al. Geotherm Energy (2021) 9:12 Page 21 of 26 Recent studies in the Molasse basin in Germany have also pointed towards the need for further stress-field characterisation (Seithel et al. 2019). Also in The Netherlands, obser - vational efforts in support of a stress-field characterisation should be intensified, as the current stress-data density impairs a more detailed assessment. The methodology presented in this paper could also be adapted and used for EGS or fractured hydrothermal systems. In these cases, a different reservoir characterization and model would need to be developed [similar to e.g. Lepillier et  al. (2019)], the con- version of the heat to electricity would need to be modelled and an economic model for the combined electricity- and heat production would need to be developed. Moreover, poroelastic effects would need to be included as their impact is much stronger in the geological settings where geothermal systems are used for electricity production. Even if more data on both the friction coefficient and stress orientation and magni - tude are available, the evaluation of all potential development options and stress condi- tions remains a major undertaking requiring considerable simulation efforts. Additional aspects such as the well positioning within the considered reservoir permeability field could further increase complexity of the required studies. Moreover, fault stability could also be assessed for faults with varying dip values in relation to both the wells and the stress field. The inclusion of mechanical and eventually chemical coupling in the sim - ulation can severely increase the computational expense (Pandey et  al. 2018). Recent improvements in computational efficiency, for example with the use of GPUs in DARTS (Khait et al. 2020), open up possibilities for a more detailed simulation, including these processes but also allowing for uncertainty quantification on multiple scales. This type of simulators will also open the pathway to studies of coupled thermo-hydro-mechanical– chemical processes. Fig. 12 Comparison of the maximum encountered pressure against fault to well distance for the range of flowrates, fault permeability and the highest and lowest injection temperature considered. The horizontal lines show the failure pressures as result of friction coefficients 0.15, 0.17, 0.19 and 0.12, 0.14, 0.16 for the injection temperature of 303.15 K and 313.15 K respectively Zaal et al. Geotherm Energy (2021) 9:12 Page 22 of 26 Fig. 13 Combined results for NPV and maximum encountered pore pressure difference as a function of FWD, flowrate and fault permeability. Both the highest and lowest injection temperature are considered. The horizontal lines show the pore pressure for failure-criterion exceedance assuming friction coefficients of 0.15, 0.17, 0.19 and 0.12, 0.14, 0.16 for the injection temperature of 303.15 K and 313.15 K respectively Conclusions This study presents a methodology to combined the assessment of the NPV and fault- stability analysis for direct-use geothermal energy production that is demonstrated in the Dutch context and is directly applicable to the WNB. THe methodology presented can be applied with minor changes to other sedimentary, direct-use systems and can be adapted to other settings with different geological and economic conditions. Reservoir simulations with a heterogeneous reservoir model representing the fluvial Delft Sand - stone considered the effect of a range of flow rates and injection temperatures. By plac - ing a vertical fault at varying distance from the fixed well positions and with different fault-permeability values the role of these conditions was evaluated. An NPV analysis was carried out using an updated economic model of the Dutch geothermal subsidy context. With simulations we resolved the pressure field inside the reservoir and close to the fault. The resulting pore pressure change was used as input for the Mohr–Coulomb failure criterion to identify the friction coefficient that results in unstable behaviour. Our main findings from these evaluations for a heterogeneous geothermal reservoir include: • both NPV and fault stability improve with a larger Fault-to-Well Distance (FWD) when fault permeability is lower than the average reservoir permeability, • fault permeability higher than the average reservoir permeability results in relatively small NPV benefits when placing the wells closer to the fault, • fault permeability at similar levels as the average reservoir permeability results in no changes in the generated NPV when choosing a different FWD, • a sealing fault limits the reservoir volume for fluid flow but heat conduction across the fault remains significant over the simulated period of 30 years, • compared to a case with a transmissive fault, a case with a sealing fault at close prox- imity to the injection well significantly alters the pressure field both between the injector and the fault but also further away from the injector inside the reservoir domain, Zaal  et al. Geotherm Energy (2021) 9:12 Page 23 of 26 • the fluid density and viscosity changes caused by a re-injection temperature decrease of 10 K (313.15 K to 303.15 K) are sufficient to cause changes in fault stability by shifting the maximum pore pressure encountered by up to 1 MPa compared to a case with higher temperature, and • the highest friction coefficient leading to exceedance of the Mohr–Coulomb failure criterion under the geological assumptions of our study was 0.17 and just below 0.14 for an injection temperature of 303.15 K and 313.15 K respectively. Such low values of friction coefficient can be encountered in clay rich fault gouges. • the presence of a transmissive fault can drastically alter the resulting shape of the cold-water plume by diverting and connecting different channelised parts of the res - ervoir. Future directions of our research include the analysis of the fault impact with variable well positions to better explore the impact of heterogeneous permeability and its inter- action with the fault hydraulic properties. Additionally, the combination of a heteroge- neous reservoir and heterogeneous fault properties will be further explored for its effects on both NPV and fault stability. We expect that the combined assessment of thermal, hydraulic, mechanical and economic aspects will become increasingly relevant as geo- thermal operations expand and the density of direct-use geothermal systems increases. A comprehensive stress-field characterisation with increased spatial density of support - ing data will be essential for a better understanding and quantification of the economic and fault-stability effects of geothermal operations. Acknowledgements The authors would like to acknowledge David Bruhn for providing feedback on the manuscript draft and Anne Pluymak- ers and Phil Vardon for their insights on geomechanics. Additionally, we would like to thank our three reviewers for their detailed, critical and constructive feedback and comments. Authors’ contributions CZ developed the fault stability model and implemented the economic model, performed the simulations in DARTS and contributed to the visualisation and the writing of the manuscript. AD conceptualised the idea for the work, supervised CZ and contributed to the development of the fault stability and economic model, preparing the visualisations, perform- ing the simulations in DARTS and writing the manuscript. FV supervised CZ and provided feedback on the DARTS simu- lations, fault stability and economic analysis, contributed to the visualisations and to writing the manuscript. All authors of this article contributed to reviewing the manuscript. All authors read and approved the final manuscript. Funding This study is partly supported by the DESTRESS project, funded by the European Union’s Horizon 2020 research and innovation program under Grant agreement ID 691728. FV was funded by the Delft Technology Fellowship of Delft University of Technology. Availability of data and materials The dataset generated and/or analysed during the current study is available from the corresponding author upon reasonable request. Declarations Competing interests The authors declare that they have no competing interests. Received: 18 December 2020 Accepted: 8 April 2021 References Altmann J, Müller B, Müller T, Heidbach O, Tingay M, Weißhardt A. Pore pressure stress coupling in 3d and consequences for reservoir stress states and fault reactivation. Geothermics. 2014;52:195–205. https:// doi. org/ 10. 1016/j. geoth ermics. 2014. 01. 004. Zaal et al. Geotherm Energy (2021) 9:12 Page 24 of 26 Bonté D, van Wees J-D, Verweij J. Subsurface temperature of the onshore Netherlands: new temperature dataset and modelling. Neth J Geosci Geol en Mijnbouw. 2012;91(4):491–515. https:// doi. org/ 10. 1017/ s0016 77460 00003 54. Buijze L, Bogert PAJ, Wassing BBT, Orlic B. Nucleation and arrest of dynamic rupture induced by reservoir depletion. J Geophys Res Solid Earth. 2019a;124(4):3620–45. https:// doi. org/ 10. 1029/ 2018j b0169 41. Buijze L, van Bijsterveldt L, Cremer H, Paap B, Veldkamp H, Wassing BB, van Wees J-D, van Yperen GC, ter Heege JH. Review of induced seismicity in geothermal systems worldwide and implications for geothermal systems in the Netherlands. Neth J Geosci. 2019b;98:E13. https:// doi. org/ 10. 1017/ njg. 2019.6. Byerlee J. Friction of rocks. Pure Appl Geophys. 1978;116(4–5):615–26. https:// doi. org/ 10. 1007/ bf008 76528. Candela T, Osinga S, Ampuero J-P, Wassing B, Pluymaekers M, Fokker PA, Wees J-D, Waal HA, Muntendam-Bos AG. Depletion-induced seismicity at the groningen gas field: Coulomb rate-and-state models including differential compaction effect. J Geophys Res Solid Earth. 2019;124(7):7081–104. https:// doi. org/ 10. 1029/ 2018j b0166 70. CBS. Aardgas en elektriciteit, gemiddelde prijzen van energieverbruikers. 2019. https:// opend ata. cbs. nl/ statl ine/#/ CBS/ nl/ datas et/ 81309 NED/ table? froms tatweb. Accessed 02 Nov 2019. Compernolle T, Welkenhuysen K, Petitclerc E, Maes D, Piessens K. The impact of policy measures on profitability and risk in geothermal energy investments. Energy Econ. 2019;84:104524. https:// doi. org/ 10. 1016/j. eneco. 2019. 104524. Daniilidis A, Doddema L, Herber R. Risk assessment of the Groningen geothermal potential: from seismic to reservoir uncertainty using a discrete parameter analysis. Geothermics. 2016;64:271–88. https:// doi. org/ 10. 1016/j. geoth ermics. 2016. 06. 014. Daniilidis A, Alpsoy B, Herber R. Impact of technical and economic uncertainties on the economic performance of a deep geothermal heat system. Renew Energy. 2017;114:805–16. https:// doi. org/ 10. 1016/j. renene. 2017. 07. 090. Daniilidis A, Khait M, Saeid S, Vokov D, Bruhn D. A high performance framework for the optimization of geothermal systems, comparing energy production and economic output. In: Proceedings world geothermal congress. 2020a;491–515. Daniilidis A, Nick HM, Bruhn DF. Interdependencies between physical, design and operational parameters for direct use geothermal heat in faulted hydrothermal reservoirs. Geothermics. 2020b;86:101806. https:// doi. org/ 10. 1016/j. geoth ermics. 2020. 101806. Daniilidis A, Nick HM, Bruhn DF. Interference between geothermal doublets across a fault under subsurface uncertainty; implications for field development and regulation. Geothermics. 2021a;91:102041. https:// doi. org/ 10. 1016/j. geoth ermics. 2021. 102041 Daniilidis A, Saeid S, Doonechaly NG. The fault plane as the main fluid pathway: Geothermal field development options under subsurface and operational uncertainty. Renewable Energy. 2021b;171:927–946. https:// doi. org/ 10. 1016/j. renene. 2021. 02. 148 Davis SD, Frohlich C. Did (or will) fluid injection cause earthquakes?—criteria for a rational assessment. Seismol Res Lett. 1993;64(3–4):207–24. https:// doi. org/ 10. 1785/ gssrl. 64.3- 4. 207. de Boer S, Bourdon E, Butterman G, et al. Evaluatie Algemeen instrumentarium geothermie. 2016. Donselaar M, Groenenberg R, Gilding D. Reservoir geology and geothermal potential of the delft sandstone member in the west Netherlands basin. In: Proceedings world geothermal congress. 2015;2015:1–9. Donselaar M. Reservoir architecture modelling for geothermal energy production—case study of the delft sandstone member, west Netherlands. In: 78th EAGE conference and exhibition 2016. EAGE Publications BV, 2016;1–5. https:// doi. org/ 10. 3997/ 2214- 4609. 20160 0596. Evans KF, Zappone A, Kraft T, Deichmann N, Moia F. A survey of the induced seismic responses to fluid injection in geothermal and CO2 reservoirs in Europe. Geothermics. 2012;41:30–54. https:// doi. org/ 10. 1016/j. geoth ermics. 2011. 08. 002. Eyerer S, Schifflechner C, Hofbauer S, Bauer W, Wieland C, Spliethoff H. Combined heat and power from hydrothermal geothermal resources in Germany: an assessment of the potential. Renew Sustain Energy Rev. 2020;120:109661. https:// doi. org/ 10. 1016/j. geoth ermics. 2011. 08. 002. Gaucher E, Schoenball M, Heidbach O, Zang A, Fokker PA, van Wees J-D, Kohl T. Induced seismicity in geothermal reser- voirs: a review of forecasting approaches. Renew Sustain Energy Rev. 2015;52:1473–90. https:// doi. org/ 10. 1016/j. rser. 2015. 08. 026. Gilding D. Heterogeneity determination of the delft subsurface for heat flow modelling. Master’s thesis, Applied Earth Sciences, Delft University of Technology. 2010. Guglielmetti L, Moscariello A, Meyer M, de Bono CN, Berthoud NA. Geothermie 2020: Exploration and development of geothermal energy in geneva. In: 81st EAGE conference and exhibition 2019 Workshop Programme. European Association of Geoscientists & Engineers. 2019;1–4. https:// doi. org/ 10. 3997/ 2214- 4609. 20190 1962. Hettema M. Analysis of mechanics of fault reactivation in depleting reservoirs. Int J Rock Mech Mining Sci. 2020;129:104290. Ikari MJ, Saffer DM, Marone C. Frictional and hydrologic properties of clay-rich fault gouge. J Geophys Res. 2009. https:// doi. org/ 10. 1029/ 2008j b0060 89. Ikari MJ, Marone C, Saffer DM. On the relation between fault strength and frictional stability. Geology. 2011;39(1):83– 6. https:// doi. org/ 10. 1130/ g31416.1. Jaeger JC, Cook NG, Zimmerman R. Fundamentals of rock mechanics. Hoboken: Wiley; 2009. Jansen JD, Singhal P, Vossepoel FC. Insights from closed-form expressions for injection- and production-induced stresses in displaced faults. J Geophys Res Solid Earth. 2019;124(7):7193–212. https:// doi. org/ 10. 1029/ 2019J B0179 32. Khait M, Voskov D. Adaptive parameterization for solving of thermal/compositional nonlinear flow and transport with buoyancy. SPE J. 2018a;23(02):522–34. https:// doi. org/ 10. 2118/ 182685- pa. Khait M, Voskov D. Operator-based linearization for efficient modeling of geothermal processes. Geothermics. 2018b;74:7–18. https:// doi. org/ 10. 1016/j. geoth ermics. 2018. 01. 012. Khait M. Delft advanced research terra simulator: general purpose reservoir simulator with operator-based lineariza- tion. 2019. http:// resol ver. tudel ft. nl/ uuid: 5f0f9 b80- a7d6- 488d- 9bd2- d68b9 d7b4b 87. Zaal  et al. Geotherm Energy (2021) 9:12 Page 25 of 26 Khait M, Voskov D, Zaydullin R. High performance framework for modelling of complex subsurface flow and trans- port applications. In: Conference proceedings, ECMOR XVII. 2020;2020:1–18. https:// www. earth doc. org/ conte nt/ papers/ 10. 3997/ 2214- 4609. 20203 5188. Knapek E, Dilger G, Hegele H. Deutschlands zukunft mit erdwärme. Tech. rep., Wirtschaftsforum Geothermie, GtV- Bundesverband Geothermie. 2013. Knoblauch TA, Trutnevyte E, Stauffacher M. Siting deep geothermal energy: Acceptance of various risk and benefit scenarios in a Swiss-German cross-national study. Energy Policy. 2019;128:807–16. https:// doi. org/ 10. 1016/j. enpol. 2019. 01. 019. Lepillier B, Daniilidis A, Gholizadeh ND, Bruna P-O, Kummerow J, Bruhn D. A fracture flow permeability and stress dependency simulation applied to multi-reservoirs, multi-production scenarios analysis. Geotherm Energy. 2019;7(1):1–16. https:// doi. org/ 10. 1186/ s40517- 019- 0141-8. Liang X, Xu T, Feng B, Jiang Z. Optimization of heat extraction strategies in fault-controlled hydro-geothermal reser- voirs. Energy. 2018;164:853–70. https:// doi. org/ 10. 1016/j. energy. 2018. 09. 043. Lukawski M, Anderson B, Augustine C, et al. Cost analysis of oil, gas, and geothermal well drilling. J Pet Sci Eng. 2014;118:1–14. Lund J, Toth A. Direct utilization of geothermal energy 2020 worldwide review. In: Proceedings world geothermal congress. 2020;2020:1–39. Markou N, Papanastasiou P. Petroleum geomechanics modelling in the eastern mediterranean basin: analysis and application of fault stress mechanics. Oil Gas Sci Technol Revue d’IFP Energies nouvelles. 2018;73:57. https:// doi. org/ 10. 2516/ ogst/ 20180 34. Mechelse E. The in-situ stress field in the Netherlands: Regional trends, local deviations and an analysis of the stress regimes in the northeast of the Netherlands. Master’s thesis, Applied Earth Sciences, Delft University of Technol- ogy. 2017. Mijnlieff HF. Introduction to the geothermal play and reservoir geology of the Netherlands. Neth J Geosci. 2020. https:// doi. org/ 10. 1017/ njg. 2020.2. Ministerie van Economische Zaken en Klimaat, Staatscourant van het Koninkrijk der Nederlanden, nr 9735. 2019. https:// zoek. offic ieleb ekend makin gen. nl/ stcrt- 2019- 9735. Accessed 05 Nov 2019. Moeck I, Kwiatek G, Zimmermann G. Slip tendency analysis, fault reactivation potential and induced seismicity in a deep geothermal reservoir. J Struct Geol. 2009;31(10):1174–82. https:// doi. org/ 10. 1016/j. jsg. 2009. 06. 012. Moeck IS. Catalog of geothermal play types based on geologic controls. Renew Sustain Energy Rev. 2014;37:867–82. https:// doi. org/ 10. 1016/j. rser. 2014. 05. 032. Morris A, Ferrill DA, Henderson DB. Slip-tendency analysis and fault reactivation. Geology. 1996;24(3):275. https:// doi. org/ 10. 1130/ 0091- 7613% 281996% 29024% 3C0275% 3Asta afr% 3E2.3. co% 3B2 Moscariello A, Guglielmetti L, Omodeo-Salé S, De Haller A, Eruteya OE, Lo H, Clerc N, Makhloufi Y, Do Couto D, De Oliveira GF, Perozzi L, De Oliveira F, Hollmuller P, Quiquerez L, Nawratil De Bono C, Martin F, Meyer M. Heat pro- duction and storage in western switzerland: advances and challenges of intense multidisciplinary geothermal exploration activities, an 8 years progress report. In: Proceedings world geothermal congress. 2020;2020:1–10. O’Sullivan MJ, Pruess K, Lippmann MJ. State of the art of geothermal reservoir simulation. Geothermics. 2001;30(4):395–429. https:// doi. org/ 10. 1016/ S0375- 6505(01) 00005-0. Osundare O, Teodoriu C, Faclone G, et al. Estimation of plugging and abandonment costs based on different EU regulations with application to geothermal wells. In: 43rd workshop on geothermal reservoir engineering. 2018;1–13. Pandey S, Vishal V, Chaudhuri A. Geothermal reservoir modeling in a coupled thermo-hydro-mechanical-chemical approach: a review. Earth-Sci Rev. 2018;185:1157–69. https:// doi. org/ 10. 1016/j. earsc irev. 2018. 09. 004. Parisio F, Vilarrasa V, Wang W, Kolditz O, Nagel T. The risks of long-term re-injection in supercritical geothermal systems. Nat Commun. 2019. https:// doi. org/ 10. 1038/ s41467- 019- 12146-0. Peaceman D. Interpretation of well-block pressures in numerical reservoir simulation with nonsquare grid blocks and anisotropic permeability. Soc Pet Eng J. 1983;23(03):531–43. https:// doi. org/ 10. 2118/ 10528- pa. Pluymaekers M, Kramers L, van Wees J-D, Kronimus A, Nelskamp S, Boxem T, Bonté D. Reservoir characterisation of aqui- fers for direct heat production: Methodology and screening of the potential reservoirs for the Netherlands. Neth J Geosci Geol en Mijnbouw. 2012;91(4):621–36. https:// doi. org/ 10. 1017/ s0016 77460 00004 1x. Rijksdienst voor Ondernemend Nederland. Regeling nationale EZ subsidies - Risico’s dekken voor Aardwarmte. Zevende en achtste openstelling: Handleiding Garantiestelling tegen het risico van misboring. 2019a. Rijksdienst voor Ondernemend Nederland. SDE+ najaar. 2019b. www. rvo. nl/ sites/ defau lt/ files/ 2019/ 10/ Broch ure% 20SDE% 20naj aar% 202019. pdf. Accessed 04 Nov 2019. Rijksoverheid. 2019. Wijzigingen vennootschapsbelasting, belastingplan. 2020. https:// www. rijks overh eid. nl/ onder wer- pen/ belas tingp lan/ belas tingw ijzig ingen- voor- onder nemers/ venno otsch apsbe lasti ng. Accessed 06 Nov 2019. Romera J. iapws.iapws97 module. 2020. https:// iapws. readt hedocs. io/ en/ latest/ iapws. iapws 97. html. Accessed 27 Apr Saeid S, Yang W, Daniilidis A, Khait M, Vokov D, Bruhn D. Lifetime and energy prediction of geothermal systems: Uncer- tainty analysis in highly heterogeneous geothermal reservoirs (Netherlands). In: Proceedings world geothermal congress. 2020;2020:1–8. Salimzadeh S, Paluszny A, Nick HM, Zimmerman RW. A three-dimensional coupled thermo-hydro-mechanical model for deformable fractured geothermal systems. Geothermics. 2018;71:212–24. https:// doi. org/ 10. 1016/j. geoth ermics. 2017. 09. 012. Schoof F, van der Hout M, van Zanten J, et al. Masterplan Aardwarmte in Nederland. Tech. rep., Platform Geothermie, DAGO, Stichting Warmtenetwerk, EBN. 2018. https:// www. ebn. nl/ wp- conte nt/ uploa ds/ 2018/ 05/ 20180 529- Maste rplan- Aardw armte- in- Neder land. pdf. Schoots K, Hammingh P. Klimaat- en Energieverkenning 2019. Tech. rep., Den Haag: Planbureau voor de Leefomgeving. 2019. https:// www. pbl. nl/ sites/ defau lt/ files/ downl oads/ pbl- 2019- klima at- en- energ iever kenni ng- 2019- 3508. pdf. Accessed 04 Oct 2019. Zaal et al. Geotherm Energy (2021) 9:12 Page 26 of 26 Seithel R, Gaucher E, Mueller B, Steiner U, Kohl T. Probability of fault reactivation in the bavarian molasse basin. Geother- mics. 2019;82:81–90. https:// doi. org/ 10. 1016/j. geoth ermics. 2019. 06. 004. Shetty S, Voskov D, Bruhn D. Numerical strategy for uncertainty quantification in low enthalpy geothermal projects. In: 43rd workshop on geothermal reservoir engineering. 2018;1–16. Soltanzadeh H, Hawkes CD. Assessing fault reactivation tendency within and surrounding porous reservoirs during fluid production or injection. Int J Rock Mech Mining Sci. 2009;46(1):1–7. Stauffacher M, Muggli N, Scolobig A, Moser C. sep, Framing deep geothermal energy in mass media: the case of Switzer - land. Technol Forecast Soc Change. 2015;98:60–70. https:// doi. org/ 10. 1016/j. techf ore. 2015. 05. 018. Takahashi M, Mizoguchi K, Kitamura K, Masuda K. Eec ff ts of clay content on the frictional strength and fluid transport property of faults. J Geophys Res. 2007. https:// doi. org/ 10. 1029/ 2006j b0046 78. Tembe S, Lockner DA, Wong T-F. Eec ff t of clay content and mineralogy on frictional sliding behavior of simulated gouges: Binary and ternary mixtures of quartz, illite, and montmorillonite. J Geophys Res. 2010. https:// doi. org/ 10. 1029/ 2009j b0063 83. TNO. Economic model ( Thermogis). 2019. https:// www. therm ogis. nl/ en/ econo mic- model. Accessed 02 Nov 2019. Troiano A, Giuseppe MGD, Troise C, Tramelli A, Natale GD. A coulomb stress model for induced seismicity distribution due to fluid injection and withdrawal in deep boreholes. Geophys J Int. 2013;195(1):504–12. https:// doi. org/ 10. 1093/ gji/ ggt229. Van Adrichem Boogaert H, Kouwe W, Stratigraphic nomenclature of the Netherlands, revision and update by Rijks Geolo- gische Dienst (RGD) and Netherlands oil and gas exploration and production association (NOGEPA). Mededelingen Rijks Geologische Dienst 50. 1993. van Balen R, van Bergen F, de Leeuw C, Pagnier H, Simmelink H, van Wees J, Verweij J. Modelling the hydrocarbon generation and migration in the west Netherlands basin, the Netherlands. Neth J Geosci Geol en Mijnbouw. 2000;79(1):29–44. Van den Bogert P. Impact of various modelling options on the onset of fault slip and fault slip response using 2-dimen- sional finite-element modelling. Restricted, Report No. SR 15. 2015. van den Bosch R, Flipse B, Vorage R. Stappenplan Winning Aardwarmte voor Glastuinbouw. Tech. rep., Kas als energie- bron. 2013. https:// www. kasal sener giebr on. nl/ conte nt/ resea rch/ 14899_ Stapp enplan_ Aardw armte_ 12- 2013. pdf. Accessed 06 Oct 2019. van ’t Spijker H, Ungemach P. Definition of Electrosubmersible Pump (ESP) design and selection workflow. 2016. https:// www. kasal sener giebr on. nl/ conte nt/ user_ upload/ Defin ition_ of_ ESP_ desig n-_ final_ versi on_ 20160 629_- public_ v2. pdf. Accessed 02 Nov 2019. Voskov DV. Operator-based linearization approach for modeling of multiphase multi-component flow in porous media. J Comput Phys. 2017;337:275–88. https:// doi. org/ 10. 1016/j. jcp. 2017. 02. 041. Voskov D. DARTS—delft advanced terra reservoir simulator. 2020. https:// darts. citg. tudel ft. nl. Accessed 09 July 2020. Vrolijk PJ, Urai JL, Kettermann M. Clay smear: review of mechanisms and applications. J Struct Geol. 2016;86:95–152. https:// doi. org/ 10. 1016/j. jsg. 2015. 09. 006. Wang Y, Khait M, Voskov D, et al. Benchmark test and sensitivity analysis for geothermal applications in the Netherlands. In: 44th workshop on geothermal reservoir engineering. 2019;1–11. Wang Y, Voskov D, Khait M, Bruhn D. An efficient numerical simulator for geothermal simulation: a benchmark study. Appl Energy. 2020;264:114693. https:// doi. org/ 10. 1016/j. apene rgy. 2020. 114693. Willems C, Nick H, Goense T, Bruhn D. The impact of reduction of doublet well spacing on the net present value and the life time of fluvial hot sedimentary aquifer doublets. Geothermics. 2017a;68:54–66. https:// doi. org/ 10. 1016/j. geoth ermics. 2017. 02. 008. Willems CJ, Nick HM, Weltje GJ, Bruhn DF. An evaluation of interferences in heat production from low enthalpy geother- mal doublets systems. Energy. 2017b;135:500–12. https:// doi. org/ 10. 1016/j. energy. 2017. 06. 129. Zang A, Oye V, Jousset P, Deichmann N, Gritto R, McGarr A, Majer E, Bruhn D. Analysis of induced seismicity in geothermal reservoirs—an overview. Geothermics. 2014;52:6–21. https:// doi. org/ 10. 1016/j. geoth ermics. 2014. 06. 005. Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Journal

Geothermal EnergySpringer Journals

Published: Apr 20, 2021

There are no references for this article.