Multiscale formulation for coupled flow-heat equations arising from single-phase flow in fractured geothermal reservoirs

Efficient heat exploitation strategies from geothermal systems demand for accurate and efficient simulation of coupled flow-heat equations on large-scale heterogeneous fractured formations. While the accuracy depends on honouring high-resolution discrete fractures and rock heterogeneities, specially avoiding excessive upscaled quantities, the efficiency can be maintained if scalable model-reduction computational frameworks are developed. Addressing both aspects, this work presents a multiscale formulation for geothermal reservoirs. To this end, the nonlinear time-dependent (transient) multiscale coarse-scale system is obtained, for both pressure and temperature unknowns, based on elliptic locally solved basis functions. These basis functions account for fine-scale heterogeneity and discrete fractures, leading to accurate and efficient simulation strategies. The flow-heat coupling is treated in a sequential implicit loop, where in each stage, the multiscale stage is complemented by an ILU(0) smoother stage to guarantee convergence to any desired accuracy. Numerical results are presented in 2D to systematically analyze the multiscale approximate solutions compared with the fine scale ones for many challenging cases, including the outcrop-based geological fractured field. These results show that the developed multiscale formulation casts a promising framework for the real-field enhanced geothermal formations.

increase the effective formation conductivity. These formations are naturally developed over large (km) length scales, while the heterogeneity of the damaged matrix needs to be resolved at fine (e.g., cm) scales. Fractures naturally add to the complexity of the mathematical formulations by introducing significant contrasts in the conductivity and geometry [1][2][3][4]. Their important role in the flow and transport of mass and energy can be properly investigated only if they are explicitly represented in the computational domain [5][6][7][8][9][10][11][12]. The embedded discrete fracture model (EDFM) [13][14][15] has been developed to resolve several geometrical challenges due to explicit treatment of the fractures. EDFM has been extended to complex scenarios in multiphase iso-thermal reservoir simulation [16][17][18], and, importantly, to geothermal systems [19]. Recently, a consistent projection-based EDFM (pEDFM) for flow has been proposed to account for all types of fracture conductivities (from flow barriers to high conductive channels) [20].
The size of final fine-scale systems describing single phase flow in fractured geothermal reservoirs, even after EDFM modeling approach, is beyond the scope of state-of-the-art commercial simulators. Upscaling these highly heterogeneous discrete quantities-in order to reduce the computational costs-would lead to inaccurate simulations, with no error control ability to the reference system. Therefore, new modeling and simulation techniques are more than ever on demand.
Multiscale finite volume methods have been developed for resolving this computational challenge by constructing coarse-scale systems based on local basis functions [21,22]. They are mainly developed for flow equations with complex fluid physics [23]. Together with their recent developments for fractured media [10,15,[24][25][26][27], they form a promising approach for real-field applications. For geothermal applications, however, the coupled flowheat equations need to be considered, which leads to additional complexities both in linear (size of the discrete system) and nonlinear (temperature-dependent coefficients and nonlinear coupling terms) aspects [28].
In this work, we propose a multiscale formulation for coupled flow-heat equations in fractured porous media, where not only the flow but also the heat equation are mapped to the coarse-scale system by using local basis functions. We investigate two different matrix-fracture coupling procedure for heat and flow basis functions, namely totally independent and semi-dependent. These two approaches differ from each other by the amount of matrixfracture coupling and number of matrix basis functions in the local system construction. The nonlinear coupling between flow and heat is treated with a sequential implicit approach.
Several challenging test cases are considered where the fractures play major role in transport of the cold water into the reservoir, and thus enhancing the production of heat. The large temperature gradients (due to slow to fast flow field) adds to the complexity of the simulations, which form a good basis to investigate the accuracy of the developed multiscale method. For all the investigated test cases, including the outcrop-based characterized formation, the basis function formulation for both the flow and heat equations are shown to be able to approximate the reference fine-scale solutions very well. Specially, the very first approximate multiscale solutions (with no iterations) are compared with the fine-scale solutions after the first Newton update, with very close agreement. Note that, even though a single-phase flow is considered here, the conservative velocity field of the Multiscale Finite Volume (MSFV) method is crucial for accurate transport of enthalpy which appears in the heat balance equation. The proposed multiscale finite volume method therefore casts a promising approach for field-scale geothermal studies.
Note that this work focuses on developing an accurate method to approximate the fine scale (fully resolved) solution. While the efficiency improvement is also part of the goal of the multiscale method, it is not extensively studied in this work and would be subject to future studies and development. However, the efficiency of the multiscale method for pressure solver in fractured reservoirs has been studied, and interested readers are referred to [24]. This paper is organized as follows. First, the mass and energy conservation equation of single-phase water system in fractured porous media are presented, together with the coupling strategy used to calculate pressure and temperature. Then, the MSFV method is introduced, and its application for pressure and temperature calculation is explained. After that, the numerical results are presented and discussed, including the simulation results of a realfield fracture geometry taken from outcrop data. Finally, the conclusions are presented.

Mass conservation equation
Mass conservation for single-phase flow in fractured porous media with embedded discrete fracture modelling (EDFM) approach reads for matrix on m ⊂ R n , and on f ⊂ R n−1 for fractures. In this work, n = 2. But these general formulations are also valid for three-dimensional domain (i.e., n = 3). Note that if gravitational effects are considered, one has to replace the pressure p with the potential (p − ρgz), where ρ, g, and z are, respectively, density, gravitational acceleration, and the coordinate along the gravity (pointing downward). Moreover, the superscripts m, f , and w indicate, respectively, the matrix, fracture, and well (source) terms. Moreover, φ stands for the porosity, λ the mobility, and q the flow rate. Note that the equation for fractures is defined in a lower dimensional space than for the matrix. The mobility is defined as λ = k/μ, where k is the absolute permeability and μ is the water viscosity. While the matrix permeability is characterized from geological inputs (possibly after proper upscaling), the permeability of the fractures can be defined (under fully developed flow assumption between parallel plates of distance (aperture) a) as k f = a 2 /12.
The well flow rates read for matrix and for fractures, where P I is the well productivity index [29], and β is the normalized well productivity index, with β m normalized with the matrix control volume V (i.e., β m = (P I λ)/V ) and β f normalized with the fracture area A (i.e., β f = (P I λ)/A). The discrete flux exchange between matrix and fractures, i.e., q mf and q f m , are modeled using EDFM approach as and Here, CI is the connectivity index between matrix and fracture, λ f −m is the effective mobility at matrix-fracture interface, and η is the normalized connectivity index, with η m normalized with the matrix control volume V (i.e., η m = CI λ f −m /V and η f normalized with fracture area A (i.e. η f = CI λ f −m /A, [24]. The discrete connectivity index CI allows for the representation of the discrete fracture element i overlapping with the matrix element j , i.e., where A i−j the surface area of fracture element i inside the element j , and d i−j is the average normal distance between the two elements [15].

Energy conservation equation
In this study, local thermal equilibrium between fluid and solid is assumed [30][31][32], i.e., the rock and the fluid have the same temperature at any given location. Under this assumption, the single-phase energy conservation equation on m ⊂ R n for matrix, and on f ⊂ R n−1 for fractures. Here, ρ r and C pr are the density and the specific heat capacity of the rock, U and h are the water specific internal energy and specific enthalpy, respectively. Also, T is the temperature and λ c is the average thermal conductivity, which is computed as λ c = φλ cw + (1 − φ)λ cr . The well source term, q * w H is defined as with q * w defined in Eqs. 3 and 4. Finally, u is the mass flow rate, which reads u = −ρλ · ∇p according to Darcy's law. From Eqs. 8 and 9, the matrix-fracture heat coupling is divided into two parts: conduction and convection. The conduction coupling terms q c are defined analogous to the matrix-fracture mass (flow) transfer, i.e., and where η c is the normalized conductive connectivity index; with η m c normalized with the matrix control volume V (i.e.
and q f m Note that the convection and conduction heat transfer are defined in a strictly conservative manner, i.e.,

Sequential implicit simulation strategy
The fluid properties depend on both pressure and temperature (see Appendix A for detailed formulations). These properties result in non-linearity of each mass and heat transfer equation, as well as the co-dependency between the coupled equations. A sequential implicit formulation is followed to treat the nonlinearly coupled mass and heat transfer equations for both fine-scale and multiscale simulations, where the pressure equation is first solved, then the total velocities are obtained and then the temperature solution in the fractured media is obtained. Notice that these equations are nonlinear functions of pressure and temperature and therefore are solved implicitly using the Newton-Raphson method.
The linearized equation for the general unknown x x x (i.e., p p p or T T T ) reads (17) or, in expanded form, ⎡ ⎢ ⎣ The superscript ν indicates the iteration stage, A A A the system matrix, and the vector f f f ν x stands for the right-hand-side terms. Note that for each unknown x x x, matrix, fracture, and well terms are present and that the system matrix and righthand-side terms depend on both p p p and T T T . As such, the complexity of the system is quite significant for real-field applications.
Algorithm 1 Sequential solver procedure Update pressure-dependent properties 4: Temperature solver: solve for T T T ν+1

MSFV method for nonlinear flow and heat equations
MSFV method approximates the fine scale solution by superposing the coarse scale solution with the basis functions as the interpolator. The approximate solution is defined as where x is a generic term denoting the unknown (i.e., p or T ), x the approximate fine scale solution, andx the coarse scale solution. The superscript * indicates the domain on which the unknowns are defined (i.e., matrix or fracture), and * • x the local basis functions in domain * coupled with domain • (i.e., matrix (m), fracture (f ), or well (w)).
Finally, N f is the number of fracture networks, N cf i the number of primal coarse cell in fracture i, and N w the number of wells (N w,inj in temperature system, which is the number of injection wells).

Multiscale grids
In the MSFV methods, two types of coarse grids are constructed and imposed on the fine scale grids. The primal coarse cells are constructed as the coarse-scale control volumes, while the dual coarse grids are overlapping coarse grids bounded by the coarse nodes (vertices) on which the local basis functions are computed (see Fig. 1). Figure 1 also illustrates the embedded fracture networks. Using EDFM benefits the multiscale implementation in that the coarsening strategy of the fracture elements is entirely independent of the matrix coarsening. Moreover, the fracture elements can be connected to any matrix cell, i.e., a vertex, edge, or an internal cell (based on dual-corsecell partition [33]).

MSFV for the flow equation
In development of an efficient multiscale method, the proper choice of basis function formulation is important. The important factors to consider for basis functions are their accuracy in representation of the underlying heterogeneity (accuracy), and their independency on the primary unknowns for adaptivity (efficiency). In this work, these aspects are being considered to formulate both pressure and temperature basis functions. The general formulation of the pressure basis function is written as is different for each coupling strategy. Equation 20 is formulated based on an equivalent incompressible system equation. This formulation is proven [34] to be the most efficient strategy (based on CPU measurements) because it eliminates the need to frequently update the local basis function, while the fully compressible coarse-scale system takes care of the global compressibility (and time-dependent) effects.
The pressure basis function formulated using Eq. 20 needs to be calculated at the beginning of the simulation, and updated only if the water mobility λ changes above a prescribed tolerance (i.e., due to the change of the water viscosity). Moreover, with regard to the matrixfracture coupling, we consider two different coupling strategies for the basis function calculation, namely the totally independent (fully decoupled) and semi-dependent (coupled) strategy. More precisely, the totally independent (decoupled) strategy formulates basis function without any coupling between matrix and fractures, while the semidependent strategy is formulated with partial coupling between matrix and fractures [24]. In the semi-dependent strategy, the matrix basis function is coupled with the fractures, but fracture basis function is decoupled from the matrix. This way, the semi-dependent approach enriches the matrix basis functions with the number of fracture course cells inside a dual-coarse domain.

Totally independent approach for basis functions
In the totally independent coupling strategy, all basis functions are calculated independent of interactions with other domains, i.e., An example of the pressure basis function calculated using this strategy is shown in Fig. 2. In Fig. 2a, it is shown that the matrix basis function is not affected by fracture existence, as well as fracture basis function not affected by matrix basis function in Fig. 2b. The basis function forms a partition of unity, meaning that the sum of all the basis function is equal to 1.

Semi-dependent approach for basis functions
In the semi-dependent coupling strategy, fracture basis function ff p is first calculated, decoupled with the matrix basis function, using These values are then used as Dirichlet boundary condition to calculate mf and setting to account for the connectivity of matrix basis function with the fracture domain. An example is shown in Fig. 3a, where ff p is plotted in the fractures and mf p is plotted in the matrix with the coupling effect clearly observed.
The matrix basis function mm p is calculated by setting therefore also accounting for fracture existence. An example is shown in Fig. 3b where the fracture basis function f m p is set to 0, and the matrix basis function mm p observing the effect of the fracture existence, as though the fractures act as flow barriers. This strategy also results in partition of unity.

Fine scale flux reconstruction
In the pressure MSFV method, one of the most important step is the fine-scale conservative flux reconstruction. In MSFV, the mass fluxes are conservative only at coarse scale. Therefore, the fine scale fluxes need to be obtained via additional reconstruction step [35]. This is especially important in multiphase flows, to accurately predict the saturation front since the fractional flow is sensitive to the flux. In geothermal simulations, conservative mass flux is also needed in the convection part of the energy balance calculation due to its velocity dependency. Therefore, it is worth revisiting the fine scale flux reconstruction in this subsection.
The mass flow rate formulation is valid at the primal coarse cell boundaries ∂ c . The conservative fine-scale flux can be reconstructed after solving locally on primal-coarse cells c , subject to the boundary condition (ρλ∇p c ) ·n n n c = (ρλ∇p ) ·n n n c (27) at ∂ c . Here,n n n c is the normal vector pointing out of the primal coarse cell boundaries, meaning that the fluxes at the coarse cell interfaces are used as Neumann boundary condition to calculate the reconstructed local pressure. The locally conservative mass flux is finally reconstructed as

MSFV for the heat equation
To exploit the efficiency of the temperature basis functions, they are formulated based on the conduction term within the whole energy balance equation. This allows for convenient implementation, as well as efficient algorithm (since basis functions are not required to be frequently updated). The general formulation of the temperature basis function can be defined as is different for each coupling approach, and is defined the same way as in pressure MSFV method (see Eqs. 21, 22, 23, and 24).
Note that the temperature basis function formulation is slightly different with pressure basis function. This is due to the fact that the well source term in the energy balance equation is defined through the enthalpy flow. Therefore, there is no explicit relation connecting the well and the matrix or fracture temperature. As such, the well basis function is omitted in Eq. 29.
The temperature basis functions depend only on thermal conductivity. And since the thermal conductivity is considered to be constant in this work, the basis functions therefore do not need to be updated frequently as they will also remain constant. In models where the thermal conductivity is considered to be non-constant, then an adaptive update of the temperature basis functions is necessary (i.e., if λ c changes above a prescribed tolerance). As will be seen in the result section, this formulation is shown to be working well to interpolate the coarse-scale temperature values to the fine scale. Note that these thermal basis functions along with the flow basis functions form the full prolongation (interpolation) operator to map between coarse and fine scale values for flow and heat.

Multiscale algebraic description and algorithm
In the algebraic formulation of the MSFV method, i.e., AMS [36,37], the multiscale procedure can be described by the prolongation P P P and the restriction R R R operators. The prolongation operator is a matrix constructed by the basis function values (interpolators) to map the coarse scale to fine scale solution. The restriction operator, on the other hand, is useful to map from fine scale to coarse scale. In finite-volume formulation, it acts as an integrator of all the fine scale fluxes, source/sink terms, as well as accumulation inside a primal coarse cell [36]. In this section, the algebraic description is explained in a generic way for both pressure and temperature calculation. More specifically, the prolongation operator reads P P P = ⎡ ⎣ P P P m P P P f P P P w ⎤ ⎦ = ⎡ ⎣ P P P mm P P P mf P P P mw P P P f m P P P ff P P P f w P P P wm P P P wf P P P ww where P P P * stores the basis functions defined on domain * , * • c,d , i.e., In the totally independent coupling strategy, the submatrices P P P mf and P P P f m are zero matrices, resulting in a sparser prolongation operator than in the semi-dependent coupling strategy, where only P P P f m is zero. Note also that P P P * w for temperature calculation has the column size of N w,inj , and P P P w * has the row size of N w,inj instead of N w .
The MSFV restriction operator is defined as and in MSFE method, the restriction is defined as the transpose of the prolongation operator, R R R F E = P P P T . Now that both operators are defined, the coarse scale system in equation is written algebraically as wherex x x ν+1/2 c is the coarse scale solution (i.e., pressure or temperature), and the superscript ν + 1/2, indicating that this stage will be complemented by a second stage smoother to be explained later.
Note that in Eq. 33, (R R RA A A ν P P P) constructs the coarse system matrix A A A ν c . The approximate fine scale solution is found as (34) or in residual form, δx δx δx ν+1/2 = P P Pδ δ δx x x ν+1/2 c = P P P(R R RA A A ν P P P) −1 R R Rr r r ν , where r r r ν is the fine-scale residual and is calculated as In each solver, both δp δp δp ν+1/2 and δT δT δT ν+1/2 are calculated first using multiscale operators (see Eq. 35), and then a 2 nd stage smoother (in this study, ILU (0)   Here, we employ 5 ILU(0) iterations per stage. This twostage multiscale procedure is repeated iteratively until the norm of residual goes below the prescribed tolerance.
The approximate fine scale solution is finally calculated as x x x ν+1 = x x x ν + δx δx δx ν+1/2 + δx δx δx ν+2/2 , where x x x ∈ {p p p , T T T }. The MS algorithms for pressure and temperature are presented in Algorithm 2.

Numerical results
In this chapter, numerical results are presented first to validate the EDFM model for coupled flow-heat equations, and then to investigate the performance of the multiscale simulation strategy for fractured reservoirs.

Test case 1: validation of EDFM
In this test case, the fine scale EDFM model is validated by comparing it to the result of the fully resolved Direct Numerical Simulation (DNS), used as a reference. The DNS result is obtained by using a very fine grid such that the fractures are captured as equi-dimensional (heterogeneous) objects [15]. EDFM, on the other hand, imposes much coarser grids and models the impact of the explicit lower-dimensional fractures by introducing fracture-matrix connectivities.
The fracture aperture is 0.0101 m, which can be fully resolved by imposing 99 × 99 DNS grid cells Fig. 4. This aperture leads to the fracture permeability of k f = 8.50 × 10 −6 m 2 . The simulation parameters are shown in Table 1. Figure 5 presents the pressure and temperature solutions obtained from EDFM and DNS simulators. Note that the EDFM solutions are obtained by imposing only 11 × 11 matrix cells and 14 fracture elements. It is clear that the EDFM solutions are in good agreement with the DNS reference ones.
The following measures are considered for the error of pressure and temperature: assuming x DNS x DNS x DNS 2 = 0, where x x x is flow rate q q q for pressure and enthalpy flux q q q H for temperature at both left and right boundary faces. The error norms for both pressure and temperature at different time steps are plotted and shown in Fig. 6. This figure also presents the EDFM error study at different times for the case when 33 × 33 EDFM grids are imposed, with 40 fracture elements. The EDFM errors are plotted against Pore Volume Injected (PVI) which is a non-dimensional time measure. The pressure errors are shown to be increasing slightly twice, and this is most likely because there are two pressure transient stages. Initially, the reservoir is hot and therefore water viscosity is lower and water density is higher, making the pressure gradient low. However, the injected water is cold and therefore increasing the pressure gradient in areas close to injection wells (which in turn decreasing the pressure farther from Fine-scale reference temperature with 99 × 99 matrix and 85 fracture elements (a) and multiscale approximate temperature solutions obtained using independent (b) and semi-dependent coupling (c) methods with 9 × 9 coarse matrix and 8 fracture grid cells at the first iteration stage before smoothing. The corresponding relative error norms (d and e) are e T 2 = 0.0343 (independent coupling) and e T 2 = 0.0198 (semi-dependent coupling). After 1 stage of smoothing these errors reduce to e T 2 = 0.0146 (independent coupling) and e T 2 = 0.0035 (semi-dependent coupling)    Fig. 13 Fracture geometry taken from outcrop data [39] (a) (b) (c) the injection wells). After the injection water gets farther into the reservoir, the pressure gradient close to the injection wells gets lower and the gradient farther in the reservoir gets higher until it reached semi-steady state. More specifically, the error is due to (1) significant difference between the grid resolutions imposed by each method and (2) the error of EDFM fracture model. Nevertheless, the two approaches are in good agreement.

Test case 2: homogeneous reservoir with a diagonal fracture
A quarter of a five-spot test case is considered in a homogeneous reservoir with a diagonal fracture. The simulation parameters are shown in Table 2. EDFM imposes 85 fracture and 99 × 99 matrix elements. The geometry of the fracture within the reservoir is shown in Fig. 7. The multiscale simulator imposes 9 × 9 coarse grids for matrix and 8 for fractures with two different coupling strategies for basis function calculation. Figures 8 and 9 show the converged solution of both fine scale reference as well as multiscale pressure and temperature. The white lines shown in the plots are the primal coarse cell boundaries. The relative error norms of the multiscale solutions are e p 2 = 2.65 × 10 −5 and e T 2 = 1.62 × 10 −5 (independent coupling), and e p 2 = 2.18 × 10 −5 and e T 2 = 1.58 × 10 −5 (semidependent coupling). It is shown that both independent and semi-dependent coupling strategies result in very good results.
The multiscale pressure and temperature solutions at the first iteration (before smoothing) are also presented in Figs. 10 and 11, respectively, to show that the multiscale provides very good approximations even with no secondstage smoother nor any other (inner and outer) iterations. These results are also compared to the reference fine scale solutions, demonstrating the accuracy of the developed multiscale formulation.
Independent coupling for pressure basis function calculation results in slightly higher error at the fracture tips, where-as expected-the interaction of matrix and fracture domain is relatively high. Note that the temperature field experiences a rapid change in the location of the fracture, due to rapid transport of cold water through the fractures. As such, the significant temperature contrast is created fairly quickly throughout the reservoir, leading to strong nonlinear time-dependent solution field. In the area near the midpoint of the fracture, the semi-dependent coupling provides a slightly lower error for temperature. This can be explained by the fact that in this area, pressure difference between matrix and fracture is lower and therefore, the heat exchange is more conduction dominated. Since the matrix-fracture coupling in the semi-dependent approach for temperature is formulated based on conduction, it leads to better approximation in this area. This is clear from results presented in Fig. 11. Nevertheless, as shown, the multiscale method can represent the complex solution field accurately, even with no smoothing stage. At the first iteration stage, the relative error norms of the multiscale solution obtained before smoothing are e p 2 = 0.0081 and e T 2 = 0.0343 (independent coupling), and e p 2 = 0.0016 and e T 2 = 0.0198 (semi-dependent coupling). After 1 stage of smoothing, the errors are reduced to e p 2 = 0.0076 and e T 2 = 0.0146 (independent coupling) and e p 2 = 0.0008 and e T 2 = 0.0035 (semi-dependent coupling). The results are summarized in Table 3.
The semi-dependent coupling strategy leads to lower errors compared to the independent coupling, especially in the area surrounding the fracture. However, it does not bring much improvements. The errors obtained using independent coupling are not significant and could be resolved with several smoothing and iterations.

Test case 3: fracture geometry from outcrop data
A quarter of a five spot test case is considered in a heterogeneous reservoir with dense and complex fracture networks taken (by applied geologists of TU Delft) from outcrop data in Brazil [39]. The base-10 logarithm of permeability and average thermal conductivity are plotted in Fig. 12, and the simulation parameters are shown in Table 4. EDFM generates 3860 fracture and 100 × 100 matrix elements. The geometry of the fractures within the reservoir is shown in Fig. 13. The multiscale simulator imposes 10 × 10 coarse grids for matrix and 386 for fractures. For this test case, all the results presented are using independent coupling method.
As shown in test case 2, the decoupled approach approximates the solution really well. Therefore, for this test case, only the results obtained using the decoupled approach are shown for conciseness. Figures 14 and 15 show the converged solution of both fine scale reference as well as multiscale pressure and temperature. The relative error norms of the multiscale solution obtained are e p 2 = 8.22 × 10 −6 and e T 2 = 7.07 × 10 −6 .
To make it more suitable for real field application and to further validate the multiscale method for this complex test case, mass and enthalpy production rate at the production  Fig. 16.
In the plot, it is shown that both multiscale and the fine scale mass and enthalpy production rate have a very good match, with e p 2 = 0.028 and e T 2 = 0.036. The error of the enthalpy production rate is slightly higher due to the fact that the enthalpy production rate is calculated based on the mass production rate, and therefore the error from the mass production rate calculation is propagated. Nevertheless, the error is still relatively low and therefore acceptable. The multiscale solutions for both pressure and temperature at the first iteration stage before smoothing, along with the fine scale reference solutions for comparison, are shown in Figs. 17 and 18. The corresponding relative error norms before smoothing are e p 2 = 0.0111 and e T 2 = 0.0467, which are relatively low for a complex model. After smoothing, the errors are further reduced to e p 2 = 0.0110 and e T 2 = 0.0180. The results of this test case are summarized in Table 5. The error reduction in pressure after the smoothing stage is not significant, while error reduction in temperature is more drastic because of the higher complexity of the heat transfer equation and the simplicity of the basis function used to calculate it. It is also clear that the temperature distribution in Fig. 18 (without smoothing) already captures the preferential flow path of the fluid. From the results obtained, it can be concluded that independent coupling method gives reasonably good approximations, even before the smoothing stage for a heterogeneous system and very dense fracture networks.

Conclusion
In this work, a multiscale method for coupled single-phase flow-heat equations in fractured reservoirs was developed. The method avoids excessive upscaling of the parameters, and honours fine-scale heterogeneity in construction of the coarse-scale system for both flow and heat equations. This is achieved by formulating flow and temperature basis functions, allowing for the accurate map between fine and coarse scale solutions. The coupling between the equations was treated by a sequential implicit framework, where both pressure and temperature systems were solved by a MSFV method. The multiscale formulation was enriched due to the presence of the fractures, with two coupling approaches for local basis functions of each solver. An EDFM approach was adapted to the framework, which allows for independent grids for matrix and fractures. This further facilitated the convenient multiscale formulation and implementation, as totally independent coarse grids were also imposed on matrix and fractures. Test cases were performed first to validate the implementation of the simulator (via comparing its results with a DNS approach), and then to systematically assess the performance of the multiscale method for heterogeneous and highly fractured media. A fracture formation from a real-field outcrop was also considered to illustrate the capacity of the algorithm in addressing complex fracture networks.
Although we employed the MSFV iterations to reach convergence in our sequential implicit framework, one can stop iterations before convergence is reached. Specially, as shown in the results, the initial multiscale approximate solutions in presence of heterogeneity and fractures were close to the fine-scale solution at the same stage of iteration. The tolerance to stop iterations of a conservative multiscale solver needs to be defined based on the influence of the solution in the overall accuracy of the coupled solutions, the stability of the time-dependent solutions, and the uncertainty within the parameters. Similar to previous studies for coupled flow and transport [40], such a study is needed for coupled P-T as a future work. Specially, in the presence of strong coupling one may consider formulating a multiscale methodology for fully implicit systems [41,42].
As for the multiscale basis functions, to exploit the maximum efficiency, the temperature basis function was formulated based on the elliptic part of the energy conservation equation (i.e., the conduction term). Numerical results showed that such an approach is well suited for the considered single-phase fluid-dynamic system, i.e., it leads to accurate results even without smoothing stage.
In this work, a robust approach for solving the coupled pressure and temperature equations in fractured heterogeneous reservoirs was developed. The results presented show promising framework for further developments for fieldscale enhanced geothermal systems. Future developments need to consider multiphase (including steam) effects for fluid and the geo-mechanical effects (including fracture activations or closures and propagation) for solid rock.
Moreover, the empirical parameters, n i , are shown in Table 6. The combination of u ws = 420 kJ/kg, C pw = 4.2 kJ/kg, and T s = 373 K was found to provide the best fitting values for internal energy calculation. More precisely, compared with the data, the density relative error norms were below 1% in most regions and 2.2 % near the critical point. Similarly, the internal energy errors were less than 6%.  [39] due to the high number of fractures in the model. The fractures are defined by two points: A and B and the x and y coordinates are given in the tables.

B.1 Fracture coordinates for test case 1
The fracture coordinates for test case 1 are listed in Table 7.

B.2 Fracture coordinates for test case 2
The fracture coordinates for test case 2 are listed in Table 8.