Performance of low-enthalpy geothermal systems: Interplay of spatially correlated heterogeneity and well-doublet spacings

(cid:129) 2600 computationally intensive simulations for geothermal doublet systems conducted. (cid:129) Discharge, well/doublet spacing, poro-perm correlation lengths and variance varied. (cid:129) A doublet spacing equal to well spacing produced consistently best performance. (cid:129) Anisotropic heterogeneity led to shorter/longer lifetime for short/long spacing. (cid:129) Su ﬃ cient lifetime for shorter well spacing than the ones conventionally designed.


Introduction
Many of the low-enthalpy deep geothermal systems are deployed in sedimentary reservoirs at depths between 2 and 2.5 km with a temperature between 70 and 90 °C [1].The most common method of geothermal energy recovery from low-enthalpy aquifers are doublet systems that utilize two wells, one for hot water production and another for cold water injection.The lifetime of the doublet (how long the doublet can produce economically), energy sweep (produced energy compared to the total amount of available energy) and energy production rate of doublets determine the performance of doublet systems.The accuracy of predictive simulation tools is essential for the successful design of doublet systems.
The sedimentary reservoirs are characterised by their lithographical, geological, structural and thermal properties.These characteristics govern the geothermal performance indicators.Various studies have focused on influence of some of the above-mentioned characteristics on doublet performances.For example, Poulsen et al. [2] investigated the impact of thermal conductivity.Several authors focused on brine viscosity and density dependence on temperature, e.g., Ma & Zheng [3] found that mean discrepancy between the simulated temperature distributions with and without considering the effects of variable density and viscosity is approximately 2.5%.Using a numerical study on density and viscosity variations with temperature, Saeid et al. [4] found that ignoring the variations leads to overestimation of the geothermal system lifetime in hot injection scenarios, and underestimation of the system lifetime in cold injection scenarios.Mottaghy et al. [5] showed that thermal and hydrogeological data are crucial to planning geothermal resource development, and numerical codes should take temperature dependence of thermal properties into account.Vogt et al. [6] demonstrated the importance of accounting for heterogeneity of rock parameters resulting in significant variations of production temperature and well pressure with time.
Several authors including [7][8][9][10] studied impacts of well spacing.Saeid et al. [7] showed that the lifetime increases linearly with the well spacing for homogeneous aquifers with no geological complexity or barriers.Willems et al. [8], also for homogeneous aquifers only, optimized well spacing so that the interference between the aquifers are minimized.Pandy & Vishal [9] showed that, for homogeneous aquifers again, at higher well spacing, the flow length/volume increased and so did the pumping power, leading to improved overall performance in heat extraction.Willems et al. [10], for homogeneous aquifers, evaluated both the possible financial advantage of well spacing reduction and its impact on doublet life time.
A recent review of Pandey et al. [11] on geothermal reservoirs coupled thermo-hydro-mechanical-chemical approaches shows that the impact of aquifer heterogeneity on performance of doublets is less studied.The presence of spatial heterogeneity inside a reservoir may induce flow channeling and reduce the volume of reservoir participating in flow fields [11].Most of the studies on performance of low enthalpy systems, as referred above, have considered homogeneous systems or simple lithographical variations.There are only few existing research works that specifically deal with heterogeneity for geothermal doublet systems [12,6,[13][14][15][16].Watanabe et al. [12] found that the most significant factors in the analysis of thermo-hydro-mechanical coupled processes in heterogeneous porous media are permeability and heat capacity.Vogt et al. [6] studied the transient temperature and pressure at the production well for 400 sets of heterogeneous realizations of the fault zone in a doublet system.The authors concluded that the distribution of porosity/permeability and thermal conductivity inside the fault zone significantly impact heat extraction rate from the reservoir.Crooijmans et al. [13] studied the impact of heterogeneity in a fluvial system for a single doublet and they have suggested a correlation for the lifetime as a function of production rate and net-to-gross values.The authors showed that at lower net-to-gross ratio the temperature drop at the production well was slower than the higher net-to-gross ratio for all sets of the heterogeneous reservoirs.Willems et al. [14] studied the impact of heterogeneity on the connectivity in geothermal systems.The authors found that the impact of heterogeneity is significant on heat production and pumping loss was less if the wells were placed parallel to the paleo flow direction.Niederau et al. [15] studied the effect of spatially heterogeneous permeability on the formation and shape of hydrothermal porous flow in the Yarragadee aquifer, Australia.Their results showed the spatially heterogeneous permeability can affect the local convection patterns.Very recently, Liu et al. [16], showed that for a doublet system with increasing correlation length, the possibility of flow channels appearing in well pair system increases, causing a short average thermal breakthrough time and a lower surface settlement around the injection well.
For heterogeneous aquifers, attention has been also mostly paid to geothermal doublet systems with fracture networks.Performance of deep geothermal doublets in fractured reservoirs have been studied either for a single fracture [17], parallel fractures [18] or complex fracture networks [19][20][21].Salimzadeh et al. [20] showed that both heat production and required pump energy is very sensitive to fracture spacing, density and connectivity.Pandey & Chaudhuri [22] studied the impact of fracture aperture heterogeneity on geothermal heat recovery in fractured aquifers and found that small correlation lengths in fracture sets do not create much variation in temperature at production well while [23] illustrated the aperture variation induced by thermal or chemical stresses can significantly influence the performance of the geothermal system.Vasilyeva et al. [19] developed multiscale model of heat and flow transport in EGS systems with varying degree of fracturisation.No systematic works, however, have been carried out on geothermal doublet systems to delineate the impact of heterogeneity on well or doublet spacing.
An important design factor for geothermal doublets is the required distance between wells in a doublet system as well as the distance between the multiple doublets.There are research studies dealing with the optimisation of low-enthalpy doublets [24,25].Chen et al. [24] coupled a complex hydrothermal simulation model and a multivariate adaptive regression spline-based surrogate model to investigate the Fig. 1.The doublets configuration and the license-boundary control region of the subsurface system used in this study.effects of geological uncertainties (fault size and geological unit permeability) on optimal well placement and control (re-injection well location, production rate) in a geothermal prospect near Superstition Mountain in Southern California, USA.The model used was homogeneous, so that in each case a constant permeability was chosen in the range of − 13.8 < log(perm) [m 2 ] < − 13.3.The authors found the optimal production rate is 30.7 kg/s and distance between production and injection well is 473 m in order to maximize the net profit after 50 years of potential geothermal extraction.Kong et al. [25] used numerical modelling, economic analysis and a homogeneous domain of a synthetic aquifer and obtained an optimal well placement of 400 m.
There are no existing work specifically focusing on well/doublet spacing for heterogeneous geothermal systems.The area of research is important because administering a sub-optimal well/doublet spacing can potentially lead to negative interferences that, in turn, may influence the utilisation efficiency of wells and doublet systems.An accurate modelling-based design that provides optimal solutions will prevent such interferences.Currently, thermal breakthrough which is the moment when the extent of re-injected cold water plume reaches the production wells is the basis to determine production lifetimes as an indicator for the temperature drop at the license boundary.However, the license areas temperature may not immediately drop to non-economic values when thermal breakthrough occurs at production wells.As a result, thermal energy is still available in economic levels and heat production could continue after thermal breakthrough as long as the average temperature at the license boundary has not experienced notable temperature drop (i.e.>1 °C).The temperature drop over the li- cense area's boundary can instead be used to determine the lifetime of the geothermal doublet systems.This can obtained utilising accurate and realistic production simulations.The temperature distribution and shape of cold water front are strongly controlled by underlying heterogeneities in the geological properties of geothermal aquifers.
To address above-mentioned knowledge gap in the design of geothermal system, in this study, we address the following questions.
1. How does well spacing or doublet spacing under different levels of heterogeneity impact performance of the geothermal system?The performance criteria include lifetime of the operation, energy or thermal sweep efficiency, energy production, and coefficients of performance.2. What is the correlation between temperature drop in the production wells and that of license boundary?3. How significant is the impact of heterogeneity on well spacing design?
To address these questions, we use a numerical simulation tool described in Section 2 alongside the geological properties, heterogeneous property generation, heat recovery scenarios as functions of well or doublet spacing and definition of performance metrics.The results and discussions are presented in Section 3 and the paper is concluded in Section 4. By answering the above-mentioned research questions, we provide a novel unprecedented understanding of interplay between "realisticallyrepresented fluvial aquifer heterogeneity and geothermal doublet  design".This is of direct relevance to subsurface energy application and geothermal heat recovery.

Methodology
The numerical simulations are carried out by ECLIPSE E300 simulator [26].Previously for geothermal heat recovery simulations (including low-enthalpy doublets or deep enhanced geothermal systems), the simulator have been successfully employed [27].Here, a multiparametric analysis is carried out on the simulation results to delineate the interplay of heterogeneity and well or doublet spacing.In the following, first we describe the methodology including a brief description of the governing equations for coupled heat transfer and flow in porous media, the geological model considered with heterogeneities of porosity and permeability fields, the simulation scenarios and the performance metrics.

Coupled heat transfer and flow in porous media
To simulate coupled heat transfer and flow in porous media, a fully implicit finite volume method utilising a two-point flux approximation scheme is employed.The flow and heat transfer solves a coupled problem integrating (i) the conservation equation for each fluid component (water) in each gridblock at each timestep (leading to the non-linear residual r fl ), and (ii) the energy conservation equation in each gridblock at each timestep (leading to the non-linear residual r e ): where m is the moles of fluid component (water) in each gridblock, V p is the pore volume, F is the net flow rate into neighbouring grid blocks, Q , where r in which ϕ is the porosity of gridblock, ρ is density of water/rock and C is the specific heat capacity of water/ rock.In Eq. ( 2), F e is the convective enthalpy flow rate into neigh- bouring gridblocks, C e is the conductive energy flow rate into neigh- bouring gridblocks, Q HL is the conductive energy flow rate to the sur- rounding rocks (heat loss), Q e is the net enthalpy flow rate into wells during the timestep [26].The net flow of water from cell i into neighbouring cells is obtained by where Γ i is the transmissibility between cells n and i which is con- structed using the harmonic mean of absolute permeability (K) of cells n and i b , w is the molar density of water, μ w is the viscosity of water, and dΦ ni is the potential difference of water phase between cells n and i.The net flow rate of energy (convection) from cell i into neighbouring cells is obtained in a similar manner: The heat conduction term for cell i is given by summing conduction between all neighbouring cells n as [26]: where Γ hi is the thermal conduction transmissibility between cells n and i constructed using the harmonic mean of the equivalent conductivity of cells n and i.The equivalent conductivity is eq w r , where λ w is the water conductivity and λ r is the rock conductivity.
In this study we have not included the heat transfer in the production and injection wells.This has been studied in detail by Saeid et al. [7] for two tubing materials for a discharge rate of 150 m 3 /h.They concluded that the heat gain and loss in the injection and production wells, respectively, after five days of operation are negligible.

Geological model and fluid properties
We assume a structurally simple, synthetic, 3D rectangle for the geological system under study with 3 km × 3 km × 500 m lengths in x, y and z directions (Fig. 1).The system comprises of an overburden (200 m thick), an aquifer (100 m thick) and an underburden (200 m thick).The system is discretized into laterally uniform 120 × 120 mesh (25 m × 25 m gridblocks).However, vertical discretization is non-uni- form.The overburden and underburden are each discretized into one layer only with 200 m thickness.Whereas the aquifer is discretized into 10 layers (each layer 10 m thick).Therefore overall we have 12 layers and 172,800 gridblocks.The top face of the structure is 2250 m deep. 1  The initial conditions are set as p init (initial pressure) equals to 200 bar at the reference depth of 2500 m, and T init (initial temperature) equals 67.5 °C at 2250 m, 75 °C at 2500 and 82.5 °C at 2750 m.The boundary conditions are fully close to flow in x-and z-directions and fully open to flow in y-direction (the direction of anisotropic heterogeneity).In order to implement open boundaries, the pore volumes of the boundary blocks in y-direction are multiplied by 1,000.This method has been previously used to emulate open flow boundary conditions or infinite domains [28][29][30].

Rock properties
Thermal conductivity of rock is 0.91 W/m/K, rock density is 2650 kg/m 3 , rock specific heat capacity is 2,000 J/kg/K and rock compressibility is × − 4.93 10 5 bar −1 at 250 bar.In this study, we neglect the dependence of rock density and heat capacity on porosity as considered by Liu et al. [16].Porosity of the overburden and underburden is set to 0.01, and permeability of the overburden and underburden is set to 0.001 mD.To generate porosity of the aquifer, a code available online [31] is used to generate correlated fields.The method utilizes the Fourier-transform of the covariance function as the power spectral density function of all realizations.Random autocorrelated fields are generated by creating random phase spectra meeting the conditions of real numbers in the physical domain.The realizations are then converted by back-transformation of the power-and phase-spectrum into the physical domain.
Correlated porosity fields are generated all with a mean of 0.17 and a variance of = σ 0.02 2 and 0.04, hereafter referred to as v1 and v2 for simplicity.The covariance matrices of porosity fields are assigned as exponential.Four different correlation lengths in lateral directions (c x denotes correlation length in x-direction and c y denotes correlation length in y-direction) are used in the process: 1 Our preliminary modelling practices refining the mesh did not yield significant change in the results of a homogeneous and several heterogeneous examples we considered for mesh sensitivity analyses.
• c x = 100 m and c y = 100 m hereafter referred to as C1 • c x = 100 m and c y = 500 m hereafter referred to as C2 • c x = 250 m and c y = 500 m hereafter referred to as C3 • c x = 500 m and c y = 500 m hereafter referred to as C4 where c x = 100 m means 4 gridblocks are correlated in x-direction.As such we generate eight fields (2 variances multiplied by 4 different correlation lengths).Eight realizations are generated for each field, totalling 64 porosity fields.Next, the permeability fields are derived from porosities through the following relationship for the fluvial Delft Sandstone aquifers in the Netherlands [32]: Figs. 2 and 3 show the profiles of porosity and permeability for one realization of each 8 pairs of correlation length and variance used to generate porosity and permeability fields, respectively.

Fluid (water) properties
Water thermal conductivity is set to 0.59 W/m/K, water specific heat is 4200 J/kg/K, water viscosity is written as a function of temperature according to IAPWS [33], and water density is input as a function of pressure and temperature through: where is the water density at reference pressure and is set to 1000 kg/m 3 , and c w is water compressibility and is set to × − 4.0 10 5 bar −1 .

Heat recovery scenarios
The geothermal heat recovery operation is carried out by two doublets using a checkboard pattern [8].The injection wells are hereafter referred to as I1 and I2, and the production wells are referred to as P1 and P2 (see Fig. 1).The distance between the injection and production wells in each doublet is L, and the distance between the two doublets is dx.The position of I1 is fixed for all of the heat recovery scenarios and it is at = x 750 m and = y 750 m.However, the positions of I2, P1 and P2 change with different L and dx values.We assign four values for well spacing, L = 500 m, 700 m, 900 m, and 1100 m.For each L, five doublet spacing, dx, are considered: Cold water is injected using two well injection rates of = Q 150 m 3 /h and 250 m 3 /h, hereafter referred to as Q1 and Q2 for simplicity.The injection rates are selected based on the actual production rates applied for the deep geothermal doublets in sedimentary aquifers [34].The maximum allowable pressure at each well is 260 bar to avoid inducing hydraulic fractures.Above this threshold, the injection rate is decreased and as a result the production rate is also decreased to equate with the total injection rate.However, assuming that for each well stimulation operations have been carried out, porosity and permeability of connection blocks at perforation depths are increased so that they are equal to the average of porosity and permeability of the perforation column.By means of this, the injection and production rates are maintained throughout the simulation time for all the simulations in this study.
The temperature of injected water is 30 °C.The operation is continued for 50 years.A robust set of timestepping criteria is employed by the simulator to ensure convergence.The detailed presentation of the timetepping strategy in thermal mode of the simulator is prohibitive to describe here and can be found in ECLIPSE Technical Description manual [26] under Thermal Features/Timestepping Criteria.Using this strategy and checking the production profiles for highest injection rate and most heterogeneous porosity/permeability, we checked that there are no significant convergence issues in our simulations.
Using 64 models of porosity and permeability fields, overall we conduct 2,560 simulations: 2 variances multiplied by 4 correlation lengths multiplied by 8 realizations for each pair of variance-correlation length multiplied by 4 L values multiplied by 5 dx values, multiplied by 2 injection rates.These parameters are summarised in Table 1.Each simulation in average takes 3,800 s to complete using a Windows Server 2012-operated HP ProLiant server with two Intel Xeon E5-2690 12-core 2.60 GHz CPU processors.A parallel code is developed in MATLAB to call 8 simulations at the same time.This parallelisation saves run time significantly.In order to make comparisons of heterogeneous systems with homogeneous systems, additional simulations are carried out with constant aquifer porosity and permeability of ϕ = 0.17

Definition of performance metrics
In order to investigate the interplay of well or doublet spacing and heterogeneity of the system, first the following parameters are defined: where q W k prd t , , is the water production rate at connection (layer) k of well W, at time t, and T ijk t is the block ijk temperature at time t.For each production well, obviously ij are fixed and only k changes from 1 to 10 (uppermost layer of the aquifer to lowermost layer).The denominator of above equation is equal to Q prd W , which is a constant.3. Converged solution for a variable α over ℓ realizations is defined as . By looking into the plot of α ℓ vs. ℓ, the con- vergence of realizations can be assessed, that is, α ℓ should converge for a particular ℓ so that for > i α ℓ, ℓ remains constant.Based on our simulations, the converged solution for 〈 〉 T W t is obtained through 6 realizations, proving that 8 realizations are sufficient for convergence.4. Life time based on production wells: LTW P1 and LTW P2 of geothermal heat recovery operation defined as the time (in years) when average 8. Coefficient of Performance defined as: where E prod is the energy produced from production wells and E pump is an estimation of pump energy losses [10], in which ε is the pump efficiency of %60, and < − > p p inj prd is calculated by using the converged solution for well bottom-hole pressures of two doublets.9. Energy Sweep is an indicator of how efficiently heat in the control volume is extracted.This is defined as: where E R is the available reservoir energy of the license area, V bi is the bulk volume of gridblock i in the license area ( × × L dx 2 2 10layers around the doublets), N b is the number of gridblocks in the license area, and E prod is the geothermal energy recovered by the doublets until the converged solution based-lifetime of the operation.

Lifetime analysis: boundary vs. production wells-based comparison
Fig. 4(a) shows LTP for various correlation lengths (C1, C2, C3 and C4), injection rates (Q1 and Q2) and variances in porosity field (v1 nd v2).The profiles are plotted versus L and for = dx L only.As a result there are 16 profiles in this figure.In this figure and the rest of figures in this manuscript, black represents Q1-v1 simulation results, green Q1-v2, red Q2-v1 and blue Q2-v2.Yellow lines represent homogeneous simulations.Some of the profiles are incomplete with respect to L. This is because the average temperature has not dropped 1 °C during 50 years of operation.From Fig. 4(a) it is clear that: • Increase in injection (Q) and variance of porosity fields (v) lead to decrease in lifetime.This is trivial as cold water front reaches to production wells earlier for higher Q and v.
• Between each fixed pairs of Q-v, there are 4 profiles for each cor-  lowest lifetimes.This point will be further analyzed and discussed in Section 3.3.
• Considering a minimum lifetime of 20 years and an average tem- perature drop of 1°C at the license boundary as the main indicators of the lifetime, the minimum well distances for Q1 and Q2 (regardless of the heterogeneity variance) are 700 m and 900 m, respectively.
• Considering a minimum lifetime of 30 years and an average tem- perature drop of 1°C at the license boundary as the main indicators of the lifetime, the minimum well distances for Q1-H, Q1-v1, Q1-v2, Q2-H, Q1-v1, and Q1-v2 (regardless of the correlation lengths) are about 650 m, 750 m, 800 m, 900 m, 950 m, and 1000 m, respectively.In order to demonstrate the impact of well or doublet spacing on lifetimes (boundary vs. production wells-based calculations), in Fig. 5 all realizations as well as converged solutions are grouped into fixed dx L / (varying L) or fixed L (varying dx).Starting from Fig. 5(a), for all realizations, LTB are plotted vs. LTW P1 .In each plot, injection rate (Q), variance (v) and L vary.The scattered markers above the 45°line, show that LTB >LTW P1 .However, there is no dominant pattern in results in order to make a conclusion that LTB is higher than LTW P1 .Moreover between different dx L / 's in Fig. 5(a), the extent of lifetimes in x-and y- axes are the same.Therefore, the variation in doublet spacing does not have a clear impact on lifetimes.Conversely, Fig. 5(c) shows that by increasing well spacing (L) the lifetimes increase and scatter markers shift towards right of the plot.
Comparing LTB with LTW P1 , the scatter markers are mostly above the 45°line, but not all of the realizations produced this trend.For the converged solutions (Fig. 5b and Fig. 5d), there is a general but not a universal trend of LTB > LTP for dx L / = 0.7 and 0.85, L = 900 m and L = 1100 m.Nevertheless, for v2 (green and blue markers), we can deduct that LTB > LTP.This is because for v2 the breakthrough of cold water is happening too early (because of heterogeneity of the system) while the licensed region still has some energy left.We discuss this in details in the next subsection.

Determining the license area's extent from temperature distribution
If LTB > LTP, the geothermal system has still energy to produce by the time that the production wells signal the temperature drop.This means that the license boundary chosen has not experienced the same level of temperature drop by the time that the production wells have.Here we show that the license boundary can be in fact calculated in a way that the boundaries experience the threshold temperature drop of 1 °C.We show that these boundaries are not necessarily the same as × dx L 2 2 which were conventionally assumed.We calculate the distribution of the actual boundaries of license area around operating wells based on a simple search algorithm.To this end, for the case of = L dx (which was shown to consistently result in op- timal well/doublet spacing solutions), and = L 500 m, 700 m and 900 m, Q1 and Q2, fixed v2 (with a general trend of LTB > LTP), and C1, C2 and C4, first we calculate LTP1, LTP5 and LTP10, where LTPn represents production lifetimes when the converged solution for production wells (〈 〉 T P t ) experiences n °C temperature drop.Then, for each realization and at LTPn, we search for a square license area around the wells so that for its boundary blocks, the average temperature drop is 1 °C.The length of this square is denoted by ′ L (= ′ dx ).Also similarly we find a boundary where the cold plume has farthest spread.That is we search for a farthest boundary away from the wells where the difference between the minimum temperature of that boundary and the initial temperature of that boundary is 1 °C.This boundary represents the extent of cold front.We denote this boundary with ″ L (= ″ dx ).Fig. 6 shows the results of this analysis.For each LTPn (LTP1, LTP5 and LTP10), we show ′ L L / and ″ L L / for C1 (averaged over all realiza- tions of C1), for C2 (averaged over all realizations of C2), and for C4 (averaged over all realizations of C4).In each subfigure there are four lines for Q1-v2 (for ′ L L / ), Q1-v2 (for ″ L L / ), Q2-v2 (for ′ L L / ) and Q2-v2 (for ″ L L / ).The following observations can be made from this figure: 1.For L < 1000 m and the well temperature drop of <10 °C, both ′ L L / and ″ L L / can be as large as 1.25 and 1.5, respectively.2. Impact of discharge (Q) on the size of the license area is not as important as how the constraint is considered for the temperature drop at the license boundary, i.e., an average 1 °C or a local 1 °C temperature.Note that obviously for a given case ″ L is always larger than ′ L .

For C4 cases ′
L and ″ L decrease by increasing the well distance (L) implying that for > L than 1000 m the license boundary can be chosen as 2 ) or slightly larger.For both C1 and C2 cases ′ L and ″ L remain more or less similar for different L values between 500 m and 900 m. 4. In most cases for the systems with a lifetime defined as 10 °C (or L is almost equal or smaller than L for a well distance of 700 m and 900 m for Q1 and Q2, respectively.These increase by 20% for the lifetime of 40 years.
Fig. 8 depicts the temperature field and the license boundaries defined based on different constraints for realization 1 of C1 and C2, and for different Q and L values.Fig. 9 shows similar results for the converged solutions.
The animations for the development of cold front for the examples shown in Fig. 8 and Fig. 9 are provided in Videos 1 to 8. The videos are zoomed in around the license areas for better visualization purposes and the stills are taken at the production time equal to 10 years.The videos represent how the swept area advances in time for different well distances and heterogeneous fields.Comparing the results of geothermal systems with L = 500 m for different Q suggests that while a minimum well distance of L = 500 m might be enough to provide a lifetime of greater than 20 years for a discharge of 150 m 3 /hr, this well distance is not proper for higher discharge values.For a discharge of 250 m 3 /hr the results suggest that a well distance of 700 m provides a lifetime of more than 25 and 35 years for LTP5 and LTP10, respectively.

Effects of correlated heterogeneity
To investigate the impact of correlation lengths on operation lifetimes, Fig. 10 shows the difference in LTP between C4 and C1 in Fig. 10  years.Such effects of correlated heterogeneity (or channelised heterogeneity) has also been studied in details in a recent publication by Lie et al. [16].
Next, we focus on 〈 〉 T W t of Sample 2 and Sample 4 (high dx but varying correlation lengths of C1 and C2, respectively corresponding to Fig. 11row b) and Fig. 11row d), and green and blue lines in Fig. 12(a)).In this instance, P2 is experiencing similar temperature drops between two cases.This is because P2 is sufficiently away from I1, so that correlated heterogeneity along I1-P2 corridor cannot lead to quick advection of cold water to P2.Nonetheless, the difference in lifetime between Sample 2 and Sample 4 is made by P1.For Sample 4, the directional correlated heterogeneity is distorting the cold water fronts of both I1 and I2 to the benefit of P1.That is, large parts of cold water fronts of I1 and I2 move alongside y-direction but in the opposite direction towards P1, so that breakthrough at P1 takes place later for Sample 4 compared to Sample 2. As a result of these complications, LTP of Sample 4 is higher than Sample 2, and this corroborates results shown in Fig. 10 The comparison is the same as comparison made between Sample 1 and Sample 3, with Sample 5 behaving similar to Sample 1 when compared to Sample 3. Due to isotropic correlated heterogeneity and proximity of doublets, the rather circular cold front did not lead to an early breakthrough of cold water front as of Sample 3. Consequently, LTP of Sample 5 is higher than Sample 3, and this corroborates results shown years.Finally, we focus on 〈 〉 T W t of Sample 4 and Sample 6 (high dx but varying correlation lengths of C2 and C4, respectively corresponding to Fig. 11row d) and Fig. 11row f), and blue and cyan lines in Fig. 12b).In contrast to comparison between Sample 4 and Sample 2, P2 (and not P1) makes the difference in determining lifetime.Anisotropic correlated heterogeneity of Sample 4 accounts for preventing the cold water front of I2 to reach P2.Consequently, LTP of Sample 4 is higher than Sample 6, and this corroborates results shown in Fig. 10(c years.The fact that either P1 Sample 2 or P2 in Sample 6 experience earlier breakthroughs than P1 or P2 in Sample 4 is random, however, it is due to isotropicity of correlated heterogeneity of Sample 2 and Sample 6 compared to Sample 4. Further to above-mentioned comparisons between different correlation length samples, the comparison of C4 with C1, and C4 with C3 is evident from discussions above and supported by Fig. 10(a) and Fig. 10(d for large dx values.While the six samples above were taken from one realization only, it is evident that the converged solution has produced the same ranking between different correlation length cases with respect to production wells-based lifetimes (LTP).
In contrast to all variations in the production wells-based lifetimes observed as a result of change in correlation lengths, Fig. 13 shows that license boundary-based lifetimes (LTB) are significantly less sensitive towards this parameter and as such can be used to reduce uncertainty.This finding is corroborated by examining LTB of six samples investigated above and shown in Fig. 11 (third column).The magnitude of difference between LTB's are less than the magnitude of difference between LTP's, for most of the comparisons.For example, (LTP -LTB Sample 4 ) = 8 y.Therefore the selected examples, clearly demonstrate that LTB is significantly less sensitive towards correlation lengths and heterogeneity of geothermal system.It is evident in Fig. 13 that the impact of the correlation length on the LTB is less than 5 years.

Coefficient of performance and energy sweep
Figs. 14(a) and 14(b) show the profiles of Coefficient of Performance (CoP) and Energy Sweep (S) calculated from Eq. ( 9) and (10) as functions of well and doublet spacing for the converged solution at LTP.To avoid crowded figures, only C1 and C2 correlation-length cases are shown.
Fig. 14(a and b) shows that 1.For each fixed well spacing, while for heterogeneous cases the optimal doublet spacing is at dx L / = 1 (in the middle), for homo- geneous cases the doublet spacing has negligible impact on CoP. 2. The highest CoP's are obtained for Q1-v2 at dx L / = 1, however, considering the profiles for other dx L / values, Q1-v1 has con- sistently led to the best CoP's.3. One clear observation is the dramatic reduction of CoP due to increase in injection rate (Q).This has also happened for homogeneous cases (compare Q1-H with Q2-H).4. Similar to the homogeneous reservoirs, for most heterogeneous cases CoP decreases as L increases for dx L / = 1. 5. The increase in correlation length mostly decreases CoP.
From Fig. 14(c and d) for S, the following observations can be made: 1. Heterogeneity in general reduces S especially for lower L values.
2. Except for the highest L = 1100 m, increasing the variance of porosity fields (v) reduces S (green and blue lines compared to black and red lines).3. Increasing the injection rate (Q) increases S. 4. It is evident that dx L / = 1 is consistently producing the maximum values for S. 5.The increase in correlation length mostly decreases S. 6.Similar to the homogeneous reservoirs, for most heterogeneous cases S decreases as L increases for dx L / = 1.
In summary, the increase in correlation length in y-direction (C2 compared to C1) will largely lead to decrease in both CoP and S, with few exceptions when > dx L / 1.The increase in variance of the porosity fields (green lines compared to black lines, and blue lines compared to red lines), largely decreases S, but the impact on CoP varies for each dx L / .These results suggest that the shorter well distance may provide higher S and CoP.While higher injection rate increases S, it has a negative impact on the CoP.
Fig. 15(a-d) show the profiles of E prod and E pump as functions of well and doublet spacing for the converged solution and the homogeneous simulations.The results show that both of these properties increase with increase in well and doublet spacings due to increasing the lifetime.However, just at dx L / = 1, the ratio of E prod over E pump produces a maximum for heterogeneous models.Also with increase in injection rate, E prod and E pump increase.However, increase in E pump is larger than increase in E prod , so that CoP decreases by increase in Q, while since E R is independent of Q, S increases with Q.The changes in CoP and S as a function of Q is illustrated in Figs.15(e) and 15(f).

Conclusions and future works
A comprehensive set of numerical simulations were carried out on synthetic homogeneous and heterogeneous low-enthalpy aquifer system, with a range of operational and physical parameters of the subsurface system.Using multiple performance metrics, the following conclusions were made: 1. Heterogeneity undermines performance of geothermal systems in terms of Life Time of the operation, 2. Spatially highly correlated heterogeneity undermines performance of geothermal systems compared to low correlated systems or randomly distributed heterogeneity, 3.An anisotropically correlated heterogeneous system performs worse than isotropically correlated heterogeneous systems for large well and doublet spacings, 4. Increase in operational flow rate can increase Energy Sweep while decrease Coefficient of Performance, 5. A doublet spacing equal to well spacing robustly leads to best performance of the operation, 6. Use of license area's boundary showed less uncertainty with respect to various operational and physical parameters, suggesting a better control criterion for the operation design.7. The difference between lifetimes based on production wells and license area's boundary is minimum for when well and doublet spacings are equal.8. Sufficient lifetime (>20 years) could be achieved for the well spacing of less than 900 m for a discharge of 250 m 3 /hr.The lifetime can be increased significantly or the well spacing can be further reduced if larger temperature drop (>1 °C) at the producers are permitted.This would however require a larger license area than the × dx L 2 2 in order to minimise the negative interference with the neighbouring geothermal projects.
Future directions include use of robust optimisation algorithms to obtain optimal well or doublet spacing for heterogeneous geothermal systems.While previously, Kong et al. [25] conducted optimization of geothermal systems using homogeneous models, the heterogeneity will have non-trivial effects on positions of wells and doublets, and consequently it will also impact optimal solutions.The well spacings can be optimised further with including the economics of the project.

Fig. 5 .
Fig. 5.The comparison between LTB and LTW from all realizations, and between LTP and LTB calculated from converged solutions: (a) LTB vs. LTW P1 per realization grouped for each dx L / (L varies in each plot), (b) LTB vs. LTP per realization grouped for each dx L / (L varies in each plot), (c) LTB vs. LTW P1 from converged solution grouped for each L (dx varies in each plot), and (d) LTB vs. LTP grouped from converged solution for each L (dx L / varies in each plot).

Fig. 6 .
Fig. 6.(a) ′ L L / and ″ L L / vs. L for LTP1, (b) for LTP5 and (c) for LTP10.From left to right: ′ L L / and ″ L L / vs. L calculated from converged solutions of C1, of C2, and of C4.

Fig. 7 .
Fig. 7. (a) ′ L L / and ″ L L / vs. L for 20 year, (b) for 30 year and (c) for 40 year.From left to right: ′ L L / and ″ L L / vs. L calculated from converged solutions of C1, of C2, and of C4.

Fig. 8 .
Fig. 8. Positions of ′ L and ″ L for selected examples of = L dx for the realization 1 and for (a) C1, LTP5, (b) C2, LTP5, (c) C1, LTPn and (d) C2, LTPn, where n is 10 °C if within 50 years 〈 〉 T P t experiences 10 °C drop, otherwise n is the temperature drop of 〈 〉 T P t at the end of simulation (50 years).From left column to right column in all rows: Q1 and = = L dx 500 m, Q2 and = = L dx 500 m, Q2 and = = L dx 700 m, and Q2 and = = L dx 900 m.In all subfigures, white, cyan and yellow dashed lines represent the × L dx 2 2 extent, ′ L extent and ″ L extent, respectively.

1 .
Licensed regions boundary: this is the shell-like lateral inner boundary blocks of the license area.The license area is taken as the × × L dx 2 2 10-layer rectangle around the wells as shown in Fig. 1. 2. Average temperature of production wells defined by

Fig. 9 . 6 . 7 .
Fig. 9. Positions of ′ L and ″ L for selected examples of = L dx for the converged solution and for (a) C1, LTP5, (b) C2, LTP5, (c) C1, LTPn and (d) C2, LTPn, where n is 10 °C if within 50 years 〈 〉 T P t experiences 10 °C drop, otherwise n is the temperature drop of 〈 〉 T P t at the end of simulation (50 years).From left column to right column in all rows: Q1 and = = L dx 500 m, Q2 and = = L dx 500 m, Q2 and = = L dx 700 m, and Q2 and = = L dx 900 m.In all subfigures, white, cyan and yellow dashed lines represent the × L dx 2 2 extent, ′ L extent and ″ L extent, respectively.

Fig. 11 .
Fig. 11.Permeability and temperature distributions (at different time: LTP, LTB and 50 years) for six selected simulation samples/realizations: row (a) C1 low dx, row (b) C1 high dx, row (c) C2 low dx, row (d) C2 high dx, row (e) C4 low dx, and row (f) C4 high dx.The first column: the permeability distribution and the locations of wells, the second column: temperature distribution at LTP (of the realization), the third column: temperature distribution at LTB (of the realization), and the fourth column: temperature distribution at the end of simulation.

Fig. 4 ( 1 . 3
Fig. 4(b) shows LTB for various correlation lengths (C1, C2, C3 and C4), injection rates (Q1 and Q2) and variances in porosity field (v1 nd v2), when = dx L. Evidently, LTB is less sensitive towards the correla- tion lengths and the profiles in each particular group of Q-v are less spread.Similar conclusions can be made for = = = dx L d x L d x L 0.7 , 0.85 , 1.15 and = dx L 1.3 as shown in Fig. 4(c) and (d).In order to demonstrate the impact of well or doublet spacing on lifetimes (boundary vs. production wells-based calculations), in Fig.5all realizations as well as converged solutions are grouped into fixed

Fig. 12 .
Fig. 12.The profiles of 〈 〉 T W t over time for six selected simulation samples: (a) C1 and C2 samples, and (b) C2 and C4 samples.
(a), C2 and C1 in Fig. 10(b), C4 and C2 in Fig. 10(c), and C4 and C3 in Fig. 10(d).All profiles are plotted with respect to dx L / for fixed L's.The following observations can be made: 1.By increasing correlation lengths in both directions from 100 m to 500 m, Fig. 10(a) shows that almost universally lifetimes decrease (negative values in y-axis).The difference between lifetimes of C4 and C1 increases as L increases.2. By increasing correlation lengths in only y-direction from 100 m to 500 m, Fig. 10(b) shows that lifetime decreases further (larger negative values in y-axis) for small doublet spacings ( = dx L 0.7 and

Fig. 15 .
Fig. 15.(a) E prod vs. dx L / for fixed L's, (b) E prod vs. L for fixed = dx L / 1, (c) E pump vs. dx L / for fixed L's calculated at LTP, (d) E pump vs. L for fixed = dx L / 1, (e) an illustration of how flow rate increase impacts CoP based on E prod and E pump dependency on Q, and (f) an illustration of how flow rate increase impacts S based on E prod and E pump dependency on Q.These are all calculated at LTP.