Analysis of conforming, non-matching, and polygonal methods for Darcy and advection-diffusion-reaction simulations in discrete fracture networks

by   Andrea Borio, et al.

The aim of this study is to compare numerical methods for the simulation of single-phase flow and transport in fractured media, described here by means of the Discrete Fracture Network (DFN) model. A Darcy problem is solved to compute the advective field, then used in a subsequent time dependent transport-diffusion-reaction problem. The numerical schemes are benchmarked in terms of flexibility in handling geometrical complexity, mass conservation and stability issues for advection dominated flow regimes. To this end, three benchmark cases have been specifically designed and are here proposed, representing some of the most critical issues encountered in DFN simulations.



There are no comments yet.


page 19

page 20

page 22

page 25

page 28

page 30

page 31


Analysis of Flux Corrected Transport Schemes for Evolutionary Convection-Diffusion-Reaction Equations

We report in this paper the analysis for the linear and nonlinear versio...

Aggregation in non-uniform systems with advection and localized source

We explore analytically and numerically agglomeration driven by advectio...

A cure for instabilities due to advection-dominance in POD solution to advection-diffusion-reaction equations

In this paper, we propose to improve the stabilized POD-ROM introduced b...

Verification benchmarks for single-phase flow in three-dimensional fractured porous media

Flow in fractured porous media occurs in the earth's subsurface, in biol...

Diffusion-redistanciation schemes for 2D and 3D constrained Willmore flow: application to the equilibrium shapes of vesicles

In this paper we present a novel algorithm for simulating geometrical fl...

An Assessment of Solvers for Algebraically Stabilized Discretizations of Convection-Diffusion-Reaction Equations

We consider flux-corrected finite element discretizations of 3D convecti...

A convergent interacting particle method and computation of KPP front speeds in chaotic flows

In this paper, we study the propagation speeds of reaction-diffusion-adv...

Code Repositories

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

The movement of liquids in the underground is heavily influenced by the presence of fractures and their relative intersections Dowd2009 ; Hardebol2015 ; Flemisch2016a . Fractures are discontinuities (here assumed planar) along which a rock has been broken, mainly due to geological movements or to artificial stimulation McClure2015 . In this work, we are considering only open structures, characterized by a geometrical aperture, that allow a liquid to flow through Erhel2009 . Possibilities are actual fractures, faults, and joints. We are thus excluding low permeable (closed/impervious) objects such as veins or dykes. For particular underground compositions (e.g., granite, shale, or sandstone) the rock permeability is several orders of magnitude smaller than the fracture permeability. It is a common choice and a reasonable approximation to ignore the rock matrix effect in the simulations and rely only on the fractures. The framework is the discrete fracture network model (DFN), where the aperture is not a geometrical constraint but a parameter in the bidimensional representation of the fractures by reduced models, see Martin2005 and the forthcoming references.

The geometrical complexity of natural fracture networks may impose difficulties in the numerical simulations, due to the presence of small intersections between fractures, different intersection configurations (e.g., Y-type or L-type), small angles of intersection, and small distances between intersection lines. In the literature, three main approaches are developed to overcome these difficulties.

The first approach considers rather standard numerical scheme for the discretization of the physical equations, and relies on robust grid generation software. The coupling conditions among fractures are imposed via Lagrange multipliers on a representation of the interfaces conforming with the computational mesh on all sides. To be specific, the edges of the computational elements have to match at the intersections. This approach, commonly known as conforming discretization, may suffer, for example, when intersections lines are very close to each other, see Hyman2014a ; Hyman2015a ; Mustapha2007a ; NGO2017 . To partially overcome these problems, a possibility is to consider the class of virtual element methods that allows grid elements of general shapes. Two different numerical methods of this type are introduced in Benedetto2014 ; Fumagalli2016a ; Benedetto2017 ; Fumagalli2017a ; BBorth .

A second possibility is to keep the intersections explicitly represented, but relaxing the conformity of the edges. This approach, named non-conforming, requires more advanced numerical schemes based on the mortar technique. In this case, we relax the actual generation of the fracture meshes which usually gives less discrete elements than the conforming approach Benedetto2016 ; Nordbotten2018 . This method may still suffer in presence of severe geometries. Also in this case the virtual element methods are an interesting option to further decrease the computational cost, see BeiraodaVeiga2014b ; BeiraoVeiga2016 .

A third family of schemes comprises the so-called non-matching discretizations. In terms of mesh generation, in this case the intersections do not place any constraints, as the fractures are meshed independently and the coupling conditions are imposed by an optimization procedure. A functional that measures the mismatch in the coupling conditions is minimized iteratively, where only the degrees of freedom involved in the intersections (the cut region) are considered

BPSc ; BPSd ; BSV ; BBBPSsupg ; BPSf . This procedure is independent on the actual numerical schemes, might take advantage of ad-hoc strategies to enrich the solution in the cut region. Extended finite elements are a successful example, see, for example, DAngelo2011 ; Formaggia2012 ; Fumagalli2012d ; Fumagalli2012bb .

The present work extends and enriches the one proposed in Scialo2017 to more complex physical phenomena. The concepts previously discussed are applied to Darcy and heat transport-and-diffusion models. The Darcy velocity is computed first and then used as an advective field for the heat equation. High quality computation or reconstruction of the Darcy velocity may significantly impact the temperature distribution. Moreover, in the heat equation the transport part usually dominates the diffusion (Péclet greater than 1) and stable or stabilized schemes are needed to avoid or limit spurious oscillations that might compromise the accuracy of the solution. In the numerical tests, we are considering several numerical schemes for the comparison to cover most of the combinations discussed and try to assess their performances. The aim of this work is thus twofold: establish a set of benchmark cases and give guidance in the development of more advanced numerical schemes to solve this problem.

The paper is organized as follow. In Section 2 the Darcy and heat models are introduced and discussed, with particular focus on the coupling conditions. Section 3 is devoted to the description of the proposed numerical schemes. Three numerical examples are presented and discussed in Section 4, comparing the performances of the considered numerical schemes. Finally, in Section 5 we draw some conclusions and suggestions for future developments.

2 Mathematical model

In this section we introduce the mathematical model used to describe the hydraulic head and Darcy velocity profiles in a discrete fracture network. Once this problem is solved, the Darcy velocity is considered as advective field to simulate the transport and diffusion of heat in DFNs.

Fractures are considered as non-overlapping planar objects, which can be connected to other fractures through intersection segments, also called traces. We consider fractures with boundary , which compose the discrete fracture network , and we denote its boundary as with outward unit normal

, defined on each fracture plane as the unit vector normal to the fracture boundary pointing outward from the fracture polygon. To keep the presentation clear, we make a distinction between external boundaries, where we will impose data, and internal interfaces, i.e. the traces, where we will couple the mathematical models on the fractures.

Given two distinct and intersecting fractures and , with , we indicate their intersection (trace) as . For simplicity, we assume that a trace is formed only by two distinct fractures, however this assumption can be relaxed. A natural order of indexes can be introduced to numerate the traces from 1 to , being the latter their cardinality. We consider also the function such that with . We have and its inverse is well defined such that where . We indicate with the union of all the traces and by the set of traces belonging to the fracture . Moreover, consider a fracture and a trace , with . naturally subdivides in two sub-regions, indicated by and , such that . To each of these sub-regions we associate an outward unit normal perpendicular to , denoted by and (with ) and a trace operator, respectively denoted by and . Generic trace operators on fracture are denoted by . An example of a simple DFN with the introduced nomenclature is given in Figure 1.

Figure 1: Representation of two fractures and intersecting in .

We present the Darcy problem in strong and weak form in Subsection 2.1, to compute the hydraulic head and Darcy velocity, whereas the heat equation is introduced along with its weak formulation and functional setting in Subsection 2.2.

2.1 The Darcy model

This section is devoted to the presentation of the mathematical models used to describe the hydraulic head field and Darcy velocity in a discrete fracture network, further details being available in BPSc ; BPSd ; Fumagalli2016a ; Scialo2017 and references therein. Unknowns and parameters restricted to a fracture are denoted by a subscript .

For clarity in the exposition, we start considering a single fracture . The Darcy model on reads: find such that

Variables, data and differential operators are defined on the tangent plane of the fracture.

is the hydraulic conductivity tensor, which is symmetric and positive definite, and

is a scalar source/sink term. In our applications, following lubrication theory Bear1972 ; Wangen2010 ; Schon2011 , we consider a scalar permeability, obtained by the cubic law: , where is the fracture aperture. Moreover, the hydraulic conductivity in (1a) is isotropic and defined as
with the fluid density, the gravity acceleration, the dynamic viscosity of the fluid.

The boundary conditions on are


where and are disjoint portions of the boundary of such that , and gives the trace on the boundary of . The data is the given pressure head at the boundary and we model as impervious. Other boundary conditions are possible, however, in order to keep the presentation as simple as possible and coherent with the examples proposed in Section 4, we consider only these, being the generalization to other boundary conditions straightforward. We assume, for a single fracture, that is not empty; in presence of intersecting fractures, instead, we can have for all the fractures but one. Equations (1a) and (1b) are well studied in the literature, refer to the aforementioned references.

Let us now consider two distinct fractures and forming an intersection . On both fractures we apply Equations (1a) and (1b) and we assume continuous coupling conditions for the hydraulic head and the normal component of the Darcy velocity at the trace. The coupling conditions between two fractures and such that , read


The case of multiple fractures follows immediately from our assumptions. The Darcy equation can be formalised in the following problem.

Problem 1 (Darcy equation on a DFN - mixed formulation).

Given the set of fractures and traces , find such that (1a)-(1b)-(1c) are satisfied.

The weak formulation of Problem 1 requires the introduction of the functional spaces on each , namely

equipped with the natural norms that make them complete. The spaces for the global network are

and , both endowed with their natural norms. The following linear forms are then introduced: , , and such that, and

We also require that , , and . The less standard functional spaces associated with are due to the low regularity of a -function at a boundary, as well as the duality pairing used instead of a scalar product. The global forms , , , and are the sums of the respective local ones, where we denoted by , the space associated to the boundary term.

Problem 2 (Darcy equation on a DFN - mixed weak formulation).

The weak formulation of Problem 1 is: find such that

Some numerical formulations require the problem written in term of pressure head alone, thus we consider also the primal formulation of Problem 1. Given a fracture the model is: find such that

Given a second distinct fracture which intersects in , the following coupling conditions are added:

Once the pressure head is computed, the Darcy velocity is given by the first equation of (1a). We can formalize also the primal formulation of the Darcy equation in the following problem.

Problem 3 (Darcy equation on a DFN - primal formulation).

Given the set of fractures and traces , find such that (2a)-(2b) are satisfied .

In this case, the weak formulation is simpler than the previous one. To keep the notation standard and since we are introducing a separate formulation of the Darcy problem, in the following we commit abuses in notation. We introduce the functional spaces

and in the special case of homogeneous head condition. The linear forms now read: and such that, ,

In this case, we ask , , and . The functional spaces on the whole network are

with their natural norms. The global forms and are the sums of the respective local ones.

Problem 4 (Darcy equation on a DFN - primal weak formulation).

The weak formulation of Problem 3 is: find such that

We note that to obtain symmetric forms in Problem 4 a lifting technique of the Dirichlet datum can be considered.

2.2 The heat equation

In this part, we present the heat equation on the DFN. Once Problem 1 is solved, the Darcy velocity can be used as advective field in the transport problem. We denote the temperature in a DFN as , and its restriction to fracture as . The heat equation on reads: given find such that

where is the end time of the simulation. Also in this case, the variables, data, and differential operators are defined on the tangent plane of the fracture. The relations to compute the physical parameters are Bear1972 ; Wangen2010 ; Schon2011 :
where is the fracture aperture, the fluid density, the fluid specific thermal capacity, the fracture effective thermal capacity, the fracture porosity, the density of the rock matrix, the specific thermal capacity of the rock matrix, the effective thermal conductivity, the effective heat transfer coefficient between fluid and rock, and, finally, is the temperature of the rock matrix, acting as external heat source/sink. When the fracture is completely open and no sediments are present, whereas values of indicate that fracture is partially filled with crystals/rocks. Equation (3a) is in conservative form, being the total flux the conserved quantity.

For the sake of brevity and clarity of explanation we do not present boundary conditions in the most general context, but conforming to the numerical tests of Section 4; generalization is however again straightforward. Recalling that is the Dirichlet portion of the boundary of fracture for the Darcy problem in Section 2.1, let us split into two parts, namely and . Thus, according to the computed Darcy velocity , on each fracture , the inflow boundary is the portion of where and conversely the outflow boundary is the portion of where , thus linking the nature of the boundary to the solution of the Darcy problem. Please note that and might be both empty for most of the fractures. Then boundary conditions on fracture are:


By (1b), on the portion both the diffusive and advective terms are zero for , so represents an impervious boundary also for the heat equation. Finally, the initial condition on is given by


We consider now two distinct fractures and with an intersection . On the temperature fulfils (3a), (3b), and (3c). The same for in the fracture . The coupling conditions on are


Also in this case, in presence of multiple fractures the generalization of the model is immediate. The heat equation can be formalized in the following problem.

Problem 5 (Heat equation on a DFN).

Given the set of fractures and traces and the Darcy velocity , find such that (3a)-(3b)-(3c)-(3d) are satisfied.

To derive the weak formulation of Problem 5 we need to introduce the Bochner space

while for the test functions we define

We need to introduce space-time forms associated to each fracture, in particular we have: , , , which are defined on , and such that

The global functional spaces are

and the global bilinear forms , , , and are the sum over all fractures of the local ones. We finally have the weak formulation of the problem.

Problem 6 (Heat equation on a DFN - weak formulation).

The weak formulation of Problem 5 is: find such that

Also in this case to obtain a symmetric problem we can use a lifting technique of the Dirichlet datum.

3 Numerical discretization

In this section, we present various discretization strategies, both well established, both unconventional, that can be used to solve the models described in Section 2. These strategies have similarities and differences that can be used to categorize them. A first point concerns the computational mesh and, in particular, how the meshing is performed at fracture intersections: it is possible to have a matching or non-matching grids. In the former case fracture grids are conforming to the intersections among fractures, while in the latter, grid elements arbitrarily cross intersections. A second issue is related to mass conservation: in computing the Darcy velocity, some schemes are locally mass conservative and some other are only globally conservative, and this property may impact the subsequent solution of the heat problem. Also, some numerical schemes are characterized by high numerical diffusivity which might impact the solution but also avoids nonphysical spurious oscillations in advection dominated flow regimes, whereas other schemes need to adopt stabilitazion techniques.

Six different approaches are considered in the present work, given as the combination of a numerical scheme for the computation of the Darcy velocity and a numerical scheme for the spatial semi-discretization of the subsequent non-stationary advection-diffusion-reaction problem (shortly denoted as Heat equation). The implicit Euler method is used, instead, in all cases, for time evolution. The approaches are listed in Table 1. The scheme tagged tpfaup is given by the combination of the Two Point Flux Approximation (tpfa) method for the Darcy problem and the TPFA with upwinding for the advection term (tpfa+upwind) for the Heat equation. Scheme mfemup uses, instead mixed finite elements (mfem) for the Darcy equation and again tpfa+upwind for the Heat equation. The method mvemup, uses instead the Virtual Element Method in mixed formulation (mvem) for the Darcy problem, on matching polygonal meshes. These schemes are implemented in PorePy, a simulation tool written in Python for fractured and deformable porous media, see Keilegavlen2017a ; Keilegavlen2017b . PorePy is freely available on GitHub along with the numerical tests proposed in Section 4. The method labeled mfemsupg is based on mixed standard finite elements for the numerical resolution of the Darcy problem and on standard finite elements (fem) with Streamline Upwind Petrov-Galerkin (supg) stabilization Franca-Frey-Hughes:1992 for the Heat equation. Finally methods femsupg and xfemsupg use a non-conventional optimization based approach for the Darcy equation and, with supg stabilization also for the Heat equation. The optimization approach can adopt different baseline discretization methods: here we consider the variants using standard finite elements (opt-fem) and extended finite elements (opt-xfem). Methods mfemsupg femsupg and xfemsupg are implemented in C++ and Matlab.

The forthcoming parts describe in more details the previous approaches, grouping them according to the coupling at the traces: matching coupling at traces in Subsection 3.1 and non-matching coupling in Subsection 3.2.

Tag Darcy eq. Heat eq. Grid
tpfaup tpfa tpfa + upwind Matching triangles
mfemup mfem tpfa + upwind Matching triangles
mvemup mvem tpfa + upwind Matching polygons
mfemsupg mfem fem + supg Matching triangles
femsupg opt-fem opt-fem + supg Non-match. triangles
xfemsupg opt-xfem opt-xfem + supg Non-match. triangles
Table 1: List of numerical schemes proposed to solve the Darcy and heat problems, with the type of meshes used.

3.1 Matching discretization at traces

Here advantages and drawbacks of a conforming discretization at the traces are discussed. As mentioned before in a conforming grid, the meshes of both the intersecting fractures match the trace with their geometry. The trace is thus entirely covered by contiguous cell edges of the two fractures, see Figure 1(a)-1(b) as an example. This approach has the clear advantage of an easy applicability to most of the existing and well established numerical methods (finite volumes, finite elements). However, in the case of complex geometries the computational cost might increase and become a severe constraint in complex fracture networks.

(a) matching triangles
(b) matching polygons obtained by coarsening
(c) non-matching triangles
Figure 2: Examples of matching meshes.

To solve the Darcy problem, we rely on different classes of numerical schemes: finite volumes, finite elements in primal and mixed formulation, and virtual elements in mixed formulation. For the finite volume class, we choose the two-point flux approximation tpfa on simplicial grids, applied to Problem 4, see Aavatsmark2007a for details. This scheme is well known in the industry field and widely used for its velocity in assembling the discrete problem and for having a narrow matrix stencil. The scheme is locally conservative and robust with respect to strong variations of the hydraulic conductivity coefficient, however it is consistent only for -orthogonal grids. Regarding the class of finite elements, we consider mixed finite elements with the pair of spaces for the Darcy velocity and piecewise constants for the hydraulic head, denoted by mfem, defined on a simplicial grid, for further references see Raviart1977 ; Roberts1991 . It is well known that is locally conservative and more robust than primal finite elements with respect to strong variations on the hydraulic conductivity coefficient. On the other hand, it gives a larger linear system, with a saddle-point structure. The mfem scheme solves the Darcy problem in the form presented in Problem 2. In some particular scenarios the regularity requirements on meshes formed by triangles are too restrictive and schemes able to handle generally shaped cells are more suitable. In these cases we rely on the new class of virtual element methods, which are variational methods where the basis functions of the discrete spaces are not prescribed a-priori, and are defined implicitly on general star-shaped elements as solutions of suitable local PDEs. See the seminal works BeiraodaVeiga2013a ; Brezzi2014 ; BeiraodaVeiga2014a ; BeiraoVeiga2016 and those related to DFN Benedetto2014 ; Benedetto2016 ; Benedetto2016b ; Benedetto2017 ; Fumagalli2016a . In our analysis we consider only the lowest order mixed (mvem) formulation, which can be viewed as a generalization mixed finite elements on generally shaped cells, solving Problem 2. Virtual element methods share many properties with their finite element counterpart: indeed, mvem are locally conservative, robust with respect to the hydraulic conductivity variation and have the same grid stencil of the . Here we use mvem on polygonal grids, obtained by coarsening a mesh originally made of triangular elements in order to reduce the number of cells required in the simulation for complex geometries, as described in Fumagalli2016a ; Fumagalli2017 .

To solve the heat equation (Problem 6), we use a tpfa scheme for the diffusive term and a weighted upwind scheme for the advective term. The advantage of this choice is that we obtain a stable scheme which respects the maximum and minimum principle, without oscillation due to high grid-Péclet numbers. However, for some applications the obtained scheme might be too diffusive. Finally, we consider also standard finite elements, with a supg stabilized discrete variational formulation, which is a globally conservative numerical scheme.

Since in all the above cases we have matching meshes, and thus the degrees of freedom on the trace of one fracture correspond to the degrees of freedom of the intersecting one, coupling conditions can be imposed strongly.

Methods tpfaup, mfemup and mfemsupg are considered here standard reference approaches of different discretization strategies (finite volume and finite element based), and are used to benchmark the behavior of the other less conventional approaches, based e.g. on polygonal or non-matching discretizations.

3.2 Non-matching discretization at traces

When dealing with huge networks, the generation of conforming meshes may require a large computational cost. Then, it is worth considering a class of methods that do not require any kind of conformity of the fracture meshes to traces, see, for example, Figure 1(c).

In Berrone2013 ; Berrone2013a ; BPSc ; BPSd ; BBoS ; BBV2019 a PDE-constrained optimization approach is proposed, based on non-conforming meshes, that can be applied both to Problem 4 and Problem 6. In this framework, the problem is rewritten as a minimization problem for a functional measuring the error in fulfilling matching conditions, constrained by local PDEs on each fracture. This approach provides not only a numerical approximation of the solution but also a directly computed approximation of the flux exchanged at traces, which is of interest for many applications. The discretization can be based on different methods: standard finite elements are the simplest choice, and the resulting scheme is denoted as opt-fem. However, as mesh elements arbitrarily cross the traces, the jump of the co-normal derivative of the solution at fracture intersections, still directly computed by the method, can not be correctly represented by non conforming finite elements. Thus the use of local extended finite elements is also considered, being at the basis of the method denoted opt-xfem. When used for advection dominated flow regimes, the supg-stabilized versions of the method are used (femsupg, xfemsupg).

The resulting numerical schemes inherit the mass-conservation properties of the local discrete formulations, thus they are globally but not locally conservative. A huge advantage of this method is that the matrices resulting from the local discrete problems can be computed in parallel, and also the solution can be computed strongly relying on parallel computing. In BSV a MPI-based parallel algorithm is proposed for assembling and solving the discrete problem using the conjugate gradient method on huge networks of fractures, whereas in BDV an implementation suitable for GPGPUs is presented.

4 Examples

In this section, we present three test cases with the aim of validating and comparing the previously introduced models. Extending the work proposed in Scialo2017 , here we analyse the various approaches for time dependent problems of advection-diffusion-reaction, where the advection velocity is computed by means of the same approach solving a diffusion problem. The key aspects of the various schemes will be highlighted and investigated, along with the impact of the lack of conservation of fluxes, both locally and through a trace, that characterizes some of the proposed approaches. For the considered test cases, both local and global quantities will be computed at different time steps, and used to assess and compare the behaviour of the various approaches, such as i) the integral mean in space of the temperature on a fracture , denoted as ; ii) the total flux mismatch on a trace, defined as the integral of the sum of the net total fluxes and entering or leaving the two fractures and meeting at trace , respectively, i.e. ,

with the unit normal vectors to , with fixed orientation on ; and iii) the averaged on the outflow boundary, , denoted as

The considered test cases are designed in order to challenge the proposed approaches with complex geometries and/or realistic models. In Subsection 4.1 we consider the effect of vanishing trace from a simple network, by a sequence of simulations. In the second example presented in Subsection 4.2 a synthetic small network of fractures is considered, to analyse the behaviour of the methods on a more general, yet simple configuration. We finally conclude with a realistic example in Subsection 4.3, where a network of 89 fractures is generated from an extrusion of an interpreted natural outcrop and physically sound values for the various parameters are used.

With the purpose of establishing a standard reference for the analysis of numerical schemes for flow in fractures, test cases in Subsections 4.1 and 4.3 are borrowed from Scialo2017 and adapted to the present context, all the geometrical data being available in FKSdata .

4.1 Vanishing trace between intersecting fractures

As a first test case, named Test Case 1, the same setting of the problem proposed in the example of (Scialo2017, , Subsection 4.3.1) is considered. In the reference, the Darcy problem was tackled, whereas here a non-stationary advection-diffusion problem for the passive scalar is solved.

In this test case the same problem is solved on different geometries. Let us consider a network composed of three fractures, named , , and , as shown in Figure 3. Fracture has a fixed position, whereas fractures , and are displaced, for each different geometry of a same distance along the -direction, such that the length of the intersection line between fractures and , denoted as , progressively reduces from the configuration shown in Figure 5, being instead fixed the intersection between , and , denoted as . In such a way 21 different configurations are obtained, with the length of the vanishing trace ranging from at configuration , as shown in Figure 3, to , at configuration . The geometries corresponding to configurations , and are shown in Figure 5. For each configuration, the Darcy problem, formulated as in 2 or 4, depending on the method, is first solved in order to compute the Darcy velocity , with a null source term. A unitary hydraulic conductivity is used for all the fractures and pressure head boundary conditions are prescribed on the bottom part of (inflow) and (outflow) equal to and , respectively (see Figure 3). On the other portions of the boundary a no-flux boundary condition is imposed. Subsequently, an advection-diffusion problem is solved, obtained from Problem 6 setting the reaction operator , and with null source. We assume a unitary coefficient , a diffusion coefficient equal to and a constant in time unitary Dirichlet boundary condition on the inflow part of the domain boundary, whereas homogeneous Neumann boundary condition are set on the rest of the boundary.

Figure 3: Geometry of Test Case 1. The inflow and outflow portions of the boundary are represented by red and blue lines, respectively.

The meshes used for the space discretization are different for conforming-triangular, conforming-polygonal and non-matching strategies, as discussed above. The mesh for the conforming schemes are generated using the Gmsh Geuzaine2009 library, and two different mesh parameters are used, corresponding to about and elements for the configuration . Clearly, as the length of trace progressively reduces from configuration to , the mesh generation tool tends to increase the number of elements in order to meet the conformity requirement without compromising mesh quality. This process inevitably leads to a large increment of the number of elements from configuration to , for each of the three initial refinement levels. The mesh for the non-matching schemes is obtained through the Triangle Shewchuk1996 software, using two mesh parameters, again corresponding to about , and elements. In this case, instead, since the mesh is independent from the traces, the number of elements is practically unaffected by the change of the geometry from configurations to (small oscillations of few elements are observed as a consequence of the change of the coordinates of fracture ). The polygonal mesh for the mvem approach is finally built gluing together the triangular elements of the conforming mesh, thus aiming at mitigating the increase of mesh elements and degrees of freedom. The number of elements for the various schemes, for each geometrical configuration is reported in Figure 4, on the leftmost column, for the coarse (top) and fine (bottom) meshes. The centre and right columns of Figure 4 report instead the number of degrees of freedom corresponding to each method, which can be used as an indication of the corresponding computational cost. We remark, however, that the computational cost is also largely affected by the approach used to solve linear systems: some of the proposed methods, e.g., have efficient parallel solvers, other take advantage of standard domain decomposition strategies to achieve computational efficiency. Also, the availability of effective preconditioners should be taken into account. Here, however, the main focus is on the analysis of the response of the methods in terms of prediction accuracy versus geometrical and model complexities typical of DFN simulations, thus we refer to the literature of each method for further details on computational efficiency issues.

In the following of this test case, for brevity, the mesh with cells on configuration will be denoted as coarse mesh and the mesh with cells on as the fine. For time discretization an equally spaced mesh with time-step equal to is considered and time-steps are performed, starting from an all zero initial condition.

Figure 4: Number of cells (left), number of DOFs for the Darcy problem (middle) and number of DOFs for the dispersion problem (right) against configuration id for the two mesh parameters of Test Case 1, coarse mesh on the top and fine mesh on the bottom grids.

The solution obtained with the scheme mfem is reported, as an example, in Figure 5 on configurations , and . On the leftmost column the pressure head distribution in the network is shown for the three geometries, whereas the remaining columns depict the temperature distribution at three time-steps corresponding to a time , and , respectively. We can notice that, as the trace between fractures and vanishes, the pressure head distribution displays a steeper gradient around the intersection, and the effective permeability of the network reduces, thus also reducing the penetration depth of the higher temperature zone in the network at fixed time frames. From a computational point of view, simulations become more and more challenging as the solution starts to display steep gradients, especially for methods built on non-conforming meshes, that are not adapted to the geometry, as shown in Figure 6.

Figure 5: Solution of Test Case 1. On the pressure head solution is shown on the leftmost column at configurations , and . On the remaining columns the solution is represented at , respectively. The colour scale spans from to .
Figure 6: Mesh for Test Case 1 on configuration for conforming (left) non-conforming (middle) and coarsened (right) schemes

For all the proposed numerical schemes, at each time-step, the following quantities are computed: the average temperature , the average temperature on the outflow portion of the boundary, and, for non-locally conservative schemes, as the optimization-based methods, the total flux mismatch at each trace.

A reference solution is computed using the mfem method on a mesh much finer than the three considered here for the simulations, counting about cells, almost independently from the configuration id, as cell size is capable of resolving the smallest length of the vanishing trace.

The plots in Figure 7 propose a comparison of average temperature for all the proposed methods against time, on three selected configuration, for all the considered meshes. In the picture, dashed black lines represent relative errors with respect to the reference solution. As expected, all the curves are in very good agreement, also on the coarsest mesh, for configuration , whereas small discrepancies appear for configuration on the coarse mesh that however disappear as the mesh is refined. Larger differences appear, for methods femsupg and xfemsupg, instead for the simulations on configuration . This is expected, since, as mentioned, when the varying trace becomes very small, methods built on non-matching meshes can not rely on the effect of mesh refinement around the vanishing trace which clearly improves representation capabilities of conforming methods. We observe that method mvemup retains good approximation capabilities despite the coarsening.

Figure 7: Average temperature on curves against time for Test Case 1 on the coarsest (first six pictures) and finest (last six pictures) for selected fracture and configurations. Columns refer to the same geometrical configuration: on column 1, on column 2 and on column 3. The average temperature curves refer to fracture .

In order to quantify what is the effect of mesh adaptation due to conformity requirements on the quality of the solution with respect to the effect of the approximations introduced by the non-matching schemes themselves, in Figure 8 (left) the outflow average temperature for the methods opt-xfem and opt-fem is reported against time, for configuration on a perturbation of the coarse mesh used for the conforming methods for configuration , overlapped to the curves of the other schemes, in transparency, on the original meshes. We observe that the curves of methods opt-fem and opt-xfem on the perturbed adapted mesh are now much closer to the curves of other conforming approaches, thus clearly showing the effect of mesh adaptation on the quality of the solution. The solution obtained with the opt-xfem approach on the perturbed mesh at is reported in Figure 8, on the right, along with a detail of the mesh around , clearly showing the non conformity of the mesh. Two aspects are to be remarked: first, the non-conforming approaches are capable of producing reasonable approximations of the solution also on the coarse mesh, which can be greatly improved refining the mesh and still using a fraction of the degrees of freedom required by the conforming approaches, see, e.g., the last two plots in the third column of Figure 7; second, non-conforming approaches allow to freely choose the refinement level of the mesh, thus allowing to efficiently use mesh adaptation strategies, only refining the mesh where required, independently of the geometrical constraints BBV2019 .

Figure 8: Outflow average temperature curve with opt-xfem and opt-fem approaches (left), and solution with opt-xfem at with a zoom of the mesh around (right), both on on the perturbed mesh of configuration .

As the non-conforming approaches opt-xfem and opt-fem are non locally conservative, Figure 9 report the value of against time on the coarsest mesh for the two traces of this example at configurations , and , from left to right, respectively, results with the opt-xfem approach are on the top, results with opt-fem on the bottom. The maximum-in-time absolute values of the total flux on traces and are reported in Table 2, computed on the finest mesh as , and . We can see that relative values of less than are obtained for all times for both methods, with the higher values corresponding to the configuration , as expected. Moreover, for the larger times, mismatch values tend to decrease or to remain constant at a fixed value. Thus this non-local-conservation has a negligible impact on the computed solution. Further, mismatch errors can be reduced by refining the mesh.

The mesh Péclet number for this problem ranges between a maximum of about to a minimum of about on the computational meshes for opt-xfem, opt-fem and mfem methods, thus a Streamline-Upwind-Petrov-Galerkin (supg) stabilization strategy was adopted. As a consequence small overshoots/undershoots in the solution are observed in supg stabilized methods, as well known in the literature, whereas the intrinsic diffusive behavior of methods using upwinding for advection prevents this kind of phenomena.

Figure 9: Total flux mismatch against time for Test Case 1 on traces (top) and (bottom) for configuration (left), (middle) and (right)
0.6072 0.4725 0.251
0.876 0.912 1.300
Table 2: Maximum value in time of total fluxes and on traces and of Test Case 1, respectively, computed with the xfemsupg method on configurations , and .

4.2 Synthetic network

In the second test, labeled Test Case 2, we consider a more complex network composed of 10 fractures with 14 traces, thus being more similar to (a portion of) realistic DFNs, still remaining simple enough to perform analyses on the solutions obtained with the considered schemes. The network is represented in Figure 10, where the inflow and outflow portions of the boundary are also marked. Also in this case the Darcy velocity is first computed solving problem 2 or 4, depending on the chosen numerical scheme, with constant hydraulic conductivity equal to one on all fractures and null source term. A unitary pressure head drop is imposed between the inflow and outflow portions of the border, all other fracture edges being insulated. The Darcy velocity is then used as an input for an advection-diffusion problem for a scalar quantity , as formulated in 6, with null source and reaction terms. A coefficient is chosen, whereas the diffusion coefficient is . A unitary Dirichlet boundary condition is prescribed on the inflow border and all other edges are insulated.

Figure 10: Geometry of Test Case 2 and two views of the network. A red line represents the inflow and a blue thick line the outflow part of the boundary, respectively.

Two meshes are used also for Test Case 2 counting about elements (coarse mesh) and elements (fine mesh), respectively. The mesh Péclet number, related to the coarse mesh for supg stabilized methods, is of about . An equally-spaced mesh is then used for time discretization with steps of length . Also in this case, the initial solution is the null function. An example of the obtained numerical solution with the tpfa method is shown in Figure 11: in the first column on the left the pressure head distribution in the network is represented, solution of the Darcy problem on the coarse (top) and fine (bottom) meshes; the remaining columns depict the solution of the dispersion problem at three selected time frames, corresponding to , and , again on the coarse (top) and fine (bottom) mesh. Coherently the heat flows from the inflow to the outflow by following a tortuous path given by the complex disposition of the fractures and traces. The solution on the coarse grid has a spreader front than on the fine grid due to the artificial diffusivity of the scheme.

Figure 11: Solution of the test case in Subsection 4.2. On the first column the pressure head solution for the two grid refinements. The other columns show at different time: 3.35, 6.25, and 12.5, respectively, for both level of refinement. In all cases the solution is rescaled in the range .

As previously, the average temperature on selected fractures, the average outflow temperature and flux mismatch at traces are used to compare and assess the approximation capabilities of the various schemes. The curves of , and are reported in Figure 12. A reference solution is computed with the mfemsupg method on a mesh counting about cells, and relative error curves with respect to this solution are shown in dashed lines in the Figure 12. We can observe that, despite the network has a larger number of fractures and fracture intersections with respect to Test Case 1, all the methods, in absence of severe geometrical features have good approximation properties, that further improve with mesh refinement. The larger differences are observed for the average outflow temperature curve related to the mvemup approach on the coarse mesh and for the average temperature curves of femsupg again on the coarse mesh. In both cases differences slightly exceed of the reference, and are reduced by mesh refinement. Concerning the mvemup the difference is caused by the coarsening process which reduces the number of elements of the original mesh to about one half. Also, the tpfa method used for advection is observed to have poor performances on irregular polygonal cells, as the ones generated by the coarsening method. Concerning the femsupg method, some discrepancies with respect to the reference are to be expected, as the approach is designed to be computationally inexpensive, nonetheless it is capable of providing satisfactory predictions.

Figure 12: Average temperature on two selected fractures (left and middle columns) and at outflow (right column) curves against time for Test Case 2. Coarse mesh on top, fine mesh at the bottom.
Figure 13: Total flux mismatch against time for all the traces of the network of Test Case 2 on the coarse mesh with methods femsupg (left) and xfemsupg (right).

Curves of the total flux mismatch at the traces are reported against time in Figure 13 for the femsupg and xfemsupg methods. In this picture values of are shown, without labels, for all the traces in the network, highlighting that, in all cases, the errors remain limited in time. Further, maximum-in-time mismatch values are lower than of the total flux, for all the traces; values of the maximum flux with respect to time on each trace are reported in Table 3, computed with the xfemsupg method on the finest mesh.

0.040 0.079 0.150 0.260 0.284 0.047 0.074
0.020 0.393 0.163 0.443 0.061 0.058 0.063
Table 3: Maximum value in time of total flux across each trace of the DFN for Test Case 2 computed with the xfemsupg method on the fine mesh.

4.3 Extruded real outcrop

In this last test case we consider a fracture network generated from an extruded outcrop, located in Western Norway. The test case is inspired by Section 4.4 of Scialo2017 . The network is composed by 89 intersecting fractures resulting in 166 traces. There are 7 non-connected fractures and with no flow boundary conditions, which will not contribute to the solution. The geometry is depicted in Figure 14. The aim of this test case is to validate the proposed numerical schemes in presence of realistic physical parameters of a real fracture network. However, we assume that all the fractures share the same values of hydraulic conductivity and heat diffusion coefficient.

Figure 14: Geometry of the test case in Subsection 4.3. On the top two screenshot of the geometry of the network. On the bottom a sketch of the network with only the largest fractures to present the boundary conditions: on the left in the case of inflow and outflow imposed on two different fractures, on the right on the same fracture. The inflow is represented in red and the outflow in blue.

We consider two distinct problems, that differ from each other from the position of the inflow and outflow boundaries. In the first case (denoted as case 1) the inflow and outflow are imposed on two different fractures, while in the second case (case 2) they belong to the same fracture. See Figure 14 where a sketch of the network is shown with the position of the inflow and outflow boundaries. On the other portions of the boundaries a no flow boundary condition is given. In both cases, we require a simulation grid with roughly 70k elements.

Fractures are immersed in granite and we assume that at the beginning of the simulation the water contained in the fractures is at (). The relations to compute the physical parameters for the simulations are the ones presented in Subsection 2.2. We assume , and . The water and rock physical parameters are reported in Table 4.

Water Rock
Dynamic viscosity
Thermal conductivity
Specific heat capacity
Heat transfer coefficient
Table 4: List of rock and water coefficients for the example in Subsection 4.3.

From these data we obtain, ,

Regarding boundary conditions, for the Darcy problem we impose a pressure head equal to at the inflow boundary and at the outflow boundary, while for the heat problem we impose () at the inflow and zero diffusive flux at the outflow. The simulation time is a year (), divided in 200 time steps.

The conforming computational mesh counts about elements, while the non-matching computational mesh has cells and are shown in Figure 15. It is possible to notice, how, in order to meet the conformity requirement, mesh elements of the conforming mesh are concentrated near the traces, whereas, the non-matching mesh has all elements of equal size evenly distributed in the network. The non matching mesh is characterized by a mesh Péclet number of about for case 1 and for case 2, reached on fracture , in both cases. The non-adapted mesh and the high mesh Péclet number make this example extremely complex for methods built on non-conforming meshes and relying on stabilization for advection dominated flow regimes.

The computed solution with the mfemup scheme is reported in Figure 16 where we display the solution of the Darcy problem on the left and the solution at the end of time evolution on the right, for both the setting of case 1, on top and case 2, at the bottom.

Figure 17 shows the curves against time of the average temperature on the inflow () and outflow () fractures, along with the average temperature on the outflow boundary. Despite the complexity of the geometry and of the model, curves appear in good agreement. For this last example no coarsening was used for the mvemup method, as the poor performances of the tpfa method on polygonal cells, already observed in Test Case 2, have a strong impact on the quality of the solution in this more complex case. The curves related to the mvemup approach on triangular meshes are almost perfectly overlapped to the curves of the mfemup method. Curves of xfemsupg and femsupg are in good agreement with those of the other methods.

Figure 15: Conforming (left) and non-matching (right) mesh for Test Case 3
Figure 16: Solution with the mfemup scheme for the Test Case 3: case 1 on the top; case 2 at the bottom. The first column shows the pressure head solution, the second column the temperature distribution at the end of time-evolution.
Figure 17: Curves of average theta on fractures (left) and and at outflow (right) against time for Test Case 3: case 1 top, case 2 bottom.

5 Conclusions

In this work, we presented a detailed comparative study of several solution strategies for single-phase flow and transport in discrete fracture networks. The proposed numerical schemes are challenged with networks of increasing geometrical complexity and with unsteady advection-reaction-diffusion problems. The characteristics of the various approaches are compared in terms of flexibility in handling geometrical complexity, local and global conservativity and stability to high Péclet numbers.

Methods based on matching grids at the traces trade simplicity in imposing coupling conditions with the lack of control on the number of mesh elements, which is actually constrained by the conformity requirement. On the other hand, non-matching and polygonal based approaches demand ad-hoc discretization strategies but allow full flexibility in meshing. Non-locally conservative schemes might be susceptible to loss in mass in particular at the intersection of fractures. The proposed examples showed that this quantity is small and, for many practical problems it is possible to conclude that it might be acceptable compared to the usual uncertainty in the model parameters. Finally, we have compared schemes that are naturally stable with respect to high Péclet number, with schemes that require a stabilization term to avoid spurious oscillations. Also in this case the obtained solutions are coherent with respect to each other. We can conclude that numerical schemes based on polygonal or non-matching meshes give a good balance in terms of computational cost and accuracy with less limitations compared to the matching cases on more standard grids.


This research has been partially supported by the MIUR project “Dipartimenti di Eccellenza 2018-2022” (CUP E11G18000350001), PRIN project ”Virtual Element Methods: Analysis and Applications” (201744KLJL_004) and by INdAM - GNCS. The second author acknowledges financial support for the ANIGMA project from the Research Council of Norway (project no. 244129/E20) through the ENERGIX program and from the “Visiting Professor Project” of Politecnico di Torino. The first and third author acknowledge financial support from Politecnico di Torino through project “Starting grant RTD”. The authors warmly thank Luisa F. Zuluaga for constructing and providing the real fracture network for the example in Subsection