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

The installed capacity of geothermal systems for direct use of heat is increasing worldwide. As their number and density is increasing, the their interaction with subsurface faults becomes more important as they could lead to safety risks from induced seismicity. Assessment and management of such risks is essential for the further development and extension of geothermal energy for heating. At the same time, the economic output 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 exceedance 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.

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 economically 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(Daniilidis et al. , 2021a. 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). Interaction of direct-use geothermal operations with pre-faults in conduction-dominated geothermal settings has shown significant effects on energy generation, system lifetime and economic output (Daniilidis et al. 2016(Daniilidis et al. , 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(Daniilidis et al. , 2020b. High fault permeability combined with the relative fault block positions can also have positive effects on energy and system lifetime 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 conduction 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). Therefore, the focus on fault stability of pre-existing faults is more relevant for these geological 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 Coulombstress changes attributed to fluid injection as the main trigger for seismic events in an Zaal et al. Geotherm Energy (2021) 9:12 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 (Parisio 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 seismicity, 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-dominated geothermal fields on economic output and fault stability for a range of fault properties and operational choices. The Delft Advanced Research Terra Simulator (DARTS) (Khait 2019) is used to perform reservoir simulations of the Delft sandstone, a heterogeneous 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 reinjection temperature of the geothermal doublet. Our final analysis compares the generated 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 outline 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 simulation 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 development 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 Netherlands 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 meandering 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 extensional 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 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 contains 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 typical 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 failure criterion. An overview of the different parts and interconnections is presented in Fig. 2.  Table 1  Zaal et al. Geotherm Energy (2021) 9:12 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; and has demonstrated fast and accurate performance for geothermal applications (Wang et al. , 2020. It is based on the framework of Operator-Based Linearization (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 c components and n p phases is given 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: in which φ is the porosity, x cj the mole fraction of component c in phase j, s j the phase saturations, ρ j the phase molar density, n p the number of phases present, n c the number of components, q j the phase rate per unit volume and v j the phase velocity.
The equation for the conservation of energy accounts for the stored energy, convection, conduction and sink and source terms: in which U j is the phase internal energy, U r the internal energy of the rock, h j the phase enthalpy, κ the thermal conductivity and ∇T the temperature gradient.
Thermal convection is controlled by Darcy's Law: in which K is the permeability tensor, k rj the relative permeability, µ j the phase viscosity, p j the vector of pressures in phase j, g j is the gravity term and ∇d is the vector of depths. 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): (1) and in which V is the control volume, q j , which is equal to q j V is a source of phase j, �ψ l is the difference in pressures between the blocks connected via interface l, T l is the difference in temperature between blocks connected via interface l, Ŵ l j , which is equal to Ŵ l k l rj µ l j is the phase transmissibility, Ŵ l is the constant geometrical part of transmissibility, Ŵ l c , equal to Ŵ l κ , equal to φ n p j=1 s l j γ l j − κ r + κ r the thermal transmissibility.
Equations 4 and 5 are approximated in time using a fully implicit method, which introduces 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 Khait and Voskov 2018a, b). During simulations the operators are evaluated based on multi-linear interpolations, which reduces computation 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 demonstrated in recent publications (Voskov 2017;Khait and Voskov 2018a, b;Wang et al. 2019, Daniilidis et al. 2020a, Saeid et al. 2020. The well index is calculated according 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): in which τ f is the shear strength of the rock (MPa), C the cohesion of the rock (MPa), µ the apparent coefficient of internal friction (−) (Byerlee 1978), σ ft the total stress normal to the failure plane (MPa), P the pore pressure (MPa) and σ ft − P the effective shear stress σ ′ ft . (4) x c jρ j s j n − �t l np j=1 x l cj ρ l j Ŵ l j �ψ l + �t np j=1 ρ p x cj q j = 0, c = 1, . . . , n c , 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 ( τ ) equals or exceeds the frictional sliding resistance, also called the failure shear stress ( τ f ) (Jaeger 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): Equation 6 gives the failure criterion describing the shear strength of the fault. In combination 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. In the reservoir simulations (see "Reservoir model, sensitivity study and input parameters" section) we assume the presence of a vertical fault in a normal faulting regime. Therefore, failure in this case does not take place on the σ 1 -σ 3 plane, but rather on the σ 2 -σ 3 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 σ 1 -σ 3 plane (a). Pore pressure increase required for a vertical fault to reach failure along the σ 2 -σ 3 plane (b) and effect of thermal stresses for a laterally confined reservoir where the thermal 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 σ 2 -σ 3 plane of a thermally stressed reservoir (d) plane that defines the angle of the vertical fault plane with respect to the horizontal components 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): in which � σ T is the thermal stress change (MPa), α t is the rock linear thermal expansion 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 implemented as outlined in previous research (Daniilidis et al. 2017(Daniilidis et al. , 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 contribute 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 drilling costs are computed as (TNO 2019): in which MD is the measured depth (m) and the doublet costs are in Euro (€). Experience 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: where C t 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 t ) in € is computed according to: where I heat is the income from the produced heat, I subsidy is the subsidy income, CapEx are the capital expenditures, OpEx annual the annual operational expenses, OpEx pumping 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 cost ), construction and unforeseen construction costs associated with preparation, drilling and installation of the equipment and drill site. Equipment includes the injector and producer 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 workover 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 maintenance 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 generated 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 : 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 subsidized 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 Economische 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 production, which is fixed, the market price of thermal energy and the full-load production 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): where q is the volumetric flow rate (m 3 /s), T the temperature difference between the produced and injected water (K), C p f the specific heat capacity of the fluid at constant pressure (J/(kgK)), ρ f the density fluid (kg/m 3 ) 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): where q is the volumetric flow rate ( m 3 /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(Daniilidis et al. , 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 heterogeneous reservoir model for the Delft area (Fig. 4) comprised of approx. 100 k elements.  Table 3 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): 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 studies Shetty et al. 2018;Daniilidis et al. 2020a). The top and bottom boundary conditions act as infinitely large domains with a fixed temperature and provide 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.
(15) K = 0.0633e 29.507φ .    For each of the cases in the sensitivity study, a simulation is performed for a combination of parameters for both reservoir conditions and operation options. These parameters 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-stability 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 temperature 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 permeability. 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 pronounced 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 permeability 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 3 /h with minor reductions in the order of 0.3 MW and after circa 7 years for the high rate of 375 m 3 /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 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 both high (375 m 3 /h) and low (m 3 /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. Contrary 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 transmissive 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 buildup 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 dissipation 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 temperature 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 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 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 3 /h) for each injection temperature depending on the fault distance to the wells in these cases (Fig. 10). We find that for the lowest flow rate of 150 m 3 /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, increasing the rate from 150 to 375 m 3 /h (a factor 2.5) increases the NPV from 5 to 22.5 M€  , d). Data shown use an injection temperature of 303.15 K and a flowrate of 9000 m 3 /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 (a factor 4.5). In the case of the highest rate of 375 m 3 /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 permeability 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 pressure (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 penetrates 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 initial stress state and the thermal stress states of the injection temperatures considered (Fig. 11a). The σ 2 and σ 3 are equally affected by the thermal stress and therefore the shapes of the pore pressure curves remain similar, with only the values of the pore pressure leading to failure is increasing. The critical angle between the maximum horizontal stress S H max (which is σ 2 ) and the fault strike shows that the angle range of 25° to 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 temperature by 10 K is almost equivalent to the effect of a decrease in flow rate of 75 m 3 /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 60 m and 90 m FWD with a rate of 375 m 3 /h and with a 60 m FWD for the 300 m 3 /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.

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 sealing 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 permeability 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 transmissive 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(Daniilidis et al. , 2020b(Daniilidis et al. , 2021aWillems et al. 2017a). For the flow rate ranges considered here and assuming current heat and electricity prices, the additional revenues from high flow rates outweigh the costs of increased pumping requirements. 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 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 suggested 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 highest 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 3 /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 initial 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 considered. 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 3 /h result in pressure levels in excess of σ 3 (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 significant mineralogy changes are incurred (Ikari et al. 2011). Other experimental studies 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 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 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. Moreover, 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 poroelastic 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 provides 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 magnitude and orientation, at a resolution higher than presently available, is crucial. This should be complemented with fault-gouge friction coefficient characterisation in the lab. Fig. 11 Pore pressure required to reach failure along the σ 2 -σ 3 plane for the initial stress state and the thermal stress states considered for a range of friction coefficients (a). The critical fault strike angle with respect to S Hmax for the different stress states for the range of friction coefficients (b) 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, observational 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 conversion 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 magnitude are available, the evaluation of all potential development options and stress conditions 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 simulation can severely increase the computational expense (Pandey et al. 2018). Recent improvements in computational efficiency, for example with the use of GPUs in DARTS , 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-mechanicalchemical processes.

Conclusions
This study presents a methodology to combined the assessment of the NPV and faultstability 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 Sandstone considered the effect of a range of flow rates and injection temperatures. By placing 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 proximity 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, 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 • 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 reservoir.
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 interaction with the fault hydraulic properties. Additionally, the combination of a heterogeneous 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 geothermal operations expand and the density of direct-use geothermal systems increases. A comprehensive stress-field characterisation with increased spatial density of supporting data will be essential for a better understanding and quantification of the economic and fault-stability effects of geothermal operations.