Regular Article
Dual Virtual Element Methods for Discrete Fracture Matrix models
Department of Mathematics, University of Bergen, 5007
Bergen, Norway
^{*} Corresponding author: alessio.fumagalli1984@gmail.com
Received:
5
December
2017
Accepted:
6
February
2019
The accurate description of fluid flow and transport in fractured porous media is of paramount importance to capture the macroscopic behavior of an oil reservoir, a geothermal system, or a CO_{2} sequestration site, to name few applications. The construction of accurate simulation models for flow in fractures is challenging due to the high ratio between a fracture’s length and width. In this paper, we present a mixeddimensional Darcy problem which can represent the pressure and Darcy velocity in all the dimensions, i.e. in the rock matrix, in the fractures, and in their intersections. Moreover, we present a mixeddimensional transport problem which, given the Darcy velocity, describes advection of a passive scalar into the fractured porous media. The approach can handle both conducting and blocking fractures. Our computational grids are created by coarsening of simplex tessellations that conform to the fracture’s surfaces. A suitable choice of the discrete approximation of the previous model, by virtual finite element and finite volume methods, allows us to simulate complex problems with a good balance of accuracy and computational cost. We illustrate the performance of our method by comparing to benchmark studies for twodimensional fractured porous media, as well as a complex threedimensional fracture geometry.
© A. Fumagalli & E. Keilegavlen, published by IFP Energies nouvelles, 2019
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
Fractures and faults can strongly influence fluid flow in a porous media, acting, depending on their permeability and porosity, as a preferential path or a barrier. As fracture aperture is several orders of magnitude smaller than any other characteristic size in the domain, fracture modeling is one of the main challenges in subsurface problems.
Geological movements, chemical reaction, or infilling processes may substantially alter the local orientation and composition of the material present in the fractures, leading to anisotropies and strong heterogeneities in both fractures and their intersections. It is thus crucial to be able to include also these phenomenological aspects in the conceptual model.
Applications where fractures can be determinant for reservoir behavior include exploitation of geothermal system, CO_{2} storage and sequestration, enhanced oil recovery, and nuclear waste disposal. In these fields the solution may dramatically depend on the presence of fractures, thus a correct derivation of suitable mathematical models and their accurate numerical solution is essential.
Our computational model is based on a Discrete Fracture Matrix (DFM) approach, e.g. [1], which aims to balance computational cost and modeling accuracy in the presence of fractures on multiple scales. Small scale fractures are incorporated in the matrix permeability by analytical or numerical upscaling techniques (e.g. see [2–4] and [5–7]), thus only macroscopic fractures and faults are considered and explicitly described. However, the computational cost associated with their equidimensional representation in the grid is still prohibitively high. Following the idea presented, among the others, in references [8–15] fractures and faults are represented as lowerdimensional objects embedded in the rock matrix. This approach is commonly denoted reduced models, hybriddimensional models, and mixeddimensional models. In this approach, fracture aperture becomes a coefficient in the equations, instead of a geometrical constraint for grid generation. The flow equations are averaged along the normal direction of the fractures, yielding a new set of reduced models that are suited for the codimensional, that is, lowerdimensional, description of fluid flow. Suitable coupling conditions are necessary to exchange information between all the objects (rock matrix, fractures, and their intersections).
The reduction of dimensionality is an essential step to overcome most of the difficulties associated to the problem. However, in the presence of several fractures that together form a complex network, more advanced numerical techniques are crucial to obtain an accurate solution in a reasonable amount of time. The main aspect is the geometrical treatment of the fracture grid with respect to the rock matrix grid. Popular approaches can be classified according to whether the grid in the rock matrix matches the fracture surfaces or not. The matching case can further be split into methods that assume full conformity between the rock matrix and fracture grid, and methods that relaxe this assumption.
In the class of conforming methods, where fracture grids are composed by faces of the rock matrix grid, several numerical schemes have been considered, including finite elements [13, 16–19], finite volumes [6, 20–23], Mimetic Finite Differences (MFD) [24, 25], and the newly introduced Virtual Finite Elements (VEM) [26–28]. For the latter see also [29–33] for the nonfractured case. All these methods highlight specific advantages for example related to local mass conservation, capability to be implemented in standard software packages, or relax some constraints on grid cell shapes to name a few.
In the class of nonconforming discretization, the main tool is a mortar coupling between the rock matrix grid and the fracture grids. The rock mesh is constrained with the position of the fractures but not strictly with the actual fracture meshes, and vice versa. Some examples are reported in [10, 27, 34] where different type of mortar variables as well as numerical schemes are considered. The mortar variables allow for a generalized formulation that unifies conforming and nonconforming approaches [35].
Finally, a fully nonmatching coupling among the grids requires adhoc solutions to establish a communication between the fractures and the rock matrix. One possibility is the class of eXtended Finite Element Methods (XFEM) [36–41], where a local enrichment with new basis functions is considered to handle the nonconformity, or the class of Embedded Discrete Fracture Matrix Methods (EDFM) [5, 41–43], where approximate formulas for fracturematrix transmissibilities are computed based on geometrical considerations.
In this paper we consider a conforming discretization with the dual virtual element approximation to simulate a mixeddimensional Darcy problem and a finite volume discretization for the mixeddimensional transport problem. The first choice is motivated by the flexibility of the VEMs with respect to the shape of the grid cells. With this choice it is possible to relax most of the difficulties related with conforming discretization, e.g. the computational cost associated to resolve, by the rock matrix grid, a complex system composed of several intersecting fractures. Moreover, this approximation ensures local mass conservation and is able to consider heterogeneous and anisotropic permeability tensors. The Darcy velocity can thus be used in the mixeddimensional transport problem. In this case we consider an upwind scheme, extended to handle the mixeddimensional nature of the problem.
The paper is structured as follows: In Section 2, we present both the physical equations and the reduced model, with the interface conditions that couple the matrixfracture system and the fracturefracture system for both the Darcy and transport problems. Section 3 deals with the weak and integral formulation of the previous physical processes. In Section 4, we present the numerical discretization of the problem with a highlight on the enrichment of the finite element spaces. In Section 5, we present some numerical experiments to assess the effectiveness of the proposed method. Finally, Section 6 is devoted to conclusion and to ongoing work.
2 Mathematical model
In this section, we introduce the mathematical models in the mixeddimensional setting, i.e. equations that couple different spatial dimensions and together describe the quantities of interest. We present two mathematical models useful for subsurface simulations: the mixeddimensional Darcy problem, presented in Section 2.1, to describe the flow and pressure, and the mixeddimensional transport problem, presented in Section 2.2, to describe the motion of a passive scalar transported by the Darcy velocity.
The motivation for the mixeddimensional formulation is the observation that the fracture aperture is orders of magnitude smaller than other characteristic sizes of the problem. Hence, a straightforward mesh construction with the two fracture surfaces as constraints would produce high numbers of cells and/or lowquality cells due to high aspect ratios or sharp angles. To avoid this geometrical constraint, fractures are represented as lower dimensional objects embedded in the rock matrix. Fracture intersections, and their intersections again, are considered as objects of even lower dimensions. To be specific, in a threedimensional domain, we consider fracture surfaces as 2d objects, the intersection between two (or more) fractures form a 1d line, and two (or more) intersection lines can meet in a 0d point. The physical processes are described via reduced models with suitable coupling conditions among the objects of different dimensions. The fracture aperture is now part of the equations and not any more a geometrical constraint. For more details on the derivation of the reduced model we refer to [8–11, 13, 26–28, 37, 39, 40, 44, 45].
In the reduced model, it is a common choice to consider the reduced scalar variables as averaged and the vector variables as integrated along cross section of the fracture. In this work we follow this approach.
In the sequel we indicate by (·,·)_{ A } the scalar product in L ^{2} (A). Moreover, the trace operator on a domain A will be indicated by ·_{ A }.
2.1 Mixeddimensional Darcy problem
Let us consider a regular domain , for N > 0, with outer boundary ∂Ω. The domain Ω is composed by a single, possibly nonconnected, equidimensional domain Ω^{ N } and several lowerdimensional domains Ω^{ d } for d < N, which are possibly nonconnected. Clearly ∪_{ d=0,…,N } Ω^{ d } = Ω and Ω^{ d } ∩ Ω^{ d }′ = ∅ for d ≠ d′. However, we indicate with Γ^{ d } = ∂Ω^{ d } ∩ Ω^{ d−1}, which is geometrically equivalent to Ω^{ d−1} but it will be more convenient to keep them separate. We have Γ = ∪_{ d=0,…,N } Γ^{ d }. See Figure 1 as an example of the subdivision.
Fig. 1 Domain subdivision for N = 2. In gray Ω^{2}, in red Ω^{1}, and in blue Ω^{0}. 
We are interested in the mixeddimensional Darcy problem which describes the pressure and velocity fields in Ω. Following the idea presented in [10], we define the pressure compound p = (p ^{0}, …, p ^{ N }) which represents the pressure in each dimension, indicated as superscript. Similarly, we introduce the velocity compound u = (0, u ^{1}, …, u ^{ N }, 0), which is the velocity field defined in Ω^{ N } or in the tangential spaces of each dimension d < N. The meaning of the two 0 s will be clarified below. In the following problem, each dimension is coupled with the lower dimension; we consider thus two normals n associated to a d < Ndimensional manifold with respect to its “two sides’’, indicated by n _{+} and n _{−}. Each disconnected object with dimension d < N when immersed in the surrounding d + 1manifold, it is assumed to be an orientable ddimensional manifold, with unique normal n lying on the tangent space of the d + 1manifold. This implies n _{+} = − n _{−}, an illustration is given in Figure 2. In the particular case of Ttype intersection only one side can be defined and the following equations should be understood in this sense. Given a side, by convention the normal vector points from the higher dimensional to the lowerdimensional domain. Moreover, we indicate by the same subscript ·_{±} the restriction of elements in Ω^{ d } defined on each side of Γ^{ d } through a suitable trace operator.
Fig. 2 Representation of the fracture normal associated with respect to each “side” (+ and −). (a) A bidimensional domain Ω^{2}, (b) The fracture is fictitiously enlarged. (c) For a 0d intersection. 
2.1.1 Fixeddimensional Darcy problem
We consider the Darcy problem for the pressure and velocity as
The notation means “on the higher dimensional object” and the means “on the lower dimensional object”, i.e. fix a dimension d the second of (1a) becomes
In the system (1a) the first equation and second equations are defined on domains of the same dimension. For lower dimensional objects, the differential operators are defined in the tangent plane of the objects; in the special 0d case these objects are void. K represents the effective permeability compound K = (1, K ^{1}, …, K ^{ N }), which is a collection of symmetric and positive defined tensors. For d < N, K ^{ d } is interpreted as laying in the tangent plane of the relevant physical object. Model (1a) allows for permeability variations both between all objects considered, and spatially within each object. The permeabilities are scaled by the “aperture” of each dimension: Denote by ϵ ^{ d } the length (N − 1), area (N − 2), or crosssectional volume (N − 3). The permeability of a lower ddimensional (d < N) object is then defined as , with the permeability of its equidimensional representation. The unit of measure of K ^{ d } is thus [m^{4−d }/s]. The source or sink term f is defined as the compound f = (f ^{0}, …, f ^{ N }), which are scalar functions defined in each dimension. Finally, the jump operator is imposed on each side of Ω^{ d−1} viewed as immersed in Ω^{ d }.
The first equation in (1a) is the generalized Darcy law in each dimension, note that for the 0d case it is simply an identity. The second equation of (1a) represents conservation of mass: The divergence term prescribes conservation of mass within the dimension, while the jump term represents the inflow/outflow from the higher dimension. Interaction with the lower dimension is incorporated into the boundary conditions, see Section 2.1.3. Due to the definition of the velocity compound, for d = N equation (1a) the inflow/outflow from the “higher dimension” is null. For d = 0 the conservation of mass equation represents the continuity of the fluxes for the 1dimensional objects involved in the intersection.
We emphasize that the mixeddimensional problem only considers direct interaction between objects one dimension apart. That is, quantities defined on Ω^{ d−1} are in communication with the quantities defined on Ω^{ d }, but not directly with quantities defined on Ω^{ d+1}.
2.1.2 Flow between dimensions
The coupling conditions between two subsequent dimensions are defined for each side as
The is the effective normal permeability, defined as . Clearly, equation (1b) is valid between dimension d and d + 1 for 0 ≤ d < N. Note that in equation (1b) it is possible to define different effective normal permeability on each side ±, however, to simplify the notation, we assume only one single value.
2.1.3 Boundary conditions
To obtain a wellposed problem, boundary conditions must be assigned on the external boundary ∂Ω and the boundary of the subdomains of each dimension, Ω^{ d }. We indicate by the (possibly empty) portion of the boundary of Ω^{ d } that intersects the boundary of Ω, for d > 0. We introduce also the portion of boundary of Ω^{ d } that does not intersect with ∂Ω, for d > 0. is further split into the internal boundaries Γ^{ d }, which are formed by the embedding of objects in Ω^{ d−1} into Ω^{ d }, and the remaining part, .
For simplicity we assume only pressure boundary conditions on , given as
(1c)with is a given pressure at the boundary for each d > 0. On the internal boundary Γ^{ d }, we impose flux conditions that represent the inflow and outflow from lower dimensions, in accordance with (1b). In practice, this couples the flux between dimensions with the pressure in both dimensions, in what can be interpreted as a Robintype condition.
On the remaining internal portion , for d > 0, we assign the socalled “tipcondition”, namely
(1d)where in the previous equation n stands for the outward unit normal of .
2.1.4 Generalized formulation
Equation (1) is the mixeddimensional Darcy problem, formulated in terms of pressure and velocity. The problem (1) may be recast into a pure pressure formulation, however the numerical scheme introduced in Section 4 considers explicitly the pressure and velocity fields as unknowns.
Following the idea presented in [10, 34, 46], we can write problem (1) in a more compact form similar to the standard Darcy formulation. To this end, we introduce the following divergence operator · between dimensions as such that
as well as a mixeddimensional gradient operator such that
Considering the Darcy velocity composed by u in Ω and u n on Γ, system (1) becomes
In the previous equation is defined accordingly to include both the effective tangential and normal permeabilities.
2.2 Mixeddimensional transport problem
Once the problem (1) is solved, the velocity field u can be used to describe the mixeddimensional transport problem. We consider the same splitting of Ω as in Section 2.1. We introduce the concentration compound as c = (0, c ^{1}, …, c ^{ N }, 0), with 0 ≤ c ≤ 1 for all d, to describe the portion of a passive tracer in a grid cell. We indicate by t the time variable and (0, T) the time interval. The proposed model extends to multiple fracture networks the model presented in [47] with multidimensional treatment of fracture intersections. Other important references on this subject are [20, 22].
2.2.1 Fixeddimensional transport problem
The model is the following, given u find c such that
(2a)where we have indicated with ∂ _{ t } c the time derivative of c, while ϕ and r are the porosity and a scalar source/sink term in each dimension, represented as compounds. Finally, ϵ represents the “fracture aperture’’ in each dimension, we have . We assume that ϵ ^{ d } > 0 for all 0 ≤ d ≤ N. For simplicity we assume that ϵ does not depend on time.
In (2a), the divergence term models the conservation of c in the current dimension d and the flow exchange with the lower dimension d – 1. The jump operator describes the flow exchange with the higher dimension d + 1. The coupling conditions, related to the definition of inflow/outflow of the lowerdimensional objects, associated to (2a) are
(2b)Also in this case, problem (2) couples the concentration defined on Ω^{ d } with Ω^{ d−1} and Ω^{ d+1}, but not directly with Ω^{ d−2} or Ω^{ d+2}.
2.2.2 Initial and boundary conditions
With (2a) we need to assign initial condition for the concentration c, as
(2c)with t _{ini} the initial time. Finally, on the inflow part of the outer boundary we assign a boundary condition for (2a)
(2d)where is the portion of such that an inflow occurs. Equation (2) is the mixeddimensional transport problem.
2.2.3 Generalized formulation
The problem (2) can be condensed, using the formalism introduce at the end of the previous subsection, as
(2abis)which is a general form of a conservative equation in the mixeddimensional setting.
3 Weak and integral formulation
In the view of introducing the numerical approximations, in this part we consider the weak formulation of (1) and the integral formulation of (2).
3.1 Weak formulation of mixeddimensional Darcy problem
We introduce the functional spaces suitable to approximate the pressure and the velocity. The domains Ω^{ d } may contain cuts or internal interfaces, in this case refer to [17, 34, 48] for a more detailed description of functional setting. Considering a fixeddimensional domain Ω^{ d } we have for the pressure , for d > 0, and . For the velocity we consider, for d > 0, the Hilbert space:
The additional condition in V ^{ d } is related to the Robintype nature of the condition (1b), see [13, 28]. The global spaces Q and V are defined as the union of the local spaces.
We introduce now the bilinear forms associated to the problem (1). We have and , defined as
with w , v ∈ V ^{ d } and q ∈ Q ^{ d }. For the coupling between dimensions, d > 0, we introduce and defined as
where and . In the previous definition we have the inverse of and n the normal of Ω^{ d−1} in the tangent space of Ω^{ d }. The global bilinear forms are defined by and by as
with w , v ∈ V and q ∈ Q. We can introduce the weak formulation of (1), which reads: find such that
(3)where the functionals are defined as
3.2 Integral formulation of mixeddimensional transport problem
In Section 4.2 we consider a finite volume approximation of problem (2), thus we introduce here its integral formulation.
Considering a suitable subdomain of Ω^{ d } called E, later will be a cell of the computational grid, we integrate equation (2) obtaining the following
(4)where I = E ∩ Ω^{ d−1} is the intersection between E and the lower dimensional domain Ω^{ d−1}, which can be an empty set.
In (4), we notice the term depending on I which describe the inflow/outflow exchange between the current and lower dimensions. Condition (2b) applies also in this case.
4 Numerical discretization
In this section we introduce the numerical scheme used to solve the mixeddimensional problems (1) and (2). In particular, for the former we consider the lowest order dual VEM approximation, while for the latter we employ a finite volume discretization with an upwind scheme and implicit Euler in time.
Let us introduce a suitable tessellation of Ω in polytopes indicated by E, such that
We denote by h _{ E } the diameter of E, by x _{ E } the center of E, by the set of faces for the cell E, and by n _{ e } with the unit normal of e pointing outward with respect to the internal part of E. The term faces indicates proper faces for the 3d grid, edges for the 2d grids, and vertexes for the 1d grids. We indicate also by the characteristic grid size and by the set of grid faces. It is important to note that the polytopes may have different spatial dimension, i.e. they are 3d or 2d or 1d or 0d objects. However, fixing a single dimension d the cells of the grid belongs to Ω^{ d }.
For simplicity, both problems (1) and (2) use the same computational grid.
In Section 4.1 we present the discretization for problem (1), while in Section 4.2 the numerical approximation for problem (2). Finally, in Section 4.3 we introduce a coarsening strategy used in the numerical examples.
4.1 Mixeddimensional dual VEM discretization
We present a numerical scheme to solve the mixeddimensional Darcy problem (1) in presence of polytopes in the grid. The resulting scheme will be locally and globally conservative, thus suitable to approximate the velocity field used in (2). We consider the dual VEM of lowest order degree, for more details on the derivation and on the analysis of the method, see [26–33].
4.1.1 Discrete function spaces
We introduce finite dimensional spaces to approximate the Darcy velocity u and the pressure p in each element of the grid. For simplicity we consider a single dimension and, if not essential, we drop the superscript d to simplify notation. Again, the differential operators are understood to be defined on the tangential space. Given an element E the local discretization space for the pressure is , where is the space of polynomial of order r on the domain E. While higherorder dual VEM approximations are possible, see [31, 32], the solution will often have low regularity caused by fractures or general permeability heterogeneities. We therefore stay with the lowest order formulation, as is common practice in the field. For the velocity we need to introduce the following space:
The shape of the functions in is not defined apriori and is implicitly defined by . The curlfree condition is necessary to uniquely define the elements in . We indicate by the velocity approximation space in the same dimensional grid and by V _{ h } the global discretization space for the velocity formed by the compound . Similarly, indicates the pressure approximation space in the same dimension d while Q _{ h } is the global discretization space. For the velocity, we impose to the faces which are not in contact with different dimensions to be single value. Otherwise, the degree of freedom is doubled and connected thorough the coupling condition (1b) to the lower dimensional object. See Figure 3 for an example.
Fig. 3 Representation of the degrees of freedom for a 2d and 1d grid. The pressure dof are represented by circles, red for the 2d grid and green for the 1d. The velocity dof are depicted by yellow diamonds for the 2d grid and purple diamonds for the 1d. The nodes of the 2d grid are moved only for visualization purpose. 
4.1.2 Approximation of bilinear forms
With the previous definition of Q _{ h } and V _{ h } it is immediate to approximate the bilinear form b for all the dimensions as well as the bilinear forms associated with the coupling conditions, α ^{ d } and β ^{ d,d−1}. The functionals F and G are discretized similarly. The bilinear forms a ^{ d } are not immediately computable with the degrees of freedom introduced, but require the definition of a suitable projection operator. First, we introduce the local space
(5)and the projection operator is defined as such that given we have for all . Note that the space in the definition of will be approximated by a (tangential) monomial basis. Where we are assuming a constant value of the permeability in the element E. For (5) to make sense, the permeability must be constant in each cell. As discussed above, the permeability will in practice often be spatially varying on all scales due to the presence of smaller scale heterogeneities. Inclusion of spatially varying permeabilities within a computational cell would require corresponding modifications of the approximation order for pressure and flux. Our preferred approach to including more permeability variations is to stay with lowest order approximations, but instead refine the grid.
With the orthogonality property it is possible to split the bilinear form a ^{ d } in a term on and one on the a ^{ d }orthogonal space of , namely
with T _{0} = I−Π_{0}. The first bilinear form is now fully computable with the velocity degrees of freedom introduced before and represents a consistency term with respect to a ^{ d } ( u , v ). The second term is not yet computable and can be replaced by a stabilization term. Following the ideas presented in [31–33], we approximate the term by,
where is the bilinear form associated with the stabilization and is a scaling parameter described in the sequel. In details, denoting by φ an element of the basis for , in our case we consider
with N _{dof} the total number of velocity degrees of freedom of the element E and is defined as
We introduce the discrete version of the bilinear form a ^{ d }, defined as
and the discrete version of the weak problem (3), which reads: find such that
(6)where . The construction of the discrete problem in term of local matrix computations is described in [28, 30].
The stabilization parameter ς is used to impose the scaling on h of the stabilization term equivalent to the consistency term. In practice we require that, for a fixed dimension d, there exists , independent from the discretization, such that
Following [28], by a scaling computation we evaluate the dependency of ς from the local grid size h _{ E } and obtain the following relation for a fix dimension d
Note that this expression is local and independent from the maximal dimension N of the problem.
4.1.3 Matrix formulation
We introduce the matrix formulation associated with the problem (6). Considering the following matrices:
where ϕ are basis for the pressure and indicates the element (i, j) in the matrix, a similar notation is used for vectors. The global problem reads for N = 3, solve the following linear system
where u _{ d } and p _{ d } represent the vectors associated with the degrees of freedom of velocity and pressure in each dimension d. We note that the u _{0} is considered only for clearness and to preserve the structure of the matrix. In practice, it is possible to remove u _{0} from the system.
4.2 Mixeddimensional FV discretization
We consider now the discretization of equation (4). For simplicity, we consider a fixed time step Δt such that T/Δt is an integer number. An implicit Euler is applied in time obtaining the semidiscrete version of (4)
where the superscript indicates the time step index, and . The approximation of the boundary integrals in the previous equation rely on an upwind scheme. We introduce a Kroneckertype delta as
Given a cell face we have
where the cells K and L share the face e. Given a cell E such that one of its faces, e, intersects one side of I, we get
where indicates a cell in the codimensional grid which is in communication with the face e. Finally, on side giving an element E of the last term of the semidiscrete problem can be approximated by
where represents the cell in the higher dimensional grid which has a face e in communication with the cell E. The previous condition applies to both the equi and codimensional coupling, see Figure 4 for an example. Note that the chosen discretization is compatible with the coupling condition (2b).
Fig. 4 Representation of the coupling between dimension for the upwind discretization scheme. (a) Between two cells in the same dimension d. (b) Between a ddimensional cell and a d − 1dimensional cell. (c) Between a d + 1dimensional cell and a ddimensional cell. 
For the matrix formulation of the transport problem, we introduce the following matrices
where is the set of all neighbor cells of E, i and j indicate generic cells and e is a face shared by i and j. For the coupling between dimensions we have
where is the set of all neighbor cells of E in the lowerdimensional grid. We obtain also
where is the set of all neighbor cells of E in the higherdimensional grid. We finally obtain the following linear system to be inverted
where c _{ d } represents the vector associated with the degrees of freedoms of the concentration in each dimension d. In the previous linear system represents the righthand side, formed as a combination of the source term and the concentration at the previous time step. We get
4.3 Coarsening strategy
The construction of conforming grids in the presence of several fractures can be a challenging task, especially in 3d. In this work, our grids are constructed by the Gmsh mesh generator [49]. In the presence of almost intersecting fractures or small fracture branches, the grid may contain a high number of simplex cells, resulting in a high computational cost. To overcome this difficulty, we exploit the ability of VEM to handle cells of arbitrary geometry. To be specific, the theory developed in [30–33] requires starshaped cells, however the study carry out in [50] shows that the VEM is able to handle also cells with cuts.
Motivated by these observations we introduce a coarsening scheme that merges small simplex cells into larger polygons or polyhedra. Starting from a given simplex grid, the algorithm computes the measure (area or volume) of the cells. Given the cell c with the smallest measure, it will be merged to one or more neighboring cells, based on their respective measure, creating a new coarse cell. The clustering stops when a cell measure threshold is reached. The cells used in discretization are thus general polytopes, formed as unions of simplexes.
The algorithm does not guarantee any regularity of the final grid, as illustrated in Figure 5.
Fig. 5 Example of the coarsening strategy adopted. (a) The computational grid is artificially forced to be finer at the tip of a fracture. (b) The resulting grid after the coarsening. Clearly the cell measures are comparable. 
We emphasize that as elements are aggregated, a new element with many faces is generated. The local matrices will be denser and larger comparing to the original grid and increase the stencil of the global matrix. Nevertheless, the coarsening algorithm will decrease the computational cost. Other strategies are possible but not considered in this work.
Finally, we note that the assumption of a cellwise constant permeability applies also in this case. Thus, merging of cells with different permeabilities requires the computation of an average permeability. Correspondingly, for highly heterogeneous media, it may be beneficial to consider also cell permeability as a parameter in the coarsening algorithm.
5 Applicative examples
We present three examples and test cases to assess the presented models and numerical schemes. The first example, presented in Section 5.1, considers an extensive validation of the mixeddimensional Darcy problem solved by VEM through a benchmark study presented in [51]. The second test case, in Section 5.2, considers the transport problem and analyzes the impact of the coarsening strategy on the results. Finally, in Section 5.3 a realistic 3d example is introduced and studied. In all the forthcoming examples we assume unitary porosity and zero source terms for the concentration and Darcy equations. The other parameters will be specified.
In all of the forthcoming examples a “Blue to Red Rainbow’’ color map is used.
The examples are part of the PorePy package, which is a simulation tool for fractured and deformable porous media written in Python. See [52] and github.com/pmgbergen/porepy for more details.
5.1 Benchmark comparison
To validate the presented model, we consider the benchmark study presented in [51], specifically cases 1 and 4. The reference solutions p ^{ref} were computed with a mimetic finite difference method [53] on a very fine grid that is able to represent the thickness of the fractures by considering the classical Darcy model. The error is thus evaluated by
for the matrix err ^{2} and for the fractures err ^{1}, where and E ^{2} is a cell in the mixeddimensional problem related to the rock matrix, respectively E ^{1} for the fracture. The errors are computed on the reference grid, p ^{1} is assumed constant on the normal direction of the fracture in the equidimensional setting.
For more detail of problem setting refer to the aforementioned work.
5.1.1 Benchmark 1: regular fracture network
The problem is inspired by Geiger et al. [54] with different boundary conditions and material properties. The domain is the unit square containing six fractures, sketched in Figure 6. The matrix permeability is taken as unitary, and the fracture aperture is set equal to 10^{−4}. No flux boundary conditions are imposed on the top and bottom of the domain, while unitary pressure is imposed on the right boundary and −1 as flux on the left boundary. The source term is set to zero. We consider two possibilities for fracture permeability: it is either highly conductive with permeability 10^{4} and has conductivity with permeability 10^{−4}. In the former case the solution obtained with the method presented previously along with the computational grid are presented in Figure 7(a), the latter in (b). We observe a good agreement between the computed and reference solutions, the latter is described in [51]. We point out that some of the elements present in Figure 7 are not convex.
Fig. 6 Domain with fractures (in red) and fracture intersections (in blue) for problem in Section 5.1.1. 
Fig. 7 Benchmark 1. (a) Pressure solution (range (1, 1.6)) with conductive fractures and computational grid. (b) Pressure solution (range (1, 3.6)) with blocking fractures and computational grid. (c) Zoom of the grid used. 
To obtain a more detailed comparison we consider two plots over line for the permeable case and one for the blocking case, shown in Figure 8. From now on, the method presented in this paper is labeled as VEM. Also in this case, the agreement between the computed and reference solutions is comparable to the other methods which are able to represent blocking fractures. The small oscillations in the line plots are related to grid effects.
Fig. 8 Benchmark 1 with conductive fractures. (a) Pressure along horizontal line at y = 0.7 with permeable fracture. (b) Pressure along vertical fracture at x = 0.5 with permeable fractures. (c) Pressure along the line (0.0, 0.1)–(0.9, 1.0) with blocking fractures. 
Finally, in both cases we consider the error decay for both the rock matrix and the system of fractures. We consider a family of three grids where the coarsening is applied in all the cases. We obtain again nonconvex elements in all the grids. Figure 9 plots for conductive and blocking fractures the error decay. The errors are comparable with those of others methods able to represent permeable and blocking fractures. In the latter case we note a stagnation of the fracture error which bounds the order of converge. This phenomenon is common also for other methods presented in the benchmark study [51].
Fig. 9 Benchmark 1, error evolution for the matrix and the fractures with conductive and blocking fractures. 
5.1.2 Benchmark 4: a realistic case
We consider now a complex system of 64 fractures from a real outcrop. The domain is Ω = (0, 700 m) × (0, 600 m), sketched in Figure 10. We consider constant rock permeability equal to 10^{−14} m^{2}, uniform fracture permeability 10^{−8} m^{2}, and fracture aperture 10^{−2} m^{2}. We impose a pressure gradient at the boundary from the left (1 013 250 Pa) to the right (0 Pa), and no flow condition on the top and bottom of the domain. The source term is set to zero. Also in this case, a triangular grid is coarsened to reduce the grid complexity. The reference grid is composed by 12 472 2d cells, 1317 1d cells, and 85 0d cells. The coarse algorithm decreases the number of 2d cells down to 4703. The total number of degrees of freedom is 19 075 for the coarse grid.
Fig. 10 Domain with fractures (in red) for problem in Section 5.1.2. 
The computed solution is represented in Figure 12(a), which matches the solution of the others method considered in [51]. Because of the complexity of the network, conforming and nonmatching methods may pose a constraint to the grid generation in particular for close fractures. However, the method presented allows general grid cells, and thus alleviate the computational cost for the simulation. For an illustration of the computational grid, we refer to Figure 11, which zooms in on a region with almost intersecting fractures.
Fig. 11 Benchmark 4. (a) Original grid composed by 12 302 2dcells and 35 153 VEM dof and coarsened grid composed by 4599 2dcells and 18 803 VEM dof. (b) A zoom on almost intersecting fractures for the original and coarsened grids, respectively. The zoom is referred to the small rectangle at position, approximately (360, 350). 
As shown in Figure 11, some of the grid cells are non starshaped or even contains cuts. For a more detailed discussion refer to [50]. Finally, to validate in more detail the computed solution we present two plots over line and compare them with the solutions obtained in the benchmark study, see Figure 12. The curves for the current method are in good agreement with the others. Again, the small oscillations that can be observed are related to grids effects.
Fig. 12 Benchmark 4. (a) Pressure solution computed with VEM. The pressure ranges in [0, 1 013 250] Pa. (b) Pressure along horizontal line at y = 500 m. (c) Pressure along vertical line at x = 625 m. 
5.2 Passive scalar transport
In this part we consider both mixeddimensional models (1) and (2) to simulate a passive scalar transport. We consider the geometry presented in Section 5.1.2 and compare the solution obtained with both grids in Figure 11 for permeable and blocking fractures. For the reference grid, the total number of degrees of freedom is 35 578. The aim of this test is to validate the quality of the Darcy velocity on the passive scalar in presence of the grid coarsening. In the following fracture aperture is constant and equal to 10^{−2} m, a pressure gradient from the right to left boundary of the domain of 3 × 10^{7} Pa, a final simulation time of 40 years.
5.2.1 Permeable fractures
We consider highly permeable fractures with tangential permeability of 5 × 10^{−6} m^{2} and normal permeability of 2.5 × 10^{−9} m^{2}. The matrix permeability is set to 2.5 × 10^{−11} m^{2}. Figure 13 compares the solutions obtained on the reference triangular grid and on the coarse grid. Moreover Figure 14(a) presents a comparison of the passive scalar production. We notice the good agreement in both the pressure and concentration fields as well as in the production curve. We can conclude that in this case the grid coarsening is not affecting the quality of the computed solutions.
Fig. 13 Permeable fractures. (a) Reference and coarse solutions for pressure and velocity, as arrows. (b) Reference and coarse solutions for concentration of the passive scalar. 
Fig. 14 Permeable fractures. (a) Comparison of passive scalar production at the outflow between the reference triangular grid and the coarse grid. (b) Temporal error decay with reference in black. 
Finally, in Figure 14(b) the temporal error decay is reported. The spatial discretization is fixed and we consider a sequence of simulation with (10, 20, 40, 80, 160, 320, 640, 1280, 2560, 5120) time steps each. The error is computed as the L ^{2} difference from a reference solution obtained with 10^{5} time steps after 10 years. A unitary error decay is achieved, coherent with the numerical scheme considered.
5.2.2 Blocking fractures
Next, we consider blocking fractures with tangential and normal permeability of 7.5 × 10^{−16} m^{2}. Matrix permeability is set to 7.5 × 10^{−11} m^{2}. Figure 15 compares the solutions obtained on the reference triangular grid and on the coarse grid. In this case we notice a pickpeak of velocity in the reference case due to small elements close to almost intersecting fractures. This may affect an explicit in time solver. Figure 16(a) presents a comparison of the passive scalar production. We notice the good agreement in both the pressure and concentration fields as well as in the production. We can conclude that also in this case the grid coarsening is not affecting the quality of the computed solutions.
Fig. 15 Blocking fractures. (a) Reference and coarse solutions for pressure and velocity, as arrows. (b) Reference and coarse solutions for concentration of the passive scalar. 
Fig. 16 Blocking fractures. (a) Comparison of passive scalar production at the outflow between the reference triangular grid and the coarse grid. (b) Temporal error decay with reference in black. 
Again, the temporal error decay is reported, see Figure 16(b). The spatial discretization is fixed and we assign (10, 20, 40, 80, 160, 320, 640, 1280, 2560, 5120) time steps, respectively. A unitary error decay is achieved, consistent with the numerical scheme considered.
5.3 Passive scalar transport on a realistic 3dnetwork
In this example we consider a realistic geometry for a geothermal system. We study a partial reconstruction of fractures from test site at SoultzsousForêts in France, for more details see [55]. The network is composed by 20 fractures represented as polygons with 10 edges each. The fracture intersections result in 33 1d objects and 4 0d objects. In this case the full model is needed to accurately simulate the fluid flow and transport on the domain. The fracture geometry is represented in Figure 17. We assume rock matrix permeability as 7.5 × 10^{−10} m^{2} and fracture permeability, in both the normal and tangential direction, equal to 5 × 10^{−5} m^{2}. The fracture aperture is ϵ ^{2} = 10^{−2} m and, for the lower dimensional objects, we consider their “aperture’’ as the square and the cube of the fracture aperture, respectively for the 1d and 0d objects. We impose a pressure from the top to bottom of the domain and no flux boundary conditions on the other sides. The reference grid is composed by 44 331 tetrahedra, 6197 triangles for the 2d grids, 151 segments for the 1d grids, and 4 pointcells for the 0d grids. After the coarsening algorithm the resulting grid is composed by 16 108 polyhedra, and for the lower dimensional objects the grids are untouched. See Figure 19(b) for an example of coarse cells. The transport simulation runs for 40 years and, for the implicit Euler scheme, we consider 100 time steps.
Fig. 17 On the left representation of the 20 fractures colored by their identification number. Pressure and velocity for both the reference and coarse grids. The pressure is scaled between 0 and 4.8 × 10^{7} Pa. 
The objective is to study the robustness of the virtual elements associated with the coarsening strategy and detect if the coarse model gives accurately enough results. The pressure solution and the concentration of the scalar passive are depicted in Figures 17 and 18, respectively, for both the reference and coarse case. Both the pressure and concentration profiles for the two grids are in good agreement. To analyze the macroscopic behavior of the resulting solution a production curve comparison is reported in Figure 19, showing a small discrepancy from the reference and the coarse production.
Fig. 18 (a) Three time steps on the reference grid for the concentration. (b) At the same steps, the concentration computed on the coarse grid. 
Fig. 19 (a) The production comparison between the reference and the coarse solution. (b) An example of five coarse cells. 
We can conclude that also in this case the adopted strategy is effective and can be applied to complex system of fractures, ensuring a good compromise between high accuracy and computational effort.
6 Conclusion
In the paper we presented two classes of mixeddimensional problems able to describe a singlephase flow and a transport of a scalar passive in fractured porous media. The latter is transported by the Darcy velocity computed by the former model. The considered mixeddimensional Darcy problem is able to represent both channel and barrier behaviors of fractures, as well as their intersections. The models thus represent a comprehensive description of this phenomena, which are key ingredients for the description of several energy application problems. The numerical schemes considered for the discretization of the Darcy equation are able to handle, supported by theoretical findings, grid cells of arbitrary geometry becoming a strong advantage when dealing with high and complex fractured porous media. Numerical results have shown the capability to apply this strategy obtaining accurate outcomes with a reasonable computational cost. The transport model can be viewed as a first attempt to introduce a full model for the description of heat exchange in a porous media, which will be a part of future investigations.
Acknowledgments
We acknowledge financial support for the ANIGMA project from the Research Council of Norway (project no. 244129/E20) through the ENERGIX program. The authors warmly thank the others components of the ANIGMA project team: Eivind Bastesen, Inga Berre, Simon John Buckley, Casey Nixon, David Peacock, Atle Rotevatn, Pål Næverlid Sævik, Luisa F. Zuluaga.
The authors wish to thank also: Runar Berge, Wietse Boon, Ivar Stefansson, Eren Uçar.
References
 Berre I., Boon W., Flemisch B., Fumagalli A., Gläser D., Keilegavlen E., Scotti A., Stefansson I., Tatomir A. (2018) Call for participation: Verification benchmarks for singlephase flow in threedimensional fractured porous media. Technical report, arXiv:1710.00556 [math.AP]. [Google Scholar]
 Mourzenko V.V., Thovert J.F., Adler P.M. (Sep 2011) Permeability of isotropic and anisotropic fracture networks, from the percolation threshold to very large densities, Phys. Rev. E 84, 036307. https://doi.org/10.1103/PhysRevE. URL https://link.aps.org/doi/10.1103/PhysRevE.84.036307. [Google Scholar]
 Sævik P.N., Berre I., Jakobsen M., Lien M. (Oct 2013) A 3d computational study of effective medium methods applied to fractured media, Transp. Porous Media, 100, 115–142. ISSN 15731634. https://doi.org/10.1007/s1124201302080. [Google Scholar]
 Ssvik P.N., Jakobsen M., Lien M., Berre I. (2014) Anisotropic effective conductivity in fractured rocks by explicit effective medium methods, Geophys. Prospect. 62, 6, 1297–1314. ISSN 13652478. https://doi.org/10.1111/13652478.12173. [CrossRef] [Google Scholar]
 Fumagalli A., Pasquale L., Zonca S., Micheletti S. (2016) An upscaling procedure for fractured reservoirs with embedded grids, Water Resour. Res. 52, 8, 6506–6525. ISSN 19447973. https://doi.org/10.1002/2015WR017729. [Google Scholar]
 KarimiFard M., Durlofsky L.J. (Oct 2016) A general gridding, discretization, and coarsening methodology for modeling flow in porous formations with discrete geological features, Adv. Water Res. 96, 354–372. ISSN 03091708. https://doi.org/10.1016/j.advwatres.2016.07.019. URL http://linkinghub.elsevier.com/retrieve/pii/S0309170816302950. [CrossRef] [Google Scholar]
 KarimiFard M., Gong B., Durlofsky L.J. (2006) Generation of coarsescale continuum flow models from detailed fracture characterizations, Water Resour. Res. 42, 10. ISSN 19447973. https://doi.org/10.1029/2006WR005015. [Google Scholar]
 Alboin C., Jaffré J., Roberts J.E., Wang X., Serres C., Chen Z., Ewing R. E., Shi Z.C. (2000) Domain decomposition for some transmission problems in flow in porous media. Numerical treatment of multiphase flows in porous media (Beijing, 1999), Springer, pp. 22–34. [CrossRef] [Google Scholar]
 Amir L., Kern M., Martin V., Roberts J.E. (2005) Décomposition de domaine et préconditionnement pour un modèle 3D en milieu poreux fracturé, in: Proceeding of JANO 8, 8^{th} Conference on Numerical Analysis and Optimization, December 2005. [Google Scholar]
 Boon W.M., Nordbotten J.M., Yotov I. (2018) Robust discretization of flow in fractured porous media, SIAM J. Numer. Anal. 56, 4, 2203–2233. https://doi.org/10.1137/17M1139102. [Google Scholar]
 Faille I., Fumagalli A., Jaffré J., Roberts J.E. (2016) Model reduction and discretization using hybrid finite volumes of flow in porous media containing faults, Comput. Geosci. 20, 2, 317–339. ISSN 15731499. https://doi.org/10.1007/s1059601695583. [Google Scholar]
 Knabner P., Roberts J.E. (2014) Mathematical analysis of a discrete fracture model coupling Darcy flow in the matrix with DarcyForchheimer flow in the fracture, ESAIM: Math. Model. Numer. Anal. 48, 1451–1472. ISSN 12903841. https://doi.org/10.1051/m2an/2014003. URL http://www.esaimm2an.org/article_S0764583X1400003X. [CrossRef] [Google Scholar]
 Martin V., Jaffré J., Roberts J.E. (2005) Modeling fractures and barriers as interfaces for flow in porous media, SIAM J. Sci. Comput. 26, 5, 1667–1691. ISSN 10648275. https://doi.org/10.1137/S1064827503429363. URL http://scitation.aip.org/getabs/servlet/GetabsServlet?prog= normal&id=SJ0CE3000026000005001667000001&idtype= cvips&gifs=yes. [Google Scholar]
 Schwenck N., Flemisch B., Helmig R., Wohlmuth B.I. (2015) Dimensionally reduced flow models in fractured porous media: crossings and boundaries, Comput. Geosci. 19, 6, 1219–1230. ISSN 14200597. https://doi.org/10.1007/s1059601595361. [Google Scholar]
 Tunc X., Faille I., Gallouet T., Cacas M.C., Havé P. (2012) A model for conductive faults with nonmatching grids, Comput. Geosci. 16, 277–296. ISSN 14200597. https://doi.org/10.1007/s105960119267x. [Google Scholar]
 Angot P. (2003) A model of fracture for elliptic problems with flux and solution jumps, C. R. Math. 337, 6, 425–430. ISSN 1631073X. https://doi.org/10.1016/S1631073X(03)003005. URL http://www.sciencedirect.com/science/article/pii/S1631073X03003005. [CrossRef] [Google Scholar]
 Angot P., Boyer F., Hubert F. (2009) Asymptotic and numerical modelling of flows in fractured porous media, M2AN Math Model. Numer. Anal. 43, 2, 239–275. ISSN 0764583X. https://doi.org/10.1051/m2an/2008052. [CrossRef] [EDP Sciences] [MathSciNet] [Google Scholar]
 Chave F.A., Di Pietro D., Formaggia L. (2017) A Hybrid HighOrder method for Darcy flows in fractured porous media, Technical report, HAL archives, 2017. URL https://hal.archivesouvertes.fr/hal01482925. [Google Scholar]
 KarimiFard M., Firoozabadi A. (2003) Numerical simulation of water injection in fractured media using the discretefracture model and the Galerkin method, SPE Reserv. Eval. Eng. 6, 2, 117–126. [CrossRef] [Google Scholar]
 Ahmed R., Edwards M.G., Lamine S., Huisman B.A.H., Pal M. (2015) Controlvolume distributed multipoint flux approximation coupled with a lowerdimensional fracture model, J. Comput. Phys. 284, 462–489. ISSN 00219991. https://doi.org/10.1016/jjcp.2014.12.047. http://www.sciencedirect.com/science/article/pii/S0021999114008705. [Google Scholar]
 Brenner K., Hennicker J., Masson R., Samier P. (September 2016) Gradient discretization of hybriddimensional Darcy flow in fractured porous media with discontinuous pressures at matrixfracture interfaces, IMA J. Numer. Anal. https://doi.org/10.1093/imanum/drw044. URL https://hal.archivesouvertes.fr/hal01192740. [Google Scholar]
 Brenner K., Hennicker J., Masson R., Samier P. (2018) Hybriddimensional modelling of twophase flow through fractured porous media with enhanced matrix fracture transmission conditions, J. Comput. Phys. 357, 100–124. ISSN 00219991. https://doi.org/10.1016/j.jcp.2017.12.003. http://www.sciencedirect.com/science/article/pii/S0021999117308781. [Google Scholar]
 Brenner K., Groza M., Guichard C., Masson R. (2015) Vertex approximate gradient scheme for hybrid dimensional twophase Darcy flows in fractured porous media, ESAIM: Math. Model. Numer. Anal. 49, 2, 303–330. https://doi.org/10.1051/m2an/2014034. [CrossRef] [Google Scholar]
 Antonietti P.F., Formaggia L., Scotti A., Verani M., Verzotti N. (2016) Mimetic finite difference approximation of flows in fractured porous media, ESAIM: M2AN 50, 3, 809–832. https://doi.org/10.1051/m2an/2015087. [CrossRef] [EDP Sciences] [Google Scholar]
 Scotti A., Formaggia L., Sottocasa F. (2017) Analysis of a mimetic finite difference approximation of flows in fractured porous media, ESAIM: M2AN. https://doi.org/10.1051/m2an/2017028. [Google Scholar]
 Benedetto M.F., Berrone S., Pieraccini S., Scialò S. (2014) The virtual element method for discrete fracture network simulations, Comput. Methods Appl. Mech. Eng. 280, 0, 135–156. ISSN 00457825. https://doi.org/10.1016/j.cma.2014.07.016. [Google Scholar]
 Benedetto M.F., Berrone S., Borio A., Pieraccini S., Scialò S. (2016) A hybrid mortar virtual element method for discrete fracture network simulations, J. Comput. Phys. 306, 148–166. ISSN 00219991. https://doi.org/10.1016/jjcp.2015.11.034. http://www.sciencedirect.com/science/article/pii/S0021999115007743. [Google Scholar]
 Fumagalli A., Keilegavlen E. (2018) Dual virtual element method for discrete fractures networks, SIAM J. Sci. Comput. 40, 1, B228–B258. https://doi.org/10.1137/16M1098231. [Google Scholar]
 da Veiga L.B., Brezzi F., Cangiani A., Manzini G., Marini L.D., Russo A. (2013) Basic principles of virtual element methods, Math. Models Methods Appl. Sci. 23, 1, 199–214. https://doi.org/10.1142/S0218202512500492. [Google Scholar]
 da Veiga L.B., Brezzi F., Marini L.D., Russo A. (2014) The hitchhiker’s guide to the virtual element method, Math. Models Methods Appl. Sci. 24, 8, 1541–1573. https://doi.org/10.1142/S021820251440003X. [Google Scholar]
 da Veiga L.B., Brezzi F., Marini L.D., Russo A. (Jun 2014) H(div) and H(curl)conforming VEM, Numer. Math. 133, 2, 303–332. ISSN 09453245. https://doi.org/10.1007/s0021101507461. [Google Scholar]
 da Veiga L.B., Brezzi F., Marini L.D., Russo A. (2016) Mixed virtual element methods for general second order elliptic problems on polygonal meshes, ESAIM: M2AN 50, 3, 727–747. https://doi.org/10.1051/m2an/2015067. [CrossRef] [EDP Sciences] [Google Scholar]
 Brezzi F., Falk R.S., Marini D.L. (2014) Basic principles of mixed virtual element methods, ESAIM: M2AN 48, 1227–1240. https://doi.org/10.1051/m2an/2013138. [CrossRef] [EDP Sciences] [Google Scholar]
 Frih N., Martin V., Roberts J.E., Saâda A. (2012) Modeling fractures as interfaces with nonmatching grids, Comput. Geosci. 16, 4, 1043–1060. ISSN 14200597. https://doi.org/10.1007/s1059601293026. [Google Scholar]
 Nordbotten J.M., Boon W., Fumagalli A., Keilegavlen E. (2018) Unified approach to discretization of flow in fractured porous media, Comput. Geosci. ISSN 15731499. https://doi.org/10.1007/s1059601897789. [Google Scholar]
 Berrone S., Pieraccini S., Scialò S. (2013) On simulations of discrete fracture network flows with an optimizationbased extended finite element method, SIAM J. Sci. Comput. 35, 2, 908–935. https://doi.org/10.1137/120882883. [Google Scholar]
 D’Angelo C., Scotti A. (2012) A mixed finite element method for Darcy flow in fractured porous media with nonmatching grids, Math. Model. Numer. Anal. 46, 2, 465–489. https://doi.org/10.1051/m2an/2011148. [CrossRef] [EDP Sciences] [Google Scholar]
 Del Pra M., Fumagalli A., Scotti A. (2017) Well posedness of fully coupled fracture/bulk Darcy flow with XFEM, SIAM J. Numer. Anal. 55, 2, 785–811. https://doi.org/10.1137/15M1022574. [Google Scholar]
 Fumagalli A., Scotti A. (2013) A numerical method for twophase flow in fractured porous media with nonmatching grids, Adv. Water Res. 62, Part C(0), 454–464. ISSN 03091708. https://doi.org/10.1016/j.advwatres.2013.04.001. URL https://www.sciencedirect.com/science/article/pii/S0309170813000523. Computational Methods in Geologic CO_{2} Sequestration. [CrossRef] [Google Scholar]
 Fumagalli A., Scotti A. (April 2014) An efficient XFEM approximation of Darcy flows in arbitrarily fractured porous media, Oil Gas Sci. Technol.  Rev. IFP Energies nouvelles 69, 4, 555–564. https://doi.org/10.2516/ogst/2013192. URL https://ogst.ifpenergiesnouvelles.fr/articles/ogst/abs/2014/04/ogst130007/ogst130007.html. [CrossRef] [Google Scholar]
 Tene M., Al Kobaisi M.S., Hajibeygi H. (2016) Multiscale projectionbased embedded discrete fracture modeling approach (famspedfm), ECMOR XIV15th European Conference on the Mathematics of Oil Recovery, August 29–1 September 2016, Beurs van Berlage, EAGE. https://doi.org/10.3997/22144609.201601890. [Google Scholar]
 Fumagalli A., Zonca S., Formaggia L. (July 2017) Advances in computation of local problems for flowbased upscaling in fractured reservoirs, Math. Comput. Simul. 137, 299–324. https://doi.org/10.1016/j.matcom.2017.01.007. URL https://www.sciencedirect.com/science/article/pii/S0378475417300320. [Google Scholar]
 Li L., Lee S.H. (2008) Efficient fieldscale simulation of black oil in a naturally fractured reservoir through discrete fracture networks and homogenized media, SPE Reserv. Eval. Eng. 11, 750–758. https://doi.org/10.2118/103901PA. [CrossRef] [Google Scholar]
 Berrone S., Pieraccini S., Scialò S. (2013) A PDEconstrained optimization formulation for discrete fracture network flows, SIAM J. Sci. Comput., 35, 2, B487–B510. https://doi.org/10.1137/120865884. [Google Scholar]
 Formaggia L., Fumagalli A., Scotti A., Ruffo P. (2014) A reduced model for Darcy’s problem in networks of fractures, ESAIM: Math. Model. Numer. Anal. 48, 1089–1116. https://doi.org/10.1051/m2an/2013132. URL https://www.esaimm2an.org/articles/m2an/abs/2014/04/m2an130132/m2an130132.html. [CrossRef] [EDP Sciences] [Google Scholar]
 Boon W.M., Nordbotten J.M., Vatne J.E. (2017) Mixeddimensional elliptic partial differential equations. Technical report, arXiv:1710.00556 [math.AP]. [Google Scholar]
 Fumagalli A., Scotti A. (2013) A reduced model for flow and transport in fractured porous media with nonmatching grids, in: Cangiani A., Davidchack R.L., Georgoulis E., Gorban A.N., Levesley J., Tretyakov M.V. (eds), Numerical mathematics and advanced applications 2011, Berlin, Heidelberg, Springer, pp. 499–507. ISBN 9783642331336. https://doi.org/10.1007/9783642331343_53. [CrossRef] [Google Scholar]
 Grisvard P. (1985) Elliptic problems in nonsmooth domains, Vol. 24, Monographs and studies in mathematics, Pitman. [Google Scholar]
 Geuzaine C., Remacle J.F. (2009) Gmsh: A 3d finite element mesh generator with builtin pre and postprocessing facilities, Int. J. Numer. Methods Eng. 79, 11, 1309–1331. ISSN 10970207. https://doi.org/10.1002/nme.2579. [Google Scholar]
 Fumagalli A. (Dec. 2018) Dual virtual element method in presence of an inclusion, Appl. Math. Lett. 86, 22–29. https://doi.org/10.1016/j.aml.2018.06.004. URL https://www.sciencedirect.com/science/article/pii/S0893965918301812. [Google Scholar]
 Flemisch B., Berre I., Boon W., Fumagalli A., Schwenck N., Scotti A., Stefansson I., Tatomir A. (January 2018) Benchmarks for singlephase flow in fractured porous media, Adv. Water Res. 111, 239–258. https://doi.org/10.1016/j.advwatres.2017.10.036. URL https://www.sciencedirect.com/science/article/pii/S0309170817300143. [CrossRef] [Google Scholar]
 Keilegavlen E., Fumagalli A., Berge R., Stefansson I., Berre I. (2017) Porepy: An open source simulation tool for flow and transport in deformable fractured rocks. Technical report, arXiv:1712.00460 [cs.CE]. URL https://arxiv.org/abs/1712.00460. [Google Scholar]
 Flemisch B., Rainer H. (2008) Numerical investigation of a mimetic finite difference method, in: Finite volumes for complex applications V – Problems and perspectives, Wiley–VCH, Germany, pp. 815–824. [Google Scholar]
 Geiger S., Dentz M., Neuweiler I. (2011) A novel multirate dualporosity model for improved simulation of fractured and multiporosity reservoirs, SPE Reservoir Characterisation and Simulation Conference and Exhibition, 9–11 October, Abu Dhabi. [Google Scholar]
 Sausse J., Dezayes C., Dorbath L., Genter A., Place J. (2010) 3d model of fracture zones at soultzsousfor Sts based on geological data, image logs, induced microseismicity and vertical seismic profiles, C. R. Geosci. 342, 7, 531–545. ISSN 16310713. doi: http://doi.org/10.1016/j.crte.2010.01.011. URL http://www.sciencedirect.com/science/article/pii/S1631071310000489. [CrossRef] [Google Scholar]
All Figures
Fig. 1 Domain subdivision for N = 2. In gray Ω^{2}, in red Ω^{1}, and in blue Ω^{0}. 

In the text 
Fig. 2 Representation of the fracture normal associated with respect to each “side” (+ and −). (a) A bidimensional domain Ω^{2}, (b) The fracture is fictitiously enlarged. (c) For a 0d intersection. 

In the text 
Fig. 3 Representation of the degrees of freedom for a 2d and 1d grid. The pressure dof are represented by circles, red for the 2d grid and green for the 1d. The velocity dof are depicted by yellow diamonds for the 2d grid and purple diamonds for the 1d. The nodes of the 2d grid are moved only for visualization purpose. 

In the text 
Fig. 4 Representation of the coupling between dimension for the upwind discretization scheme. (a) Between two cells in the same dimension d. (b) Between a ddimensional cell and a d − 1dimensional cell. (c) Between a d + 1dimensional cell and a ddimensional cell. 

In the text 
Fig. 5 Example of the coarsening strategy adopted. (a) The computational grid is artificially forced to be finer at the tip of a fracture. (b) The resulting grid after the coarsening. Clearly the cell measures are comparable. 

In the text 
Fig. 6 Domain with fractures (in red) and fracture intersections (in blue) for problem in Section 5.1.1. 

In the text 
Fig. 7 Benchmark 1. (a) Pressure solution (range (1, 1.6)) with conductive fractures and computational grid. (b) Pressure solution (range (1, 3.6)) with blocking fractures and computational grid. (c) Zoom of the grid used. 

In the text 
Fig. 8 Benchmark 1 with conductive fractures. (a) Pressure along horizontal line at y = 0.7 with permeable fracture. (b) Pressure along vertical fracture at x = 0.5 with permeable fractures. (c) Pressure along the line (0.0, 0.1)–(0.9, 1.0) with blocking fractures. 

In the text 
Fig. 9 Benchmark 1, error evolution for the matrix and the fractures with conductive and blocking fractures. 

In the text 
Fig. 10 Domain with fractures (in red) for problem in Section 5.1.2. 

In the text 
Fig. 11 Benchmark 4. (a) Original grid composed by 12 302 2dcells and 35 153 VEM dof and coarsened grid composed by 4599 2dcells and 18 803 VEM dof. (b) A zoom on almost intersecting fractures for the original and coarsened grids, respectively. The zoom is referred to the small rectangle at position, approximately (360, 350). 

In the text 
Fig. 12 Benchmark 4. (a) Pressure solution computed with VEM. The pressure ranges in [0, 1 013 250] Pa. (b) Pressure along horizontal line at y = 500 m. (c) Pressure along vertical line at x = 625 m. 

In the text 
Fig. 13 Permeable fractures. (a) Reference and coarse solutions for pressure and velocity, as arrows. (b) Reference and coarse solutions for concentration of the passive scalar. 

In the text 
Fig. 14 Permeable fractures. (a) Comparison of passive scalar production at the outflow between the reference triangular grid and the coarse grid. (b) Temporal error decay with reference in black. 

In the text 
Fig. 15 Blocking fractures. (a) Reference and coarse solutions for pressure and velocity, as arrows. (b) Reference and coarse solutions for concentration of the passive scalar. 

In the text 
Fig. 16 Blocking fractures. (a) Comparison of passive scalar production at the outflow between the reference triangular grid and the coarse grid. (b) Temporal error decay with reference in black. 

In the text 
Fig. 17 On the left representation of the 20 fractures colored by their identification number. Pressure and velocity for both the reference and coarse grids. The pressure is scaled between 0 and 4.8 × 10^{7} Pa. 

In the text 
Fig. 18 (a) Three time steps on the reference grid for the concentration. (b) At the same steps, the concentration computed on the coarse grid. 

In the text 
Fig. 19 (a) The production comparison between the reference and the coarse solution. (b) An example of five coarse cells. 

In the text 