Influential factors on the development of a low-enthalpy geothermal reservoir: A sensitivity study of a realistic field

A realistic deep low-enthalpy geothermal reservoir based on real data with high detail and complicated sedimentary structure is utilized to perform sensitivity analyses of the geological features influencing reservoir properties. We perform simulations using the Delft Advanced Research Terra Simulator (DARTS). Compelling numerical performance of DARTS makes it suitable for handling a large ensemble of models including efficient sensitivity and uncertainty analyses. The major finding is that shale facies, generally ignored in hydrocarbon reservoir simulations, can significantly extend the predictive lifetime of geothermal reservoirs exploited by deep well doublets. It is important to accurately account for the shale facies in the simulation, though with an additional computational overhead. The overburden layers can improve doublet performance, but the impact depends on reservoir heterogeneity. In addition, heterogeneity will also divert the flow path with even a minor shift in the well placement. The discharge rate, an essential parameter of geothermal operation strategy, inversely corresponds to the doublet lifetime but positively correlates with the energy production for studied parameter ranges. Low sensitivity of doublet lifetime to vertical-horizontal permeability ratio and permeability-porosity correlation is observed. All these systematic findings for a realistic geothermal field with characterization at unprecedented level of detail can help to provide a general guideline for forward simulation and farther improve the profitability of geothermal energy production in realistic deep geothermal reservoirs through computer-assisted modeling and optimization. © 2021 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).

Influential factors on the development of a low-enthalpy geothermal reservoir: A sensitivity study of a realistic field

Introduction
In the context of global warming and the consequent strife to control greenhouse gas (GHG) emissions, deep geothermal energy in the subsurface, long overlooked in many parts of the world, has started to attract increased attention [1], since GHG emissions during its life cycle are low [2]. Besides, the worldwide distribution [3] and astonishing subsurface reserves [4] make deep geothermal energy promising and potentially competitive to other low-carbon energy resources. The global installed geothermal energy capacity has reached 13 GWe by the end of 2018 [5] and increases rapidly. In particular cases, the development of geothermal energy can take advantage of depleted oil-gas fields [6e8] to reduce early-stage investment, e.g., for drilling a pair of deep geothermal wells or synergy potential with oil exploitation [9] to maximize the Net Present Value (NPV).
Although deep geothermal energy is an attractive resource, the uncertainty in size and properties of subsurface geothermal reservoirs directly impacts the technical and economic feasibility of geothermal projects [10]. The lifetime of deep geothermal wells and heat recovery rate of the reservoir often vary a lot with both reservoir parameters and operational management [11].
Shale facies are generally regarded as flow barriers due to their low permeability. In isothermal simulations, the shale facies are often removed from the reservoir characterization to reduce the model size. However, for geothermal reservoirs, this simplified treatment will introduce inaccuracy to the simulation results. Since heat transport depends on both thermal convective and conductive effects, the impermeable shale facies still impact the heat transfer process via conduction. With process-based realizations, Crooijmans et al. [12] found that shale facies play an important role in doublet performance when the net-to-gross ratio is less than 50%. But the utilized reservoir parameters are based on conceptual estimation. How the effects of the shale facies will play out in realistic geothermal reservoirs, where detailed geological information is necessarily constrained to the measurements (e.g., log and core data), remains an open question. Once the realistic impact of shale is clarified, more advanced or optimized operational strategies can therefore be proposed to improve the development performance. In this study, the impact of shale facies to lowenthalpy geothermal development is sufficiently considered and studied with the utilized realistic geological model.
Besides the sealing effect to trap heat and fluid [13], the burden layers in geothermal reservoirs can also provide thermal recharge to the payzone [14]. However, inconsistent treatment and recognition of burden layers still exist in the previous studies. In several studies, the burden layers are either taken as non-conductive (i.e., thermal insulation [15]) or found to have a weak impact on doublet performance (e.g., with synthetic process-based heterogeneous models [16]). In contrast, Willems et al. [17] and de Bruijn et al. [18] found the burden layers have a positive impact on heat production by charging the target aquifer based on a simplified homogeneous model. In reality, however, heterogeneity is a ubiquitous feature of subsurface systems [19e21]. It becomes more meaningful and representative to investigate the influence of burden layers in realistic contexts. One point to mention is that sufficient characteristic resolution is also important for the investigation. Oversimplification or coarser representation of burden layers will weaken, or even ignore their influence on geothermal doublet performance. Taking all of the above into account, the sensitivity of the geothermal production to overburdern layers with the presence of realistic heterogeneity is studied in this work.
Reasonable operational strategies are also essential for the success of a geothermal project with the existence of subsurface uncertainties. For example, the discharge rate directly impacts the thermal propagation and energy output [22,23]. Saeid et al. [24] found a higher discharge rate results in shorter doublet lifetime, sharper thermal front and less energy extraction. In addition, well interference in a multi-well system can influence doublet lifetime by changing the development of the cold water plume. To maximize the doublet lifetime, an optimized doublet configuration and well spacing are generally needed [17,25]. Daniilidis et al. [26] performed systematic sensitivity analyses of doublet performance to operational parameters and doublet configuration within a faulted reservoir, highlighting the importance of the well placement strategy for profitability of geothermal projects. Furthermore, Shetty et al. [16] concluded that a variation of up to 50% in lifetime can be achieved by a variation of less than 50m in doublet location after a set of systematic studies on the influence of doublet position on project lifetime. However, these findings are based on conceptual synthetic models. The doublet performance and reservoir feedback to different operational strategies with realistic geological structures can be very different, therefore special interests exist and new findings are expected in realistic reservoirs, like the one used in this study.
To analyse and investigate the influence of different parameters on geothermal energy production, numerical simulation has been widely employed for sensitivity analysis and optimization. Several simulators have been used for geothermal applications [27]. However, the physical complexity and large size of realistic models often challenge conventional simulation techniques. To be more specific, the complicated physical processes (e.g., multi-phase flow [28,29] or multi-component reactive transport [30]) encountered in geothermal applications require robust, flexible, and efficient solutions. In addition, large models with multi-million computational cells [31,32] are often needed to characterize and predict the behavior of a geothermal reservoir, which slows down the simulation process. Furthermore, to quantify uncertainties and optimize development strategies [33], a large ensemble of models are necessary to cover the wide range of parameter settings, which requires high-performance and reliability of forward simulations.
In this study, we utilize the Delft Advanced Research Terra Simulator (DARTS), developed for modelling of various industrial applications [34]. DARTS provides a unified framework for simulating various types of flow and transport processes for both forward and inverse problems in porous media [35], including the isothermal formulation mostly used for hydrocarbon reservoirs [36] and a thermal-compositional simulation used for geothermal applications. To solve highly nonlinear problems, the Operator-Based Linearization (OBL) approach [37] is employed in DARTS. The OBL approach has been proposed for multi-phase flow and transport applications and aims to improve the simulation performance [38]. For spatial and temporal discretization, the finitevolume fully-implicit method in combination with a two-point flux approximation on unstructured grids is implemented in DARTS. Besides conventional temporal and spatial discretizations, our framework also utilizes discretization in physical space with the OBL approach. Owing to this approach, DARTS can achieve fast simulation, and the capacity of it to simulate geothermal applications has been verified with systematic benchmark tests [39]. The integrated framework of DARTS [34] helps to implement and investigate the designed scenarios on sensitivity analyses, for example, deactivating and activating specific regions of the reservoir.
DARTS is therefore ideally suited as a tool to address the primary objective of this work: analysing the influence of different parameters on the lifetime of realistic geothermal projects. In this study, a realistic geothermal reservoir based on real data with an unprecedented detailed geological characterization is utilized in the simulation study. The impact of shale facies on the production of a realistic geothermal field is clearly identified. The sensitivity study demonstrates that the presence of realistic heterogeneity generates large variations in geothermal production while varying different parameters (e.g., overburden layers, well placement, etc.), which has not been revealed in previous studies. These new findings suggest that an accurate uncertainty analysis for low-enthalpy geothermal applications should be performed based on highly detailed reservoir models to account for a realistic variability in energy production.
The rest of this paper is organized as follows. First, we introduce the governing equations and correlations used to model the system. Then, a set of forward simulations are deployed and sensitivity analyses of different parameters (e.g., facies, overburden layer, discharge rate, well placement, etc.) to doublet performance is evaluated. The results of this study introduce the basis for a detailed discussion of these influential factors, which can provide guidance and instruction for future geothermal field operations around the world.

Methodology
The target geothermal reservoir is at low-enthalpy conditions. For low-enthalpy geothermal systems, liquid water is generally selected as the main working fluid with several additional chemicals. Without loss of generality, only reservoir brine with fixed composition is included in the studied system. The governing equations and nonlinear formulations describe the geothermal system as follows: v vt where Eq.
The different quantities in Eqs. (1) and (2) are defined as follows: C The fluid internal energy per volume is expressed as: C The rock internal energy per volume is expressed as: , which can be selected based on user's preference, e.g., the freezing point of the working fluid.
C The rock is compressible, which is reflected by the change of porosity through: where 4 0 is the initial porosity, c r is the rock compressibility [bars À1 ] and p ref is the reference pressure [bars], which can be selected based on user's preference, e.g., atmospheric pressure.
C In addition, Darcy's law is used to describe the fluid flow in the reservoir, A finite-volume approach is taken in spatial discretization on a general unstructured grid, and backward Euler approximation is applied for time discretization.
The discretized mass conservation equation is The discretized energy conservation equation is Here V is the control volume [m 3 ], j is the phase potential [bars], q f ¼q f V is the fluid source/sink term [m 3 /s], G l is the convective transmissibility which approximates the divergence term in the computational grids, G l c is the conductive transmissibility which is expressed with the help of geometrical shape coefficient (G l g ) as: The set of equations above is solved in fully-coupled fully-implicit manner in the DARTS framework utilizing the OBL approach [37]. In DARTS, the molar nonlinear formulation [28,40] is taken, where pressure and enthalpy are selected as primary variables for geothermal simulation. The Newton-Raphson approach is adopted to linearize the nonlinear formulation. The linear system of equations on each nonlinear iteration can be expressed in the following form: where J and r are the Jacobian matrix and residual defined at the kth nonlinear iteration respectively, u denotes the physical state (in this study, pressure and enthalpy) of all grid blocks.
The linearization approach requires an ensemble of the Jacobian matrix that contains the derivatives of different terms of the discretized governing equations with respect to primary variables. Conventional linearization includes accurate evaluation of fluid and rock properties and requires the solution of a highly nonlinear system of equations using the chain rule and inverse theorem [41]. Therefore, a nonlinear solver often performs extra iterations to capture small variations in properties, which are sometimes negligible because of the uncertainties in property evaluation and numerical nature of their representation.
In the OBL approach, a mesh is introduced in the space of nonlinear unknowns, and conventional property evaluations are only performed at the vertices of the mesh, namely the "supporting points". Various physical properties involved in discrete approximation of governing equations and dependent on state variables (nonlinear unknowns) are combined in state operators. Continuous physical description between supporting points is introduced through multilinear interpolation of operator values calculated in mesh points. The nonlinear physics of the original problem is represented by a piece-wise linear approximation, which often helps the nonlinear solver to converge faster during the simulation depending on the introduced physical mesh resolution [38]. A more detailed description of the OBL approach can be found in Ref. [37].

Sensitivity analysis using a geothermal reservoir case study
The development of a geothermal project can be impacted by a lot of factors from various aspects. In this section, we perform sensitivity analyses using a realistic geological model of a highly heterogeneous low-enthalpy geothermal field. Simulation studies performed on such a highly detailed geological model are not yet common for geothermal applications. Many studies conducted for hydrocarbon reservoirs are often misleading due to the specifics of energy flow. The influence of reservoir geological (e.g., facies, overburden layer, etc.) and operational (e.g., well placement, discharge rate, etc.) parameters on the doublet lifetime is analyzed below.

Reservoir description
The research area is located in the West Netherlands Basin (WNB), which is an inverted rift basin in the Netherlands. Sediments of this basin range in age from Jurassic to recent time and are overlying Triassic and older sediments. The Upper Jurassic and Lower Cretaceous start with the continental sediments of the Nieuwerkerk Formation. These sediments were deposited in subsiding half-grabens, while adjacent highs were subjected to erosion. In the Nieuwerkerk Formation, the Delft Sandstone has in general good reservoir properties [11,17,25].
"The Delft Sandstone is interpreted to be deposited as stacked distributary-channel deposits in a lower coastal plain setting resulting in massive sandstone sequences." [42] The thickness of the Delft Sandstone is influenced by the syn-rift deposition of the sediments and therefore it is of variable width with a thickness up to 130 m. The sandstone consists of fineto coarse-grained sand bodies, and their lateral continuity is difficult to predict. The Berkel Sandstone Member and Berkel Sands-Claystone Member have a shallow marine depositional setting. The facies range from upper shoreface to lower shoreface of a coastal-barrier system. Fig. 1 and Fig. 2 show the modelled porosity distribution of the target reservoir.
Owing to the detailed exploration data and recorded geological statistics from previous hydrocarbon extraction in the region, higher-fidelity data interpretation and geological modeling are enabled for the target reservoir in comparison with what is typically available and used for modelling low-enthalpy geothermal reservoirs. Furthermore, based on completeness of the fundamental data sets, it will be much easier to implement different realizations for sensitivity analysis in DARTS, which can be observed in the following investigations.

Numerical model description
The model extends by 9000 m Â 4200 m in the horizontal dimension. Two doublets were placed in the reservoir and operated with constant rate control. Since it is difficult to predict the lateral continuity, the reservoir boundary condition is defined as no-flow. Table 1 provides the parameter settings for the base case. As is shown, a permeability contrast by five orders of magnitude is present in the reservoir. With the water cycling in the system, we try to observe the evolution of the thermal plume in the reservoir and production temperature of both doublets by varying reservoir parameters and operational schemes. The project lifetime in this study has been defined as a drop of the production temperature to 335 K.

Sensitivity analysis
The sensitivity analysis is performed below for various reservoir parameters and operational management schemes.

Shale layers
The model includes intersecting sandstone and shale facies. The facies distribution (Fig. 3) corresponds to circa 0.8 million grid blocks for the sandstone (Fig. 4a) and 2.4 million blocks for shale facies (Fig. 4b). To evaluate the impact of the shales on thermal breakthrough, two scenarios are simulated: including and excluding the shale grid blocks. The production temperature of each doublet is depicted in Fig. 5.
The decline of production temperature becomes mild when the shale layers are included in the model, and therefore the lifetime of both doublets is significantly extended. This observation demonstrates that the shale facies should be sufficiently considered in the simulation and prediction of real geothermal projects. In addition, to what extent the shale facies will impact the lifetime of a geothermal project can hardly be quantified lacking forward numerical simulations with a detailed geological characterization. The presence of the shale layers in the simulation also allows the use of higher discharge rates that result in higher energy production for an equivalent system lifetime.
A considerable difference in the shape of the thermal plume is evident with and without shale layers (Fig. 6). The cold front propagation is more restricted due to the conductive thermal recharge when the shales are included. As the injected cold water transports through the sandstone layers, it is re-heated, extracting energy from the sandstone layers. As time evolves, a temperature gradient is built up between the sandstone bodies and the neighboring shale layers with the shales providing thermal recharge by heat conduction. The spatial intercalation of the sandstone and shale facies increases the contact area between them and amplifies the effect of the thermal recharge. The difference in the shape of the thermal plume becomes more pronounced with longer simulation time. The large impact of shale facies on thermal propagation also reveals the importance of geothermal modeling at sufficient geological resolution.

Overburden layer
To investigate the impact of overburden layer on heat production, ghost impermeable layers are placed on top of the reservoir to mimic the overburden. The underburden layers are contained in the geological model; therefore their influence is inherent and not considered as a factor here. No convective flow is assumed inside the overburden layer because of its low permeability. The overburden affects the thermal recharge in the reservoir only through heat conduction. For an accurate representation of overburden in the model, we used a detailed conformal grid resolving the scale of reservoir layers. The production temperature of each doublet in presence and absence of the overburden layer is shown in Fig. 7.
The lifetime of two doublets displays different sensitivity to inclusion of the overburden layer: the lifetime of doublet 1 is extended (Fig. 7a) while the lifetime of doublet 2 is not impacted at all (Fig. 7b). The reservoir heterogeneity plays an essential role in the extent that overburden layers can take effect. The highpermeable region in the top layers around the production well of doublet 1 (see Fig. 2) benefits from the thermal recharge of the overburden, dalaying the cold front. With cold front propagation, the temperature difference is built up between the overburden and permeable layers. Since heat conduction is proportional to the temperature gradient and contact area, the impact of heat conduction on heat propagation increases during the development process. The heat conduction of the overburden layer recharges the cold water and delays the thermal breakthrough on the top layers. Therefore, compared to the case without overburden, the lifetime is extended. This observation on the positive effect of overburden layers on heat production is similar to previous studies [17,18]. Doublet 2 mainly targets the Delft Sandstone (bottom reservoir layers). The production well of doublet 2 does not pass highly permeable layers on top of the reservoir, as is displayed in Fig. 2. The overburden layer is regarded as "non-conductive" in this case since no temperature gradient is actually built up between the overburden and its adjacent layers. Therefore, the thermal recharge from the overburden layer does not affect thermal breakthrough of doublet 2. Notice that the different conclusions on the impact of overburden layers on geothermal production can only be made based on detailed modeling of a geothermal field.

Well placement
Local changes in well placement can impact the flow path and sweep area of the doublet, consequently influencing doublet performance and lifetime. Shetty et al. [16] already examined the impact of the well location on the production temperature of a doublet based on a synthetic geological model and found a large variation in lifetime with such changes. As shown in Shetty et al. [16], the fluid streamlines and flow paths can change significantly, although the doublet was only re-located locally (less than 50 m) around the reference position. In this section, we follow a similar workflow and investigate the influence of the well placement on lifetime in the realistic geothermal reservoir.
The schematic for well placement is depicted in Fig. 8. Besides the reference locations, neighboring grid cells are also considered as the possible doublet locations, which corresponds to approximately 50 m variation in space. In each sensitivity study, both wells of one doublet are simultaneously moved in the lateral direction. For each simulation, only the position of one doublet is changed, while the other doublet remains at its initial location. The simulation results are displayed in Fig. 9 for each doublet.
As shown in Fig. 9, the well location can significantly influence       the doublet lifetime of the geothermal reservoir. For doublet 1 (Fig. 9a), the lower and upper limits of doublet lifetime are reached at the South (S) and the North-East (NE) locations respectively. The well location selected at the North-East can extend the doublet lifetime by around 10 years compared to the South position which accounts for around 20% variation in the lifetime. For doublet 2 (Fig. 9b), the lower and upper limits of doublet lifetime are achieved at the East (E) and the South-West (SW) locations respectively. Approximately 24 years of lifetime of the doublet can be gained when its location is changed from the worst to the best case (only by circa 50 m), which is more than a 25% variation. This confirms the conclusion of Shetty et al. [16], which was just based on an unconstrained conceptual geological model. The finding demonstrates the necessity to perform uncertainty analyses based on the existing geological knowledge leading to a more robust evaluation and reduction of the risk of earlier thermal breakthrough, while the heat recovery of geothermal fields can be robustly optimized.

Discharge rate
The development scheme directly impacts the lifetime and heat production of geothermal reservoirs. Besides, heat production can also be affected by doublet interference. In this study, both doublets are operated with flow rate control, which is the operation scheme usually applied in industrial practice. Here, only the injection and production rate of doublet 1 are tuned while the production temperature and energy production of both doublets are recorded to evaluate heat production from doublet 1 and its interference to doublet 2.
The production temperature of doublet 1 varies a lot with the rate changes (Fig. 10a). A shorter doublet lifetime is observed with the increase of discharge rate, which happenes because larger discharge rates accelerate the propagation of the cold front resulting in an earlier thermal breakthrough. The lifetime of doublet 1 is extended by 43 years with a 3 times reduction in injection and production rate. It is obvious that the lower rate gives slower propagation of the cold front and elongates the lifetime. However, Fig. 10b shows that the total energy production decreases with the rate reduction, although the lifetime is extended which is happening due to a reduction in swept region of the reservoir.
Notice that the input energy needed to maintain a higher production rate is also expected to be larger. The required pressure differences for these 3 production schemes are 22.4, 44.7 and 67.2 bars for production rates varying from 2500 m 3 /d to 7500 m 3 / d. A quick evaluation of the required pump energy can be expressed as the function of pressure difference (DP) and discharge rate (Q): where h denotes the pump efficiency. Detailed economical analysis (including factors related to NPV) is hence required before a final decision about the preferred discharge rate can be made [43e45]. The production temperature of doublet 2 (Fig. 11a) does not change with the variation in the rate of doublet 1, which demonstrates that production in doublet 2 is not impacted by doublet 1. The energy production in Fig. 11b also states that doublet 2 is not influenced by varying the production rate of doublet 1. This can be  explained by the fact that the two doublets target different parts of the reservoir in the neighbourhood of each other: doublet 1 targets the upper parts of the reservoir, while doublet 2 more focuses on the bottom reservoir layers, which can be seen by their different responses to overburden layers. This is an important observation for a typical sedimentary environment of low-enthalpy geothermal systems. The two doublets located next to each other can be operated independently.

Vertical-horizontal permeability ratio
The permeability ratio between vertical and horizontal directions varies with the geological and sedimentary process and can influence the direction of fluid flow among layers, which in turn influences the thermal propagation in the geothermal reservoir. In this section, the influence of vertical and horizontal permeability ratio to production temperature is investigated. Three permeability ratios are selected: k v /k h ¼ 0.01, 0.1 and 1 which indicates changes in vertical permeability while horizontal permeability remains fixed. The production temperature of both doublets is depicted in Fig. 12.
The lifetime of doublet 1 is slightly influenced by the permeability ratio (Fig. 12a). With the increase in permeability ratio, the lifetime is slightly extended, which means more heat is swept. In contrast, the inter-layer communication is reduced and the diversion of flow between layers is restricted for lower permeability ratios. Therefore, a lower permeability ratio will lead to earlier thermal breakthrough, and the total heat sweep efficiency of the reservoir decreases.
In Fig. 12b, the lifetime of doublet 2 seems insensitive to the permeability ratio. Due to the complexity of geological formations in the reservoir, this doublet may be drilled through multiple impermeable layers, which already separate the fluid and heat flow in the vertical direction as barriers. Therefore, even if the permeability ratio is adjusted for different simulation runs, the influence of permeability ratio to flow does not take effect.
In addition, the flow along the horizontal (x) direction is much larger than that along the vertical (z) direction for permeable sandstone layers (Fig. 13). When a smaller vertical-horizontal permeability ratio is selected, the connectivity in vertical direction will become worse, which will give rise to even weaker flow between inter-layers. Therefore, the variation of permeability ratio does not significantly impact the production temperature even when high-resolution geological realism is involved in the model.

Permeability-porosity correlations
The permeability-porosity correlation is another source of uncertainty in the geothermal reservoir. Here, we use two different permeability-porosity (k À 4) relations (Fig. 14) based on the correlations from TNO-AGE [46] for the typical sandstones of the WNB.
As displayed in Fig. 15a, the production temperature of doublet 1 is only insignificantly affected by the change of correlation. As doublet 1 mainly targets the top Berkel Sandstone, the preferential channels inside the control area of doublet 1 are associated with high porosity (Fig. 2). Two k À 4 correlations generate similar permeability distribution in the high porosity values which is quite typical for sandstones. Therefore, the convective flow of doublet 1 only slightly varies between the two scenarios using different correlations.
However, a more pronounced difference is observed in the production temperature of doublet 2 (in Fig. 15b), where the lifetime is noticeably extended with correlation 1. This can be explained by the existence of an evenly distributed permeability  field within the control area of doublet 2. With correlation 1, the number of flow paths in doublet 2 is increased with an improved k À 4 correlation (compared to correlation 2). In this scenario, the injected cold water sweeps a larger area of the reservoir and the thermal breakthrough is delayed. This difference can only be observed with modeling at high-resolution geological accuracy.

Discussion and conclusions
In this study, an extensive investigation is performed based on a realistic heterogeneous reservoir model at geological-scale resolution utilizing the open-source Delft Advanced Research Terra Simulator (DARTS). Detailed sensitivity analysis in this work shows the importance of full uncertainty quantification for meaningful optimization of geothermal development.
C Shale layers should always be considered when modeling fluid flow in geothermal reservoirs. The energy stored in shales can be exploited through heat conduction with the neighboring high-permeable layers. Besides, the intersecting distribution of shales and sandstones is also essential for the efficient exploitation of energy stored in shale layers. The inclusion of the shale layers in the simulation can have one of two results: extend the predicted system lifetime due to thermal recharge (more than 70% extension for the lifetime of one of the doublets in this study), or enable the use of higher rate controls and therefore higher energy production for a similar system lifetime. C Overburden can recharge its adjacent high-permeable layers by heat conduction. As observed in this study, the influence of overburden layers to doublet lifetime is different at two doublets locations inside the same geological formation (i.e., more than 50% extension for the lifetime for one doublet while almost no extension for another). Due to the limited working distance, heat conduction only impacts several layers close to it. Besides, the influence of the overburden to production temperature closely correlates to the contribution of adjacent layers to the total fluid production (i.e., the reservoir heterogeneity). This finding attributes to the realistic high-resolution geological representation, which is hard to achieve in conceptual models unconstrained to real observations. C Due to the reservoir heterogeneity, the flow path varies significantly with local well positioning. In the studied realistic reservoir model, a synchronous minor location shift of both wells by approximately 50 m horizontally leads to a variation of lifetime of more than 25%. This indicates the   necessity to perform an uncertainty analysis of well placement based on detailed geological characterization. A striking extension in the lifetime of geothermal projects may be expected after the optimization of well positions in the presence of uncertainties. C Higher discharge rates yield faster thermal breakthrough and shorten the doublet lifetime. On the other hand, the cumulative energy production increases with the production rate, which indicates that a longer lifetime does not mean better heat recovery. Larger rates can extract more heat to the surface in a shorter time, which is attractive for industrial operators. Certainly, the energy investment to achieve a larger rate is higher as well. The optimized development scheme should take into account energy cost, financial investment and operational scheme comprehensively to maximize the NPV. C The interference between well doublets in realistic geothermal production scenarios is not necessarily directly related to their relative position. The performance of doublet 2 in our study is not influenced by the rate variation in doublet 1, reflecting the weak interference between these two doublets in the range of selected operation strategies. These findings support current technological practice with doublet-based development of low-enthalpy geothermal reservoirs. However, only a detailed uncertainty quantification can directly evaluate the risk of reduction in the energy production due to geothermal doublets interference. C The influence of vertical-horizontal permeability ratio on the temperature profile is relatively minor in the studied scenario for both doublets, which is related to the property distribution of the reservoir. However, the permeability ratio is still a factor worth checking for the development of a geothermal reservoir because of its influence to the propagation of cold front and heat sweep efficiency. C The in-situ porosity distribution is essential for the k À 4 correlation and has a noticeable effect on the doublet lifetime. We demonstrate that the influence of k À 4 correlation to doublet lifetime varies within the same model. If a preferential channelized porosity distribution exists between injection and production well, the variation of k À 4 correlation will not largely impact the heat production with fixed rate control. It is important to mention that the given correlations are based on the same facies distribution. The influence of k À 4 correlation can vary to different extent with different facies interpretation in the geological model. C Thanks to the high performance of DARTS framework, the set of simulation runs was completed in a limited computational time. Besides, the flexibility of DARTS makes it convenient to adjust different parameters in the sensitivity study, such as activation and deactivation of shale layers and the overburden, based on the realistic data set. This provides a generic and robust framework for the sensitivity analysis of practical geothermal projects.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.