Heterogeneity Preserving Upscaling for Heat Transport in Fractured Geothermal Reservoirs

by   Anna Nissen, et al.
University of Bergen

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.



There are no comments yet.


page 6

page 9

page 10

page 11

page 13

page 14

page 16

page 17


Projection-based resolved interface mixed-dimension method for embedded tubular network systems

We present a flexible discretization technique for computational models ...

Streamline-Based Simulation of Carbon Dioxide Sequestration in Saline Aquifers

Subsurface sequestration of CO2 has received attention from the global s...

An Online Generalized Multiscale finite element method for heat and mass transfer problem with artificial ground freezing

In this paper, we present an Online Generalized Multiscale Finite Elemen...

An Asymptotically Compatible Formulation for Local-to-Nonlocal Coupling Problems without Overlapping Regions

In this paper we design and analyze an explicit partitioned procedure fo...

Energetic Stable Discretization for Non-Isothermal Electrokinetics Model

We propose an edge averaged finite element(EAFE) discretization to solve...

Combined Newton-Raphson and Streamlines-Upwind Petrov-Galerkin iterations for nano-particles transport in buoyancy driven flow

The present study deals with the finite element discretization of nanofl...

A Characteristic Mapping Method for Transport on the Sphere

A semi-Lagrangian method for the solution of the transport equation on a...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

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,


in combination with Darcy’s law,


is combined into an elliptic equation for the pressure, ,


where K

is the permeability tensor,

is the viscosity, which is assumed to be constant, v is the Darcy flux field and contains any sources or sinks. The relation between hydraulic aperture, , and fracture permeability, , is modelled as


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


where subscripts and denote fluid and rock, respectively, and denotes effective. Here, is the density, is the specific heat capacity, and 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, , given by


In addition, is the energy source term, and 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 , defined as


where is a representative length scale, typically the distance between injector and producer wells, is an average velocity, and is the average heat conductivity. The flow injection rate generally determines the flow rate through the fracture network and subsequently the highest value of 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 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 be a set of fine grid cells, = {, where is the total number of fine cells. Integrating equation (3) over each grid cell and using the divergence theorem gives


where is the boundary of grid cell and n is the outward normal of . Similarly, equation (6) becomes


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



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


where is a vector containing cell-centered temperature values, and approximate the advective and the conductive terms in eq. (9), respectively, and approximate the right hand side. , , and 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:

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

  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. 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. 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.

Figure 1: Construction of the coarse grid. Fine scale grid with fractures explicitly represented (top left). TOF indicator and outlined coarse grid based on TOF indicator (top right). Distance indicator and outlined coarse grid based on distance indicator (bottom left). Final coarse grid, where the TOF coarse grid and the distance coarse grid have been combined (bottom right).

3.2 Coarse Scale Discretization

Let be a set of coarse grid cells, = {, , where 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 .

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,


where denotes a face between fine grid cells and , and a face between coarse grid cells and .

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, , is related to the corresponding fine scale discretization matrix, , by


where R is a restriction operator represented by a matrix of size , and P is a prolongation matrix of size . The restriction operator, R, maps the fine scale discretization to the coarse scale. We use a finite volume restriction matrix, with elements given by


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


where the columns , , 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 contain some fine cells that are included in . 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


i.e. . We note that is a partition of unity, and forms a (piecewise constant) discretization on the coarse grid.

Now, successive iterations for basis function are given by


where , is a relaxation parameter, and 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 as the discrete energy of the basis function for coarse grid cell , by


Oscillations in the basis function will increase the discrete energy, and we therefore terminate the iterations for basis function if for some iteration we have


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 . While the termination criterion (19) does not guarantee monotone basis functions, the approach has been sufficient for the fracture networks considered here. In addition to the termination condition (19) 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 iterations, 24 () of the basis functions are still updating, and they are therefore not included in the figure.


In addition to non-monotonic basis functions, excessive smoothing may also result in eigenvalues of

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 such that , . If the diagonal elements of 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, , is symmetric, since in this case , and thus . 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 . That is, the coarse-scale conduction retains a dependency on the fine scale grid size , rather than on a representative coarse grid size .

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.

Figure 2: Initial (constant) basis functions for three coarse cells from the coarse grid in figure 1 are shown in the top figures. The corresponding smoothed basis functions after termination of iterations are shown in the bottom figures after 9, 50, and 9 iterations, respectively. The boundaries of the interaction regions for the basis functions are depicted with black lines in the bottom figures.
Figure 3: Histogram over the iteration count for basis functions in the simulations in section 4.2. There are 1412 of a total of 1852 basis functions for which iterations are terminated early. 440 basis functions are still updating when the smoothing is terminated based on the residual, and therefore not included in the histogram.

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, kJ/(mK) and J(mKs), respectively. For the volumetric heat capacity of the fluid and viscosity, kJ(mK), and cP are used, respectively. The porosity in the matrix is set to . The reservoir domain is m 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 is given by


where refers to the grid cell index. The relative error over space in the energy at time , , is then computed as


where is the energy vector for the upscaled solution projected on the fine grid, given by , and is the energy vector for a reference solution computed on the fine grid. 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 . 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 grid cells and the fine grids contain , , , and grid cells, respectively, including one-dimensional fracture cells. Water with temperature 20°C is injected at a rate of 0.1 in the middle of the domain (at ), and extracted where each fracture meets the exterior boundaries. The total simulation time is 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 . 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 , rather than the coarse resolution .

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 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 rather than , the number of iterations should scale with . 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.

Figure 4: Reference temperature profile from fine scale simulation in the refinement study presented in section 4.1. The fine grid contains 413449 grid cells, including fracture cells.
Figure 5: Upscaled temperature profiles using piecewise constant basis functions (left) and smoothed basis functions (right), respectively, for the simulations in section 4.1. The fine scale underlying grids have a total of and grid cells, respectively, depicted from top to bottom. The coarse grid contains grid cells.
No. fine cells Coarsening factor (no. of iterations)
6889 11
26569 44
104329 174
413449 688
Table 1: Accuracy with respect to coarsening factor for the refinement study in section 4.1. and denote the energy error over space for constant and smoothed basis functions, respectively.
Figure 6: One dimensional temperature profiles at for a successive grid refinement, with constant basis functions (left) and smoothed basis functions (right), respectively, for the refinement study in section 4.1. The coarse grid contains grid cells for all simulations, the coarsening factors are for each case given in the figures.

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 . 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 m to 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 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 grid cells and an upscaled grid with grid cells, i.e the coarsening factor is . 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 and 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.

[width=7cm,height=7cm]Fractures_runEX3_1point5_linewidth.png InjectionProduction

Figure 7: Fractures in the discretization for the simulations in section 4.2 (left) and the field with heat Péclet number per grid cell (right). The heat Péclet number plot is in logarithmic scale with base . Injection and production locations are indicated by arrows in the fracture network plot.
Figure 8: The coarse grid partition used for the numerical experiments in section 4.2. The coarse grid has 1852 grid cells and is generated from a fine grid with 103893 grid cells.
Figure 9: Time sequence of temperature profiles for reference solution (top) and upscaled solution with 1852 grid cells with smoothed basis (middle) and constant basis (bottom). The temperature profiles are plotted after and years, respectively, from left to right. Up to 111 iterations are used to smoothen the basis functions. The total relative energy errors at years are for the upscaled method with smoothed basis and for the upscaled method with constant basis, both compared to the fine scale solution. The upscaling factor is 56.
Figure 10: Production temperature as a function of time for the simulations in section 4.2. Temperature profiles are shown in figure 9.

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 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 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 grid cells and the upscaled grid grid cells, i.e. the coarsening factor is . The smoothening of basis functions is terminated after iterations. The injection rate is 5 and the injection temperature is 20°C.

Figure 12 shows the fine scale and upscaled temperature profiles at and 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 energy errors over space changes over time, as shown in figure 14. As time progresses, the difference between the 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.

[width=7cm,height=7cm]Fractures_runEX4_1point5_linewidth.png InjectionProduction

Figure 11: Fractures in the discretization for the simulations in section 4.2 (left) and the field with heat Péclet number per grid cell (right). The heat Péclet number plot is in logarithmic scale with base . Injection and production locations are indicated by arrows in the fracture network plot.
Figure 12: Fine scale and upscaled temperature profiles after , and years, respectively, shown from left to right. The fine scale simulation uses 355104 grid cells and the upscaled simulations uses a grid with with 25616 grid cells, i.e. the upscaling factor is . The fine scale temperatures are shown in the top figures, temperatures with smoothed basis functions in the middle, and temperatures using constant basis functions in the bottom figures, respectively. The fracture network used in the simulations is shown in figure 11.
Figure 13: Production temperature as a function of time for the simulations in section 4.3. Temperature profiles are shown in figure 12.
Figure 14: Reletive energy error over space as a function of time for upscaled simulations in section 4.3, for a constant basis and a smoothed basis, respectively.

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.

This work is supported by VISTA project no. 6357.


  • (1) J.E. Aarnes, and Y. Efendiev, An adaptive multiscale method for simulation of fluid flow in heterogeneous porous media, Multiscale Model. Simul. 5(3), pp. 918–939 (2006).
  • (2) J.E. Aarnes, V. L. Hauge, and Y. Efendiev, Coarsening of three-dimensional structured and unstructured grids for subsurface flow, Adv. Water Resour. 30(11), pp. 2177–2193 (2007).
  • (3) I. Aavatsmark, An introduction to multipoint flux approximations for quadrilateral grids, Comput. Geosci. 6(3), pp 405–432 (2002).
  • (4) P. M. Adler and J.-F. Thovert, Fractures and Fracture Networks, 430 p., Kluwer Academic Publishers, Dordrecht, the Netherlands (1999).
  • (5) P. M. Adler, J.-F. Thovert, J.-F. and V. V. Mourzenko, Fractured Porous Media, 184 p., Oxford University Press, United Kingdom (2012).
  • (6) K. Aziz, A. Settari, Petroleum Reservoir Simulation, 476 p., Science Publishers Ltd., London, United Kingdom (1979).
  • (7) G. I. Barenblatt, I. P. Zheltov, and I. N. Kochina, Basic concepts in the theory of seepage of homogeneous liquids in fissured rocks, J. Appl. Math. Mech. 24(5), pp. 1286–1303 (1960).
  • (8) M. Bai, D. Elsworth and J.-C. Roegiers, Multiporosity/multipermeability approach to the simulation of naturally fractured reservoirs, Water Resour. Res., 29(6), pp. 1621–1633 (1993).
  • (9) C.C. Barton, Fractal analysis of scaling and spatial clustering of fractures, Fractals in the Earth Sciences, Springer US, pp. 141–178 (1995).
  • (10) B. Berkowitz, Characterizing flow and transport in fractured geological reservoirs: A review, Adv. Water Resour., 25, pp. 861–884 (2002).
  • (11) E. Bonnet, O. Bour, N. E. Odling, P. Davy, I. Main, P. Cowie and B. Berkowitz, Scaling of fracture systems in geological media, Rev. Geophys. 39(3), pp. 347–383 (2001).
  • (12) D. Bruel and M. C. Cacas, Numerical modeling technique: contribution to the Soultz HDR project, In Geothermal Energy in Europe -The Soultz Hot Dry Rock Project, Gordon and Breach Science Publishers, pp. 267–279 (1992).
  • (13) M. C. Cacas, E. Ledoux, G. de Marsily, A. Barbreau, P. Calmels, B. Gaillard and R. Margritta, Modeling fracture flow with a stochastic discrete fracture network: Calibration and validation 2. The transport model Water Resour. Res., 26(3), pp. 491–500 (1990).
  • (14) P. Cheng, Heat transfer in geothermal systems, Adv. Heat Transfer 14, pp. 1–105 (1979).
  • (15) C. D’Angelo, and A. Scotti, A mixed finite element method for Darcy flow in fractured porous media with non-matching grids, ESAIM-Math. Model. Num., 46(2), pp. 465–489 (2012).
  • (16) P. Dietrich, R. Helmig, M. Sauter, H. Hötzl, J. Köngeter and G. Teutsch, Flow and Transport in Fractured Porous Media, 447 p., Springer-Verlag Berlin Heidelberg, (2005).
  • (17) A. Fox, P. La Pointe, J. Hermanson, and J. Öhman,Statistical geological discrete fracture network model: Forsmark modelling stage 2.2, Stockholm: SKB, 2007.
  • (18) H. Hajibeygi, G. Bonfigli, M. A. Hesse, and P. Jenny, Iterative multiscale finite-volume method, J. comput. Phys. 227, pp. 8604–8621 (2008).
  • (19) V. L. Hauge, K.-A. Lie and J. R. Natvig, Flow-based coarsening for multiscale simulation of transport in porous media, Comput. Geosci., 16(2), pp. 391–408 (2012).
  • (20) V. Hauge Multiscale methods and flow-based gridding for flow and transport in porous media, PhD Thesis, NTNU (2010).
  • (21) K. Hayashi, J. Willis-Richards, R. J. Hopkirk and Y. Niibori, Numerical models of HDR geothermal reservoirs – a review of current thinking and progress, Geothermics, 28(4-5), pp. 507–518 (1999).
  • (22) H. Hajibeygi, D. Karvounis and P. Jenny, A hierarchical fracture model for the iterative multiscale finite volume method, J. Comput. Phys., 230(24), pp. 8729–8743 (2011).
  • (23) R. Holm, R. Kaufmann, B.-O. Heimsund, E. Øian and M. S. Espedal, Meshing of domains with complex internal geometries, Numer. Linear Algebra Appl., 13, pp. 717–731 (2006).
  • (24) P. Jenny, S. H. Lee and H. A. Tchelepi, Multi-scale finite-volume method for elliptic problems in subsurface flow simulation, J. Comput. Phys., 187(1), pp. 47–67 (2003).
  • (25) M. Karimi-Fard, L. Durlofsky, and K. Aziz, An efficient discrete-fracture model applicable for general-purpose reservoir simulators, SPE Journal, pp. 227–236 (2004).
  • (26) M. Karimi-Fard, B. Gong, and L. J. Durlofsky, Generation of coarse-scale continuum flow models from detailed fracture characterizations, Water Resour. Res. 42(10), W10423 (2006).
  • (27) M. Karimi-Fard and L. J. DurlofskyA general gridding, discretization, and coarsening methodology for modeling flow in porous formations with discrete geological features, Adv. Water Resour., 96, pp. 354–372 (2016).
  • (28) E. Keilegavlen, J. M. Nordbotten, Inexact linear solvers for control volume discretizations in porous media, Comp. Geosci. 19(1), pp. 159–176 (2015).
  • (29) K.-A. Lie, S. Krogstad, I. S. Ligaarden, J. R. Natvig, H. M, Nilsen, B. Skaflestad, Open-source MATLAB implementation of consistent discretisations on complex grids, Comput. Geosci., 16(2), pp. 297–322 (2012).
  • (30) V. Martin, J. Jaffré, and J. E. Roberts, Modeling fractures and barriers as interfaces for flow in porous media, SIAM J. Sci. Comput., 26(5), pp. 1667–1691 (2005).
  • (31) O. Møyner and K.-A. Lie, A multiscale restriction-smoothed basis method for high contrast porous media represented on unstructured grids, J. Comput. Phys., 304, pp. 46–71 (2016).
  • (32) R. Munier, Statistical analysis of fracture data, adapted for modelling Discrete Fracture Networks-Version 2 Stockholm: SKB, 2004.
  • (33) H. Mustapha and K. Mustapha, A new approach to simulating flow in discrete fracture networks with an optimized mesh SIAM J. on Sci. Comput., 29(4), pp. 1439–1459 (2007).
  • (34) J. M. Nordbotten and P. E. Bjørstad, On the relationship between the multiscale finite volume method and domain decomposition preconditioners, Comput. Geosci., 12(3), pp. 367–376 (2008).
  • (35) K. Pruess and T. Narasimhan, On fluid reserves and the production of superheated steam from fractured, vapor-dominated geothermal reservoirs, J. Geophys. Res., 87(B11), pp. 9329–9339 (1982).
  • (36) K. Pruess, Heat transfer in fractured geothermal reservoirs with boiling, Water Resour Res., 19(1), pp. 201–208 (1983).
  • (37) K. Pruess and T. N. Narasimhan, A practical method for modeling fluid and heat flow in fractured porous media, Soc. Petrol. Eng. J., 25(1), pp. 14–26 (1985).
  • (38) K. Pruess and Y. Wu, A new semi-analytical method for numerical simulation of fluid and heat flow in fractured reservoirs, SPE Adv. Tech. Ser., 1(2), pp. 63–72 (1993).
  • (39) V. Reichenberger, H. Jakobs, P. Bastian, and R. HelmigA mixed-dimensional finite volume method for two-phase flow in fractured porous media, Adv. Water Resour., 29(7), pp. 1020-1036 (2006).
  • (40) T. H. Sandve, I. Berre, and J. M. Nordbotten, An efficient multi-point flux approximation method for Discrete Fracture-Matrix simulations, J. Comput. Phys., 231(9), pp. 3784–3800 (2012).
  • (41) T. H. Sandve, E. Keilegavlen, and J. M. Nordbotten, Physics-based preconditioners for flow in fractured porous media, Water Resour. Res., 50(2), pp. 1357–1373 (2014).
  • (42) S. Shah, O. Møyner, M. Tene, K.-A. Lie and H. Hajibeygi, The multiscale restriction smoothed basis method for fractured porous media (F-MsRSB), J. Comput. Phys., 318, pp. 36–57 (2016).
  • (43) A. Silberhorn-Hemminger, Modellierung von Kluftaquifersystemen: Geostatistische Analyse und deterministisch-stochastische Kluftgenerierung, Dissertation, Institut für Wasserbau, Universität Stuttgart, 2002. ISBN: 3-933761-17-4.
  • (44) M. Tene, M. S. Al Kobaisi and H. Hajibeygi, Algebraic multiscale solver for flow in heterogeneous fractured porous media, In: Proceedings of SPE Reservoir Simulation Symposium (2015).
  • (45) A. Toselli and O. B. Widlund, Domain decomposition methods: algorithms and theory, 450 p., Springer Berlin Heidelberg, (2005).
  • (46) P. M. Adler and J.-F. Thovert, Fractures and Fracture Networks, 430 p., Kluwer Academic Publishers, Dordrecht, the Netherlands (1999).
  • (47) J. van Lent, R. Scheichl, I. G. Graham, Energy-minimizing coarse spaces for two-level Schwarz methods for multiscale PDE, Numer. Linear Algebra Appl., 16(10), pp. 775-799 (2009).
  • (48) P. Vanek, J. Mandel, M. Brezina, Algebraic multigrid by smoothed aggregation for second and fourth order elliptic problems, Computing, 56(3), pp. 179-196 (1996).
  • (49) J. E. Warren and P. J. Root, The behavior of naturally fractured reservoirs, Soc. Petrol. Eng. J., 3(3), pp. 245–255 (1963).
  • (50) K. Watanabe and H. Takahashi, Fractal geometry characterization of geothermal fracture networks, J. Geophys. Res., 100(B1), pp. 521–528 (1995).
  • (51) H. Zhou and H. A. Tchelepi, Operator-based multiscale method for compressible flow, Soc. Petrol. Eng. J., 13(2), pp. 267–273 (2008).