Benchmarks for single-phase flow in fractured porous media Bernd Flemisch1, Inga Berre2, Wietse Boon2, Alessio Fumagalli2, Nicolas Schwenck1, Anna Scotti3, Ivar Stefansson2, and Alexandru Tatomir4 1Department of Hydromechanics and Modelling of Hydrosystems, University of Stuttgart, 7 1 Pfaffenwaldring 61, 70569 Stuttgart, Germany, [email protected] 0 2 2Department of Mathematics, University of Bergen, All´egaten 41, 5007 Bergen, Norway n a 3Laboratory for Modeling and Scientific Computing MOX, Politecnico di Milano, J 5 p.za Leonardo da Vinci 32, 20133 Milano, Italy ] 4Department of Applied Geology, Geosciences Center, University of Go¨ttingen, A N Goldschmidtstrasse 3, 37077 Go¨ttingen, Germany . h t January 9, 2017 a m [ 1 This paper presents several test cases intended to be benchmarks for numer- v 6 ical schemes for single-phase fluid flow in fractured porous media. A number 9 of solution strategies are compared, including a vertex and a cell-centered fi- 4 1 nite volume method, a non-conforming embedded discrete fracture model, a 0 primal and a dual extended finite element formulation, and a mortar discrete . 1 fracture model. The proposed benchmarks test the schemes by increasing the 0 7 difficulties in terms of network geometry, e.g. intersecting fractures, and physi- 1 cal parameters, e.g. low and high fracture-matrix permeability ratio as well as : v heterogeneous fracture permeabilities. For each problem, the results presented i X by the participants are the number of unknowns, the approximation errors in r the porous matrix and in the fractures with respect to a reference solution, and a the sparsity and condition number of the discretized linear system. All data and meshes used in this study are publicly available for further comparisons. 1 Introduction In porous-media flow applications, the domains of interest often contain geometrically anisotropic inclusions and strongly discontinuous material coefficients that can span sev- eral orders of magnitude. If the size of these heterogeneities is small in normal direction compared to the tangential directions, these features are called fractures. Fractures can 1 act both as conduits and barriers and affect flow patterns severely. Target applications concerning fractured porous-media systems in earth sciences include groundwater resource management, renewable energy storage, recovery of petroleum resources, radioactive waste reposition, coal bed methane migration in mines, and geothermal energy production. The analysis and prediction of flow in fractured porous media systems are important for all the aforementioned applications. Many different conceptual and numerical models of flow in fractured porous-media systems can be found in the literature. Even though fractured porous-media systems have been of interest to modelers for a long time, they still present challenges for simulators. During the last 70 years, different modeling approaches have been developed and gradually improved. Comprehensive reviews can be found in Berkowitz (2002); Dietrich et al. (2005); Hoteit and Firoozabadi (2008); Neumann (2005); Sahimi(2011);SinghalandGupta(2010). Roughly, thefracturedporousmediasystemsare classified in two broad categories: discrete fracture-matrix (DFM) models and continuum fracture models. Within this paper, we will only consider DFM models. The DFM models consider flow occurring in both the fracture network and the surround- ing rock matrix. They account explicitly for the effects of individual fractures on the fluid flow. An efficient way to represent fractures in DFMs is the hybrid-dimensional approach, e.g. Helmig (1997); Firoozabadi and Monteagudo (2004); Karimi-Fard et al. (2004); Mar- tin et al. (2005); Reichenberger et al. (2006), where fractures in the geometrical domain are discretized with elements of co-dimension one with respect to the dimension of the surrounding matrix, such as one-dimensional elements in two-dimensional settings. The aforementioned classical DFM approaches all rely on matching fracture and matrix grids in the sense that a fracture element coincides geometrically with co-dimension-one mesh enti- ties, i.e. faces of matrix grid elements. In addition to the classical models, several so-called non-conforming DFM models have been developed in recent years, such as EDFM (Moinfar et al., 2014a; Hajibeygi et al., 2011), XFEM-based approaches (D’Angelo and Scotti, 2012; Schwenck et al., 2015; Huang et al., 2011), or mortar-type methods (Frih et al., 2012). Benchmarking represents a methodology for verifying, testing and comparing the model- ing tools. Various codes have been developed by academic institutions or companies based on different conceptual, mathematical, and numerical models. Even though benchmarking studies are increasing in all fields of engineering and workshops have been organized around specific problems (e.g. Class et al. (2009)), there are still only a limited number of studies of this type in the field of geoscience. Some are related to a specific application and are flexibleastohowtheproblemismodeledintermsofassumptionsregardingthephysicsand the selection of the domain, see Dahle et al. (2010); Nordbotten et al. (2012); Caers (2013); Kolditz et al. (2015). Others (De Dreuzy et al., 2013; Caers, 2013), like ours, focus on the comparison of numerical schemes. One of the common requirements when selecting the test problems for comparing numerical schemes is that they allow the examination of the capabilities of each of the compared methods. Therefore, our benchmark study proposes a set of problems starting from simple geometries and then gradually increasing the geo- metrical complexity. The test problems are specifically selected to make clear distinctions between the different methods. The main focus of this work is to use existing and new computational benchmarks for fluid flow in fractured porous media to compare several DFM-based numerical schemes in a systematic way. We would also like to invite the scientific community to follow up on this study and evaluate further methods by means of the proposed benchmarks. In order to facilitate this, the paper is accompanied by grid and result files in the form of a Git repository at https://git.iws.uni-stuttgart.de/benchmarks/fracture-flow. 2 The remainder of this paper is organized as follows. In Section 2, we formulate the model problem in terms of the partial differential equation to be solved. The participating DFM models are described in Section 3. The central Section 4 proposes the benchmarks and compares the results of the different methods. Finally, Section 5 concludes with a summary and outlook. 2 The model problem We are considering an incompressible single-phase flow through a porous medium, assumed to be described by Darcy’s law, resulting in the governing sytem of equations u = −Kgradp, (1a) divu = q, (1b) in an open bounded domain D ⊂ RN, subject to boundary conditions p = p on ∂D , (1c) D D u·n = q on ∂D , (1d) N N with ∂D = ∂D ∪∂D and ∂D ∩∂D = ∅. In equations (1), u denotes the macroscopic D N D N fluid velocity whereas K and p stand for absolute permeability and pressure. Let us assume that D contains several fractures, that all together constitute a single domain Γ of spatial dimension N such that Γ ⊂ D, which is a possibly unconnected, open subset of D. The surrounding porous rock, namely, the remaining part of D, is called Ω = D \ Γ. Assuming that the fracture aperture ε at each point of Γ is small compared to other characteristic dimensions of the fractures, the full-dimensional domain Γ can be reduced to the (N−1)-dimensional fracture network γ. This reduction involves modeling choices resulting in different hybrid-dimensional problem formulations that form the basis for the methods presented in the following section. 3 Participating discretization methods Within this section, the discretization methods participating in this benchmark study are described. The purpose of this article is the comparison of well-known, established and/or at least published methods. Therefore, only the most significant aspects of each method are summarized. We do not show a comparison against analytical solutions here. The analysis of the methods and theoretical results such as proofs of optimal convergence can be found in the corresponding references. A summary of all participating methods is provided in Table 1. In the sequel, we will denote with d.o.f. the degrees of freedom associated to a specificmethod. Weindicatealsothetypeofconformityrequiredtothecomputationalgrid with respect to the fractures and the assumption that the pressure is considered continuous across the fractures. 3.1 Vertex-centered, continuous-pressure, conforming lower-dimensional DFM (Box-DFM) The lower-dimensional representation of fractures allows easier mesh generation in com- parison to the equi-dimensional approach, as it circumvents the appearance of very small 3 method d.o.f. frac-dim conforming p-cont. Box-DFM p (vert) dim-1 yes yes CC-DFM p (elem) dim-1 yes no EDFM p (elem) dim-1 no yes mortar-DFM p (elem), u (faces) dim-1 geometrically no P-XFEM p (vert) dim-1 no no D-XFEM p (elem), u (faces) dim-1 no no MFD p (faces) dim geometrically no Table 1: Participating discretization methods. elements when discretizing the interior of the fracture (i.e. , within the fracture width). The conforming mesh generation algorithm honors the geometrical characteristics of the fracture system. Conform meshing implies that the fractures are discretized with a set of line elements (in a 2D domain) that are also the edges of the triangular finite elements. The spatial discretization in Box-DFM is performed with the Box method, a vertex- centered finite-volume method proposed in, e.g. Helmig (1997) which combines the advan- tages of finite element and finite volume grids, allowing unstructured grids and guaran- teeing a locally conservative scheme (Reichenberger et al. (2006)). Figure 1 illustrates a two-dimensional representation of the dual-grid with two finite elements E and E sharing 1 5 the same edge (ij ) that represents a lower-dimensional fracture with the aperture ε . The 1 ij main characteristic in terms of the fractured system is that the pressure is required to be continuous, in particular in those vertices whose control volumes overlap both fracture and matrix regions. The Box-DFM method used for this paper is implemented in the open-source numerical simulator DuMux. A detailed description of the conceptual, mathematical and numerical model and code implementation is published in Tatomir (2012). The Box-DFM simulation code used for the benchmark studies is publicly available under https://git.iws.uni- stuttgart.de/dumux-pub/Flemisch2016a.git. 3.2 Cell-centered, discontinuous-pressure, conforming DFM (CC-DFM) The control volume finite difference method uses a two-point flux approximation (TPFA) based on the cell-center pressure values for the evaluation of the face fluxes, and is a widely applied and standard method for simulation of flow in porous media. The domain is partitioned with fractures coinciding with the interior faces between matrix cells just as describedinSection3.1. Thefluxoverthefacebetweenmatrixcellsiandj isapproximated by u = T (p −p ), (2) ij ij i j where p and p are the pressures in the neighboring cells and T is the face transmissibility, i j ij computed as the harmonic average of the two half transmissibilities corresponding to the face and the two cells. The half transmissibility of cell-face pair i is in turn given as A n(cid:62)K α = i i i ·d , (3) i d(cid:62)d i i i where A and n are the area and unit normal vector of the face, K is the permeability i i i assigned to the cell and d is the distance vector from cell center to face centroid. i 4 Figure 1: Conceptual representation of the Box-DFM method: (left-hand side) The dual finite element and finite volume mesh from which the control volume B around i nodei iscreated. Nodeiissurroundedbynodes{j , j , j , j , j }, wheresegment 1 2 3 4 5 ij representsbothafractureandasharedFEedge; (right-handside)Sub-control 1 volume (SCV) bE1 in element E has barycenter G and the mid-points of the i 1 1 edges ij and ij are M , respectively M . The SCV face fE1 is the segment 1 2 ij1 ij2 ij1 G M which contains the integration point xE1 where the normal vector nE1 is 1 ij1 ij1 ij1 applied. In addition to the unknowns given at the centroids of the matrix cells, unknowns are associated to the centroids of the fracture cells. The fracture cells are associated with apertures, which multiplied with the length give the volume of these cells. The aperture is also used to construct hybrid faces for the matrix-fracture interfaces. These faces, parallel to the fracture but displaced half an aperture to either side, enable us to compute the half transmissibilities between the fracture cell and the matrix cells on the two sides. These faces are indicated by the dashed blue lines in Figure 2, where the computational domain is superimposed on the geometrical grid. The result is a hybrid grid with fractures which are lower dimensional in the grid, but equidimensional in the computational domain at the cost of a small matrix volume error corresponding to the overlap of the matrix cells with the fracture cells. Following the method proposed by Karimi-Fard et al. (2004), the intermediate fracture intersection cell drawn with dashed red lines in Figure 2 is removed, leading to direct coupling of the fracture cells neighbor to the intersection. The purpose of this is both to obtain a smaller condition number and to avoid severe time-step restrictions associated with small cells in transport simulations. To each new face between cell i and j, face transmissibilities are assigned, calculated using the star delta transformation as described in Karimi-Fard et al. (2004): α α i j T = , (4) ij n (cid:80) α k k=1 withndenotingthenumberoffracturecellsmeetingattheintersection. Asthiselimination disregards all information on the permeability of the intersection, it should be used with caution in cases of crossing fractures of different permeability. We encounter this feature in 4.3, and include results both with and without the elimination for one of the test cases presented in that section. InspiredbytheCC-DFMmethodbyKarimi-Fardetal.(2004)presentedabove, amethod 5 (a) (b) Figure 2: (a)Conceptualdecompositionofthedomainaccordingtoelementdimensionwith the matrix depicted in black, fractures in blue and their intersections in red. (b) The computational domain of the CC-DFM. Dashed lines are faces of the fracture cells. based on the multi-point flux approximation has also been developed Sandve et al. (2012). The MPFA variant of the method reduces errors associated with the TPFA approach for grids that are not close to K-orthogonal, and avoids errors related to the splitting of the fluxes in the star-delta transformation. We refer to Sandve et al. (2012) for a thorough comparison of the TPFA and MPFA CC-DFM approaches. The implementation of both methods is available in the open-source Matlab Reservoir Simulation Toolbox http://www. sintef.no/projectweb/mrst/. 3.3 Continuous-pressure, non-conforming embedded DFM (EDFM) Recently, non-conforming methods for the treatment of lower-dimensional fractures have been developed, for example in Moinfar et al. (2014a, 2011); Hajibeygi et al. (2011), to avoid the time-consuming construction of complex matrix grids which explicitly represent the fractures. They are mostly used in the context of single and multi-phase flow simula- tions for petroleum engineering applications and require the normal fracture permeability to be orders of magnitude higher than the matrix permeability, as in the case of enhanced reservoir exploitation and fractures stimulation. In this field of applications corner-point grids are normally employed to describe the geological layers, e.g. different rock type, of the reservoir. An adaptation of such computational grids to the fractures could be unaffordable for real cases. The numerical method belongs to the family of two-point schemes, where a one-to-one connection between the degrees of freedom is considered through the trans- missibility concept (Eymard et al. (2000)). References on the embedded discrete fracture method (EDFM) can be found, for example, in Li and Lee (2008); Panfili et al. (2013); Moinfar et al. (2014b); Panfili and Cominelli (2014); de Araujo Cavalcante Filho et al. (2015); Fumagalli et al. (2016). In practice, the meshes of the fractures are generated on top of the rock grid so that each rock cell cut by fractures contains exactly one fracture cell per fracture. Intersections between fractures are computed without affecting the creation of the grids of fractures and rock and used to compute approximate transmissibilities between different fracture 6 cells. See Figure 3 as an example. A degree of freedom that represents a pressure or Figure 3: Example of meshes, for both fractures and rock matrix, suited for EDFM. The rock matrix is considered as a background mesh. Each fracture cell is represented by two blue dots and the green dots are the non-matching intersection among fractures. a saturation value is assigned to each matrix cell and to each fracture cell. This means that transmissibilities between matrix and fracture cells, as well as those between different fracture cells, need to be computed. We compute the transmissibility between a fracture cell and a matrix cell T and the half-transmissibility T between two intersecting fracture fm i cells (related to the fracture i) through the following approximate expressions n(cid:62)K·n k ε f f i i T = A and T = s . fm i d d f,m i,s Here A is the measure of the fracture cell in the current rock cell, n is the normal of the f fracture cell and d is an average distance between the fracture cell and the matrix cell, f,m see Li and Lee (2008). For the fracture-fracture transmissibility, s indicates the measure of the intersecting segment, k the scalar permeability of the fracture, ε the aperture and d i i i,s istheaveragedistancebetweenthefracturecellandtheintersectingsegment. Thestandard harmonic average is considered to compute the transmissibility between the two fracture cells. Standard formulae for fracture-fracture as well as matrix-matrix transmissibilities are computed by means of a two-point flux approximation. It is worth to notice that the recent extension of EDFM called Projection-based EDFM (pEDFM), proposed in Tene et al. (2016), is also able to handle low permeable fractures. Finally, even if the proposed benchmark cases are two-dimensional the method can be extended to three dimensions without any additional constraints. 3.4 Cell-centered, discontinuous-pressure, geometrically-conforming mortar DFM (mortar-DFM) The key concept behind the mortar-DFM, as described more thoroughly in Boon et al. (2016), is the idea that fractures can be considered as interfaces between different sub- domains. This has been explored previously by Martin et al. (2005); Frih et al. (2012), amongothers. Inthiscontext, itisinterestingtoconsiderdomaindecompositiontechniques such as the mortar method to model flow through the fractured porous medium. The mortar method is generally used to couple equations in different sub-domains by introducing a so-called mortar variable, defined on the interface. In case of modeling 7 fracture flow, a well-explored choice of the mortar variable is the fracture pressure (Martin et al., 2005). The method considered here, however, uses as the mortar variable the flux, which leads to a stronger sense of mass conservation for flows between the matrix and fractures. One of the main advantages of the close relationship to mortar methods is the capability to handle non-matching grids. In particular, two sub-domains bordering a fracture can be meshed independently on both sides, as illustrated in Figure 4. The Figure 4: The mortar-DFM allows for non-matching grids along fracture interfaces. Frac- ture and matrix flows are coupled using a mortar variable, defined on a coarser grid (green dots). difficulty in mesh generation is then relieved significantly since only the geometry of the fractures needs to be respected. Byconstruction, themortar-DFMisapplicabletoproblemsinarbitrarydimensions. The governing equations in the matrix and the fractures (as well as fracture intersections in 3D) are identical and thus all fractures, intersections and tips are handled in a unified manner. Consequently, although only two-dimensional problems are considered in this case study, the discretization scheme is not at all limited to the presented benchmark problems and 3D cases can easily be considered. With the use of mixed finite elements, mass is conserved locally in the matrix, fractures, and fracture intersections. The flux u and pressure p are modeled using the lowest-order Raviart-Thomas elements and piecewise constants, i.e. RT −P . Additionally, the mor- 0 0 tar variable is given by piecewise constants on a separately generated, lower-dimensional, mortar grid. This grid matches with the surrounding grids in case of matching grids and is coarser otherwise (Boon et al., 2016). The resulting mixed finite element formulation is a saddle-point problem, which may be challenging to solve numerically. To relieve this, the flux variables may be eliminated through hybridization, which leads to a less expensive scheme containing solely the cell-center pressures. Two implementations of the method have been developed, both of which are used in this benchmark study. The first version, implemented in MATLAB, has the capability of handling non-matching grids along fractures for two-dimensional problems. It is well- suited for simpler geometries, containing relatively few fractures, such as those considered in Benchmarks 1-3. The second version has been implemented for 3D problems and higher- order spaces on matching grids using the open-source finite element library FEniCS (Logg et al., 2012). This version is more efficient for complex cases such as Benchmark 4. 8 3.5 Discontinuous-pressure, non-conforming primal XFEM (P-XFEM) The primal XFEM method participating in this benchmarking study is described in detail in Schwenck (2015), see also Flemisch et al. (2016); Schwenck et al. (2015). The method is based on the hybrid-dimensional problem formulation investigated in Martin et al. (2005), where conditions for the coupling between fracture and matrix are derived: {{u ·n}} = k /ε p (5a) m γ f,n m γ (cid:16)(cid:74) (cid:75) (cid:17) ξ u ·n = k /ε {{p }} −p (5b) 0 m γ f,n m γ f (cid:74) (cid:75) Here, the subscripts “m” and “f” indicate matrix and fracture quantities, while {{·}} and γ · denote the average and the jump of a matrix quantity over the fracture γ, respectively. γ (cid:74) (cid:75)The coupling conditions (5) can be used to define a source term for the fracture flow problem,whiletheyyieldaninterfaceproblemforthematrixdomain. Forthediscretization of this interface problem, the methodology presented in Hansbo and Hansbo (2002) is used, which amounts to applying the eXtended Finite Element Method (XFEM). Together with an independent standard discretization of the lower-dimensional fracture problem, this yields a hybrid-dimensional, non-conforming primal XFEM-based method. The XFEM space is built enriching the standard Lagrangian P (or Q for quads) finite-element spaces, 1 1 whose degrees of freedom are located at the vertices of the full-dimensional grid of the matrix Ω and the lower-dimensional grid of the fracture γ. A representative example of Figure 5: Example of meshes, for both fractures and rock matrix, suited for P-XFEM. The fracture grid vertices are indicated by the blue dots. matrix and fracture grids is illustrated in Figure 5. Unlike the EDFM method, see Figure 3, the fracture grid vertices can be placed arbitrarily without taking into account the matrix grid. On the other hand, the method requires matching fracture branch grids in the form of vertices placed at the fracture intersections. In particular, special care has to be taken of intersecting and immersed fractures (Schwenck et al., 2015). The method is implemented on top of the DUNE framework Bastian et al. (2008) and the discretization module DUNE-PDELab Bastian et al. (2010). For the enrichment of the finite-element spaces in the context of XFEM, the modules DUNE-Multidomain and DUNE-MultidomaingridareemployedMu¨thing(2015). ThesimulationcodefortheXFEM approach and for the benchmarks studied here is publicly available under https://git. iws.uni-stuttgart.de/dumux-pub/Flemisch2016a.git. 9 3.6 Discontinuous-pressure, non-conforming dual XFEM (D-XFEM) The dual XFEM method participating in his benchmark is based on D’Angelo and Scotti (2012). The method, originally derived for a domain cut by one fracture, was further developedinFormaggiaetal.(2014),FumagalliandScotti(2014)toaccountforintersecting fractures with different permeabilities. The same equations and coupling conditions as for the primal XFEM are used, but in a dual formulation where Darcy law and mass conservation give rise to a saddle-point problem for the fluid mean velocity and pressure, bothinthefractureandinthesurroundingmedium. Moreover, unlikethepreviousmethod, this method employs triangular/tetrahedral grids. The usual lowest order RT −P pair for 0 0 velocityandpressureisenrichedfollowingHansboandHansbo(2002)intheelementsofthe porous medium cut by a fracture, or in the elements of a fracture at the intersection with other fractures. Indeed, triangular/tetrahedral grids are arbitrarily cut by triangulated lines/surfaces in 2D and 3D respectively. These surfaces can, in turn, intersect each other in a non-conforming way, as shown in Figure 6. Figure 6: A portion of the grid cut by two fractures: in the two dimensional case they can split the elements in two (grey), three (yellow), or four (red) independent parts, where the restrictions of the basis functions are defined. The fracture grids are irrespective of the bulk grid and of each other, i.e. the intersection point i is not p a point of the grid. In the current implementation of the method no special enrichment is added in the bulk elements containing the fracture tips. Instead, fractures are artificially extended up to the boundary of the domain, and in the extension we prescribe the same permeability of the surrounding porous medium to obtain a “virtual” fracture with no effects on the flow. The method has been implemented on the basis of the Getfem++ library, http:// download.gna.org/getfem/html/homepage/, whichprovidessupportforthecomputation of the intersections and the quadrature on sub-elements thanks to an interface with QHull, http://www.qhull.org/. 3.7 Reference Solutions calculated with mimetic finite differences (MFD) The reference solutions are computed on very fine grids that discretize both matrix and fractures by full-dimensional triangular or quadrilateral elements. A mimetic finite differ- 10