© A. Fumagalli & E. Keilegavlen, published by IFP Energies nouvelles, 2019

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₂ 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 equi-dimensional 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 lower-dimensional objects embedded in the rock matrix. This approach is commonly denoted reduced models, hybrid-dimensional models, and mixed-dimensional 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 co-dimensional, that is, lower-dimensional, 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 non-fractured 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 non-conforming 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 non-conforming approaches [35].

Finally, a fully non-matching coupling among the grids requires ad-hoc 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 non-conformity, or the class of Embedded Discrete Fracture Matrix Methods (EDFM) [5, 41–43], where approximate formulas for fracture-matrix transmissibilities are computed based on geometrical considerations.

In this paper we consider a conforming discretization with the dual virtual element approximation to simulate a mixed-dimensional Darcy problem and a finite volume discretization for the mixed-dimensional 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 mixed-dimensional transport problem. In this case we consider an upwind scheme, extended to handle the mixed-dimensional 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 matrix-fracture system and the fracture-fracture 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 mixed-dimensional 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 mixed-dimensional Darcy problem, presented in Section 2.1, to describe the flow and pressure, and the mixed-dimensional transport problem, presented in Section 2.2, to describe the motion of a passive scalar transported by the Darcy velocity.

The motivation for the mixed-dimensional 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 low-quality 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 three-dimensional 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 ² (A). Moreover, the trace operator on a domain A will be indicated by ·| _A .

2.1 Mixed-dimensional Darcy problem

Let us consider a regular domain , for N > 0, with outer boundary ∂Ω. The domain Ω is composed by a single, possibly non-connected, equi-dimensional domain Ω ^N and several lower-dimensional domains Ω ^d for d < N, which are possibly non-connected. 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.

We are interested in the mixed-dimensional Darcy problem which describes the pressure and velocity fields in Ω. Following the idea presented in [10], we define the pressure compound p = (p ⁰, …, p ^N ) which represents the pressure in each dimension, indicated as superscript. Similarly, we introduce the velocity compound u  = (0, u ¹, …, 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 < N-dimensional 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 + 1-manifold, it is assumed to be an orientable d-dimensional manifold, with unique normal n lying on the tangent space of the d + 1-manifold. This implies n ₊ = − n ₋, an illustration is given in Figure 2. In the particular case of T-type 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 lower-dimensional 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.

2.1.1 Fixed-dimensional Darcy problem

We consider the Darcy problem for the pressure and velocity as

(1a)

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 ¹, …, 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 cross-sectional volume (N − 3). The permeability of a lower d-dimensional (d < N) object is then defined as , with the permeability of its equi-dimensional 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 ⁰, …, 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 1-dimensional objects involved in the intersection.

We emphasize that the mixed-dimensional 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

(1b)

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 well-posed problem, boundary conditions must be assigned on the external boundary ∂Ω and the boundary of the sub-domains 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 Robin-type condition.

On the remaining internal portion , for d > 0, we assign the so-called “tip-condition”, namely

(1d)where in the previous equation n stands for the outward unit normal of .

2.1.4 Generalized formulation

Equation (1) is the mixed-dimensional 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 mixed-dimensional gradient operator such that

Considering the Darcy velocity composed by u in Ω and u n on Γ, system (1) becomes

(1a, 1b-bis)

In the previous equation is defined accordingly to include both the effective tangential and normal permeabilities.

2.2 Mixed-dimensional transport problem

Once the problem (1) is solved, the velocity field u can be used to describe the mixed-dimensional transport problem. We consider the same splitting of Ω as in Section 2.1. We introduce the concentration compound as c = (0, c ¹, …, 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 multi-dimensional treatment of fracture intersections. Other important references on this subject are [20, 22].

2.2.1 Fixed-dimensional 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 lower-dimensional 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 mixed-dimensional 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

(2a-bis)which is a general form of a conservative equation in the mixed-dimensional 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 mixed-dimensional 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 fixed-dimensional 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 Robin-type 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 mixed-dimensional 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 sub-domain 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 mixed-dimensional 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 Mixed-dimensional dual VEM discretization

We present a numerical scheme to solve the mixed-dimensional 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 higher-order 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 a-priori and is implicitly defined by . The curl-free 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.

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 ₀ = I−Π₀. 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:

and vectors

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 ₀ is considered only for clearness and to preserve the structure of the matrix. In practice, it is possible to remove u ₀ from the system.

4.2 Mixed-dimensional 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 semi-discrete 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 Kronecker-type 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 co-dimensional grid which is in communication with the face e. Finally, on side giving an element E of the last term of the semi-discrete 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 co-dimensional coupling, see Figure 4 for an example. Note that the chosen discretization is compatible with the coupling condition (2b).

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 lower-dimensional grid. We obtain also

where is the set of all neighbor cells of E in the higher-dimensional 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 right-hand 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 star-shaped 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.

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 cell-wise 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 mixed-dimensional 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 ² and for the fractures err ¹, where and E ² is a cell in the mixed-dimensional problem related to the rock matrix, respectively E ¹ for the fracture. The errors are computed on the reference grid, p ¹ is assumed constant on the normal direction of the fracture in the equi-dimensional 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⁻⁴. 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⁴ and has conductivity with permeability 10⁻⁴. 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.

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.

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 non-convex 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].

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⁻¹⁴ m², uniform fracture permeability 10⁻⁸ m², and fracture aperture 10⁻² m². 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.

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

As shown in Figure 11, some of the grid cells are non star-shaped 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.

5.2 Passive scalar transport

In this part we consider both mixed-dimensional 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⁻² m, a pressure gradient from the right to left boundary of the domain of 3 × 10⁷ Pa, a final simulation time of 40 years.

5.2.1 Permeable fractures

We consider highly permeable fractures with tangential permeability of 5 × 10⁻⁶ m² and normal permeability of 2.5 × 10⁻⁹ m². The matrix permeability is set to 2.5 × 10⁻¹¹ m². 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.

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 ² difference from a reference solution obtained with 10⁵ 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⁻¹⁶ m². Matrix permeability is set to 7.5 × 10⁻¹¹ m². 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.

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 3d-network

In this example we consider a realistic geometry for a geothermal system. We study a partial reconstruction of fractures from test site at Soultz-sous-Forê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⁻¹⁰ m² and fracture permeability, in both the normal and tangential direction, equal to 5 × 10⁻⁵ m². The fracture aperture is ϵ ² = 10⁻² 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 point-cells 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.

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.

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 mixed-dimensional problems able to describe a single-phase 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 mixed-dimensional 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

All Figures