Comparison of the response to geometrical complexity of methods for unstationary simulations in discrete fracture networks with conforming, polygonal, and non-matching grids

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, two benchmark cases, along with an additional one from a previous work, have been specifically designed and are here proposed and investigated, representing some of the most critical issues encountered in DFN simulations.


Introduction
The movement of liquids in the underground is heavily influenced by the presence of fractures and their relative intersections [25,27,28,30,41]. Fractures are discontinuities (here assumed planar) along which a rock has been broken, mainly due to geological movements or to artificial stimulation [48]. In this work, we are considering only open structures, characterized by a geometrical aperture, Members  that allow a liquid to flow through [28]. 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 [47] 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, called conforming-mesh discretization, may suffer, for example, when intersections lines are very close to each other; see [26,42,43,49,50]. To partially overcome these problems, a possibility is to consider the class of polytopal methods, like mimetic finite difference and polytopic discontinuous Galerkin (DG) that allow grid elements of general shapes. Numerical methods of this type are introduced in [2,3,7,9,10,29,34,35]. A recent work [58] extends the use of the enriched Galerkin (EG) approach [46] to mixed-dimensional problems in porofractured domains.
A second possibility is to keep the intersections explicitly represented, but relaxing the conformity of the edges. This approach, named non-conforming mesh discretization, 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 fewer discrete elements than the conforming mesh approach [6,22,23,[51][52][53]. This method may still suffer in the presence of severe geometries. Also, in this case, the virtual element methods are an interesting option to further decrease the computational cost; see [60,62].
A third family of schemes comprises the so-called nonmatching-mesh 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 [5,[16][17][18][19]. This procedure is independent on the actual numerical schemes, which 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, [24,31,38,39]. An alternative, based on the cut finite element method, is proposed in [21].
The present work extends and enriches the one proposed in [37] 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 number 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, a subset of the previously cited papers that are in the expertise of this work authors, for the comparison to cover most of the combinations discussed and try to assess their performance. 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 follows. 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.

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 an advective field to simulate the transport and diffusion of heat in DFNs.
Fractures are considered non-overlapping planar polygons, which can be connected to other fractures through intersection segments, also called traces. We consider N Ω fractures Ω i ⊂ R 3 with boundary ∂Ω i , which compose the discrete fracture network Ω = ∪ N Ω i=1 Ω i , and we denote its boundary as ∂Ω with outward unit normal n ext , defined on each fracture plane as the unit vector normal to the fracture boundary pointing outward from the fracture polygon.
Given two distinct and intersecting fractures Ω i and Ω j , with i = j , we denote their intersection (trace) by Γ k = Ω i ∩ Ω j . 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 Γ k from 1 to N Γ , the latter being their cardinality. We consider also the function t : {1, . . . , We denote by Γ = ∪ N Γ k=1 Γ k the union of all the traces and by Γ Ω i the set of traces belonging to the fracture Ω i . Moreover, consider a fracture Ω i and a trace Γ k , with k ∈ Γ Ω i . Γ k naturally subdivides Ω i in two sub-regions, indicated by To each of these sub-regions, we associate an outward unit normal perpendicular to Γ k , denoted by n i,+ and n i,− (with n i,+ = −n i,− ) and a trace operator, respectively denoted by tr k i,+ and tr k i,− . Generic trace operators on fracture Ω i are denoted by tr i . An example of a simple DFN with the introduced nomenclature is given in Fig. 1.
We present the Darcy problem in strong and weak form in Section 2.1, to compute the hydraulic head and Darcy

The Darcy model
This section is devoted to a brief presentation of the mathematical models used to describe the hydraulic head field h and Darcy velocity u: further details being available in [16,17,34,37] and references therein. Unknowns and parameters restricted to a fracture Ω i are denoted by a subscript i.
For clarity in the exposition, we start considering a single fracture Ω i . The Darcy model on Ω i reads: find (u i , h i ) such that: Variables, data, and differential operators are defined on the tangent plane of the fracture. K i is the effective hydraulic conductivity tensor, which is symmetric and positive definite, and f i is a scalar source/sink term. In our applications, following lubrication theory [4,56,63], we consider the intrinsic fracture permeability, obtained by the cubic law, scaled by the aperture: k i = 2 i /12, where i is the fracture aperture. Moreover, the effective hydraulic conductivity in Eq. (1a) is isotropic and defined as: with ρ w the fluid density, g the gravity acceleration, and μ the dynamic viscosity of the fluid. The boundary conditions on Ω i are: where ∂Ω i,D and ∂Ω i,N are a partition of ∂Ω i and tr ∂ i is the trace operator on ∂Ω i . The boundary conditions considered here are chosen to be coherent with the examples in Section 4. The generalization to other boundary conditions is straightforward. We assume that ∂Ω i,D = ∅ for at least one i. Equations (1a) and (1b) are well studied in the literature; refer to the aforementioned references.
Let us now consider two distinct fractures Ω i and Ω j forming an intersection Γ k . On both fractures (1a) and (1b) are considered 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 Ω i and Ω j such that t (i, j) = k, read: l∈{i,j } tr k l,+ u l · n l,+ + tr k l,− u l · n l,− = 0 The case of multiple fractures follows immediately from our assumptions. The Darcy equation can be formalized in the following problem.
It is known that solving Problem 1 is equivalent to solve one of the following two weak problems, corresponding to the mixed and primal formulations, respectively. Problem 2 (Darcy equation on a DFN-mixed weak formulation) Let us define the functional space V as Then, the problem reads: Problem 3 (Darcy equation on a DFN-primal weak formulation) Let us define the functional space V as Then the problem reads: find h ∈ V such that

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 Ω i as θ i . The heat equation on Ω i reads: given where T ∈ R 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 [4,56,63]: where i is the fracture aperture, ρ w the fluid density, c w the fluid specific thermal capacity, c e,i = φ i ρ w c w + (1 − φ i )ρ m c m the fracture effective thermal capacity, φ i the fracture porosity, ρ m the density of the rock matrix, c m the specific thermal capacity of the rock matrix, λ e,i = λ φ i w λ 1−φ i m the effective thermal conductivity, γ e,i the effective heat transfer coefficient between fluid and rock, and, finally,θ i is the temperature of the rock matrix, acting as external heat source/sink.
For brevity, we present only boundary conditions conforming to the numerical tests of Section 4. Recalling that ∂Ω i,D is the Dirichlet portion of the boundary of fracture Ω i for the Darcy problem in Section 2.1, let us split ∂Ω i,D into two parts, namely ∂Ω inflow Note that, by Eq. (1b), tr ∂ i (D i ∇θ i +u i θ)·n ext = 0 on ∂Ω i,N . The coupling conditions between two fractures Ω i and The heat equation can be formalized in the following problem. We consider the following primal weak formulation of problem 4.

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 unphysical spurious oscillations in advection dominated flow regimes, whereas other schemes need to adopt stabilization 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 impact of different time advancing schemes is also of interest but is not considered here, the focus being the geometrical complexity and thus the spatial discretization techniques. 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 [44,45]. 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 [32] 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 detail the previous approaches, grouping them according to the coupling at the traces: matching coupling at traces in Section 3.1 and non-matching coupling in Section 3.2.

Matching discretization at traces
Here, advantages and drawbacks of a conforming mesh 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 Fig. 2a and 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 cost for mesh generation might increase and become a severe constraint in complex fracture networks.
To solve the Darcy problem, we rely on different classes of numerical schemes: finite volumes, finite elements in primal and mixed formulations, 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 3; see [1] 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 effective hydraulic conductivity coefficient; however, it is consistent only for K-orthogonal grids. Regarding the class of finite elements, we consider mixed finite elements with the pair of spaces RT 0 − P 0 for the Darcy velocity and piecewise constants for the hydraulic head, denoted by MFEM, defined on a simplicial grid; for further references, see [54,55]. It is well known that RT 0 − P 0 is locally conservative and  [20,59,61,62] and those related to DFN [6][7][8][9]34]. In our analysis, we consider only the lowest order mixed (MVEM) formulation, which can be viewed as a generalization RT 0 − P 0 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 effective hydraulic conductivity variation and have the same grid stencil of the RT 0 − P 0 . 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 [33,34].
To solve the heat equation (Problem 5), 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 P 1 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 mesh discretizations.

Non-matching mesh 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, Fig. 2c.
In [11,12,[14][15][16][17], a PDE-constrained optimization approach is proposed, based on non-conforming meshes, that can be applied both to Problems 3 and 5. 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 P 1 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 nonconforming P 1 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 [19], 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 [13] an implementation suitable for GPGPUs is presented.

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 [37], here we analyze 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 behavior 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 Φ Γ,i and Φ Γ,j entering or leaving the two fractures Ω i and Ω j meeting at trace Γ , respectively, i.e., δΦ Γ = |Φ Γ,i + Φ Γ,j |: with n iΓ the unit normal vectors to Γ , with fixed orientation on Ω i ; and (iii) the averaged θ on the outflow boundary, ∂Ω outflow D , denoted as: Furthermore, we denote the norm of total flux at each trace Γ on Ω i as: The considered test cases are designed in order to challenge the proposed approaches with complex geometries and/or realistic models. In Section 4.1, we consider the effect of a vanishing trace from a simple network, by a sequence of simulations. In the second example presented in Section 4.2, a synthetic small network of 10 fractures is considered, to analyze the behavior of the methods on a more general yet simple configuration. We finally conclude with a realistic example in Section 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 Sections 4.1 and 4.3 are borrowed from [37] and adapted to the present context, all the geometrical data being available in [36].

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 [37, 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 Ω l , Ω r , and Ω c , as shown in Fig. 3.
Fracture Ω l has a fixed position, whereas fractures Ω r and Ω c are displaced, for each different geometry of a same distance along the z-direction, such that the length of the intersection line between fractures Ω l and Ω r , denoted as Γ 0 , progressively reduces from the configuration shown in Fig. 5, being instead fixed at the intersection between Ω r and Ω c , denoted as Γ 1 . In such a way, 21 different configurations are obtained, with the length of the vanishing trace Γ 0 ranging from 0.6 at configuration C0, as shown in Fig. 3, to 0.01, at configuration C20. The geometries corresponding to configurations C0, C10, and C20 are shown in Fig. 5. For each configuration, the Darcy problem, formulated as in (2) or (3), depending on the method, is first solved in order to compute the Darcy velocity u, with a null source term. A unitary effective hydraulic conductivity K is used for all the fractures and pressure head boundary conditions are prescribed on the bottom part of ∂Ω l (inflow) and ∂Ω r (outflow) equal to 1 and 0, respectively (see Fig. 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 5 setting the reaction operator r(·, ·) ≡ 0, and with null source. We assume a unitary coefficient ζ , a diffusion coefficient D equal to 10 −4 , and a constant in time unitary Dirichlet boundary condition on the inflow part of the domain boundary, whereas homogeneous Neumann boundary conditions are set on the rest of the boundary.
The meshes used for the space discretization are different for conforming-triangular, conforming-polygonal, and nonmatching gridding strategies, as discussed above. The mesh for the conforming mesh schemes is generated using the Gmsh [40] library, and two different mesh parameters are used, corresponding to about 10 3 and 10 4 elements for the configuration C0. Clearly, as the length of trace Γ 0 progressively reduces from configuration C0 to C20, 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 C0 to C20, for each of the three initial refinement levels. The mesh for the non-matching mesh schemes is obtained through the Triangle [57] software, using two mesh parameters, again corresponding to about 10 3 and 10 4 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 C0 to C20 (small oscillations of few elements are observed as a consequence of the change of the coordinates of fracture Ω l ). 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 Fig. 4, on the leftmost column, for the coarse (top) and fine (bottom) meshes. The center and right columns of Fig. 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,  others 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 10 3 cells on configuration C0 will be denoted as coarse mesh and the mesh with 10 4 cells on C0 as the fine. For time discretization an equally spaced mesh with time-step equal to 0.05 is considered and 300 time-steps are performed, starting from an all zero initial condition.
The solution obtained with the scheme MFEMUP is reported, as an example, in Fig. 5 on configurations C0, C10, and C20. 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 t = 1.25, t = 2.5, and t = 5, respectively. We can notice that, as the trace between fractures Ω l and Ω r vanishes, the pressure head distribution displays a steeper gradient around the intersection, and the effective permeability of the network reduces, thus also reducing the Fig. 6 Mesh for test case 1 on configuration C20 for conforming (left) non-conforming (middle) and coarsened (right) schemes Fig. 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: C0 on column 1, C10 on column 2, and C20 on column 3. The average temperature curves refer to fracture Ω m 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 Fig. 6.
For all the proposed numerical schemes, at each timestep, the following quantities are computed: the average Fig. 8 Outflow average temperature curve with XFEMSUPG and FEMSUPG approaches (left), and solution with XFEMSUPG at t = 15 with a zoom of the mesh around Γ 0 (right), both on the perturbed mesh of configuration C20 temperature θ Ω , the average temperature θ outflow 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 MFEMUP method on a mesh much finer than the three considered here for the simulations, counting about 6.3 × 10 5 cells, almost independently from the configuration id, as cell size is capable of resolving the smallest length of the vanishing trace.
The plots in Fig. 7 propose a comparison of average temperature for all the proposed methods against time, on three selected configurations, 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 C0, whereas small discrepancies appear for configuration C10 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 C20. 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.
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   Fig. 8 (left) the outflow average temperature for the methods XFEMSUPG and FEMSUPG is reported against time, for configuration C20 on a perturbation of the coarse mesh used for the conforming methods for configuration C20, overlapped to the curves of the other schemes, in transparency, on the original meshes. We observe that the curves of methods FEMSUPG and XFEMSUPG 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 XFEMSUPG approach on the perturbed mesh at t = 15 is reported in Fig. 8, on the right, along with a detail of the mesh around Γ 0 , 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 Fig. 7; second, nonconforming-mesh 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 [12]. As the non-conforming mesh approaches XFEMSUPG and FEMSUPG are non locally conservative, Fig. 9 reports the value of δΦ Γ against time on the coarsest mesh for the two traces of this example at configurations C0, C10, and C20, from left to right, respectively, results with the XFEMSUPG approach are on the top, results with FEMSUPG on the bottom. The maximum-in-time total flux norms on traces Γ 0 and Γ 1 are reported in Table 2, computed on the finest mesh asΦ Γ 0 = max t Φ Γ 1 ,c +Φ Γ 1 ,r , used to compute relative flux mismatches. We can see that relative values of less than 1% are obtained for all times for both methods, with the higher values corresponding to the configuration C20, 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. Furthermore, mismatch errors can be reduced by refining the mesh.
The mesh Péclet number for this problem ranges between a maximum of about 6 × 10 2 to a minimum of about 100 on the computational meshes for XFEMSUPG, FEMSUPG, and MFEMSUPG methods, thus a streamline upwind Petrov-Galerkin (SUPG) stabilization strategy was adopted. As a consequence, small overshots/undershots 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.

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 Fig. 10, where the inflow and outflow portions of the boundary are also marked. Also, in this case, the Darcy velocity u is first computed solving Problem 2 or 3, depending on the chosen numerical scheme, with constant effective hydraulic conductivity equal to 1 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 advectiondiffusion problem for a scalar quantity θ , as formulated in (5), with null source and reaction terms. A coefficient ζ = 1 is chosen, whereas the diffusion coefficient is D = 10 −4 . A unitary Dirichlet boundary condition is prescribed on the inflow border and all other edges are insulated.
Two meshes are used also for test case 2 counting about 10 3 elements (coarse mesh) and 4 × 10 4 elements (fine mesh), respectively. The mesh Péclet number, related to the coarse mesh for SUPG-stabilized methods, is of about 100. An equally spaced mesh is then used for time discretization with 500 steps of length 0.05. 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 Fig. 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 t = 3.35, t = 6.25, and t = 12.5, again on the coarse (top) and fine (bottom) meshes. 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. As previously, the average temperature θ Ω on selected fractures, the average outflow temperature θ outflow , and flux mismatch δΦ Γ at traces are used to compare and assess the approximation capabilities of the various schemes. The curves of θ Ω 1 , θ Ω 3 , and θ outflow are reported in Fig. 12. A reference solution is computed with the MFEMSUPG method on a mesh counting about 2 × 10 4 cells, and relative error curves with respect to this solution are shown in dashed lines in Fig. 12. We can observe that, despite the network having a larger number of fractures and fracture intersections with respect to test case 1, all the methods, in the 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 10% 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.
Curves of the total flux mismatch at the traces are reported against time in Fig. 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. Furthermore, maximum-in-time mismatch values are lower than 1% of the maximum in time total flux norm, for   Table 3, computed with the XFEMSUPG method on the finest mesh.

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 [37]. The network is composed by N Ω = 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 Fig. 14. The aim of this test case is to validate the proposed numerical schemes in the presence of realistic physical parameters of a real fracture network. However, we assume that all the fractures share the same values of effective hydraulic conductivity and heat diffusion coefficient.
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 Fig. 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 we have thermal equilibrium, the water contained in the fractures is at θ i = 353.15 K (80 • C), for all i = 1, . . . , N Ω . The value of θ i , which represents the constant temperature of the rock matrix, is at the same value,θ i = θ i , for all i = 1, . . . , N Ω . The relations to compute the physical parameters for the simulations are the ones presented in Section 2.2. We assume i = 2 mm ∀i, and φ i = 0.95 ∀i. The water and rock physical parameters are reported in Table 4.
From these data, we obtain: ∀i = 1, . . . , N Ω ,   The conforming computational mesh counts about 7 × 10 4 elements, while the non-matching computational mesh has 2 × 10 4 cells and is shown in Fig. 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 3 × 10 4 for case 1 and 6.6 × 10 4 for case 2, reached on fracture Ω 0 , in both cases. The nonadapted mesh and the high mesh Péclet numbers 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 Fig. 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 (Ω 0 ) and outflow (Ω 1 ) 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. Flux mismatch values have been computed for these methods for those traces having a non-negligible maximum-in-time total flux norm, i.e., traces Γ k withΦ Γ k greater than 1% of max kΦΓ k , i.e., the maximum total flux norm through all the traces in the network. Indeed, in this more complex example, some traces are excluded from the main flow path, and thus have an almost zero net flux. The computed relative mismatch values are lower than 1% also for this example.

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 advectionreaction-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. The main focus of the present work is the response of the considered numerical schemes to geometrical complexity and thus other relevant aspects such as the impact of different time integration schemes, or of strongly heterogeneous fracture properties, such as fracture transmissivity, and also a more detailed discussion on computational cost, have not been addressed.
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. The obtained results show that the presence of traces much smaller than the characteristic size of the network is a major source of complexity. Methods based on non-matching meshes and polygonal meshes are both capable of producing reliable approximations, also in time-evolving simulations. Lack of local conservation at traces appears to have an acceptable impact on the quality of the solution and on the measured global quantities, compared with other sources of error and to the usual uncertainty in model parameters. This is a non-trivial result, mainly for time-evolving simulations, where an accumulation of local errors could have a strong impact on the solution. Peclét-related instabilities are effectively cured by SUPG stabilization, providing solutions comparable with those obtained by methods intrinsically stable thanks to numerical diffusivity.
Non-matching mesh discretizations are thus suggested to provide reliable predictions in complex geometrical configurations, and benefit from mesh adaptivity to recover high accuracy levels. The VEM-based approach based on the coarsening of a conforming mesh allows to mitigate the impact of handling large meshes for complex geometries, and is promising of being effective on complex geometrical configurations in combination with more advanced discretization strategies for unstationary advection-reaction simulations.