Comparison and Application of non-Conforming Mesh Models for Flow in Fractured Porous Media using dual {L}agrange multipliers
Patrick Zulian, Philipp Schädle, Liudmila Karagyaur, Maria Nestola
CC OMPARISON AND A PPLICATION OF NON -C ONFORMING M ESH M ODELS FOR F LOW IN F RACTURED P OROUS M EDIA USINGDUAL L AGRANGE MULTIPLIERS
A P
REPRINT
Patrick Zulian
Institute of Computational Science, USI Lugano, 6904 Lugano, Switzerland [email protected]
Philipp Schädle
Institute of Geophysics, Department of Earth Sciences, ETH Zürich, 8092 Zürich, Switzerland. [email protected]
Liudmila Karagyaur
Institute of Computational Science, USI Lugano, 6904 Lugano, Switzerland [email protected]
Maria G. C. Nestola
Institute of Computational Science, USI Lugano, 6904 Lugano, Switzerland.Institute of Geochemistry and Petrology, Department of Earth Sciences, ETH Zürich, 8092 Zürich, Switzerland [email protected] , [email protected] August 17, 2020 A BSTRACT
Geological settings with reservoir characteristics include fractures with different material and geo-metrical properties. Hence, numerical simulations in applied geophysics demands for computationalframeworks which efficiently allow to integrate various fracture geometries in a porous mediummatrix. This study presents a modeling approach for single-phase flow in fractured porous mediaand its application to different types of non-conforming mesh models. We propose a combinationof the Lagrange multiplier method with variational transfer to allow for complex non-conforminggeometries as well as hybrid- and equi-dimensional models and discretizations of flow throughfractured porous media. The variational transfer is based on the L -projection and enables an accurateand highly efficient parallel projection of fields between non-conforming meshes (e.g., betweenfracture and porous matrix domain).We present the different techniques as a unified mathematical framework with a practical perspective.By means of numerical examples we discuss both, performance and applicability of the particularstrategies. Comparisons of finite element simulation results to widely adopted 2D benchmark casesshow good agreement and the dual Lagrange multiplier spaces show good performance. In anextension to 3D fracture networks, we first provide complementary results to a recently developedbenchmark case, before we explore a complex scenario which leverages the different types of fracturemeshes. Complex and highly conductive fracture networks are found more suitable in combinationwith embedded hybrid-dimensional fractures. However, thick and blocking fractures are betterapproximated by equi-dimensional embedded fractures and the equi-dimensional mortar method,respectively. a r X i v : . [ phy s i c s . c o m p - ph ] A ug on-Conforming Mesh Models for Flow in Fractured Porous Media using Lagrange multipliers Fluid flow through fractured porous media is a crucial process in the context of numerous subsurface applications, e.g.groundwater management, geothermal energy utilization, CO sequestration, hazardous waste storage, and enhancedoil and gas recovery [1–8]. Often, fluid flow velocities in fractures and in the porous medium matrix range over manyorders of magnitudes. Therefore, single fractures and networks of fractures largely govern the characteristics of fluidtransport in fracture-dominated porous media. More specifically, the detailed fracture geometry of each individualfracture has a significant influence on the fluid flow in a fracture-dominated system. While the aperture width of afracture can vary over several orders of magnitudes, the fracture can additionally be filled with an infilling materialexhibiting permeability values ranging over many orders of magnitude. As a consequence, these two aspects lead tostrong differences in the transmissivity and therefore the fractures ability to permit fluid transport. Thereby, the presenceand permeability of an infilling material determines if a single fracture acts as a conduit or a barrier for fluid flow [9–12].A detailed description of fluid flow through fractured porous media therefore requires comprehensive knowledge aboutthe hydraulic properties of each individual fracture. Such properties are very difficult to obtain in the field and adetailed deterministic description of such systems is not possible [13]. Consequently, stochastic investigations areuniversally conducted to describe fractured media and to account for associated uncertainties [9, 14]. Due to the largenumber of forward simulations required by stochastic studies, highly efficient and accurate numerical methods andmesh generation approaches are needed [15–19].To numerically model fractured systems, two method classes are widely used, continuum models and discrete fracturemodels. In the class of continuum models, fractures and porous-medium matrix are represented by separate continuawithin the same geometric mesh [20–23]. Effective flow properties are obtained by upscaling and information betweenthe continua needs to be transferred. In the class of discrete fracture models, fractures are represented as discretedomains in a numerical mesh [24, 25]. Here, two concepts are distinguished where the porous-medium matrix is eitherrepresented in discrete-fracture-matrix models (DFM) or neglected in discrete-fracture-networks models (DFN) [9, 26].In classic DFMs, fracture and porous-medium matrix domains are explicitly meshed and conforming at the domaininterface, i.e. conforming geometry and discretization [27–30]. Due to the large complexity of fracture networks, meshgeneration for such conforming DFMs can be very challenging and time consuming [31–33]. Such challenges occur inthe matrix mesh domain as well as the fracture mesh domain. Due to the large length-to-width ratio of many fractures itis common to use lower-dimensional elements to represent the fracture domain, e.g. [34–37]. This avoids elements withlarge aspect ratios in the fracture mesh and thus improves numerical performance. However, fractures with considerablylarge aperture width are not ideally represented by lower-dimensional elements [38]. Those fractures might be lessnumerous but require to generate fracture domain meshes which are equi-dimensional to the porous medium matrixdomain. Further challenges are posed by the matrix mesh generation around the fractures. This is particularly difficultwhere fractures are close to each other or intersect with very low angles. Such configurations can lead to elementswith large aspect ratios or non-physical connections, which requires fine tuning of the meshes to improve performanceand stability of the solution [38]. It is important to bear in mind that this is even more challenging in 3D and makesstochastic studies with DFMs very challenging. These challenges in model generation combined with the requirementsand geometrical complexities of DFMs motivated a large number of method developments and improvements, andyielding this to be an active research field.The drawbacks associated with the classic discrete-domain approaches have triggered research on, and the developmentof, numerical methods that allow to use individual meshes for the fracture and porous medium matrix domains. Suchmethods rely on the concept of non-conforming discretizations, while the meshes might be geometrically conformingor non-conforming. Methods with non-conforming discretization but conforming geometries require the element facetsof the fracture domain to align with the neighboring element facets of the porous-medium matrix domain withoutcoinciding, e.g. mortar methods [39–44] and discontinuos Galerking methods [45]. In contrast, fully non-conformingmethods, i.e. non-conforming discretization and non-conforming geometry, require no geometrical relationship betweenthe fracture and the porous-medium matrix domains. These methods exist for finite volume schemes, e.g. (p)EDFM[46–50] and for finite element schemes, e.g. extended finite element methods (XFEM) [51], or continuous Galerkinmethod where fractures are represented as Dirac functions [52]. XFEM approaches exist with primal formulations [53–55] and dual formulations [56, 57]. An exposition of different coupling techniques in discrete fracture networks isprovided in [58], where the focus is on handling complex fracture networks without accounting the effects of thesurrounding rock. Recently, Köppel et al. [59] presented an alternative fully non-conforming finite element formulationfor which Schädle et al. [60] demonstrated its applicability in 3D. Often, non-conforming mesh methods allow foran automated mesh generation and model setup process, which enables stochastic studies with a large number ofdifferent fracture network geometries. In particular, methods that handle fracture and matrix meshes separately (i.e.non-conforming geometries), significantly reduce the work required for preparing the simulation geometries. Berreet al. [61] provide a good overview of existing conceptual and discretization methods and discuss the differences inconforming and non-conforming methods.The mutually non-conforming discretizations of non-conforming mesh approaches require the use of coupling tech-2on-Conforming Mesh Models for Flow in Fractured Porous Media using Lagrange multipliersniques and information transfer between the domains. Lagrange multipliers are a common tool to couple systems ofequations and they are widely used throughout various fields [62, 63]. More specifically, fictitious domain [64] andmortar [65] methods are discretization techniques based on the Lagrange multiplier [66, 67]. In the context of DFMs,Lagrange multipliers have thus far been applied for both, fully non-conforming and non-conforming discretizationswith conforming geometries. Frih et al. [68] and Boon et al. [69] presented approaches where the porous-matrixdomain is split along the fracture planes. The resulting sub-domains are then meshed independently and glued togetherusing the mortar method. Alternatively, discontinuous Galerkin methods are used to handle non-conforming interfacesbetween any pair of elements, hence providing flexibility at a finer granularity with respect to the mortar method [45].A fully non-conforming approach using Lagrange multipliers to couple fracture and porous matrix discretizations isstudied by Köppel et al. [59] and Schädle et al. [60]. Particularly in [60], the L -projection with different discreteLagrange multipliers are used to transfer information between the fracture and porous-medium matrix domains in avariationally consistent way. Both, Köppel et al. [59] and Schädle et al. [60] studied flow through DFMs with fracturesof codimension one. While these hybrid-dimensional DFMs are widely applied and very efficient with respect tomeshing and numerical performance they are less suited for fractures with large aperture widths [38]. For such cases,equi-dimensional descriptions of the porous-medium matrix and the fractures provide more accurate results. However,this requires a volume-to-volume coupling between the fracture and the matrix discretizations. Consequently, the L -projection needs to be constructed considering volumetric polyhedral intersections. Volumetric coupling combinedwith variational transfer has been studied for several applications related to fluid-structure interaction [62, 70]. Notethat, as for the hybrid-dimensional case, local mass conservation is generally not guaranteed and jumping pressurecoefficients can not be represented properly with our continuous Galerkin approach and Lagrange finite element spaces.One underlying assumptions of the Lagrange multiplier approach is a continuity of the solution across the fractures,which prevents them to act as a flow barriers. However, some fractures are populated with an infilling material withvery low permeability and therefore hinder fluid flow. To model such fractures, equi-dimensional sub-domains withvery low permeability values need to be employed and then coupled to the porous-medium matrix using the mortarmethod. Furthermore, Schädle et al. [60] found that the solution is less accurate in areas with fracture intersections andboundaries, which is more enhanced by steep pressure gradients located in these areas. Local adaptive mesh refinementwould improve accuracy in these cases. Schädle et al. [60] employed dual-Lagrange multipliers [66, 71], which haveshown to have a positive impact on the condition number and further allow to construct a symmetric positive-definitesystem of equations. Solving such systems is much more convenient than solving saddle-point problems, whichare indefinite systems and are typically harder to solve. These types of systems can be solved with a wide range ofmethods (e.g., conjugate gradient method) and preconditioners (e.g., multigrid [72]), and enable to perform large scalecomputations in a convenient manner. Furthermore, the dual-Lagrange mulitplier space facilitates combining differentcoupling strategies, i.e. fully non-conforming hybrid-dimensional and equi-dimensional as well as mortar coupling, ina simple and unified way.This paper presents a unified framework based on the method of Lagrange multipliers, which combines embeddeddiscretization methods with non-conforming domain decomposition approaches. Embedded discretization methods aredesigned to couple overlapping meshes which are mutually non-conforming and can have non-matching geometricfeatures. Here, the finite element discretization of the matrix is coupled with any number of fracture discretizations,which can be hybrid-dimensional and equi-dimensional. Non-conforming domain decomposition techniques, such asmortar, allow us to split the domain into multiple sub-domains, discretize them independently, and couple them at theirinterfaces.This combination allows to employ each technique to scenarios where it is most suited, i.e. blocking fractures, fractureswith large apertures, or fractures with large aspect ratio in networks with many fractures. By choosing the dualLagrange multipliers space for each of the aforementioned coupling approaches, the arising systems of equations areeasily combined and condensed, as mentioned earlier, into a unique symmetric positive-definite matrix. We showhow non-conforming adaptive mesh refinement is combined with the variational transfer, and we employ it to controlthe error as well as to reduce the number of degrees of freedom in the arising system of equations. We study theapproaches both in isolation and combined. In particular, the accuracy and performance of the presented frameworks isdemonstrated by comparison to standard benchmark cases in 2D and 3D as well as a realistic scenario which combinesthe different approaches.In Section 2 we describe the overall methodology. We illustrate the unified formulation of the flow model (Section 2.1),its variational formulation (Section 2.2), the finite element discretization (Section 2.3), the necessary steps for couplingthe different discretizations with dual Lagrange multipliers (Section 2.5), and the construction of the algebraic systemof equations. In Section 2.6, we show how non-conforming mesh refinement can be integrated within the couplingframework, followed by some specific details about the implementation in Section 2.7. Numerical investigations andexperiments are illustrated and discussed in Section 3. Finally, a conclusion of our findings and future developments areprovided in Section 4. 3on-Conforming Mesh Models for Flow in Fractured Porous Media using Lagrange multipliers The method of Lagrange multipliers allows to discretize flow problems for porous media with two main types ofnon-conformity. First, the matrix is split into sub-domains which can be discretized independently then glued togetherusing the mortar method [71]. Not only the sub-domains can differ in terms of permeability, but their interface canrepresent fractures. Second, the fractures are represented as separate bodies embedded in the matrix. Such fracturesare described by either lower-dimensional manifolds or equi-dimensional manifolds. In this section, we present aunified framework to describe the different geometric representations depicted in Figure 1, and the related discretizationtechniques. Ω Ω ξ γ γ γ n γ n γ − n γ Figure 1: Two-dimensional example with an embedded lower-dimensional fracture γ , a lower dimensional (line)fracture γ at the interface ξ , , between sub-domains Ω and Ω of the matrix, and an embedded equi-dimensional(polygon) fracture γ . The dashed lines represent an nonexistent spacing which is employed only for illustrativepurposes. Let Ω ⊂ R d , d ∈ { , } be the matrix domain with the following decomposition into N sub-domains Ω = N (cid:91) i =1 Ω i , where Ω i ∩ Ω j = ∅ , i (cid:54) = j . If Ω i and Ω j are connected, hence Ω i ∩ Ω j (cid:54) = ∅ , their interface is described by ξ ij = ∂ Ω i ∩ ∂ Ω j ∩ Ω . With Ξ = { ξ ij } .Let γ ⊂ Ω be a manifold of dimension d or d − describing the fracture domains with the following decomposition γ = N γ (cid:91) k =1 γ k . When required we distinguish the dimension of the manifold γ k , with γ d − k we have a lower-dimensional fracture, andwith γ dk an equi-dimensional fracture. If the interface ξ = ξ ij is interpreted as a lower-dimensional fracture we employthe short-hand notation γ ξ = γ d − ξ .Steady state fluid flow in the matrix Ω is governed by ∇ · ( − K ∇ p ) − λ = f in Ω , (1)with p = p on ∂ Ω D , where p is the pressure, K ∈ R d × d is a bounded symmetric positive definite permeability tensor, f is a sink/source term, p is a given pressure on the boundary ∂ Ω of the domain of interest Ω .Flow in the fracture-network γ is described by ∇ · ( − K γ ∇ p γ ) + λ = f γ in γ, (2)with p γ = p γ on ∂γ D . 4on-Conforming Mesh Models for Flow in Fractured Porous Media using Lagrange multipliersThe fluid exchange between Ω and γ is given by the Lagrange multiplier λ ∈ Λ( γ ) . The function spaces V Ω and V γ aredefined by V Ω = H (Ω) , W Ω = H (Ω) ,V γ = H ( γ ) , W γ = H ( γ ) ,V = V Ω × V γ , W = W Ω × W γ , where H is the Sobolev space of weakly differentiable functions, and H ⊂ H its restriction to functions vanishingat the boundary. The Lagrange multiplier space is defined as the dual of W with the product Λ = N γd − (cid:89) k =1 H − ( γ d − k ) × N γd (cid:89) k =1 H − ( γ dk ) , N γ = N γ d − + N γ d . Note that the Lagrange multiplier is extended by zero outside γ . For each γ k ∈ γ the corresponding Lagrange multiplieris denoted with λ k ∈ Λ k = Λ( γ k ) .Depending on the type of fracture the pressure term and the Lagrange multiplier have slightly different meanings. Forthe lower-dimensional fracture, i.e., γ d − k the pressure term p γ k represents the average pressure across the fracture withtangential permeability K γ k , and the Lagrange multiplier λ k = [ K ∇ p · n γ k ] = ( K ∇ p − K γ k ∇ p γ k ) · n γ k p = p ( x ) , x ∈ γ d − k represents the jump of the fluid pressure gradient in normal direction n γ k with respect to the fracture surface γ k .For the embedded equi-dimensional fracture the Lagrange multiplier can be thought as a reactive force field introducedin order to ensure the continuity of the pressure.For a more compact notation, the aperture of lower dimensional fractures is neglected and considered in the permeabilitytensor. With ( · , · ) Ω and ( · , · ) γ we denote the L -inner product over Ω and γ , respectively. The variational formulation is foundby multiplying (1) and (2) by test functions and integrating over the domains Ω and γ , using integration by parts. Hence,the weak form of the coupled system of equations is given as follows: find ( p, p γ ) ∈ V and λ ∈ Λ , such that ( K ∇ p, ∇ q ) Ω + ( K γ ∇ p γ , ∇ q γ ) γ − ( λ, q − q γ ) γ = ( f, q ) Ω + ( f γ , q γ ) γ ∀ ( q, q γ ) ∈ W, (3)and the weak equality condition ( p − p γ , µ ) γ ∀ µ ∈ Λ , (4)are satisfied. The variational formulation introduced in Section 2.2 is discretized using the finite element method. Depending onthe settings different meshes M i = M Ω i and M γ k are used for the different sub-domains of the porous matrix Ω i , i = 1 , . . . N , and for the fractures γ k , k = 1 . . . N γ respectively. The presented techniques and their implementationallow for an arbitrary choice of M λ k , which is the mesh associated with the Lagrange multiplier, however since werestrict ourselves to a particular choice of multiplier space, we set M λ k = M γ k .The method allows for a wide variety of elements for each of the meshes. For a manifold with dimension d we employeither Lagrange elements P k , or tensor-product elements Q k of order k ∈ { , } W h,α = { w ∈ W ( α ) : ∀ E ∈ M α ,w | E ∈ (cid:40) P k if E is a simplex Q k if E is a hyper-cuboid (cid:41) } ,α ∈ { Ω , γ d , γ d − } Λ h,β β ∈ { γ d , γ d − } , (5)5on-Conforming Mesh Models for Flow in Fractured Porous Media using Lagrange multiplierswhere Λ h is a discrete Lagrange multiplier space. Let { ϕ i } i ∈ J be a basis of W h, Ω , { θ j } j ∈ J γ a basis of W h,γ , and { ψ k } k ∈ J γ a basis of Λ h , where J and J γ ⊂ N are index sets of the node-sets of their respective meshes M and M γ .Writing the functions q ∈ W h, Ω , q γ ∈ W h,γ , and µ ∈ Λ h,γ in terms of their respective bases and coefficients, they read q = (cid:80) i ∈ J q i ϕ i , q γ = (cid:80) j ∈ J γ q jγ θ j , and µ = (cid:80) k ∈ J γ µ kγ ψ k .The dual shape functions ψ j ∈ Λ h ( γ ) are constructed in such a way that they satisfy the bi-orthogonality condition [71]: ( θ i , ψ j ) γ h = δ ij ( θ i , γ h ∀ i, j ∈ J γ , (6)and integral positivity ( ψ j , γ h > . (7)Note that (7) is naturally satisfied for first order finite elements. For second order elements we follow the constructiondescribed in [66, 73].After reformulating the variational problem (3) as a set of point-wise equations, the discrete problem for the porousmatrix reads: (cid:88) i ∈ J p i ( K ∇ ϕ i , ∇ ϕ j ) Ω h − (cid:88) k ∈ J γ λ k ( ψ k , ϕ j ) γ h = ( f, ϕ j ) Ω h , ∀ j ∈ J, (8)which translates to the linear system Ap − B T λ = f .The fracture equations result in, (cid:88) i ∈ J p i,γ ( K γ ∇ θ i , ∇ θ j ) γ h + (cid:88) k ∈ J γ λ k ( ψ k , θ j ) γ h = ( f γ , θ j ) γ h , ∀ j ∈ J γ , (9)which translates to the linear system A γ p γ + D T λ = f γ . The weak-equality condition (4) results in, − (cid:32)(cid:88) i ∈ J p i ( ϕ i , ψ j ) Ω h − (cid:88) k ∈ J p k,γ ( θ k , ψ j ) γ h (cid:33) = 0 , ∀ j ∈ J λ , (10)which translates to the linear system − Bp + Dp γ = . The discretization of the complete problem results insaddle-point system as follows: (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) A 0 − B T γ D T − B D 0 (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) pp γ λ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ff γ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , (11)However, the trivially invertible matrix D enables us to perform block Gaussian elimination and obtain the followingstatically condensed system [74] ( A + T T A γ T ) p = f + T T f γ , (12)where T = D − B . Once the system is solved for p , the solution for the fracture network can be computed by p γ = Tp .The resulting system matrix is symmetric positive definite which allows us to adopt optimal solution strategies such asMultigrid methods [72]. The presented methods handles non-conforming meshes in two stages. The first stage involves the domain decompositionof the porous matrix. For instance, a matrix domain Ω can be split into the two domains Ω and Ω with interface ξ , ,this interface could be interpreted as a fracture. However, in discrete settings (Fig. 2 (a)) it could also be a convenientway to handle different resolutions for the meshes M and M , and ensure the continuity of the solution at ξ , using the standard mortar approach introduced in [65]. In practice, we need to assign the standard master and slaverole, for instance, we assign the master role to ξ and slave role to ξ . Once the transfer operator is assembled theporous-medium-matrix system is condensed. Note that in the resulting system of equations the degrees of freedomassociated with the slave discretization are eliminated after static condensation (e.g., by replacing the related rows of thematrix with the identity and the related right-hand side value with zero). Once this is achieved we go to the next stage.The second stage involves the embedded case (Figure 2(b)) where we compute the transfer operator between W h and W h,γ and condense the resulting system as described in Section 2.3.6on-Conforming Mesh Models for Flow in Fractured Porous Media using Lagrange multipliers ℳ ℳ ξ ℳℳ γ (a) Non conforming interface ξ , between mesh M and mesh M . (b) Meshes for fracture M γ and porous matrix M for the embedded scenario. Figure 2: Non-conformity with matching geometry (a) and non-matching geometry (b).
The computation of the so called mortar integrals, which are the integral terms in (8), (9), and (10), associated withthe Lagrange multiplier require special handling. In fact, quadrature formulas have to be generated in the intersectionbetween elements of the matrix and the fracture. The computation of intersections differs depending on which type offracture is considered. Table 1 lists, for each type of coupling, the roles and the intersection types. Here, the master roleis given to the domain covering completely the slave domain. The slave role is given to the domain with which weassociate the Lagrange multiplier space.For the coupling at an interface ξ , either for non-conforming domain decomposition or for an interface fracture γ ξ ,the coupling is performed on a common surface description. This operation requires intersecting oriented planarpolygonal elements in the case of a three-dimensional problem, or intersecting oriented line elements in the case oftwo-dimensional problem.For the embedded scenario the polytopal element of the matrix are intersected with the lower- or equi- dimensionalpolytopal elements of the fracture. The case where we have M ⊂ R d and M γ being a ( d − -dimensional manifoldmesh, requires particular handling when the fracture elements are aligned with the surface of the matrix elements. Insuch case it is likely that some intersections might be computed twice, hence these duplicate intersection are detectedand removed.The ( d − -dimensional dimensional fractures represented with the mortar method require a careful set-up. In fact,in the case of intersecting fractures, multiple sub-domains (more than two) might intersect at one point or edge. Themesh primitives, i.e. edges and nodes, generating such intersections require a specific handling. The first approachconsists of ignoring the side elements that are incident to the aforementioned mesh primitives when defining thediscrete Lagrange multiplier as in Krause et al. [75], hence introducing discontinuities of the solution at these interfaces.The second approach, which would require a set-up similar to Farah et al. [76], consists of defining one-to-manyrelationships for the intersecting primitives. Here, one master primitive has to be determined and continuity is eitherenforced using interpolation for intersecting nodes, or with a weak equality condition for intersecting edges. Thethird approach, consists of explicitly defining the entire surface mesh for which the discrete Lagrange multiplier isconstructed. These complications with the mortar method are not present in the equi-dimensional case if the discretefractures are represented with a unique conforming mesh, since the role of master and slave can be trivially assigned toporous-matrix and fracture respectively.The computation of the intersections listed in Table 1 is performed with suitable variants of the Sutherland-Hodgmanclipping algorithm [77]. Once the intersection is computed, if required, this intersection is meshed into a simplicialcomplex so that we can map quadrature rules to each simplex and integrate exactly. Let us recall the definition of element E from Section 2.3. With ∂E we denote the boundary of element E and with ¯ E its closure. A mesh is said to be conforming if ¯ E i ∩ ¯ E j , i (cid:54) = j is a common vertex, edge, face, or ∅ . A node is said to behanging if it lies on the interior of an edge or face of another element. A mesh containing at least one hanging node iscalled non-conforming. 7on-Conforming Mesh Models for Flow in Fractured Porous Media using Lagrange multipliers Matrix (master) Fracture (slave) Intersection type Ω ⊂ R γ polyhedron-polyhedron Ω ⊂ R γ polyhedron-polygon ∂ Ω i ∩ ξ i,j ⊂ R ∂ Ω j ∩ ξ i,j ⊂ R polygon-polygon (oriented) Ω ⊂ R γ polygon-polygon Ω ⊂ R γ polygon-segment ∂ Ω i ∩ ξ i,j ⊂ R ∂ Ω j ∩ ξ i,j ⊂ R segment-segment (oriented)Table 1: Intersection for different coupling types. The standard master and slave roles used in the mortar literature areassociated with matrix and fracture, respectively. c ij p i p j E l E l +1 E l +1 E l +1 E l +1 Figure 3: Bi-linear elements E k with different levels of refinement k ∈ { l, l + 1 } . With p i , p j we denote the parentnodes that are generated on level l , with c ij = ( p i + p j ) / we denote the child node generated on level l + 1 by splittingthe edge ( i, j ) . Depending on the finite element discretization the splitting is done in different (multiple) locations. Theprolongation operation at node c ij for this particular case would simply be u ( c ij ) = ( u ( p i ) + u ( p j )) / .Non-conforming adaptive mesh refinement has the advantage that can be applied to any type of element in a ratherstraight-forward and independent manner. However, once the elements marked by the error estimator are refined, theresulting mesh might have hanging nodes, as shown in Figure 3.As a consequence, continuity of the solution is to be enforced either by employing variational restriction or discontinuosGalerkin methods. Here, we consider the variational restriction technique which is thoroughly explained in ˇCervenýet al. [78] also for high-order discretizations.Let R i , i ∈ { m, s } be a suitable restriction operator that splits contributions of each hanging node to its adjacentnodes, where m stands for master and s for slave. For combining adaptivity with DFMs and the method of Lagrangemultipliers, we consider the constrained spaces arising from the refinement and variational restriction, which requiresus to perform some slight modifications to the final steps of the assembly procedure for the transfer operator, in any ofthe cases we mentioned in previous sections. We recall the definitions of the coupling matrix B and mass-matrix D from Section 2.3, and define the modified transfer operator T R = ( R Ts DR s ) − ( R Ts BR m ) . The operator T R allows us to transfer between constrained spaces, however for also setting the values in the hangingnodes we apply the prolongation operator P s = R Ts as follows T = P s T R . This small modification allows to use T as in the standard case without any special treatment as described in Section 2.2. The routines described in this paper are implemented within the open-source software library
Utopia [79]. In this work,
Utopia uses libMesh [80] for the finite element discretization,
MOONoLith [81] for the intersection detection, and
PETSc [82] for the linear algebra calculations. The software developed for this contribution is used by means of a JSON8on-Conforming Mesh Models for Flow in Fractured Porous Media using Lagrange multipliers(JavaScript Object Notation) input file where any number of mesh files can be linked to the simulation and coupledtogether automatically for creating complex networks of fractures.
First, the focus is on a comparison of results obtained with the equi-dimensional embedded technique and the mortarmethod to a specific selection of commonly used 2D and 3D benchmarks [38, 83]. Here, we use non-conforming meshrefinement in proximity of the fractures for maintaining the size of the mesh small while achieving small deviationsfrom the reference solution. Second, we show how employing an adaptive mesh refinement, with gradient-recoverybased error estimator [84], allows us to refine only where the error in the solution is estimated to be higher. Finally,we present a complex 3D experiment inspired by realistic scenarios where we conveniently mix all the techniques wecovered in this article. In each of the following sections we discuss the practicalities of the different techniques and howto combine them.For compactness, in the following sections we use the abbreviation ED for equi-dimensional and HD for hybrid-dimensional. We report exclusively the error of the solution associated with the matrix discretization. This is donefor avoiding redundant information, since the solution for the fracture is just the L -projection of the solution for thematrix. The 2D settings allows us to provide a simpler and clearer presentation of the numerical results. Hence, with a selectionof 2D benchmarks from Flemisch et al. [38], we complement and extend the previous contribution presented in Schädleet al. [60]. We verify the embedded ED approach and our implementation of the mortar method with the dual Lagrangemultiplier. We show how adaptive mesh refinement allows us to solve problems with a smaller number of degrees offreedom while achieving the desired accuracy in the solution. For all 2D cases we use the reference solution proposedin Flemisch et al. [38] which is computed the mimetic finite difference method. (a) Embedded/conductive (b) Mortar/blocking
Figure 4:
Benchmark 1 in Flemisch et al. [38]. Pressure solution for regular fracture network with six equi-dimensionalfractures for conductive fractures with the equi-dimensional embedded method (a) and blocking fractures with theequi-dimensional mortar method (b).We consider the same settings and reference solutions used in Flemisch et al. [38]
Benchmark 1 , with both conductiveand blocking fractures, as shown in Fig. 4. Both settings have the same square domain
Ω = [0 , and boundaryconditions. We impose Dirichlet conditions on the right boundary, where the pressure is set to constant value 1. Weimpose no-flow conditions on the bottom and top sides. Permeability of the matrix is uniform K = I , where I is theidentity matrix. We distinguish conductive and blocking scenarios for the fracture permeability. For the conductivescenario, in the HD case the permeability tensor is described as K γ = (cid:15) I · , where (cid:15) = 10 − is the fractureaperture, whereas for the ED case we have K γ = I · . For the blocking scenario, we only have the ED case with K γ = I · − . 9on-Conforming Mesh Models for Flow in Fractured Porous Media using Lagrange multipliers (a) Line AA (cid:48) (b) Line BB (cid:48) Figure 5:
Benchmark 1 in Flemisch et al. [38]: Pressure profiles along the lines AA (cid:48) (a) and BB (cid:48) (b) with a zoom intothe area with the largest deviation. Blue lines indicate Emebedded-ED, dashed red lines indicate Mortar-ED.Table 2:
Benchmark case 1 in Flemisch et al. [38]. For each method we report number of elements in the matrix ( ) and fracture ( ), number of degrees of freedom ( d.o.f. ), normalized number of non zero entries ( nnz/size ),condition number ( (cid:107) · (cid:107) - cond .) and error in the matrix ( err m ) with respect to the reference solution. Method (cid:107) · (cid:107) -cond. err m Embedded-ED 10 656 triangles 8648 5853 1.3e-3 3.7e6 3.3e-5Embedded-ED 47 142 triangles 8648 25 992 2.8e-4 1.8e8 1.9e-6Embedded-ED 96 762 triangles 8648 53 379 1.4e-4 1.2e9 2.8e-7Mortar-ED 12 791 triangles 34 592 32 935 2.3e-4 5.2e10 6.4e-8Mortar-ED 36 331 triangles 34 592 45 172 1.7e-4 8.8e10 2.8e-8Mortar-ED (blocking) 922 triangles 34 592 26 520 2.7e-4 9.4e7 3.3e-8Embedded-HD 1089 quads 112 1156 9.7e-3 3.0e4 9.7e-3Embedded-HD 16 641 quads 448 16 900 5.6e-4 1.6e6 2.5e-3Embedded-HD 66 049 quads 896 66 564 1.4e-4 1.3e7 1.3e-3
A particular emphasis is placed on equi-dimensional fractures and the comparison between non-conforming embed-ded/immersed fractures and geometrically conforming fractures which are glued together with the matrix using themortar method. For the hybrid-dimensional embedded case we consider the results of Schädle et al. [60] using themethods of dual-Lagrange multipliers and static condensation. In Table 2, we report results for three different resolutionsfor each of the three strategies. In the following paragraphs we illustrate the different set-ups, results, and limitations ofthe embedded and mortar methodologies.The employed embedded techniques enforce the continuity of the solution at the intersection of the matrix and fracturemeshes. Hence, steep pressure jumps and barriers can not be represented. Consequently, for the embedded case werestrict our study to conductive fractures. In particular, we study the equi-dimensional embedded technique for 3, 5, and6 levels of non-conforming mesh refinement (Section 2.6) performed exclusively in proximity of the fractures. Thisrefinement patterns are generated automatically using the variational transfer algorithm for marking the elements of thematrix that are intersecting with the fractures. On Fig. 5, we plot the solution over the lines AA (cid:48) (y = 0.75) and BB (cid:48) (x= 0.5) and it can be observed that even for low resolutions, the solutions are in agreement with the reference resultsof Flemisch et al. [38]. From Table 2, it can be observed that with this particular set-up we reach an error in the order of − . Note that the number of elements in the fracture network do not influence the number of degrees of freedom,hence they do not count for the computational cost of solving the linear system but only for the set-up phase which istypically cheaper.For the mortar method based experiment we are required to represent the fracture explicitly in the matrix mesh. For thisparticular scenario, the fracture is modelled as an equi-dimensional geometry which is meshed independently fromthe matrix. This allows us to refine all the different sub-domains in a completely independent manner and glue themtogether using the mortar method. We manually refine the matrix around the two fractures crossing at the center of thedomain. This is done for having a higher resolution in the region of interest of the benchmark. From Table 2, it can beobserved that even for low mesh resolutions of the matrix the error reaches an order of − for both, conductive andblocking fractures. However, the mesh resolution of the fractures is very high from the start. Despite this fact, in ourexperiments the actual degrees of freedom are only associated with one layer of nodes in the middle of the fracture,while the ones at the interface are eliminated by means of the mortar constraints as mention in Section 2.4.10on-Conforming Mesh Models for Flow in Fractured Porous Media using Lagrange multipliersWe can observe that the embedded methodologies generate linear systems with comparable condition numbers andnumber of degrees of freedom. Whereas, the mortar-ED method gives rise to larger systems with larger conditionnumbers. Only with the mortar-ED we are able to solve the blocking scenario 4(b), however already for this simpleexperiment the mesh set-up is more laborious due to the matching geometry. Figure 6:
Benchmark 2 in Flemisch et al. [38]. Embedded-ED: spatial distribution of the piezometric head [m]. Solutionprofiles are compared along the line AA (cid:48) with coordinates y = − [m] (marked in green). Figure 7: Embedded-ED adaptive mesh refinement. The meshes have been obtained by performing first an adaptiverefinement in the region where fractures are located, then by using the gradient-recovery based error estimator.We test the capability of the adaptive mesh refinement by considering the hydrocoin benchmark [85] for two differenttechniques: the embedded-ED and the mortar-ED. The matrix consists of a rectangular box with a broken line locatedon the top, whereas the fracture network consists of two oblique lines. As shown in Fig. 6, we prescribe the piezometric11on-Conforming Mesh Models for Flow in Fractured Porous Media using Lagrange multipliersTable 3: Uniform refinement for the Embedded-ED vs Mortar-ED cases. Here err m is the error computed with respectto the reference solution. UR refers to the number of uniform refinements. Method (cid:107) · (cid:107) -cond. err m UR Embedded-ED 272 triangles 3 480 159 4.3e-2 2.9e9 5.1e-4 1Embedded-ED 1 088 triangles 3 480 589 1.2e-2 1.1e10 8.5e-5 2Embedded-ED 4 352 triangles 3 480 2 265 3.1e-2 4.3e10 1.7e-5 3Embedded-ED 17 408 triangles 3 480 8 881 8.1e-4 1.8e11 4.0e-6 4Embedded-ED 69 632 triangles 3 480 35 169 2.1e-4 6.8e11 1.9e-6 5Embedded-ED 278 528 triangles 3 480 139 969 5.5e-5 2.7e12 1.3e-6 6Mortar-ED 1 116 triangles 220 246 9.7e-3 3.6e6 2.1e-4 1Mortar-ED 4 464 triangles 880 2 973 2.5e-3 7.6e6 1.9e-4 2Mortar-ED 17 856 triangles 3520 11 285 6.4e-4 2.7e7 2.0e-4 3Mortar-ED 71 424 triangles 14 080 43 941 1.7e-4 1.1e8 2.0e-4 4 P i e z o m e t r i c h e a d [ m ] Reference1 UR2 UR3 UR4 UR5 UR6 UR 0 500 1000 1500Arc length [m]107.5110.0112.5115.0117.5 P i e z o m e t r i c h e a d [ m ] Reference1 AR2 AR3 AR4 AR5 AR6 AR (a) Uniform Refinement (UR) (b) Adaptive Refinement (AR)
Figure 8: Pressure along line AA (cid:48) for the uniform (A) and the adaptive (B) Embedded-ED test cases. Here, UR andAR refer to different steps of uniform and adaptive refinements performed on the initial matrix mesh. In particular,the number of AR are obtained combining a gradient recovery strategy with an adaptive refinement performed on theoverlapping region between the matrix and the fracture zone. More details can be found in tables 4 and 3.head on the Dirichlet boundary on the top and Neumann no-flow on the remaining sides of the rectangular box. Thepermeability is set equal to K γ = 10 − [m/s] for the fractures and K = 10 − [m/s] for the matrix. Fracture aperture isabout (cid:15) = 10 . for the left fracture and about (cid:15) = 15 for the right fracture.We use a gradient recovery strategy for the a-posteriori error estimation to guide the adaptive refinement (AR) of thematrix and the fracture meshes, and compare the numerical results with those computed employing uniform refinement(UR). To this aim, we estimate the accuracy of each numerical simulation by computing the error with respect to thereference solution obtained by Flemisch et al. [38] using a mimetic finite difference (MFD) method on a very fine meshwith triangles and dofs.In Tables 3 and 4 we report the results related to the embedded-ED technique for UR and AR test cases, respectively.In particular, we specify the number of elements in the matrix mesh ( ) and in the fracture mesh ( ), theTable 4: Adaptive refinement for the Embedded-ED and Mortar-ED cases. Here, err g is the threshold used for thegradient recovery strategy. AR refers to the number of adaptive refinements, whereas the number within parenthesesrefers to the steps of adaptive refinement performed on the overlapping region between the matrix and the fracturemeshes. Method (cid:107) · (cid:107) -cond. err m g Embedded-ED 611 triangles 3 480 352 2.0e-2 1.5e10 8.9e-5 0 (2) -Embedded-ED 1 544 triangles 3 480 874 7.8e-3 5.3e10 1.6e-5 1 (2) 5Embedded-ED 2 378 triangles 3 480 1354 4.5e-3 8.1e10 7.2e-5 2 (2) 5Embedded-ED 7 913 triangles 3 480 4371 1.6e-3 3.5e11 2.9e-6 4 (2) 1Embedded-ED 9 147 triangles 3 480 4991 1.5e-3 2.8e11 2.0e-6 4 (4) 1Embedded-ED 21 599 triangles 3 480 11658 1.02e-3 1.06e12 1.4e-6 4 (6) 1Mortar-ED 2 127 triangles 202 1375 5.9e-3 3.7e11 2.1e-4 1 (1) 5Mortar-ED 8004 triangles 332 2791 1.1e-3 1.9e12 1.9e-4 2 (2) 5Mortar-ED 19 515 triangles 442 10 600 7.1e-4 1.2e12 1.9e-4 3 (2) 5Mortar-ED 28 692 triangles 457 15 284 4.8e-4 1.7e12 1.9e-4 4 (2) 5 P i e z o m e t r i c h e a d [ m ] Reference1 UR2 UR3 UR4 UR 0 500 1000 1500Arc length [m]107.5110.0112.5115.0117.5 P i e z o m e t r i c h e a d [ m ] Reference1 AR2 AR3 AR4 AR (a) Uniform Refinement (UR) (b) Adaptive Refinement (AR)
Figure 9: Pressure along line AA (cid:48) for the uniform (A) and the adaptive (B) Mortar-ED test cases. Again, UR and ARrefer to the steps of uniform and adaptive refinement performed on the initial matrix mesh. We use the MFD solutionproposed in Flemisch et al. [38] as a reference.number of degrees of freedoms ( d.o.f ), the density of non zeros entries ( nnz/size ), the condition number ( (cid:107) · (cid:107) -cond ),and the error computed with respect to the reference solution ( err m ). In Table 3 we also specify the number of uniformrefinements ( UR ), whereas in Table 4 we report the number of adaptive refinements ( AR ), and the error-thresholdused for the gradient recovery strategy ( err g ). We point out that both the UR and the AR strategies are only performedon the matrix mesh.In Fig. 7 we show the meshes obtained for all the six AR test cases. The AR based on the gradient recovery strategy iscombined with an AR perform ed on the overlapping region between the matrix and the fracture network. The numberof AR performed in the overlapping zone are reported in Table 4 and specified within parentheses.From the results summerized in Tables 3 and 4, we observe that the error err m progressively reduces by increasing thenumber of refinements. However, the adaptive strategy allows accurate results with coarser meshes when compared tothe uniform technique. While a uniform refined mesh with more than
200 000 elements ( ) is needed to get anerror ( err m ) close to − , an adaptively refined mesh with a total number of elements ten times smaller ensures thesame accuracy.Fig. 8 shows the piezometric head obtained over line AA (cid:48) for the UR (a) and the AR (b) test cases. One may note largediscrepancies in the area of the left fracture, especially for the results referred to the uniform refined meshes. In thisregard, we found that the a-posteriori error err g was higher close the sharp top boundary and in the zone occupied bythe left fracture. Indeed, Fig. 7 (3)-(5) show that the AR strategy produces a mesh size reduction in such regions. Hence,all the AR test cases reveal a better agreement with the reference solution when compared to the UR scenarios.For completeness, the numerical results obtained for the test case AR are presented in Fig. 6. Here, the colour refersto the spatial distribution of the piezometric head, whereas the contour lines are in black.We perform the same analysis for the mortar-ED scenario and collect all the numerical results in Tables 3 and 4. Again,the use of AR allows us to achieve an accuracy comparable to the UR test cases with less elements in both the matrixmesh ( ), and the fracture mesh ( ).In Fig. 9 we observe that the piezometric head profile computed along line AA (cid:48) converges to the reference solution(MFD) by increasing the number of refinements for all the numerical simulations. We point out that, while theembedded-ED test cases reveal large discrepancies with respect to the reference solution in the region of the left fracture,only little differences are observed for the mortar-ED scenarios. Indeed, the mortar approach requires a matchinggeometry at the interface between fracture and matrix, and consequently more accurate results are achieved even withcoarser meshes.
We consider the same settings and reference solution used in Flemisch et al. [38]
Benchmark 4 . For this benchmark,a more realistic fracture network is taken into consideration. Here, the size of the domain is 700 m x 600 m with afracture network of 64 fractures divided in 13 connected groups. The matrix permeability is K = I · − m All thefractures have the same permeability − m and the same aperture (cid:15) = 10 − m. Hence, we have K γ = (cid:15) I · − for the HD case, and K γ = I · − for the ED one. There are no-flow boundary conditions on top and bottom of thedomain. A pressure of Pa is imposed on the left boundary and of 0 Pa on the right boundary. Due to the rather13on-Conforming Mesh Models for Flow in Fractured Porous Media using Lagrange multipliersTable 5:
Benchmark 4 in Flemisch et al. [38]. Mesh information, characteristics of system matrix and error in the matrixfor the
2D Realistic fracture network benchmark in Embedded-ED and Embedded-HD cases.
Method (cid:107) · (cid:107) -cond. Embedded-ED 1405 quads 199 996 1487 1.9e-4 2.0e15Embedded-ED 39 085 quads 199 996 41 773 2.3e-4 2.2e17Embedded-ED 94 513 quads 199 996 102 486 9.1e-5 5.8e17Embedded-HD 4200 quads 1024 4331 3.8e-3 1.3e16Embedded-HD 67 200 quads 4096 67 721 1.6e-4 2.1e17 complex geometry of the fracture network, the mortar method is impractical, hence it is not used here. We exclusivelypresent results for the equi-dimensional embedded method. The mesh of the matrix has been automatically refinedat the intersection with the fracture, for each experiment with 3, 6, and 7 levels of non-conforming mesh refinement,respectively. It can be observed in Table 5 the condition numbers have comparable magnitudes for both HD and EDversions. The condition number becomes large with the extremely varying mesh size resulting from the heavily focusedrefinement around the fractures.As can be observed in Fig. 10, the results obtained are similar to those obtained by other methods in the field although noreference solution is available. It can be noticed that the equi-dimensional variant is closer to the
Box method comparedto most other methods, including the hybrid dimensional results presented in [60]. From a practical perspective theequi-dimensional technique is slightly more complex since the fractures are extruded in normal direction, althoughautomatically, thus requiring many more elements. P r e ss u r e [ M P a ] BoxTPFAEDFMFlux-MortarEmbedded-HDEmbedded-ED P r e ss u r e [ M P a ] BoxTPFAEDFMFlux-MortarEmbedded-HDEmbedded-ED (a) Pressure solution (b) Line AA (cid:48) (c) Line BB (cid:48)
Figure 10:
Benchmark 4 in Flemisch et al. [38]. Pressure solution for real fracture network with 64 fractures (a).Pressure profiles along the lines (b) AA (cid:48) nd (c) BB (cid:48) for both Embedded-HD and Embedded-ED techniques.
Flow through fractured porous media is largely governed by 3D effects. Therefore, this section presents an applicationand evaluation of the dual Lagrange multiplier methods in 3D. First, results obtained with three Lagrange multipliermethods (embedded HD, embedded ED, mortar HD) are compared to results of methods presented in the benchmarkstudy Berre et al. [83]. More complex benchmark cases in 3D are studied for the embedded HD method as part ofthe aforementioned benchmark study. The embedded HD method is preferably used for those geometrically complexcases as it eases meshing of the fracture networks and the porous matrix mesh can be chosen regular. Fracture networkmesh generation for the embedded ED method is very challenging, which is particularly relevant at the fractureintersection. Furthermore, all 3D benchmark cases in Berre et al. [83] are designed for very small fracture aperturesand therefore equi-dimensional fracture meshes result in poor mesh quality. It is also important to bear in mind thatequi-dimensional fracture meshes are mainly necessary for large aperture values and the benchmark cases yield noreason to use equi-dimensional meshes. Furthermore, the mortar method is less suited for complex fracture geometriesdue to the complexity of both, setting up the mesh and dealing with over-constrained scenarios for the mortar conditions.In a final realistic scenario the strength of each method is demonstrated and they are applied in a combined scenario asit could be typical for fractured systems. 14on-Conforming Mesh Models for Flow in Fractured Porous Media using Lagrange multipliers In this section three Lagrange multiplier methods are compared to benchmark
Case 1: Single Fracture presentedby Berre et al. [83], which already includes results obtained with the hybrid embedded method, presented in Schädleet al. [60]. Nevertheless, for completeness and for better comparison of the methods presented in this study the resultsof the hybrid embedded approach are again explicitly presented. To ease comparison, the results presented in thebenchmark study are summarized and only the mean of all results as well as the standard deviation are plotted. Thecomparison of vastly different discretization methods is enabled by interpolating each solution to evenly spacedpoints along a line. Following this, the mean and standard deviation are computed at each of these points. It is importantto note that the mean of all results is not necessarily the correct solution and just provides a measure of comparison tothe results presented in the aforementioned benchmark study. For the detailed benchmark results the interested reader isreferred to Berre et al. [83]. Finally, improvements of the accuracy can be demonstrated for the embedded ED approachby adaptive mesh refinement.Figure 11 shows the model domain of benchmark
Case 1: Single Fracture which is adapted from Zielke et al. [86]and Barlag et al. [87]. This case consists of a single fracture intersecting a matrix cube with
100 m edge length. Thelowest
10 m thick layer of the matrix block (Matrix ) has an increased permeability. Fluid injection occurs at a Dirichletboundary condition (BC ) located above the fracture at the upper most
10 m thick layer of the matrix block. A Dirichletboundary condition (BC ) which acts as the outflow, is located at Matrix . The injection pressure is fixed at andthe production pressure at . The permeability K in the Matrix is I · − [m/s] and in the Matrix I · − [m/s].The fracture has a permeability K γ of I γ · − [m/s] with an aperture (cid:15) of − [m]. Matrix 1 K =1x10 -6 I [m/s]
100 mA A' B C p D = m B C p D = m
100 m10 m10 m80 m10 m10 m10 m70 m Matrix 2 K =1x10 -5 I [m/s] Fracture K =1x10 -3 I [m/s]a=1x10 -2 [m] Figure 11: Model domain (outlines) of benchmark
Case 1: Single Fracture in Berre et al. [83] with Matrix 1 and 2 anda single fracture intersecting Matrix 1. Boundary 1 (BC ) is the inflow and boundary 2 (BC ) the outflow. Results arecompared along the line AA (cid:48) through the matrix domains.The results obtained with the hybrid- and equi-dimensional embedded method and the hybrid mortar method arepresented in Figure 12. In the upper part of the figure, three different mesh sizes ( ∼ ∼ ∼ methods presented in Berre et al. [83]. It is important to keep in mind that the hybrid-dimensionalembedded results are also part of the benchmark results. The pressure solution for all methods is compared along theline AA (cid:48) (see Fig. 11).Overall, the results of the methods presented here and the methods presented in Berre et al. [83] show good convergencetowards a common solution. Particularly for the very coarse mesh of only ∼
1k cells in the matrix domains the Lagrangemultiplier methods show some deviations. These deviations are more pronounced for the embedded HD method at thefirst half of line AA (cid:48) and for the embedded ED and the mortar method at the second half. The deviations are partiallydue to the fact that most of the other methods presented in the benchmark study represent the fracture geometry in thematrix explicitly in matching mesh geometries. Embedded techniques typically require a resolution that is roughly twiceas high as fitted mesh techniques to achieve the same accuracy. Already for the case with ∼
10k cells the embeddedHD and ED methods show very similar results. For the coarsest case neither of the methods is preferable as they allshow deviations in different regions. Furthermore, the solution of the mortar method shows a kink at ∼
115 m for allmesh sizes. As mentioned above, this is due to over-constrained dofs, which poses significant technical challenges to be15on-Conforming Mesh Models for Flow in Fractured Porous Media using Lagrange multipliersautomated in 3D.As shown in 3.1.2 adaptive mesh refinement allows to obtain more accurate results while reducing the number of matrixelements. Here the matrix mesh for the embedded ED approach is adaptively refined with up to two refinement steps,starting at ∼
500 cells by progressively reducing the threshold of error adopted for the gradient recovery strategy. Forbrevity, we refer to the error threshold as err g and to the number of adaptive refinement steps as AR. Thus, we employ err g = 15 and AR = 1 for the coarsest mesh, err g = 5 and AR = 2 for the middle mesh, and err g = 0 . and AR = 2 for the finest test case. One may note that the the error threshold, err g , is progressively reduced to increase the accuracyof the numerical results. The matrix meshes resulting from the adaptive mesh refinement have ∼ ∼ ∼ amr has some more cells than the compared uniform mesh, the pressure solutionalong parts of the line AA (cid:48) is closer to the mean of all benchmark methods. For the two finer meshes, an improvementin accuracy is clearly visible. Taken together, all Lagrange multiplier methods presented here show good convergenceand match well with other methods presented in the benchmark study. P r e ss u r e [ m ]
1k cells
10k cells
Mean allEmbedded HDEmbedded EDMortar HDStd. dev. all0 50 100 150Arc length [m]1234 P r e ss u r e [ m ] Embedded ED
Mean all1k cells Uniform1.3k cells AMRStd. dev. all 0 50 100 150Arc length [m]
Embedded ED
Mean all10k cells Uniform9k cells AMRStd. dev. all 0 50 100 150Arc length [m]
Embedded ED
Mean all100k cells Uniform28k cells AMRStd. dev. all
Figure 12: Pressure results along the line AA (cid:48) for ∼ ∼ ∼ ∼ ∼ ∼
28k for adaptive mesh refinement (lower figure). Results are shown with uniformrefinement for the hybrid-dimensional (HD) and equi-dimensional (ED) embedded method and the hybrid mortarmethod and with adaptive mesh refinement for the equi-dimensional embedded method. The black line shows the meanresults of all benchmark methods and the grey range shows the standard deviation, respectively [83].
Geological setting with fractures and faults often require to represent such features with permeability and apertureranging over several orders of magnitude. Therefore, their combined representation in numerical models is crucial for acomplete description of geological settings. The setup of present realistic scenario is loosely based on the geologicalsetting at the Grimsel Test Site, (GTS) [88]. However, it is important to note that the steady-state flow field, ascomputed here, is difficult to achieve in experiments conducted in such laboratories with a very low permeability rockmatrix. Nevertheless, the given setup allows to demonstrate the different strength of the particular methods describedin this study, i.e. complex fracture networks, fractures with large aperture widths, and blocking fractures. Even more16on-Conforming Mesh Models for Flow in Fractured Porous Media using Lagrange multipliersimportantly, this realistic scenario shows the integration of several coupling strategies in a single joint frameworkmethod.Figure 13 shows the model domain with a complex fracture network (blue) located between two large features withtwo intersecting fractures each. One of these features acts as blocking fractures (yellow) with low permeability values,the fractures in the other feature have large apertures (magenta). Furthermore, two boreholes are drilled into the rockdomain with one of them ending in the large aperture fractures and the other one in the rock domain. The upper rightcorner of the square domain is not modeled as it acts as an access tunnel with atmospheric pressure. Consequently, thepressure at the outer boundary is fixed by a Dirichlet boundary condition of p D = ( x − y + 100) 0 . (green). Thisboundary condition results in a pressure of at the access tunnel location and at the lower right corner. Atthe top and bottom ( z -direction) of the domain no-flow boundary conditions are applied.In the numerical model the fractures and fracture network are represented by three different methods described above.The complex fracture network (blue) is meshed with lower-dimensional manifolds and embedded in the matrix meshand has been initially presented by Schädle et al. [60]. In this study the fracture aperture and permeability and thematrix permeability are chosen so that the fractures and the matrix contribute similarly to the overall flow. In contrast,the rock matrix in the present experiment holds a low permeability and the fracture properties are chosen to dominatethe flow. The fracture radius distribution in the network follows a power law, with truncations at . and
10 m .The fractures are circular, randomly oriented, and distributed in a cubical area of the model domain with a sidelength of
25 m . The fracture aperture is ∼ − m and the permeability chosen to be ∼ − m . Furthermore, thetwo fractures with large apertures (magenta) are represented by equi-dimensional domains embedded in the matrixdomain. One of these fractures is intersected by a borehole which is explicitly meshed as a sub-domain. The twofractures are circular with an aperture width of × − m and radii of
25 m . Further, the infilling material of thesefractures is assumed to have a permeability of × − m . In the borehole domain a forcing function of . × − is applied in a volume of × − m , acting as an injection borehole with a pressure of ∼ . . To representfractures with low permeability, acting as blocking fractures, two fractures (yellow) are described by equi-dimensionaldomains coupled to the matrix mesh by the mortar method. These fractures are rectangular with a side length of
35 m and permeability of × − m . Generally, the shape of the fractures might of any shape, e.g. circular or rectangular.The second borehole acts as a sink with a fixed atmospheric pressure. This borehole is described by a mortar couplingof a small domain with the dimensions of the borehole and a fixed pressure of . Finally, the matrix permeabilityis × − m . p D = ( x - y + ) . xy z Figure 13: Model domain of the complex case with a combined application of all presented methods. The upper rightcorner represents a ventilated access tunnel, as it is typical in deep underground laboratories. Fractures computed withthe hybrid embedded method are depicted in blue and fractures computed with the equi-dimensional embedded methodin pink. The yellow domains, with low permeability, are equi-dimensional and coupled by the mortar method. Thegreen planes show the outer boundaries and the black lines indicate the injection and production boreholes.Figure 14 shows the pressure distribution across three planes intersecting the model domain parallel to the xy -plane.Throughout all planes the low permeability fractures act as discontinuities for pressure. Moreover, for plane (a) aclear discontinuity can be observed across the two fracture cross sections. With the high aperture and permeability ofthe complex fracture network and the low permeability of the rock matrix the pressure across the fracture network is17on-Conforming Mesh Models for Flow in Fractured Porous Media using Lagrange multipliersequilibrated, connecting the upper right with the lower left part of the domain. With the injection borehole in one of thetwo equi-dimensional fractures, these fractures and the surrounding rock matrix area subject to the largest pressurevalues. In plane (b) the production borehole locally reduces the pressure to atmospheric pressure. However, due to thelow permeability of the rock matrix the gradient around this borehole is very steep and the influence on the overallsolution is limited. Additionally, the fixed pressure boundary condition at the outer boundary forces the pressure to steepgradients close to the embedded ED fractures with high injection pressure. These steep gradients result from boundaryeffects and for a representative study of such a geological setting the domain would have to be extended in theseareas. However, the goal of this study is to demonstrate the application of Lagrange multipliers for different couplingstrategies, spanning from equi-dimensional to lower dimensional models, and from embedded to mortar techniques.Ultimately allowing for more flexibility in the treatment of the fracture and matrix configuration. Furthermore, suchsteep gradients are generally difficult to resolve in numerical models, thus adding further complexity to this test case.In summary, this realistic scenario demonstrates the strength of each method and the ability of the presented unifiedframework to combine all of these approaches while yielding smooth pressure results in geologically complex settings.More specifically, the hybrid-dimensional embedded approach eases meshing of complex fracture networks of fractureswith small aperture widths, the equi-dimensional embedded approach allows to consider fractures with large aperturewidths, and the equi-dimensional mortar approach enables to consider blocking fractures. P r e ss u r e [ M P a ] (a) Plane at z = − (b) Plane at z = 0 m (c) Plane at z = 5 m Figure 14: Pressure solution of realistic scenario case at three planes intersecting the domain at − , , and in z -direction. This study expands on previous works based on the application of the Lagrange multiplier method to computesingle-phase fluid flow problems in fracture dominated porous media. In particular, we employ the finite elementmethod in combination with an L -projection operator to couple different types of non-conforming meshes. Thenon-conformity might arise at the interface of independently meshed sub-domains, at hanging nodes resulting fromadaptive mesh refinement, or by combining multiple overlapping meshes. Furthermore, fractures are either describedby equi-dimensional or hybrid-dimensional domains. The applied Lagrange multiplier is discretized using dual basisfunctions, which provide two main advantages. First, the number of degrees of freedom is reduced to the ones of thebackground mesh representing the porous-matrix. Second, the arising symmetric-positive-definite linear systems areconvenient work with.Here, we present a unified framework covering all coupling techniques mentioned above. The different mesh solutionsare compared with state-of-the-art benchmark cases in 2D and 3D, the numerical performance is studied, and usecases are presented in isolation and combination. Overall, the results suggest that the presented tool-set is capableof computing fluid-flow through complex and heterogeneous rock formations in a robust and convenient way. It isimportant to bear in mind that realistic scenarios of fracture dominated rock formations may include fractures withgeometric and physical properties ranging over many orders of magnitude. Therefore, the presented single jointframework allows to deeply exploit non-conforming hybrid- and equi-dimensional fracture models and efficientlycombine these models. By using the dual Lagrange multiplier we are able to combine multiple complex fracturenetworks without changing the size of the algebraic system arising from the porous-medium matrix, although we havemore non-zero entries associated with the coupled degrees of freedom. It is also worth to point out that the conditioningof the system is not worsened, as it is the case when using other types of Lagrange multipliers [60].18on-Conforming Mesh Models for Flow in Fractured Porous Media using Lagrange multipliersWith the integration of adaptive mesh refinement in the solution process, the error in the solution can be controlled in anautomated way either by means of an error estimator or by pre-defining areas of interest for refinement. This allowed tocomplement the discussion in Schädle et al. [60] about the necessity of having a finer mesh around fractures and theirtips and intersections.The study of a realistic geological setting, inspired by the Grimsel Test Site, demonstrates the advantage of theunified framework to represent vastly different fracture geometries and properties within a single numerical model.This allows to expand on existing numerical studies of such systems by the opportunity to include a large range offracture representations. Furthermore, the highly efficient and convenient tools combined with the eased meshing ofnon-conforming meshes enables stochastic studies for a wide range of fracture dominated systems.Further investigations based on this unified framework would benefit by focusing on mass conservation properties ofthe finite element discretization and application to transport problems. Acknowledgment
P.Z., M.G.C.N., and L.K. thank the SCCER-SoE and SCCER-FURIES programs, and the PASC project FASTER:Forecasting and Assessing Seismicity and Thermal Evolution in geothermal Reservoirs. P.S. thanks the Werner SiemensFoundation for their endowment of the Geothermal Energy and Geofluids group at the Institute of Geophysics, ETHZurich.
Authorship statement
P.Z. lead the drafting of the manuscript, implemented most methods and numerical tools, developed parts of theconceptual models, and lead parts of the numerical experiments. P.S. contributed drafting the manuscript, lead thedevelopment of the 3D conceptual models, their validation, and presentation of the results. L.K. produced, collected,and prepared most of the 2D benchmark results. M.G.C.N contributed drafting the manuscript, implemented theadaptive refinement strategy and its integration with the variational transfer, lead the numerical experiments related tothe adaptive refinement, and contributed to the 2D and 3D numerical experiments.
Conflict of interest
The authors declare that they have no conflict of interest.
Computer Code Availability
All methods and routines, used for this study, are implemented with the open-source software library
Utopia [79].
Utopia ’s lead developer is author Patrick Zulian at USI Lugano, Switzerland. The contact address and e-mail of PatrickZulian are as follows:Institute of Computational ScienceUniversità della Svizzera italiana (USI - University of Lugano)Via Giuseppe Buffi 13CH-6900 [email protected]
Utopia was first available in 2016, the programming language is
C++ and it can be accessed through a git repository ora docker container on: https://bitbucket.org/zulianp/utopia (approx. 40 MB of uncompressed data), https://hub.docker.com/r/utopiadev/utopia .The software dependencies are as follows:PETSc ( ),must be compiled with
MUMPS enabledlibMesh for the FE module ( https://github.com/libMesh )There are no hardware requirements given by
Utopia . Potential hardware or software requirements of the underlyinglibraries libMesh and
PETSc are not stated here. 19on-Conforming Mesh Models for Flow in Fractured Porous Media using Lagrange multipliers
References [1] Alain Pochon, Jean-Pierre Tripet, Ronald Kozel, Benjamin Meylan, Michael Sinreich, and François Zwahlen.Groundwater protection in fractured media: a vulnerability-based approach for delineating protection zones inswitzerland.
Hydrogeology Journal , 16(7), 2008. doi: 10.1007/s10040-008-0323-0.[2] T. Read, O. Bour, V. Bense, T. Le Borgne, P. Goderniaux, M.V. Klepikova, R. Hochreutener, N. Lavenant, andV. Boschero. Characterizing groundwater flow and heat transport in fractured rock using fiber-optic distributedtemperature sensing.
Geophysical Research Letters , 40(10):2055–2059, 2013. doi: 10.1002/grl.50397.[3] Jefferson W Tester, B Anderson, A Batchelor, D Blackwell, Ronald DiPippo, E Drake, John Garnish, B Livesay,Michal C Moore, Kenneth Nichols, et al. The future of geothermal energy: Impact of enhanced geothermalsystems (egs) on the united states in the 21st century.
Massachusetts Institute of Technology , 209, 2006.[4] Mark W. McClure and Roland N. Horne. Correlations between formation properties and induced seismicityduring high pressure injection into granitic rock.
Engineering Geology , 175:74–80, 2014. ISSN 0013-7952. doi:http://dx.doi.org/10.1016/j.enggeo.2014.03.015.[5] Clare E. Bond, Ruth Wightman, and Philip S. Ringrose. The influence of fracture anisotropy on co2 flow.
Geophysical Research Letters , 40(7):1284–1289, 2003. doi: 10.1002/grl.50313. URL https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1002/grl.50313 .[6] E. Bonnet, O. Bour, N. E. Odling, P. Davy, I. Main, P. Cowie, and B. Berkowitz. Scaling of fracture systems ingeological media.
Reviews of Geophysics , 39(3):347–383, 8 2001. ISSN 8755-1209. doi: 10.1029/1999RG000074.URL http:https://doi.org/10.1029/1999RG000074 .[7] Anders Rasmuson and Ivars Neretnieks. Radionuclide transport in fast channels in crystalline rock.
WaterResources Research , 22(8):1247–1256, 1986.[8] Florian Amann, Valentin Gischig, Keith Evans, Joseph Doetsch, Reza Jalali, Benoît Valley, Hannes Krietsch,Nathan Dutler, Linus Villiger, Bernard Brixel, et al. The seismo-hydromechanical behavior during deep geothermalreservoir stimulations: open questions tackled in a decameter-scale in situ stimulation experiment.
Solid Earth , 9(1):115–137, 2018.[9] Jean-Raynald de Dreuzy, Yves Méheust, and Géraldine Pichot. Influence of fracture scale heterogeneity on theflow properties of three-dimensional discrete fracture networks (dfn).
Journal of Geophysical Research: SolidEarth , 117(B11), 2012.[10] R.W. Zimmerman, S. Kumar, and G.S. Bodvarsson. Lubrication theory analysis of the permeability of rough-walled fractures.
International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts ,28(4):325–331, 1991. ISSN 0148-9062. doi: https://doi.org/10.1016/0148-9062(91)90597-F. URL .[11] Anozie Ebigbo, Philipp S Lang, Adriana Paluszny, and Robert W Zimmerman. Inclusion-based effective mediummodels for the permeability of a 3d fractured rock mass.
Transport in Porous Media , 113(1):137–158, 2016.[12] Daniel Vogler, Florian Amann, Peter Bayer, and Derek Elsworth. Permeability evolution in natural fracturessubject to cyclic loading and gouge formation.
Rock Mechanics and Rock Engineering , 49:3463–3479, 2016. doi:https://doi.org/10.1007/s00603-016-1022-0.[13] Shlomo P. Neuman.
Stochastic approach to subsurface flow and transport: a view to the future , pages 231–241.International Hydrology Series. Cambridge University Press, 1997. doi: 10.1017/CBO9780511600081.016.[14] Brian Berkowitz. Characterizing flow and transport in fractured geological media: A review.
Advances in WaterResources , 25(8):861–884, 2002. doi: 10.1016/S0309-1708(02)00042-8.[15] M. C. Cacas, E. Ledoux, G. Marsily, B. Tillie, A. Barbreau, E. Durand, B. Feuga, and P. Peaudecerf. Modelingfracture flow with a stochastic discrete fracture network: calibration and validation: 1. the flow model.
WaterResources Research , 26(3):479–489, 3 1990. ISSN 0043-1397. doi: 10.1029/WR026i003p00479. URL https://doi.org/10.1029/WR026i003p00479 .[16] Alex Hobé, Daniel Vogler, Martin P. Seybold, Anozie Ebigbo, Randolph R. Settgast, and Martin O. Saar.Estimating fluid flow rates through fracture networks using combinatorial optimization.
Advances in WaterResources , 122:85–97, 2018. ISSN 0309-1708. doi: https://doi.org/10.1016/j.advwatres.2018.10.002. URL .[17] Shlomo P Neuman. Trends, prospects and challenges in quantifying flow and transport through fractured rocks.
Hydrogeology Journal , 13(1):124–147, 2005.[18] Jean-Raynald de Dreuzy, Géraldine Pichot, Baptiste Poirriez, and Jocelyne Erhel. Synthetic benchmark formodeling flow in 3d fractured media.
Computers & Geosciences , 50:59–71, 2013.20on-Conforming Mesh Models for Flow in Fractured Porous Media using Lagrange multipliers[19] Benoît Dessirier, Chin-Fu Tsang, and Auli Niemi. A new scripting library for modeling flow and transport infractured rock with channel networks.
Computers & Geosciences , 111:181–189, 2018.[20] J.E. Warren and P.J. Root. The Behavior of Naturally Fractured Reservoirs.
Society of Petroleum Engi-neers Journal , 3(3):245–255, sep 1963. ISSN 0197-7520. doi: 10.2118/426-PA. URL .[21] GI Barenblatt, Iu P Zheltov, and IN Kochina. Basic concepts in the theory of seepage of homogeneous liquids infissured rocks strata.
Journal of applied mathematics and mechanics , 24(5):1286–1303, 1960.[22] Hossein Kazemi. Pressure transient analysis of naturally fractured reservoirs with uniform fracture distribution.
Society of petroleum engineers Journal , 9(04):451–462, 1969.[23] H Kazemi, LS Merrill Jr, KL Porterfield, PR Zeman, et al. Numerical simulation of water-oil flow in naturallyfractured reservoirs.
Society of Petroleum Engineers Journal , 16(06):317–326, 1976.[24] Jahan Noorishad and Mohsen Mehran. An upstream finite element method for solution of transient transportequation in fractured porous media.
Water Resources Research , 18(3):588–596, 6 1982. ISSN 0043-1397. doi:10.1029/WR018i003p00588. URL https://doi.org/10.1029/WR018i003p00588 .[25] R. G. Baca, R. C. Arnett, and D. W. Langford. Modelling fluid flow in fractured-porous rock masses by finite-element techniques.
International Journal for Numerical Methods in Fluids , 4(4):337–348, 4 1984. ISSN1097-0363. doi: 10.1002/fld.1650040404. URL https://doi.org/10.1002/fld.1650040404 .[26] Jeffrey D Hyman, Satish Karra, Nataliia Makedonska, Carl W Gable, Scott L Painter, and Hari S Viswanathan.dfnworks: A discrete fracture network framework for modeling subsurface flow and transport.
Computers &Geosciences , 84:10–19, 2015.[27] Bernd Flemisch, Melanie Darcis, K Erbertseder, B Faigle, A Lauser, Klaus Mosthaf, S Müthing, Philipp Nuske,A Tatomir, M Wolff, and Helmig Rainer. Dumux: Dune for multi- { phase, component, scale, physics,. . . } flowand transport in porous media. Advances in Water Resources , 34(9):1102–1112, 2011. doi: 10.1016/j.advwatres.2011.03.007.[28] Konstantin Lipnikov, Gianmarco Manzini, and Mikhail Shashkov. Mimetic finite difference method.
Journal ofComputational Physics , 257:1163–1227, 2014. doi: 10.1016/j.jcp.2013.07.031.[29] I-Hsien Lee and Chuen-Fa Ni. Fracture-based modeling of complex flow and co2 migration in three-dimensionalfractured rocks.
Computers & geosciences , 81:64–77, 2015.[30] I-Hsien Lee, Chuen-Fa Ni, Fang-Pang Lin, Chi-Ping Lin, and Chien-Chung Ke. Stochastic modeling of flow andconservative transport in three-dimensional discrete fracture networks.
Hydrology and Earth System Sciences , 23(1):19–34, 2019. doi: 10.5194/hess-23-19-2019.[31] Mauro Cacace and Guido Blöcher. Meshit—a software for three dimensional volumetric meshing of complexfaulted reservoirs.
Environmental Earth Sciences , 74(6):5191–5209, 2015.[32] Randi Holm, Roland Kaufmann, Bjørn-Ove Heimsund, Erlend Øian, and Magne S Espedal. Meshing of domainswith complex internal geometries.
Numerical Linear Algebra with Applications , 13(9):717–731, 2006.[33] Daniela Blessent, René Therrien, and Kerry MacQuarrie. Coupling geological and numerical models to simulategroundwater flow and contaminant transport in fractured media.
Computers & Geosciences , 35(9):1897–1906,2009.[34] M. Karimi-Fard, L.J. Durlofsky, and K. Aziz. An efficient discrete-fracture model applicable for general-purposereservoir simulators.
SPE Journal , 9(2):227–236, 2004.[35] II Bogdanov, VV Mourzenko, J-F Thovert, and PM Adler. Two-phase flow through fractured porous media.
Physical Review E , 68(2):026703, 2003.[36] JEP Monteagudo and Abbas Firoozabadi. Control-volume method for numerical simulation of two-phaseimmiscible flow in two-and three-dimensional discrete-fractured media.
Water resources research , 40(7), 2004.[37] Rainer Helmig et al.
Multiphase flow and transport processes in the subsurface: a contribution to the modeling ofhydrosystems.
Springer-Verlag, 1997.[38] Bernd Flemisch, Inga Berre, Wietse Boon, Alessio Fumagalli, Nicolas Schwenck, Anna Scotti, Ivar Stefansson,and Alexandru Tatomir. Benchmarks for single-phase flow in fractured porous media.
Advances in WaterResources , 111:239–258, 2018. ISSN 0309-1708. doi: https://doi.org/10.1016/j.advwatres.2017.10.036. URL .21on-Conforming Mesh Models for Flow in Fractured Porous Media using Lagrange multipliers[39] Eirik Keilegavlen, Runar Berge, Alessio Fumagalli, Michele Starnoni, Ivar Stefansson, Jhabriel Varela, and IngaBerre. Porepy: An open-source software for simulation of multiphysics processes in fractured porous media. arXiv preprint arXiv:1908.09869 , 2019.[40] Jan Martin Nordbotten, Wietse M Boon, Alessio Fumagalli, and Eirik Keilegavlen. Unified approach todiscretization of flow in fractured porous media.
Computational Geosciences , 23(2):225–237, 2019. doi:10.1007/s10596-018-9778-9.[41] Philippe Devloo, Wenchao Teng, and Chen-Song Zhang. Multiscale hybrid-mixed finite element method for flowsimulation in fractured porous media.
Computer Modeling in Engineering & Sciences , 119(1):145–163, 2019.doi: 10.32604/cmes.2019.04812.[42] Omar Duran, Philippe R B Devloo, Sônia M Gomes, and Frédéric Valentin. A multiscale hybrid method for darcy’sproblems using mixed finite element local solvers.
Computer Methods in Applied Mechanics and Engineering ,354:213–244, 2019. doi: 10.1016/j.cma.2019.05.013.[43] Konstantin Brenner, Mayya Groza, Cindy Guichard, Gilles Lebeau, and Roland Masson. Gradient discretizationof hybrid dimensional darcy flows in fractured porous media.
Numerische Mathematik , 134(3):569–609, 2016.doi: 10.1007/s00211-015-0782-x.[44] Konstantin Brenner, Julian Hennicker, Roland Masson, and Pierre Samier. Gradient discretization of hybrid-dimensional darcy flow in fractured porous media with discontinuous pressures at matrix–fracture interfaces.
IMAJournal of Numerical Analysis , 37(3):1551–1585, 2016. doi: 10.1093/imanum/drw044.[45] Chiara Facciolà, Paola Francesca Antonietti, and Marco Verani. Mixed-primal discontinuous galerkin approx-imation of flows in fractured porous media on polygonal and polyhedral grids.
PAMM , 19(1):e201900117,2019.[46] Liyong Li and Seong H Lee. Efficient field-scale simulation of black oil in a naturally fractured reservoir throughdiscrete fracture networks and homogenized media.
SPE Reservoir Evaluation & Engineering , 11(04):750–758,2008.[47] Hadi Hajibeygi, Dimitrios C. Karvounis, and Patrick Jenny. A hierarchical fracture model for the iterativemultiscale finite volume method.
Journal of Computational Physics , 230(24):8729–8743, 2011. ISSN 00219991.doi: 10.1016/j.jcp.2011.08.021.[48] Matei ¸Tene, Sebastian BM Bosma, Mohammed Saad Al Kobaisi, and Hadi Hajibeygi. Projection-based embeddeddiscrete fracture model (pedfm).
Advances in Water Resources , 105:205–216, 2017.[49] Ali Moinfar, Abdoljalil Varavei, Kamy Sepehrnoori, and Russell T. Johns. Development of an Efficient EmbeddedDiscrete Fracture Model for 3D Compositional Reservoir Simulation in Fractured Reservoirs.
SPE Journal , 19(02):289–303, apr 2014. ISSN 1086-055X. doi: 10.2118/154246-PA. URL .[50] Kirill D Nikitin and Ruslan M Yanbarisov. Monotone embedded discrete fractures method for flows in porousmedia.
Journal of Computational and Applied Mathematics , 364:112353, 2020. doi: 10.1016/j.cam.2019.112353.[51] Bernd Flemisch, Alessio Fumagalli, and Anna Scotti. A Review of the XFEM-Based Approximation ofFlow in Fractured Porous Media. In
Advances in Discretization Methods , pages 47–76. Springer Interna-tional Publishing, 2016. doi: 10.1007/978-3-319-41246-7_3. URL http://link.springer.com/10.1007/978-3-319-41246-7{_}3 .[52] Ziyao Xu and Yang Yang. The hybrid dimensional representation of permeability tensor: A reinterpretation ofthe discrete fracture model and its extension on nonconforming meshes.
Journal of Computational Physics , 415:109523, 2020. doi: https://doi.org/10.1016/j.jcp.2020.109523.[53] D Capatina, R Luce, H El-Otmany, and N Barrau. Nitsche’s extended finite element method for a fracture modelin porous media.
Applicable Analysis , 95(10):2224–2242, 2016.[54] Hao Huang, Ted A Long, Jing Wan, and William P Brown. On the use of enriched finite element method to modelsubsurface features in porous media flow problems.
Computational Geosciences , 15(4):721–736, 2011.[55] Nicolas Schwenck, Bernd Flemisch, Rainer Helmig, and Barbara I. Wohlmuth. Dimensionally reduced flowmodels in fractured porous media: crossings and boundaries.
Computational Geosciences , 19(6):1219–1230,dec 2015. ISSN 1420-0597. doi: 10.1007/s10596-015-9536-1. URL http://link.springer.com/10.1007/s10596-015-9536-1 .[56] Carlo D’Angelo and Anna Scotti. A mixed finite element method for darcy flow in fractured porous media withnon-matching grids.
ESAIM: Mathematical Modelling and Numerical Analysis , 46(2):465–489, 2012.22on-Conforming Mesh Models for Flow in Fractured Porous Media using Lagrange multipliers[57] Alessio Fumagalli and Anna Scotti. A numerical method for two-phase flow in fractured porous media withnon-matching grids.
Advances in Water Resources , 62:454–464, 2013.[58] Alessio Fumagalli, Eirik Keilegavlen, and Stefano Scialò. Conforming, non-conforming and non-matchingdiscretization couplings in discrete fracture network simulations.
Journal of Computational Physics , 376:694–712,2019.[59] M. Köppel, V. Martin, J. Jaffré, and Jean E Roberts. A lagrange multiplier method for a discrete fracture modelfor flow in porous media.
Computational Geosciences , 23(2):239–253, 2018.[60] Philipp Schädle, Patrick Zulian, Daniel Vogler, Sthavishtha R. Bhopalam, Maria G.C. Nestola, Anozie Ebigbo,Rolf Krause, and Martin O. Saar. 3d non-conforming mesh model for flow in fractured porous media using lagrangemultipliers.
Computers & Geosciences , 132:42 – 55, 2019. ISSN 0098-3004. doi: https://doi.org/10.1016/j.cageo.2019.06.014. URL .[61] Inga Berre, Florian Doster, and Eirik Keilegavlen. Flow in fractured porous media: A review of conceptual modelsand discretization approaches.
Transport in Porous Media , pages 1–22, 2018.[62] Maria Giuseppina Chiara Nestola, Barna Becsek, Hadi Zolfaghari, Patrick Zulian, Dario De Marinis, Rolf Krause,and Dominik Obrist. An immersed boundary method for fluid-structure interaction based on variational transfer.
Journal of Computational Physics , 398:108884, 2019. ISSN 0021-9991. doi: https://doi.org/10.1016/j.jcp.2019.108884. URL .[63] Sarah Osborn, Patrick Zulian, Thomas Benson, Umberto Villa, Rolf Krause, and Panayot S. Vassilevski. Scalablehierarchical pde sampler for generating spatially correlated random fields using nonmatching meshes.
NumericalLinear Algebra with Applications , 25(3):e2146, 2018. doi: 10.1002/nla.2146. URL https://onlinelibrary.wiley.com/doi/abs/10.1002/nla.2146 .[64] Roland Glowinski, Tsorng-Whay Pan, and Jacques Periaux. A fictitious domain method for dirichlet problem andapplications.
Computer Methods in Applied Mechanics and Engineering , 111(3-4):283–303, 1994.[65] Christine Bernardi, Yvon Maday, and Anthony T Patera. Domain decomposition by the mortar element method.In
Asymptotic and numerical methods for partial differential equations with critical parameters , pages 269–286.Springer, 1993.[66] Alexander Popp, Barbara I Wohlmuth, Michael W Gee, and Wolfgang A Wall. Dual quadratic mortar finiteelement methods for 3d finite deformation contact.
SIAM Journal on Scientific Computing , 34(4):B421–B446,2012.[67] Cyrill Von Planta, Daniel Vogler, Xiaoqing Chen, Maria GC Nestola, Martin O Saar, and Rolf Krause. Simulationof hydro-mechanically coupled processes in rough rock fractures using an immersed boundary method andvariational transfer operators.
Computational Geosciences , 23(5):1125–1140, 2019.[68] Najla Frih, Vincent Martin, Jean Elizabeth Roberts, and Ali Saâda. Modeling fractures as interfaces withnonmatching grids.
Computational Geosciences , 16(4):1043–1060, sep 2012. doi: 10.1007/s10596-012-9302-6.URL http://link.springer.com/10.1007/s10596-012-9302-6 .[69] Wietse M Boon, Jan M Nordbotten, and Ivan Yotov. Robust discretization of flow in fractured porous media.
SIAM Journal on Numerical Analysis , 56(4):2203–2233, 2018.[70] C. Hesch, A.J. Gil, A. Arranz Carre no, J. Bonet, and P. Betsch. A mortar approach for fluid–structure interactionproblems: Immersed strategies for deformable and rigid bodies.
Computer Methods in Applied Mechanics andEngineering , 278:853 – 882, 2014. ISSN 0045-7825. doi: https://doi.org/10.1016/j.cma.2014.06.004. URL .[71] Barbara I Wohlmuth. A mortar finite element method using dual spaces for the lagrange multiplier.
SIAM journalon numerical analysis , 38(3):989–1012, 2000.[72] William L Briggs and Steve F McCormick.
A multigrid tutorial , volume 72. Siam, 2000.[73] Bishnu P Lamichhane, Rob P Stevenson, and Barbara I Wohlmuth. Higher order mortar finite element methods in3d with dual lagrange multiplier bases.
Numerische Mathematik , 102(1):93–121, 2005.[74] Mario Paz and William Leigh. Static condensation and substructuring. In
Integrated Matrix Analysis of Structures ,pages 239–260. Springer, 2001.[75] Dorian Krause, Thomas Dickopf, Mark Potse, and Rolf Krause. Towards a large-scale scalable adaptive heartmodel using shallow tree meshes.
Journal of Computational Physics , 298:79 – 94, 2015. ISSN 0021-9991. doi:https://doi.org/10.1016/j.jcp.2015.05.005. URL . 23on-Conforming Mesh Models for Flow in Fractured Porous Media using Lagrange multipliers[76] P. Farah, W. A. Wall, and A. Popp. A mortar finite element approach for point, line, and surface contact.
International Journal for Numerical Methods in Engineering , 114(3):255–291, 2018. doi: 10.1002/nme.5743.URL https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.5743 .[77] Ivan E Sutherland and Gary W Hodgman. Reentrant polygon clipping.
Communications of the ACM , 17(1):32–42,1974.[78] Jakub ˇCervený, Veselin Dobrev, and Tzanio Kolev. Non-conforming mesh refinement for high-order finiteelements, 2019.[79] Patrick Zulian, Alena Kopaniˇcáková, Maria Chiara Giuseppina Nestola, Andreas Fink, Nur Fadel, Victor Magri,Teseo Schneider, and Eric Botter. Utopia: A C++ embedded domain specific language for scientific computing.Git repository. https://bitbucket.org/zulianp/utopia, 2016. URL https://bitbucket.org/zulianp/utopia .[80] B. S. Kirk, J. W. Peterson, R. H. Stogner, and G. F. Carey. libMesh : A C++ Library for Parallel Adaptive MeshRefinement/Coarsening Simulations.
Engineering with Computers , 22(3–4):237–254, 2006.[81] Patrick Zulian. ParMOONoLith: parallel intersection detection and automatic load-balancing library. Gitrepository. https://bitbucket.org/zulianp/par_moonolith, 2016. URL https://bitbucket.org/zulianp/par_moonolith .[82] Satish Balay, William D. Gropp, Lois Curfman McInnes, and Barry F. Smith. Efficient management of parallelismin object oriented numerical software libraries. In E. Arge, A. M. Bruaset, and H. P. Langtangen, editors,
ModernSoftware Tools in Scientific Computing , pages 163–202. Birkhäuser Press, 1997.[83] Inga Berre, Wietse M Boon, Bernd Flemisch, Allessio Fumagalli, Dennis Gläser, Eirik Keilegavlen, AnnaScotti, Ivar Stefansson, Alexandru Tatomir, Konstantin Brenner, Samuel Burbulla, Philippe Devloo, Omar Duran,Marco Favino, Julian Hennicker, I-Hsien Lee, Konstantin Lipnikov, Roland Masson, Klaus Mosthaf, MariaGiuseppina Chiara Nestola, Chuen-Fa Ni, Kirill Nikitin, Philipp Schädle, Daniil Svyatskiy, Ruslan Yanbarisov,and Patrick Zulian. Verification benchmarks for single-phase flow in three-dimensional fractured porous media. arXiv preprint arXiv:2002.07005 , 2020.[84] Ningning Yan and Aihui Zhou. Gradient recovery type a posteriori error estimates for finite element approximationson irregular meshes.
Computer methods in applied mechanics and engineering , 190(32-33):4289–4299, 2001.[85] Thomas J Nicholson, TJ McCartin, Paul A Davis, and Walt Beyeler. Nrc experiences in hydrocoin: An internationalproject for studying ground-water flow modeling strategies. In
Waste management’87: Waste isolation in the US,technical programs, and public education . 1987.[86] W Zielke, R Helmig, K P Krohn, H Shao, and J Wollrath. Discrete modelling of transport processes in fracturedporous rock. In . International Society for Rock Mechanics and Rock Engineering, 1991.[87] C Barlag, R Hinkelmann, R Helmig, and W Zielke. Adaptive methods for modelling transport processes infractured subsurface systems. In , volume 284, 1998.[88] F. Amann, V. Gischig, K. Evans, J. Doetsch, R. Jalali, B. Valley, H. Krietsch, N. Dutler, L. Villiger, B. Brixel,M. Klepikova, A. Kittilä, C. Madonna, S. Wiemer, M. O. Saar, S. Loew, T. Driesner, H. Maurer, and D. Giardini.The seismo-hydromechanical behavior during deep geothermal reservoir stimulations: open questions tackled in adecameter-scale in situ stimulation experiment.