Heterogeneity Preserving Upscaling for Heat Transport in Fractured Geothermal Reservoirs

TL;DR

This paper introduces a novel upscaling method for simulating heat transport in fractured geothermal reservoirs, effectively capturing heterogeneity and complex fracture networks while reducing computational costs.

Contribution

The paper presents a new heterogeneity-preserving upscaling technique that explicitly models fractures in coarse grids for geothermal heat transport simulations.

Findings

Method performs well for large upscaling ratios.

Applicable to highly complex fracture networks.

Reduces computational effort while maintaining accuracy.

Abstract

In simulation of fluid injection in fractured geothermal reservoirs, the characteristics of the physical processes are severely affected by the local occurence of connected fractures. To resolve these structurally dominated processes, there is a need to develop discretization strategies that also limit computational effort. In this paper we present an upscaling methodology for geothermal heat transport with fractures represented explicitly in the computational grid. The heat transport is modeled by an advection-conduction equation for the temperature, and solved on a highly irregular coarse grid that preserves the fracture heterogeneity. The upscaling is based on different strategies for the advective term and the conductive term, respectively. The coarse scale advective term is constructed from sums of fine scale fluxes, whereas the coarse scale conductive term is constructed based on numerically computed basis functions. The method naturally incorporates a coupling between the matrix and the fractures via the discretization, so that explicit transfer terms that couple solution variables in the fractures and the matrix are avoided. Numerical results show that the upscaling methodology performs well, in particular for large upscaling ratios, and that it is applicable also to highly complex fracture networks.

Tables

Table 1. Table 1: Accuracy with respect to coarsening factor for the refinement study in section 4.1 . 𝒆 𝑪 ​ 𝑩 subscript 𝒆 𝑪 𝑩 e_{CB} and 𝒆 𝑺 ​ 𝑩 subscript 𝒆 𝑺 𝑩 e_{SB} denote the ℓ 𝟐 subscript bold-ℓ 2 \ell_{2} energy error over space for constant and smoothed basis functions, respectively.

No. fine cells	Coarsening factor $𝑪𝑭$	${𝒆}_{𝑪𝑩}$	${𝒆}_{𝑺𝑩}$ (no. of iterations)
6889	11	$1.06\cdot {\mathrm{𝟏𝟎}}^{-\mathrm{𝟐}}$	$1.06\cdot {\mathrm{𝟏𝟎}}^{-\mathrm{𝟐}}\left(\mathrm{𝟏}\right)$
26569	44	$3.00\cdot {\mathrm{𝟏𝟎}}^{-\mathrm{𝟐}}$	$1.63\cdot {\mathrm{𝟏𝟎}}^{-\mathrm{𝟐}}\left(\mathrm{𝟓}\right)$
104329	174	$5.99\cdot {\mathrm{𝟏𝟎}}^{-\mathrm{𝟐}}$	$1.90\cdot {\mathrm{𝟏𝟎}}^{-\mathrm{𝟐}}\left(\mathrm{𝟏𝟓}\right)$
413449	688	$8.98\cdot {\mathrm{𝟏𝟎}}^{-\mathrm{𝟐}}$	$2.15\cdot {\mathrm{𝟏𝟎}}^{-\mathrm{𝟐}}\left(\mathrm{𝟑𝟎}\right)$

Equations

$$\nabla \cdot \text{v}=q,$$

$$\text{v}=-\frac{\text{K}}{\mu }\nabla p,$$

$$-\nabla \cdot \left(\frac{\text{K}}{\mu }\nabla p\right)=q,$$

$$k_f=\frac{a^2}{12}.$$

$$(\rho c_p{)}_{eff}=\varphi (\rho c_p{)}_f+(1-\varphi )(\rho c_p{)}_r,$$

$$(\rho c_p{)}_{eff}\frac{\partial T}{\partial t}+(\rho c_p{)}_f\nabla \cdot (\text{v}T)-\nabla \cdot (C\nabla T)=q_e.$$

$$Pe=\frac{L\cdot u\cdot (\rho c_p{)}_{eff}}{\alpha },$$

$${\int }_{\partial {\omega }_i}\text{v}\cdot \text{n}dS={\int }_{{\omega }_i}qdV,i=1,\cdots N_f,$$

$$(\rho c_p{)}_{eff}{\int }_{{\omega }_i}\frac{\partial T}{\partial t}dV-(\rho c_p{)}_f{\int }_{\partial {\omega }_i}(\text{v}T)\cdot \text{n}dS$$

$$+{\int }_{\partial {\omega }_i}(C\nabla T)\cdot \text{n}dS={\int }_{{\omega }_i}{q}_{e}dV,i=1,\cdots N_f.$$

$${\text{A}}_f{\text{p}}_f=\text{q},$$

$$\frac{d{\text{T}}_f}{dt}+({\text{A}}_{f,conv}+{\text{A}}_{f,cond}){\text{T}}_f={\text{q}}_e,$$

$${\int }_{\partial {\Omega }_{\ell }}(\text{v}T)\cdot \text{n}dS\approx$$

$$\sum _{k\ne \ell }\max \left(T_{\ell }\sum _{{\gamma }_{ij}\in {\Gamma }_{k\ell }}{\text{v}}_{ij}-T_k\sum _{{\gamma }_{ij}\in {\Gamma }_{k\ell }}{\text{v}}_{ij}\right),$$

$${\text{A}}_{c,cond}=\text{R}{\text{A}}_{f,cond}\text{P},$$

$$\displaystyle\textbf{R}_{i,j}=\left\{\begin{tabular}[]{cc}1&if $\omega_{j}\in\Omega_{i}$,\\ 0&otherwise,\end{tabular}\right.$$

$$\text{P}=\left[P_1{P}_2\cdots P_{N_c}\right],$$

$$\displaystyle\textbf{P}_{i,j}^{(0)}=\left\{\begin{tabular}[]{cc}1&if $\omega_{i}\in\Omega_{j}$,\\ 0&otherwise,\end{tabular}\right.$$

$$P_i^{(n+1)}=P_i^{(n)}-\omega {\text{D}}^{-1}{\text{A}}_{f,cond}{P}_i^{(n)},$$

$$E_i=P_i^T{\text{A}}_{f,cond}{P}_i.$$

$$E_i^{(n+1)}>E_i^{(n)},$$

$$e_i(t)=\left[\varphi (\rho c_p{)}_f+(1-\varphi )(\rho c_p{)}_r\right]V_i{T}_i(t)$$

$$=(\rho c_p{)}_{eff}{V}_i{T}_i(t),$$

$$\epsilon (t)=\frac{\sqrt{{\sum }_{i=1}^{N_f}|{\tilde{e}}_{c,i}(t)-e_{ref,i}(t){|}^2}}{\sqrt{{\sum }_{i=1}^{N_f}|e_{ref,i}(t){|}^2}},$$

Full text

††journal:

1 Introduction

Geothermal reservoirs are typically situated in igneous rocks, where the permeability of the reservoir mainly is governed by discrete fractures. The fractures occur on a range of scales, where large-scale connected structures tend to completely dictate flow paths. Flow and transport processes are strongly dominated by these structural heterogeneities. While large-scale fractures dictate preferential fluid pathways, fine-scale fractures contribute significantly to the subsurface heat exchange from the rock to the brine. The size of the fractures in a fractured reservoir typically follow a power law distribution, which means that at any scale of investigation, flow will be controlled by single fractures that cannot be represented by an upscaled permeability BonnetEtAl2001 . Due to the lack of scale separation in fractured formations, particular care must be taken to appropriately account for local flow and transport characteristics.

Modeling approaches for fractured reservoirs can be classified into three main categories based on their spatial representation in the considered formation DietrichEtAl2005 : single continuum models, multi-continuum models, and discrete fracture-matrix models. Single continuum models average spatial properties of the rock over grid cells that are much larger than the narrow width of the structural heterogeneities (fractures), assuming that the concept of representative elementary volume is valid. An advantage of this approach is that standard reservoir simulators can be applied. In multi-continuum models, the different characteristics of flow in different structural components are preserved for each region of rock; see, e.g., Berkowitz2002 and references therein. Each structural component, e.g. fractures on different scales and matrix, is modeled by a representative continuum, and interacts with the other continua comprising the same region. In this way integral transport behavior can be better captured compared to when using a single continuum model. In contrast, discrete fracture network (DFN) models represent fractures explicitly and solely considers flow in the network of explicily represented fractures; see, e.g., WatanabeTakahashi1995 ; BruelCacas1992 ; HayashiEtAl1999 . Discrete fracture-matrix (DFM) models conceptually represent a combination of DFN models and multi-continuum models, with explicitly represented fractures between regions of continuity (e.g., DietrichEtAl2005 ; ReichenbergerEtAl2006 ; MartinEtAl2005 ), and is the approach we consider here. The advantage of DFM models is that they allow explicit representation of fractures that cannot be represented by an equivalent continuum at the chosen scale in simulation models, while at the same time account for (upscaled) permeability due to small scale fractures and pores in the regions surrounding the fractures. Therefore, the DFM models account for the lack of scale separation by modelling fractures that dominate flow characteristics at a certain scale explicitly, while at the same time accounting for the region outside these fractures being permeable. This is of particular importance in geothermal heat transport, where the permeability of the region surrounding the fractures with dominating flow plays an important role for heat transport.

Discretization methods for DFM models are challenging due to the large difference in scales between fractured and non-fractured regions. Co-dimension one approaches treat fractures as lower dimensional objects in the geometrical domain, and associate the removed dimension with a hydraulic aperture in the computations. This simplifies modelling at the cost of a small volumetric error, see e.g. Karimi-Fard2004 ; SandveEtAl . The explicit representation of fractures in DFM models necessitate large flexibility in terms of grid generation. Frameworks that slightly modifies the fracture network to align with a background grid have been developed to avoid excessive grid refinement close to e.g. fracture intersections MusthaphaMusthapha2007 ; Karimi-FardDurlofsky2016 . In addition to gridding and discretization challenges, numerical representation using DFM models typically leads to a high number of degrees of freedom compared to using continuum models. Efficient numerical methods are therefore important.

In this paper we present a heterogeneity-preserving upscaling method for advective-conductive heat transport in fractured reservoirs. Our main contribution is two-fold: First, we introduce an upgridding method that is tailored to drainage of energy from the matrix to the fracture system. Second, we study discretization schemes for heat transport on the resulting coarse grid. Following Hauge2012 , the coarse grid used for heat transport is based on a fine-scale velocity field, which is assumed to be known. The coarse cells are constructed by merging fine-scale cells with similar flow properties and distance from the fracture network. The resulting coarse grid is well suited to represent both the slow drainage of heat from matrix to the fracture network, and the fast transport within the network. We emphasize that the user interacts with the coarsening algorithm only via a few meta-parameters.

Depending on the geometry of the fracture network, the coarse grid can have cells with highly irregular shapes, including non-convexity and high aspect ratios. While the advective term can be discretized from the known fine-scale velocity field in the manner of Hauge2012 , the complex grid poses challenges for the discretization of the heat conduction term. To overcome this, we use ideas from multiscale methods developed for elliptic and parabolic equations. Multiscale methods have been used extensively in porous media applications to efficiently discretize the pressure equation in a way that accurately represents fine-scale heterogeneities, see, e.g., Jenny2003 ; ZhouTchelepi2008 ; Hajibeygi2011 . Due to the emphasis on material heterogeneities, the majority of the work on multiscale methods applied to porous media problems has with a few exceptions been focused on structured grids. The framework presented in Moyner2016 , on the other hand, makes few assumptions with regards to grid geometry, and uses an algebraic procedure to compute multiscale basis functions. This makes the framework in Moyner2016 well suited for general coarse grids, and we therefore use it as a basis for discretization of the coarse-scale conductive term. However, the highly irregular shape of the coarse cells makes it necessary to introduce modifications to avoid numerical instabilities. We also study a non-conforming discretization of the coarse advection term, i.e. based on piecewise constant basis functions. This second approach can be expected to have inferior approximation properties, but it does provide a simple and robust alternative to the multiscale-based discretizations.

To study the performance of our methods we consider fracture network ranging from a Cartesian structure to a highly complex network that consists of a large number of stochastic fractures. We then study the accuracy of the methods under varying upscaling ratios and fluid injection rates. The numerical examples show that both discretizations capture the main features of the transport, compared to a fine-scale simulation. As expected, the multiscale based method gives the most accurate predictions, in particular in regions where conduction dominates.

The outline of the paper is as follows. In section 2 we describe the governing equations for flow and transport in geothermal reservoirs and the discretization on the fine scale. The framework for coarse scale discretization is presented in section 3. Numerical experiments are presented in section 4 and the paper ends with conclusions and a discussion in section 5.

2 Fine Scale Model

In this section we present the governing equations for flow and heat transport. We also derive the fine scale discretization that later will be the starting point for our upscaling methodology.

2.1 Flow and Heat Transport Equations

We consider single phase incompressible flow. Mass conservation of the fluid,

[TABLE]

in combination with Darcy’s law,

[TABLE]

is combined into an elliptic equation for the pressure, $p$,

[TABLE]

where K is the permeability tensor, $\mu$ is the viscosity, which is assumed to be constant, v is the Darcy flux field and $q$ contains any sources or sinks. The relation between hydraulic aperture, $a$, and fracture permeability, $k_{f}$, is modelled as

[TABLE]

The effective volumetric heat capacity for rock and fluid combined is given by

[TABLE]

where subscripts $f$ and $r$ denote fluid and rock, respectively, and $eff$ denotes effective. Here, $\rho$ is the density, $c_{p}$ is the specific heat capacity, and $\phi$ is the porosity. We assume that the specific heat capacities of fluid and rock varies slowly with time. Assuming local thermal equilibrium for a single phase flow, the energy equation that describes heat tranfer through the reservoir can be simplified to a linear advection-conduction equation for the temperature, $T$, given by

[TABLE]

In addition, $q_{e}$ is the energy source term, and $C$ is the effective thermal conductivity for the porous medium saturated with fluid at local thermal equilibrium. For simplicity, we will assume a constant value for conductivity. The assumptions made above are reasonable for water-filled geothermal reservoirs, including a wide range of temperatures for enhanced geothermal systems (EGS) situated at large depths. In deep reservoirs water remains a compressed liquid for high temperatures due to high hydrostatic pressures Cheng79 . We obtain a decoupled system of equations (3) and (6) that together with appropriate boundary conditions can be solved sequentially in numerical simulations. The coupling is one-way via the Darcy flux, v, due to the assumptions of an incompressible fluid with constant viscosity, and the pressure equation is only solved once, at the start of the simulation.

The ratio between heat advection and conduction is commonly characterized by the heat Péclet number $Pe$, defined as

[TABLE]

where $L$ is a representative length scale, typically the distance between injector and producer wells, $u$ is an average velocity, and $\alpha$ is the average heat conductivity. The flow injection rate generally determines the flow rate through the fracture network and subsequently the highest value of $Pe$ in the reservoir. Depending on e.g. the matrix-fracture permeability contrast, the heat advection in the matrix can be considerably slower compared to that in the fractures leading to a second response from advection through the matrix Pruess1983 . As an effect the values of $Pe$ can, depending on the local flow field, vary by orders of magnitude in different regions of the reservoir.

2.2 Fine Scale Discretization

We use an unstructured computational grid, constructed as a triangulation constrained to the fractures, with a slightly modified version of the algorithm presented in Holm2006 . Equations (3) and (6) are discretized using a finite volume discretization. Let $\omega$ be a set of fine grid cells, $\omega$ = {$\omega_{i}\},i=1,\cdots,N_{f}$, where $N_{f}$ is the total number of fine cells. Integrating equation (3) over each grid cell $\omega_{i}$ and using the divergence theorem gives

[TABLE]

where $\partial\omega_{i}$ is the boundary of grid cell $\omega_{i}$ and n is the outward normal of $\partial\omega_{i}$. Similarly, equation (6) becomes

[TABLE]

A two-point flux approximation is used to discretize both equation (8) and the conductive part of equation (9). The fractures are incorporated by the hybrid discretization developed in Karimi-Fard2004 , where the permeability heterogeneity is preserved by the introduction of cells along the fractures. The advective fluxes in equation (9) are approximated using standard upstream weighting AzizSettari1970 .

The discretization of equation (3) gives rise to a linear system on the form

[TABLE]

where $\textbf{p}_{f}$ is a vector containing cell-centered pressure values. Equation (10) can be solved either via a direct solver or by using an appropriate iterative solver, see, e.g., SandveWRR . Similarly, after discretizing equation (9) in space, the semi-discrete temperature equation reads

[TABLE]

where $\textbf{T}_{f}$ is a vector containing cell-centered temperature values, $\textbf{A}_{f,conv}$ and $\textbf{A}_{f,cond}$ approximate the advective and the conductive terms in eq. (9), respectively, and $\textbf{q}_{e}$ approximate the right hand side. $\textbf{A}_{f,conv}$, $\textbf{A}_{f,cond}$, and $\textbf{q}_{e}$ are all scaled with the inverse of the effective heat capacity.

For time propagation we use an implicit second order accurate backward differentiation formula (BDF2), which has favorable stability properties, with time steps sufficiently small for temporal diffusion not to pollute the simulation accuracy.

3 Heat Transport Upscaling

The fine scale discretization with explicitly represented fractures in the computational grid is often too computationally expensive to use directly in the heat transport simulations. As an alternative approach, the fine scale discretization can be used as the basis for numerical upscaling, where most of the computations are carried out on a coarser grid. In this section, we introduce an upscaling method for simulation of heat transport in fractured media. The upscaled transport model requires the construction of a coarse grid, together with a discretization of equation (9) on that grid. We assume that a fine-scale velocity field is available, so that fine scale heterogeneities in the flow model can be incorporated into the coarse scale transport discretization. If the coupling between flow and transport is assumed to be weak, a single fine scale pressure solve may be sufficient, in which case the pressure solve will not constitute a large part of the overall simulation cost. Alternatively, iterative solvers or multiscale methods tailored for fractured media can be applied to effectively produce conservative approximations to the velocity field on the fine scale, see for example SandveEtAl ; Shah2016 ; Tene2015 . Because of the assumption of incompressible flow, we use a direct fine scale pressure solver in this work. Since no additional errors from an approximate velocity field are introduced, errors due to the coarse scale transport solver can be isolated, and the performance of the method can be thoroughly investigated.

3.1 Coarse Scale Grid Construction

To design adequate coarse grids, insight into the behavior of transport processes in fractured media is needed. Global transport is often dominated by flow in the fractures, while heat stored in the matrix is transported into the fracture network by a combination of matrix flow and conduction. This suggests two criteria for the construction of coarse grids. First, the coarse grid should honor the structure of the flux field to adapt to large scale connectivity, i.e. large scale fractures. Specifically, we have found it useful to apply the flow-based indicator framework introduced by Aarnes et al. AarnesHaugeEfendiev2007 , and to use time-of-flight (TOF) as indicator function for the coarse grid. Second, the assumption that heat drains from the matrix to nearby fractures suggests that the coarse cell geometry in the matrix should reflect the distance to nearby fractures. We achieve this by classifying fine scale cells located in the matrix according to their geometric distance from the closest fracture and refining the coarse grid according to the distance measure. It is instructive to draw the parallel between this second criterion and the so-called multiple interactive continuum (MINC) models Pruess1985 . Both models use distance from fractures as a proxy for the time it takes for heat to drain towards the fractures. A major difference is that since our starting point is a fine-scale grid, estimation of interaction between coarse cells is relatively simple even for general fracture geometries. Moreover, our approach can readily be combined with heat advection.

We refer to the mapping from the fine grid to the coarse grid as a coarse grid partition, i.e. the description of which fine grid cells constitute each coarse grid cell. The different steps in the construction of the coarse grid are visualized in figure 1 and briefly described in the following steps:

Generate an initial coarse grid partition of a fine grid based on a coarsening indicator, e.g., based on TOF. 2. 2.

Refine the initial coarse grid using a distance based grid partition. This partition is defined by a predefined coarse grid resolution parameter and a distance indicator (measuring the distance between the cell center of a fine cell and the closest fracture). 3. 3.

Hybrid coarse cells (consisting of both fine fracture cells and fine matrix cells) are separated so that all coarse cells either contain only fine fracture cells or fine matrix cells. 4. 4.

Coarse cells are ensured to be (or split into) connected components.

The coarse cells can be merged and refined in an iterative loop in order to obtain a coarse grid with preferred properties Hauge2012 . The focus of this work is not to obtain an optimal coarse grid, but rather to design a coarse scale discretization that works well also for highly irregular coarse grids. We remark that the coarse grid in figure 1 gives a fair representation of the grid cell shapes we have encountered. For transport problems dominated by fracture flow, the number of grid cells can be vastly reduced and still capture large scale behavior for a coarse grid constructed based on flow indicators compared to when a more structured coarse grid is used. However, with a more irregular coarse grid the discretization on the coarse scale needs to be carefully constructed. This is described in the following sections.

3.2 Coarse Scale Discretization

Let $\Omega$ be a set of coarse grid cells, $\Omega$ = {$\Omega_{\ell}\}$, $\ell=1,\cdots,N_{c}$, where $N_{c}$ is the total number of coarse cells. The coarse scale finite volume discretization takes the same form as the fine scale discretization in equation (9), with integration over a coarse cell $\Omega_{\ell}$.

In the following, we detail the coarse scale discretization of the advective and conductive terms, respectively. As we will see, the former can be treated relatively easily, while the latter requires more attention.

3.2.1 Coarse Scale Advective Term

The treatment of the coarse scale advective term assumes that fine scale fluxes are known. Following Hauge2012 , a control volume discretization of the advective term is obtained by integrating the (known) fine scale fluxes to net fluxes on the coarse scale,

[TABLE]

where $\gamma_{ij}$ denotes a face between fine grid cells $i$ and $j$, and $\Gamma_{k\ell}$ a face between coarse grid cells $k$ and $\ell$.

3.2.2 Coarse Scale Conductive Term

For the coarse scale conductive term we take advantage of the fact that the same structure is found in elliptic and parabolic pressure equations. In porous media applications, coarse scale discretizations of elliptic equations with strongly variable coefficients (related to heterogeneous and anisotropic media) has received much attention the last decade in the context of multiscale methods Jenny2003 ; AarnesEfendiev2006 ; ZhouTchelepi2008 ; Hajibeygi2008 . Due to the strong hyperbolic component of the heat transport, a projection of the full temperature profile from the coarse to the fine scale is highly challenging, and we therefore cannot directly apply a multiscale framework to the advection-conduction equation. Nevertheless, we construct the upscaled discretization of the conductive term inspired by multiscale metods.

The coarse scale finite volume discretization matrix for the conductive term, $\textbf{A}_{c,cond}$, is related to the corresponding fine scale discretization matrix, $\textbf{A}_{f,cond}$, by

[TABLE]

where R is a restriction operator represented by a matrix of size $N_{c}\times N_{f}$, and P is a prolongation matrix of size $N_{f}\times N_{c}$. The restriction operator, R, maps the fine scale discretization to the coarse scale. We use a finite volume restriction matrix, with elements given by

[TABLE]

see ZhouTchelepi2008 . P can be expressed as a column matrix,

[TABLE]

where the columns $P_{i}$, $i=1,\cdots,N_{c}$, contain prolongation operators or basis functions.

A coarse scale discretization on the form (13) is often used for multiscale methods for elliptic equations with variable coefficients, see for instance ZhouTchelepi2008 ; Moyner2016 . In that case, the basis functions can be combined to directly map the coarse scale solution variables to a fine scale approximation. For elliptic equations that stems from strongly heterogeneous and anisotropic porous media problems, a main difficulty with respect to multiscale methods is to define basis functions that give a numerical solution with desired properties. More generally, the construction of the prolongation matrix P is a long-standing issue in the field of coarse discretizations, and it is also closely related to the formulation of coarse spaces in the domain decomposition and multigrid literature Toselli2005 .

It is a common feature of many multiscale methods and numerical upscaling methods that they are most easily realized on relatively structured grids. As exemplified in figure 1, the grid coarsening strategy presented in Section 3.1 can give rise to highly irregular coarse grids. To define a discretization of the conductive term on these grids based on geometric arguments is very challenging. Instead we consider an algebraic construction of the basis functions in the coarse discretization presented in Moyner2016 , which makes the method much more amenable to complex coarse grids. The algebraic construction is closely related to ideas in smoothed aggregation algebraic multigrid methods Vanek1996 ; Nordbotten2008 . In the following subsections, we describe how the methodology in Moyner2016 is applied and modified to obtain a suitable discretization of the conductive term in our advection-conduction equation.

3.2.3 Interaction Regions

The region of support for a basis function is defined in terms of interaction regions. These regions resemble discretiztion stencils for multipoint flux approximations Aavatsmark2002 in the sense that all coarse cells that share a vertex with coarse cell $i$ contain some fine cells that are included in $P_{i}$. Conceptually, we need to consider three different types of interactions:

1. matrix cells neighboring matrix cells,
2. matrix cells neighboring fracture cells, and
3. fracture cells neighboring fracture or matrix cells. Interactions regions for these three types are illustrated in figure 2, together with the corresponding coarse cells and basis functions for these cells, respectively. The interaction regions are constructed by the algorithm in Moyner2016 , with some adaptions to account for the presence of fracture cells. Specifically, we do not allow for interaction between coarse matrix cells that are separated by a fracture cell. This is due to the assumption of advective-dominated heat transport in the fractures.

3.2.4 Algebraic Smoothing

Having defined the support of the basis functions in terms of the interaction regions, we go on to explicitly compute the basis functions. We use the algebraic construction from Moyner2016 , described briefly below. The initial description of the elements in the prolongation matrix P is

[TABLE]

i.e. $\textbf{P}^{(0)}=\textbf{R}^{T}$. We note that $\textbf{P}^{(0)}$ is a partition of unity, and forms a $P_{0}$ (piecewise constant) discretization on the coarse grid.

Now, successive iterations for basis function $P_{i}=\textbf{P}_{\cdot,i}$ are given by

[TABLE]

where $\textbf{D}=diag(\textbf{A}_{f,cond})$, $\omega\in(0,1]$ is a relaxation parameter, and $n$ denotes the iteration number. The effect of the iterations is to smoothen the basis functions. The support of a basis function is limited to its interaction region, this is enforced after each iteration by simply truncating values that fall outside the interaction region. However, the truncated basis functions will generally not form a partition of unity. As a remedy, the functions are rescaled cell-wise to ensure that they sum to one within each fine scale cell Moyner2016 .

3.2.5 Non-Oscillatory Basis Functions

The above construction of basis functions works well for coarse cells with a reasonably regular geometry. However, we have observed that when the coarse grid becomes highly irregular, as illustrated in figure 1, the basis functions may have a non-monotonic behavior, which consequently lead to oscillations in the numerical solution. To see why the non-monotonic profiles occur, we note that for cells with large aspect ratios or with complex shapes, interaction regions associated with neighboring cells are not necessarily in contact with the cell center. Thus the rescaling to preserve a partition of unity, which is carried out for each interaction region individually, can cause the basis function not to be monotonically decreasing from the cell center. In this and the following subsection we discuss two approaches to mitigate this effect.

As a partial remedy to the oscillations we limit the number of iterations in the construction of basis functions by minimizing an energy functional vanLent2009 . Specifically, we define $E_{i}$ as the discrete energy of the basis function for coarse grid cell $i$, $\Psi_{i}$ by

[TABLE]

Oscillations in the basis function will increase the discrete energy, and we therefore terminate the iterations for basis function $P_{i}$ if for some iteration $n$ we have

[TABLE]

while continuing to update basis functions that have not been terminated. The termination freezes the basis functions in all fine-scale cells contained within the interaction region of coarse cell $i$. While the termination criterion (23) does not guarantee monotone basis functions, the approach has been sufficient for the fracture networks considered here. In addition to the termination condition (23) for each basis function, global criteria for terminating the iterations are based on the size of the residual and by a user given total number of iterations. Figure 3 shows the distribution of basis functions terminated at different numbers of iterations for the fracture network considered in section 4.2. When the smoothening is terminated based on the residual after $111$ iterations, 24$\%$ ($440/1852$) of the basis functions are still updating, and they are therefore not included in the figure.

Remark In addition to non-monotonic basis functions, excessive smoothing may also result in eigenvalues of $\textbf{A}_{c,cond}$ with negative real parts. This consequently leads to blow up of the numerical solution as time propagates. If needed, this can be avoided by introducing an additional criteria on the diagonal elements of $\textbf{A}_{c,cond}$ such that $\textbf{A}_{c,cond,ii}>0$, $i=1,\cdots,N_{c}$. If the diagonal elements of $\textbf{A}_{c,cond}$ are strictly positive, then by the Gershgorin circle theorem the eigenvalues have real parts that are strictly positive.

3.2.6 Piecewise Constant Basis Functions

As noted above, oscillations in the basis functions can to a large degree be avoided by reducing the number of smoothing iterations. It is of particular interest to investigate the limit option of taking no iterations at all, and thus represent the coarse-scale temperature with piecewise constant basis functions. This approach avoids the oscillations, and as the interaction regions are no longer required, the method is simpler to implement. Moreover, it leads to a symmetric discretization of the coarse conduction term, assuming that the fine scale discretization matrix, $\textbf{A}_{f,cond}$, is symmetric, since in this case $\textbf{P}=\textbf{P}^{(0)}=\textbf{R}^{T}$, and thus $\textbf{A}_{c,cond}=\textbf{R}\textbf{A}_{f,cond}\textbf{R}^{T}$. However, the piecewise constant representation of temperature implies that temperature gradients can only be present on the boundary between the coarse-scale cells, where they are represented by the fine-scale discretization $\textbf{A}_{f,cond}$. That is, the coarse-scale conduction retains a dependency on the fine scale grid size $h$, rather than on a representative coarse grid size $H$.

As the above considerations show, the constant basis function will therefore lead to inferior approximations, in particular for large upscaling ratios. The effect will be most pronounced in regions where conduction is the dominant transport mechanism. Nevertheless, the approach is interesting, in that it is considerably simpler and more numerically robust than the smoothed basis functions.

3.2.7 Relation to Other Methods

It is of interest to discuss differences with related methods for coarse discretizations in fractured media. Recently, Tene et al. Tene2015 and Shah. et al Shah2016 developed multiscale methods for solving the pressure equation in fractured media. In both these works the fracture networks and the matrix are separated into two independent coarse grids, and exchange between fractures and matrix are modelled via a Peaceman-type coupling term that is proportial to the difference between the matrix and fracture pressures Hajibeygi2011 . While this approach can be advantageous for the pressure equation with high permeability contrast between fracture and matrix, it is less relevant for heat conduction that has no abrubt jumps in the parameter field.

4 Numerical Experiments

In this section we consider three different test cases to illustrate the applicability and performance of the upscaling framework. Numerical experiments are run both with a smoothed basis and with constant basis functions, respectively. First, we present a refinement study for a structured fracture network on an underlying Cartesian fine scale grid to investigate the practical importance of sufficiently smooth basis functions. We follow by investigating the performance of the numerical method for two more complex fracture networks using unstructured triangular fine scale grids.

In the numerical simulations we use standard constant values for granite for the volumetric heat capacity of the rock and the thermal conductivity, $(\rho c_{p})_{r}=2170$ kJ/(m3K) and ${C}=2.1$ J$/$(mKs), respectively. For the volumetric heat capacity of the fluid and viscosity, $(\rho c_{p})_{f}=4180$ kJ$/$(m3K), and $\mu=1$ cP are used, respectively. The porosity in the matrix is set to $\phi=0.001$. The reservoir domain is $\Omega=[0,1000]\times[0,1000]$ m2 for all simulations. Initially the temperature in the reservoir is set to 100°C, and then water with temperature 20°C is injected somewhere in the reservoir. The precise injection and production setups as well as the fracture networks are described for each case in subsections 4.1, 4.2, and 4.3, respcetively.

The total internal energy per grid cell at time $t$ is given by

[TABLE]

where $i$ refers to the grid cell index. The relative $\ell_{2}$ error over space in the energy at time $t$, $\varepsilon(t)$, is then computed as

[TABLE]

where $\tilde{e}_{c}$ is the energy vector for the upscaled solution projected on the fine grid, given by $\tilde{e}_{c}=\textbf{R}^{T}e_{c}$, and $e_{ref}$ is the energy vector for a reference solution computed on the fine grid. $N_{f}$ refers to the total number of fine grid cells.

4.1 Accuracy with Respect to Coarsening Ratio

We conduct a refinement study using a fine underlying cartesian grid to better understand the practical importance of the iterative framework. The fracture network consists of six fractures that are aligned with the grid, three in each direction. First, a coarse mesh is constructed based on a partitioning of the fine grid. Because of the structured nature of this fracture network we base the initial coarse grid partition solely on the distance to nearest fracture. Using a flow based partition leads to similar results. The upscaling ratio, or coarsening factor between the coarse and the fine grid is defined as $CF=N_{f}/N_{c}$. The coarsening factor is consecutively increased by refining the fine grid by a factor two in each spatial direction, while the same coarse grid is used. The coarse grid contains $601$ grid cells and the fine grids contain $6889$, $26569$, $104329$, and $413449$ grid cells, respectively, including one-dimensional fracture cells. Water with temperature 20°C is injected at a rate of 0.1 $dm^{2}/s$ in the middle of the domain (at $(500,500)$), and extracted where each fracture meets the exterior boundaries. The total simulation time is $T=30$ years.

The reference temperature for the finest grid is shown in figure 4. Upscaled temperature profiles for the two largest refinement factors are shown in figure 5 for

smoothed basis functions and for constant basis functions, respectively. Figure 6 shows one dimensional cross-sections of the temperature profiles at vertical value $y=375$. As seen from the figures, when using the smoothed basis functions, the upscaled method preserves the main features of the fine-scale solution, even for high upscaling ratios. The constant basis functions overestimate the diffusion, and the heat transfer from the matrix to the fractures is not well captured. The error increases with the coarsening factor, as can also be seen in the energy error given in table 1. This is consistent with the discussion in section 3.2.6, in that the gradient is approximated with respect to the fine-scale grid size $h$, rather than the coarse resolution $H$.

To investigate the impact of smoothing iterations we have run a number of simulations varying the maximal number of iterations, without terminating the smoothing of any basis function early. The upscaled simulations in this section are therefore run with a fixed number of iterations, given in table 1. When the coarsening factor $CF$ is larger more iterations are needed in order to reduce the error in the temperature equation. This is not surprising since for each iteration the support of a basis function for a coarse cell extends approximately one fine grid cell in each spatial direction. Thus, for the temperature gradient to be taken with respect to $H$ rather than $h$, the number of iterations should scale with $H/h$. However, at a certain number of iterations, depending on the refinement ratio, the error typically starts to increase with more interations. This is due to that the explicit representation of fractures leads both to a non-uniform coarse grid structure and to irregular shapes of the interaction regions for each coarse cell. What we have observed is that at a certain point more iterations introduce oscillations in the the basis functions. To ensure the construction of (overall) monotone basis functions, we enforce termination of the iterative procedure for each basis function through the criteria of energy reduction described in section 3.2.5 in the numerical experiments presented in the following sections.

4.2 Investigation of Heat Transport Characteristics

The previous example showed that the coarse discretization, with grids adapted to the fracture network, combined with smoothed basis functions gives a good approximation of the fine-scale temperature profiles. When the conduction term is discretized with constant basis functions, heat conduction is overestimated, and the approximation quality may suffer. How severe this effect is naturally depends on the relative importance of heat advection and conduction, as quantified by the heat Péclet number.

In this set of numerical experiments we consider a more realistic setup and investigate how variations of the heat Péclet number throughout the reservoir affect the numerical results. We consider a cell based heat Péclet number, defined for each cell by equation (7), where fine scale velocities averaged on each grid cell are used for $u$. The fracture network used in the simulations in this section are generated using the stochastic fracture generator Frac3D SilberhornHemminger2002 , except from a few of the largest fractures with assumed known positions. The fracture network is shown in figure 7, and apertures range from $a=10^{-3}$ m to $a=0.1$ m for the largest fractures. In practice, the largest of these values corresponds to fracture corridors and fault zones rather than individual features, but we will for simplicity refer to them as fractures. Note in particular that the large scale fractures do not form one single connected network. The average heat Péclet numbers per fine grid cell are also shown in figure 7. Note that the cell based heat Péclet numbers vary by orders of magnitude in different regions of the reservoir due to the fracture-matrix heterogeneity. From figure 7 we observe that, as expected, advection is the dominant process in the fractures. Since the fractures do not form a continuous pathway between inlet and outlet, the fluid needs to flow through the matrix. This is reflected in relatively high flow rates in the matrix in the middle of the domain, whereas the flow is smaller towards the domain boundary. Based on the previous example, we expect the largest differences between the smoothed and constant basis functions in regions where the heat Péclet numbers attain low values.

In the simulations 1 $dm^{2}/s$ of water with temperature 20°C is injected into the domain through a fracture near the bottom left corner and leaves through a fracture near the upper right corner. Exact locations for injection and production are shown in figure 7. We use two computational grids, one fine reference grid with $103893$ grid cells and an upscaled grid with $1852$ grid cells, i.e the coarsening factor is $CF=56$. The grid partitioning for the coarse grid is shown in figure 8, to illustrate the nonuniform grid structure. A time sequence of temperature profiles for simulations using the fine grid and the upscaled grid with both smoothed and constant basis functions are shown in figure 9 at times $T=0.5,1$ and $5$ years, respectively. For shorter times, advection is dominating and both upscaled simulations capture fine scale behavior well, with some additional numerical diffusion due to the coarse spatial discretization. For longer times the difference between the smoothed basis and the constant basis is more visible. The smoothed basis performs better compared to the constant basis with respect to capturing heat conduction in the matrix, leading to a less diffused temperature profile. As expected, the differences are largest towards the domain boundary, where the flow rates are small, and conduction is the dominant transport process. The production temperatures for both the fine scale simulation and the upscaled simulations are shown in figure 10. Both upscaled methods compare well to the fine scale reference simulation, and with smoothed basis functions the upscaled production temperature is closer to that of the reference solution compared to using the constant basis.

4.3 A Highly Complex Fracture Network

As a final demonstration of the performance of the coarse discretization framework we apply the method to a very challenging fracture network, shown in figure 11. The fractures are generated with Frac3D SilberhornHemminger2002 , as one set of fractures with high fracture density and the aperture is set to $a=0.01$ m. All fractures are represented both in the pressure grid and in the transport grid to demonstrate the applicability of the method. For efficiency, disconnected fractures can be excluded in the explicit representation and upscaled into the porous medium. We consider a constant background permeability of $10^{4}$ mD approximating the upscaled fractured media. The fracture permeability is modeled by equation (4), leading to a contrast between fracture permeability and background permeability of almost six orders of magnitude. In figure 11 the average heat Péclet numbers per fine grid cell are shown in addition to the fracture network, as well as production and injection locations. We observe that due to the high fracture density, there are only small matrix regions with low flow rates. The fine grid contains $355104$ grid cells and the upscaled grid $25616$ grid cells, i.e. the coarsening factor is $CF=14$. The smoothening of basis functions is terminated after $45$ iterations. The injection rate is 5 $dm^{2}/s$ and the injection temperature is 20°C.

Figure 12 shows the fine scale and upscaled temperature profiles at $T=0.2,0.6$ and $1$ years, respectively. Both methods preserve the fracture heterogeneity very well. Over longer times, it is clear from figure 12 that the smoothed basis adds less diffusion compared to the constant basis. Figure 13 shows the production temperature profiles, and we see that even for this highly complex fracture network, the upscaled methods perform very well, with a similar asymptotic behavior for both methods. This is also reflected by how the $\ell_{2}$ energy errors over space changes over time, as shown in figure 14. As time progresses, the difference between the $\ell_{2}$ energy errors for a smoothed basis and a constant basis increase in favor for the smoothed basis. Nevertheless, the numerical experiments clearly show that both methods can handle complex fracture networks.

5 Concluding Remarks

In this paper, we have presented an upscaling methodology for heat transfer in geothermal fractured reservoirs, modeled by an advection-conduction equation for the temperature. The methodology uses flow- and structure-based upgridding strategies well suited for advection-driven processes. The upgridding is combined with a discretization of the coarse scale conductive term based on an algebraic construction of basis functions.

We present two discretizations of the coarse conductive term: One based on piecewise constant functions, and a second formed by a the carefully constructed smoothing procedure of the constant functions. While the former has a simple construction, the latter leads to a superior represenatation of heat transfer via conduction from matrix to fractures. We demonstrate how to terminate the smoothing iterations to avoid oscillations in the numerical solution.

Numerical simulations show that our upscaled framework produces accurate results compared to fine scale numerical solutions for relatively large upscaling ratios, i.e. at a significantly lower cost compared to fine scale simulations. We also demontrate the applicability of the upscaling framework to a highly complex fracture network, and show that the structural heterogeneity of the fractured reservoir can be preserved.

Both the upgridding and the construction of basis functions are based on algebraic approaches, and in principle extensions to three dimensions should therefore be feasible. However, the upgridding procedure will likely produce no less complex shapes of the coarse cells, and the construction of smoothed basis functions may therefore require further investigation to be robust in 3D. On the other hand, constant basis functions should be readily applicable, and as the numerical examples show will also give satisfactory accuracy in many cases.

In total, we believe that the present framework has the potential to significantly facilitate field-scale simulation of heat recovery from fractured rocks.

Acknowledgements.

This work is supported by VISTA project no. 6357.