Delft University of Technology Adaptive Dynamic Multilevel Simulation of Fractured Geothermal Reservoirs

Article history: Received 17 January 2020 Received in revised form 11 May 2020 Accepted 12 May 2020 Available online 22 May 2020


Introduction
Geothermal energy resources are becoming more important in the transition towards renewable energy. Due to low carbon footprints [1][2][3][4] and high potential of sustainability, the demand for geothermal energy has been rising and is expected to increase further within the next decades. Accurate scalable modeling of fluid and heat transport in geothermal reservoirs is of paramount importance in order to fulfil both scientific and societal expectations on successful field development plans [5,6].
While maintaining field scalability, simulation models should accurately capture high contrasts in both mass and heat transport properties inside the porous media. Moreover, a large number of reservoirs contain complex geological networks of fractures and faults with a broad range of conductivities. Such extreme contrasts have significant impact on flow and heat patterns. Therefore, high fidelity representation of the physical phenomena within the heterogeneous fractured reservoirs is crucial [7]. However, while achieving this task, the computational science community faces a variety of complex challenges, as the models demand high-resolution grids to be imposed on the entire reservoir (in orders of kilometres) [8]. Non-linear behavior of the system due to strong mass-heat coupling results in poor stability and convergence. In case of multi-phase flow (e.g., high-enthalpy systems) these issues become more severe [9]. The geo-mechanical processes (i.e., elastic and plastic deformation) [10][11][12], reactive transport (e.g., geo-chemical interaction between the substances) [13][14][15] and compositional changes are other notable challenges that the system might face. Expectedly, the presence of fractures and faults [16][17][18] with wide heterogeneity contrasts compared to the rest of the formation significantly increases the computational complexity. These challenges introduce high demands for developing advanced simulation methods that are able to provide efficiency (i.e., applicable to real field-scale problems) while maintaining accuracy at the desired level. To address these challenges, various advanced simulation schemes have been developed.
Discrete Fracture Models (DFM) aim to represent fractures as conforming lower dimensional domains by discretizing them at the matrix cell interfaces [19][20][21]. On the other hand, the Embedded Discrete Fracture Model (EDFM) makes use of non-conforming grids to discretize fractures as lower dimensional objects fully independent from the matrix grid structure [22][23][24][25][26]. Conservative flux transfer terms are introduced between overlapping matrix and fracture elements. EDFM is applicable to provide acceptable solutions for highly conductive fractures. Even for high-conductive fractures, it leads to very small, yet unphysical, flux leakage. More importantly, EDFM cannot accurately represent flow barriers (e.g., sealing fractures and impermeable faults). In order to resolve this limitation, projection-based EDFM (pEDFM) was introduced [27] where consistent connectivities between matrix and fractures were developed. Therefore, pEDFM can be applied to fractured porous media with generic range of conductivity contrast between fracture and the hosting rock. Even with pEDFM independent grid, the size of the reservoir and density of the fractures make any fine-scale simulation approach impractical.
Multiscale Finite-Volume (MsFV) [28][29][30][31][32][33] and Multiscale Finite-Element (MsFE) [34][35][36] methods have been devised to accomplish efficient simulations in presence of heterogeneity, by solving coarse-scale systems, while avoiding upscaling, yet preserving fine-scale heterogeneities [37,38]. Dynamic Local Grid Refinement (DLGR) [39][40][41][42][43][44][45][46][47] is another method, which employs fine-scale resolution in parts of the domain where needed (i.e., moving fronts). Algebraic dynamic multilevel (ADM) method has been developed to address the multilevel multiscale challenges on equations with elliptic, parabolic and hyperbolic systems of PDEs [48]. By extending the grid refinement strategy and taking the advantages of multilevel multiscale basis functions, the ADM method is able to algebraically map the fine-scale discrete systems (either fully implicit or sequential) to a dynamic multilevel system for all unknowns. Once the system is solved at ADM resolution, it is prolonged back to fine-scale resolution. In case of fully implicit (FIM) formulation, as the system is solved at ADM resolution for all unknowns, reconstruction of conservative flux field is not required. This method has been further developed for multiphase flow in heterogeneous formations with compositional and capillary effects [49], and fractured reservoirs using EDFM scheme [52,50]. Recent developments also include multiscale finite element method for flow-heat transfer modeling using EDFM approach [37,38]. However, such method has not yet been developed for coupled mass-heat transport in fractured porous media with pEDFM method.
In this article, an Algebraic Dynamic Multilevel (ADM) method for simulation of non-isothermal single-phase flow in fractured porous media is presented. Due to the complex non-linearity of the system (strong mass-heat coupling), two sets of equations (mass and energy conservation) are solved using the FIM approach. In the energy balance equation, conduction terms are described for both fluid and rock. Projection-based EDFM (pEDFM) is employed in order to explicitly represent fractures and to provide independent gridding of matrix and fractures. Here, the applicability of the pEDFM implementation [27] has been extended to a fully generic 3D geometry where it allows for including fractures (or flow barriers) with any orientation on the Cartesian coordinate system. This is crucial for practical applications. While maintaining the geometrical flexibility of EDFM, the matrix-matrix and fracture-matrix connectivities are adjusted to account for projection of fracture plates on the matrix transmissibilities. This allows for consistent modeling of fractures with all conductivity (from high to low) compared with that of the matrix. The fine-scale discrete system including both matrix and fracture network is obtained for two main unknowns, namely, pressure and temperature. We note that local thermal equilibrium is assumed, meaning that the temperatures in fluid and rock are considered to be equal. This assumption can be recognized as safe for the majority (or perhaps all) of the geo-science applications [51] due to large contact area between the liquid and solid phases. Once the fine-scale system is discretized using pEDFM, it is mapped into a dynamic multilevel resolution grid (i.e., ADM grid). The ADM grid is constructed based on a hierarchically nested multilevel multiscale grid resolution, where locally computed multiscale basis functions are employed at each coarsening level [34,28,30]. These basis functions are computed only at the beginning of the simulation as a pre-processing step. This provides higher computational efficiency, while capturing the fine-scale heterogeneous properties accurately.
In order to construct the multilevel multiscale operators for the geothermal fractured system, different types of basis functions are computed for pressure and temperature unknowns. For pressure, as a global variable, multiscale basis functions are computed. However, temperature resembles a local variable for which constant basis functions are used. This distinction in computation of basis functions is preserved on a hierarchical pattern for all the coarsening levels across the domain. To capture the effect of permeable and impermeable fractures, all basis functions are fully coupled, meaning that extra basis functions are computed for fractures as well. The construction of all basis functions are done fully algebraically. In this process, the multilevel prolongation operator is assembled so that the partition of unity holds [53,54]. The multilevel finite-volume restriction operator ensures mass conservation inside the domain. At each time-step, the ADM map is obtained by using a front tracking technique to detect the region where the solution gradient is higher than a user-defined threshold. These sub-domains of sharp gradient are resolved at fine-scale level. Moreover, regions near the source terms are always kept at fine-scale resolution to maintain accuracy. The rest of the domain is resolved at hierarchically nested coarser resolution. This ADM map is able to adapt itself to the moving front dynamically and efficiently, while capturing the physics of interest in both flow and heat transport with an accuracy set by the user-defined threshold. The Newton linearization scheme is employed to find nonlinear solutions iteratively. Once the ADM solution is found, it is prolonged back to the fine scale using prolongation operators.
Various homogeneous and heterogeneous test cases are considered and their numerical results are presented to demonstrate the accuracy and applicability of the devised method. Furthermore, its sensitivity towards fractures with different ranges of conductivity is studied.
This article is organized as follows. The governing equations and FIM approach are described in Section 1. The fine-scale discretization using pEDFM is covered in Section 2. The ADM method as the simulation strategy for fractured geothermal reservoirs is presented in Section 3. The test cases and their numerical results are shown in Section 4. The paper is concluded in Section 5.

Governing equations
In this section, the two main governing equations are described. Here, the temperatures of fluid and rock are assumed to be identical.

Mass balance
Mass balance equation for non-isothermal single-phase fluid flow in porous media with n frac discrete embedded fractures reads ∂ ∂t for the rock matrix (m) and ∂ ∂t .., n frac }, (2) for the lower dimensional fracture ( f i ). Here, the main unknown is the pressure p. Moreover, φ and K are the porosity and permeability of each medium. K is a tensor to account for a generic anisotropic case.
The Peaceman well model is used to obtain the well source terms for matrix and fractures Here, P I is the well productivity index and * is the effective mobility ( = K /μ) between the well and the penetrated grid cell in the medium. V and A are the control volume and control area used in the discrete system for matrix m and fracture f i respectively. The flux exchange terms Q mf i , Q f i m (matrix-fracture connectivities) and Q f i f j (fracture-fracture connectivities) are written as: (7) where C I denotes the connectivity index between each two non-neighboring elements. As an example, the connectivity index between i-th matrix element and j-th fracture element is calculated as where A ij is the area fraction of fracture cell j overlapping with matrix cell i and d ij is the average distance between these cells.

Energy balance
Assuming local equilibrium, energy balance on the entire domain reads (8) in the rock matrix (m) and .., n frac }, (9) in the lower dimensional fracture ( f i ). Here, beside pressure p, the second main unknown is T as the temperature in both fluid and solid rock. (ρ U ) ef f is the effective property defined as where U f and U r denote the specific internal energy in fluid and rock respectively. Distinctly, (ρ U ) m Moreover, H f is the specific fluid enthalpy. The three mentioned terms can be expressed as non-linear functions of pressure and temperature. ef f is the effective thermal conductivity of the saturated rock defined as Here, f and r are the thermal conductivities in fluid and rock, respectively. The subscripts f and r indicate fluid and solid rock. Note that, m Lastly, R mf i and R f i m are the conductive heat flux exchange between matrix m and the overlapping fracture f i . R f i f j denotes the conductive heat flux exchange from j-th fracture to the i-th fracture where the intersection occurs. Similar to mass flux exchange, the conductive flux exchange terms are non-zero only for the existing matrix-fracture overlaps or fracture-fracture intersections.
hold as well to honor the conservation of energy.
To obtain the conductive heat flux exchanges, i.e., R mf i , R f i m (matrix-fracture connectivities) and R f i f j (fracture-fracture connectivities), the embedded discrete scheme is used, i.e., where * ef f is obtained as harmonically-averaged property between the two non-neighboring elements. C I is identical to the connectivity index used in Eq. (5). The independent matrix 3D grid and the fractures 2D grids are shown. Note that each domain has its own grid and that any fracture orientation can be considered.

Correlations
The following correlations are used to compute certain fluid and rock properties.

Fluid viscosity:
Viscosity as a function of temperature is calculated via the following correlation [55], Fluid density: The density of fluid is a function of pressure and temperature, defined as [51] ρ f (P , where P S = 1 bar. c w (T ) and ρ f s (T ) are computed from the following empirical correlations [8,56] Fluid Entalphy: Fluid enthalpy is obtained via [51] H f (P , where u ws = 420000 J/kg.

Fine-scale discrete system using pEDFM
The coupled system of non-linear equations described above, i.e., Eqs. (1)-(9) with two main unknowns (i.e., p, T ) is discretized spatially using two-point-flux-approximation (TPFA) finite-volume scheme and temporally using the backward (implicit) Euler scheme [57]. The discretization results in set of structured grids for a three-dimensional (3D) porous medium with 2D fracture planes included that are generated independently. A visual example of the gridding structure is shown in Fig. 1. The obtained non-linear FIM system is solved using the Newton-Raphson iterative method [57].

Mass and heat fluxes
The mass flux exchange between each two neighboring control volumes i and j (inside one medium) using TPFA scheme can be written as ij denotes the transmissibility between grid cells i and j. A ij is the interface area between the two grid cells, d ij is the distance between their cell centers and K H ij is the harmonic average of the permeabilities set at the interface of grid cells i and j. The properties with superscript * are obtained using the upwind scheme. Note that density ρ and viscosity μ are functions of pressure and temperature.
The mass flux between each two non-neighboring grid cells i and j (either matrix-fracture or fracture-fracture connectivities) is obtained via EDFM formulation. This flux between a matrix (denoted as m) cell i and a fracture (denoted as f ) cell j is defined as where the transmissibility T mf ij = K H ij C I ij . K H ij is the harmonically averaged permeability between the overlapping matrix and fracture elements, and C I ij is the connectivity index between them. The EDFM models the matrix-fracture connectivity where A mf ij is the area fraction of fracture cell j overlapping with matrix cell i (see Fig. 2, left) and d ij is the average distance between these cells [24].
The mass flux exchange between two intersecting fracture elements i (on fracture f ) and j (on fracture g) is given as where T f g ij is the transmissibility between the two cells obtained via a lower dimensional connectivity index formulation. Mind that the intersection of 2D fracture plates is a line and for 1D fracture line-segments it is a point. Fig. 2 (right) shows a visualization example of intersection between two non-neighboring 2D fracture elements. The intersection forms a line segment I ij (shown in red color) with the average distances from the intersection segment of d Thus, T f g ij is obtained using an harmonic-average formulation, i.e., where, e.g., C I f iI ij represents the connectivity index between the 2D fracture element i and the 1D intersection line segment I ij , as shown in Fig. 2.
The convective heat flux exchange between the neighboring control volumes i and j reads Here, H * f is the enthalpy of fluid determined at the interface between grid cells i and j. One can conclude F ij = H * f F ij . Similarly, the convective heat flux exchange between non-neighboring elements can be obtained via multiplication of their mass flux exchange (F ij ) by the effective enthalpy (H * f ) determined at the intersection of two overlapping elements.
The conductive heat flux exchange between each two neighboring cells i and j is written as where, T cond,ij = The conductive heat flux exchange between non-neighboring elements i and j is obtained as for matrix-fracture connectivities and for fracture-fracture connectivities ij is the harmonic average of the effective conductivity between the overlapping matrix and fracture elements. A similar procedure as in Eqs (18) will be used to obtain the fracture-fracture connectivity indices.
At each time-step the fine-scale discrete mass balance equation reads for element i in matrix m and for element i in fracture f h . N m and N f k denote the number of elements in matrix m and fracture f k respectively. N n is the number of neighboring grid cells.
Similarly, the fine-scale discrete form of the energy balance equation at each time-step is written as: The residuals for the other three equations (i.e., the mass balance in fractures (r As the residual at time-step n + 1 (r n+1 ) is a non-linear function of pressure and temperature unknowns at time-step n + 1 (i.e., p n+1 and T n+1 ), the Newton-Raphson iteration method is applied to solve the system at each time-step. In other words, Here, the superscript ν denotes the iteration index. At each Newton iteration, the linearized system of equation is written as J ν δx ν+1 = −r ν . In this system, J ν is the Jacobian matrix including all derivatives. Each block J e α contains the derivatives of equation e with respect to unknown α, i.e., J e α = ∂r e /∂α (containing the list of, ∂r M B In order to achieve a converged solution, the following conditions have to be assured: Here, each threshold ( x ) is a user-defined tolerance set initially as input at the beginning of the simulation.
To improve the convergence, at the beginning of each time-step, the mass balance and energy balance equations are solved as decoupled systems in a sequential manner. In other words, first, decoupled linear system of mass balance equation (32) is solved and a vector of pressure updates (i.e., δ p) is obtained. After updating the pressure, all dependent properties and derivatives are re-calculated. Similarly, in the next step, the decoupled linear system of energy balance equation (33) is solved and a vector of temperature updates is acquired. Consequently, the temperature unknowns are updated and all dependent parameters and derivatives are recalculated. After these steps, the fully coupled linear system (30) is solved iteratively.
The extension to pEDFM is explained in section 2.2.2.
Obtaining the solution of this linearized system of equations is computationally challenging due to the strong coupling of mass and heat transfer, especially for real field-scale applications. In section 3 the simulation framework of this work is presented with which resolving the computational challenge is aimed using the ADM method.

pEDFM implementation
To correct for EDFM limitation on fractures with generic conductivity due to parallel transmissibilities, the matrix-matrix, fracture-matrix and fracture-fracture connectivities are modified in the overlapping regions. The mentioned modifications eliminate the parallel transmissibilities, such that pEDFM is applicable to any conductivity contrast between the matrix and fractures. Initially, all connectivities between two neighboring matrix cells that are disconnected due to the overlapping fractures are detected. Due to geometrical algorithm devised during the development of this method, a continuous projection path (visible in Fig. 3 as solid lines in light blue color) is automatically obtained on the interfaces as such it disconnects the neighboring connections letting the flux occur only on one route (i.e., through matrix-fracture-matrix).
Let fracture element f overlap matrix grid cell i (as shown in Fig. 3) with an area fraction A if . A set of projections is defined on the interface between the overlapped matrix grid cell i and its neighboring grid cells (in orange) that are affected by the crossing (i.e., j and k). Please note that in the 3D dimensional case, there will be three projections. For each dimension (i.e., x, y and z) the projection area fractions are obtained via: where γ is the angle between the fracture element and the interface connecting the matrix grid cell i and the neighboring grid cell in the corresponding dimension (shown in Fig. 3 on the zoomed-in section and highlighted in red color as A if ⊥x and A if ⊥y ). New transmissibilities are defined to connect the fracture element f to each non-neighboring matrix grid cells (i.e., j and k in the 2D example shown in Fig. 3): where, d i e f is the average distance between the fracture element f and matrix grid cell i e . i e f is the effective fluid mobility between these cells. As a result of the new transmissibilities, the connectivity between the matrix grid cell i and its corresponding neighboring cells is modified: x e i e f , X ∈ {x, y, z}.
For simplicity of the implementation, the modified transmissibilities are obtained by multiplication of coefficient α as a fraction of the projected cross-section area, and the cross-section area of the corresponding interface. Please note that for all overlapping fracture elements (except for the boundaries of the fractures), the projection will cover the whole interface.
Therefore, α is 1.0 for most of the cases, resulting in zero transmissibility between the matrix grid cells (i.e., T ii e = 0), thus removing the parallel transmissibilities [27].

Solution strategy
At the beginning of every Newton iteration, the fine-scale fully-implicit linearized system of equations from Eq. (30) is used to algebraically assemble a reduced linear system on a dynamic multilevel grid resolution which is defined at the beginning of each time-step. At the first stage (i.e., obtaining the dynamic multilevel grid), sets of  (37) for matrix m and fracture f i . The coarsening ratios γ l are constant on the hierarchical set of coarsening levels l, but both properties are defined independently for matrix and each fracture, providing practical flexibility. Combination of grid cells at different resolutions results in construction of the ADM grid. To map from the fine-scale resolution to this dynamic multilevel resolution, ADM employs a series of restriction (R) and prolongation (P) operators, which are constructed and computed only at the beginning of the simulation (see the offline stage of Fig. 4). The ADM system is written as Note that δx 0 represents the reference fine-scale solution. The ADM Restriction (R l−1 l ) and prolongation (P l l−1 ) operators are obtained from the static multilevel multiscale restriction (R l−1 l ) and prolongation (P l l−1 ) operators respectively. The static prolongation operator P l l−1 is constructed by assembling the locally computed basis functions (see section 3.2) at level l and is written as Note that P l l−1 consists of two main diagonal blocks for the two unknowns (i.e., pressure p and temperature T ). Within each main block, the off-diagonal terms refer to the coupling terms between matrix and fractures. The static restriction operator R l−1 To assure mass conservation, a finite-volume restriction operator is used. Therefore, R l−1 l elements contain either 1 or 0 where For the test cases that are presented in this article, different interpolators are employed for each of the variables.

Multilevel multiscale basis functions
The However, by choosing the fully-coupled approach, each medium will influence the local solution (i.e., the basis functions) of the overlapping domain, hence the off-diagonal blocks contain non-zero elements as well [54,50]. The prolongation operator is constructed fully algebraically [54] and it includes the basis functions at all coarsening levels. To compute the basis functions, first the fine-scale pressure system of equations is assembled. Consequently, the coarse grids are imposed on the fine-scale computational grid cells for matrix and fractures. The approach to construct the hierarchically nested primal and dual coarse grids is flexible in this work. One option is to impose the coarse grids by determining identical primal coarse grids at each coarsening level, namely each primal coarse grid at level l has identical size of γ l = γ l x × γ l y × γ l z . In this case the coarse nodes are not located on the corners or edges of the global domain. Fig. 5 shows an example of the primal and coarse grids of this option for a 2D domain of 75 × 75 grid including 3 fractures. The other option is to include the coarse nodes on the corners and edges of the global domain, resulting in primal coarse grids with different sizes. Fig. 6 illustrates the coarse grid construction of this approach. Please note that the level at which each of the media (either matrix or fractures) is coarsened can be defined independently.
Once the coarse grids are constructed, the basis functions are calculated algebraically for each coarsening level. The multiscale procedure converts the two-point flux approximation (TPFA) connectivities of the fine-scale system (i.e., A 0 ) to multi-point flux approximation (MPFA) scheme at coarse-scale systems for each coarsening level l (i.e., A l ). Therefore, prior to construction of coarse-scale system at level l, the system at level l − 1 is reduced to the TPFA scheme, namely  Here, the T P F A_Operator removes the MPFA components from the entries of the coarse-scale system matrix. This procedure is performed independently on each medium [50]. Fig. 7 shows coupled basis functions for a matrix coarse node (left) and two fracture coarse nodes (right).

Selection of the grid resolution
The fine-scale grid resolution is mapped into the ADM grid resolution using the ADM operators. This dynamic grid resolution is obtained by a user-defined grid selection criterion as a front tracking technique. At each time-step n, the ADM grid is selected based on the solution at time-step n − 1 (explicit scheme) using a threshold set as the input parameter for the grid selection criterion. The condition for which the criterion is considered, compares the difference of the values of a main unknown (in this case, the temperature) between neighboring cells and the grid cells in a coarse node (spatial gradient of the unknown) [48]. Let   The coarse grid block I is refined from coarsening level l to resolution l − 1 if the condition holds. Here, N refers to all coarse grid cells neighboring coarse grid cell I (at coarsening level l). Fig. 8 shows an example of a 3D domain including one fracture with two coarsening levels over matrix and one coarsening level over fracture. Please note the independency of coarsening strategy and grid selection procedure for matrix and fracture. It can be realized that the moving front of the flooding fluid will be captured with higher resolution as it corresponds to sharper gradients. Moreover, the grid cells around the wells are always kept at fine-scale resolution to guarantee accurate capture of source term fluxes.

Numerical results
Numerical results are presented in this section. First, the pEDFM model is validated. Thereafter, ADM results are compared against fine-scale simulations, used as a reference.

Test case 1: validation of pEDFM
This test case demonstrates the validation of pEDFM implementation. For this purpose, a 2D homogeneous domain of 3 m × 3 m is considered (see Fig. 9). Two fractures intersect in the middle of the domain forming a cross-shape. Both fractures are 1.5 [m] long with aperture of 8 × 10 −3 [m]. In terms of fractures conductivity, two scenarios are taken into account. In the first scenario, fractures are a f = 10 8 times more permeable than the rock matrix (K m = 10 −14 m 2 and Fig. 9. Visualization of test case 1 geometry, a 2D domain with a cross-shaped fracture network at the center.

Table 1
Input parameters of test case 1 for fluid and rock properties.

Property Value
Rock thermal conductivity ( r ) Fluid thermal conductivity ( f ) n e g l i g i b l e Rock density (ρ r )  Table 1 shows the input parameters of this test case.
To obtain the reference solution (denoted as DNS), a 375 × 375 grid resolution is imposed on the domain. This allows for fully resolving the flow inside the fracture. In this case fractures are defined as channels along the middle row and column of the discretized domain. EDFM and pEDFM simulations are run with four different matrix gridding resolution of 75 × 75, 45 × 45 and 25 × 25. At all runs, the fractures aperture is identical and set to be a f = 8 × 10 −3 [m] which is approximately equal to the matrix grid cell size of DNS discretization. The size of grid cells inside the fractures have similar dimensions as matrix grid cells. In cases with highly conductive fractures, the simulation is run for 6 [hours] and in the cases with low-permeability fractures, the simulation is run for 12 [hours]. The results are provided at 50 isochronal intervals.
In the first scenario, fractures are considered to be highly conductive (see Table 1 for more details). The results of the simulations as plots of pressure and temperature solutions are shown in the Fig. 10 and 11.
The plots in Fig. 12 present the pressure and temperature errors as functions of the simulation time. The error of variable x (i.e., pressure or temperature), denoted as e x , is calculated via At the center of the domain where the fractures intersect, a difference in the temperature distribution between the DNS, EDFM and pEDFM methods can be observed. Please note that the difference arises from the discretization approaches. In DNS method, fractures are actually channels with aperture of one grid cell. As a result, those matrix grid cells are flooded with the cold fluid. This is not the case for EDFM and pEDFM.
In the second scenario, both fractures are set to be 10 8 times less permeable than the matrix rock. Figs. 13 and 14 show the simulation results of this scenario as pressure and temperature solutions respectively. Fig. 15 shows the pressure and temperature errors versus the simulation time for the second scenario. As the EDFM is incapable of capturing the low permeable fractures and the results are not accurate and representative, the errors are only shown for pEDFM.

Test case 2: ADM with a 2D homogeneous fractured reservoir
In this test case, the accuracy of the ADM method for mass-heat transport in a homogeneous domain containing fractures with mixed conductivities is studied. For this purpose, a 2D fractured 100 m ×100 m domain with 30 fractures is considered. The length of each fracture is different but the size of their aperture is identical and set to a f = 5 · 10 −3 m. A 136 × 136    The coarsening level for matrix is l m = 2 and for fractures is l f = 1. Over the entire computational grids, the coarsening ratio is γ = 3 × 3. The ADM simulations are run for three different thresholds of coarsening criterion, namely, where, x is either pressure or temperature and the subscripts ADM and FS denote ADM and fine-scale, respectively. The simulation is run for 1000 [days] and the results are printed at 100 isochronal intervals.  Fig. 18 provides more details of the error analysis for test case 2. On the top row pressure and temperature errors ( Fig. 18a and 18b) versus simulation time-steps are shown. Fig. 18d plots the percentage of active grid cells used during the simulation. Lastly, averaged errors and average active grid cells are plotted over the entire simulation for all three thresholds.
Please note that the regions around the wells are always kept at fine-scale resolution and as a result (depending on the size of the domain) a noticeable portion of active grid cells are imposed. In this test case, approximately 5 percent of total grid cells are kept at fine-scale due to near-well refinement.
Please note that in this test case, intersection between permeable fractures and impermeable fractures occur. By default, the impermeable fractures are considered as "non-sealing" flow barriers, meaning that the crossing of the highly permeable   Fig. 18a and 18b show ADM error of pressure and temperature over simulation time-steps. The errors are calculated based on Eq. (47). Fig. 18c shows the percentage of active grid cells used during simulation versus the time-steps. Lastly, Fig. 18d demonstrates the averaged values of the mentioned parameters over the entire simulation time for three thresholds. fractures cause flow leakage. This effect can be seen in the Fig. 16 where the pressure gradient across the permeable fracture stays negligible on both sides of the impermeable fracture.

Test case 3: ADM with a 2D heterogeneous fractured reservoir
This test case, assesses the ADM performance on a 2D heterogeneous and fractured domain. Similar to the previous test case, a 2D domain with dimensions of 100 m × 100 m is considered. A network of 30 fractures with identical geometry as in test case 2 exists. The rock matrix is discretized into 136 × 136 fine-scale grid cells and 291 grids are imposed on fracture network (in total 18516 grid cells). The permeability of the matrix ranges from K m min = 1.2 × 10 −15 m 2 to K m max = 1.2 ×10 −12 m 2 . All fractures are highly conductive with permeability of K f = 10 −8 m 2 . The well pattern used in this test case is identical to test case 2. Moreover, the same coarsening strategy is used in this test case. The input parameters of Table 2 holds for this test case Fig. 19 and 20 respectively. Similar to the previous test case, the errors (e) and the percentage of active grid cells ( AGC ) are indicated under the ADM plots. Fig. 21 shows the error analysis for test case 3. Pressure and temperature errors versus simulation time-steps are illustrated on the top row ( Fig. 21a and 21b). Fig. 21c shows the percentage of active grid cells during the simulation. Average errors and active grid cells over the entire simulation time can be seen for all four ADM thresholds on Fig. 21d.

Test case 4: ADM with a 3D homogeneous fractured reservoir
In this test case, a 3D 100 m × 100 m × 40 m containing 15 lower dimensional fractures with different geometrical properties is considered. An 81 × 81 × 27 Cartesian grid is imposed on rock matrix and fracture network is discretized in to 1377 grid cells (summation of 178524 grid cells). The rock matrix has permeability of K m = 10 −14 m 2 . Fracture network consists of both highly conductive fractures with permeability of K f = 10 −6 m 2 and flow barriers with permeability of K f = 10 −22 m 2 . Fig. 23a shows the structure of fracture network for this 3D test case. Two injection wells exist on the bottom left and top left boundaries with pressure of p inj = 5 · 10 7 Pa. Similarly, two production wells are located at the bottom right and top right boundaries with pressure of p prod = 1 · 10 7 Pa. All wells are vertical and perforate the entire thickness of the reservoir. Simulation is run for 200 days and the results are printed at 100 isochronal intervals. Two coarsening levels with coarsening ratio of γ m = 3 × 3 × 3 is considered for matrix and one coarsening level with coarsening ratio of γ f = 3 × 3 for fractures.

Conclusions
An algebraic dynamic multilevel (ADM) method for fully-implicit simulation of single-phase thermal flow in fractured heterogeneous porous media using pEDFM is presented. First, the pEDFM model [27] is extended to account for 3D fracture geometries. The fine-scale fully implicit system resulting from geothermal-pEDFM model is then mapped into a dynamic multilevel (i.e., ADM) grid, defined for matrix and fractures independently. Such map is obtained using sequences of ADM prolongation and restriction operators which are assembled from the static multilevel multiscale operators. The static operators are calculated only at the beginning of the simulation after the construction of the coarse grids for both matrix and fractures on all coarsening levels. The prolongation operator is constructed based on locally computed basis functions. Important to note is that the matrix-fracture coupling is taken into account for multiscale basis functions. The ADM grid selection strategy with front tracking technique employs the fine-scale grid cells only where and when needed.
Fine-scale simulation for a 2D geothermal system allows to investigate the pEDFM and EDFM accuracy, both compared to the direct numerical solution (DNS). Both highly conductive and flow-barrier fractures are considered. It is shown that in the presence of highly conductive fractures, pEDFM and EDFM performed accurate, though EDFM still allows for small unphysical leakage. However, EDFM fails to capture the low conductive fractures, while pEDFM shows a good representation of DNS with satisfactory accuracy.
Moreover, numerical results are obtained for 2D and 3D geothermal reservoirs to investigate ADM performance. The ADM results on these test cases are compared to those obtained from fine-scale simulations. The sensitivity of ADM to different grid coarsening criteria is also studied. It is observed that ADM is capable of providing accurate simulations by employing only a fraction of the fine-scale grid cells in the sub-domains where it is needed. Due to the rarefaction of the temperature distribution (highly diffused temperature at the front), more fine-scale grid cells are used at the temperature front depending on the threshold. Additionally, ADM provides an algebraic framework which brings a scalable simulation method for thermal fluid flows. One can expect that by increasing the size of the domain, the percentage of active grid cells used during simulation reduces. Therefore, ADM introduces a promising simulation strategy for real-field geothermal reservoir simulations.
All software developments of this work is made available open source at https://gitlab .com /DARSim2simulator.

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.