3D non-conforming mesh model for flow in fractured porous media using Lagrange multipliers

TL;DR

This paper introduces a novel 3D modeling approach for simulating flow in fractured porous media using non-conforming meshes, combining Lagrange multipliers with an $L^2$-projection transfer for accurate, parallel variable coupling.

Contribution

It presents a new method integrating Lagrange multipliers and $L^2$-projection for efficient, accurate simulation of flow in complex 3D fracture networks with non-conforming meshes.

Findings

Good agreement with 2D benchmarks

Validated on 3D fracture networks against conforming mesh results

Capable of modeling realistic 3D fracture networks

Abstract

This work presents a modeling approach for single-phase flow in 3D fractured porous media with non-conforming meshes. To this end, a Lagrange multiplier method is combined with a parallel $L^2$-projection variational transfer approach. This Lagrange multiplier method enables the use of non-conforming meshes and depicts the variable coupling between fracture and matrix domain. The $L^2$-projection variational transfer allows general, accurate, and parallel projection of variables between non-conforming meshes (i.e. between fracture and matrix domain). Comparisons of simulations with 2D benchmarks show good agreement, and the method is further validated on 3D fracture networks by comparing it to results from conforming mesh simulations which were used as a reference. Application to realistic fracture networks with hundreds of fractures is demonstrated. Mesh size and mesh convergence are investigated for benchmark cases and 3D fracture network applications. Results demonstrate that the Lagrange multiplier method, in combination with the $L^2$-projection method, is capable of modeling single-phase flow through realistic 3D fracture networks.

Equations

$$\nabla \cdot (-\mathbf{K}\nabla p)-\lambda =f$$

$$p=0$$

$${\nabla }_{\gamma }\cdot (-{\mathbf{K}}_{\gamma }{\nabla }_{\gamma }{p}_{\gamma })+\lambda =f_{\gamma }$$

$$p_{\gamma }=0$$

$$V_{\Omega }=H_0^1(\Omega ),$$

$$\mathbf{V}=V_{\Omega }\times V_{\gamma },$$

$${\int }_{\Omega }\mathbf{K}\nabla p\cdot \nabla q+{\int }_{\gamma }{\mathbf{K}}_{\gamma }{\nabla }_{\gamma }{p}_{\gamma }\cdot {\nabla }_{\gamma }{q}_{\gamma }-$$

$${\int }_{\gamma }\lambda (q-q_{\gamma })={\int }_{\Omega }fq+{\int }_{\gamma }{f}_{\gamma }{q}_{\gamma },$$

$${\int }_{\gamma }(p-p_{\gamma })\mu =0,$$

$$V_{h,\Omega }=$$

$$q{|}_M\in \left\{\begin{array}{rl}{\mathbb{P}}^2(M)& \text{if }M\text{ is: a triangle, pyramid, tetrahedron, or prism}\\ {\mathbb{Q}}^{2,2}(M)& \text{if }M\text{ is: a quadrilateral or hexahedron}\end{array}\right\},$$

$$V_{h,\gamma }=$$

$$q_{\gamma }{|}_m\in \left\{\begin{array}{rl}{\mathbb{P}}^1(m)& \text{if }m\text{ is: a triangle or line segment}\\ {\mathbb{Q}}^{1,1}(m)& \text{if }m\text{ is: a quadrilateral}\end{array}\right\},$$

$${\Lambda }_h=$$

$${\text{err}}_n^2=$$

$${\text{RMS}}_{\mathcal{M},(\gamma )}=$$

Full text

Comparisons of simulations with 2D benchmarks show good agreement, and the method is further validated on 3D fracture networks by comparing it to results from conforming mesh simulations which were used as a reference. Application to realistic fracture networks with hundreds of fractures is demonstrated. Mesh size and mesh convergence are investigated for benchmark cases and 3D fracture network applications. Results demonstrate that the Lagrange multiplier method, in combination with the $L^{2}$-projection method, is capable of modeling single-phase flow through realistic 3D fracture networks.

Keywords: Embedded discrete fracture model, Flow in 3D fractured porous media, Finite element method, Non-conforming grids

1 Introduction

Fractured rock formations in the subsurface are of crucial importance in a variety of reservoir applications, such as geothermal energy extraction, CO2 sequestration, nuclear waste storage, and unconventional oil and gas recovery [51, 37, 11, 12, 48, 2]. As fluid flow velocities in fractures are often magnitudes higher than in the rock matrix, individual fractures as well as fracture networks commonly govern the overall fluid transport characteristics of the entire fracture-dominated porous medium. Here, the geometric fracture configuration and the hydraulic properties of individual fractures, such as fracture permeability fields, largely determine, where preferential fluid flow may occur, as fractures with particularly high or low permeability can act as flow conduits or ”bottlenecks”, respectively [16, 56, 17]. Furthermore, Ahkami et al. [1] describe and visualize experimentally that the permeability of the porous-medium matrix influences fluid flow in the fractures of a fractured porous medium.

An in-depth understanding of processes in fractured rock masses thus requires knowledge about the hydraulic parameters of each fracture in a fracture network. As these parameters are notoriously difficult to obtain in the subsurface, reliance on stochastic investigations are often required, where hundreds or more system realizations typically have to be performed to assess uncertainties [7, 15, 28, 41]. To facilitate solving large numbers of numerical simulations of fluid flow through fracture-network systems in three dimensions (3D), highly efficient and accurate numerical methods and mesh generation approaches are required.

Generally, in numerical models, two method classes are used when representing fractured porous media, i.e. fractures embedded in a porous-medium matrix. Fractures are either represented by a continuum approach [55, 6, 33, 32] or as discrete domains in a numerical mesh [42, 4]. In the continuum approach, the fractures and the porous-medium matrix share the same geometric mesh with separate continua. The respective flow properties are obtained by upscaling and information needs to be transfered between the continua. In contrast, the classic discrete-domain approach explicitly meshes both fractures and porous media, with the two meshes conforming at the boundaries of the domains (i.e. conforming numerical mesh). Since fracture configurations in fracture networks can be arbitrarily complex, mesh generation for discrete fracture networks (DFN) [30, 43, 44, 22] or discrete fracture models (DFM) [14, 29, 9] with the background matrix can be very difficult and time consuming. Due to the large length-to-width ratio of fractures and the need to overcome very small elements in the fracture domain, fractures may be represented by lower-dimensional elements, e.g. [31, 10, 39, 26]. Still, generating smooth matrix meshes on, and around, fracture intersections and fracture tips can lead to very small elements and significant increases in the number of degrees of freedom. In contrast, areas in the (porous) medium that are void of fractures might contain elements with large edge lengths, leading to large differences in element size, compared to elements close to fractures. The strong influence of fractures on fluid mass and energy transfer processes in a wide range of applications, as mentioned above, and the associated difficulties in model generation, have rendered related method improvements an area of active research.

The aforementioned shortcomings of classic discrete-domain approaches have led to an increased focus on the development of numerical methods that allow the use of independent meshes for the fracture domain and the matrix domain. Such non-conforming mesh approaches might be based on mortar methods [21, 13], where the mesh for the discrete fractures and the matrix are required to align geometrically, but consist of independent discretizations. Larger flexibility is offered by methods that handle fracture and matrix meshes separately (i.e. no aligned geometries). Such methods exist for finite volume schemes, e.g. (p)EDFM [25, 50, 38] and for XFEM-based approaches [19, and references therein] for finite elements. A review of existing mathematical and conceptual models for flow in fractured porous media is given by Berre et al. [8]. Recently, Köppel et al. [35] proposed a Lagrange multiplier method for a non-conforming finite element formulation. This method enables the use of independent meshes for fractures and matrices by applying variational transfer between the two mesh domains. The variational transfer allows projection of variables between the fractures and the matrix domain. With the geometric mapping between the fractures and the matrix domain established, the Lagrange multiplier accounts for the variable coupling at the domain boundaries. In contrast to XFEM methods, this approach does not enrich the finite element space locally, so that the pressure across the fractures is assumed to be continuous and the fractures have a higher permeability than the matrix. Köppel et al. [35] developed a general Lagrange multiplier method for 2D or 3D model domains. They further show the uniqueness of the solution for the primal formulation of the continuous problem. Due to the large technical complexities in 3D, they focused on the implementation and verification in 2D.

However, flow through fractured rock formations is governed by 3D effects, as strong heterogeneities affect flow properties in the fracture and porous-medium domains [52, 16, 54, 53]. To apply the Lagrange multiplier method to 3D fracture networks, the variational transfer between fractures and matrix for non-conforming methods needs to be implemented both accurately and with parallel processing capabilities. To this end, Krause and Zulian [36] developed a general, accurate and parallel variational transfer approach, which requires no prior knowledge about the relationship between the two meshes. More specifically, an $L^{2}$-projection variational transfer operator has been shown [27] to provide better approximations than interpolations. Previous applications of this $L^{2}$-projection include fluid–structure interaction (FSI) problems and mechanical contact of rough fractures [40, 46, 47, 45]. Typically, these applications are solved in equi-dimensional domains. To address the representation of fractures by lower-dimensional manifolds (i.e. surface elements), the $L^{2}$-projection algorithm has been extended to surface–volume interactions. This enables transfer of information between surface elements and volume elements for fractures and matrices, respectively. Building on the work of Krause and Zulian [36] and Köppel et al. [35], the aim of this work is to demonstrate their methods’ applicability to steady-state, single-phase fluid flow in 3D fractured porous media.

This paper presents an application of the Lagrange multiplier method in combination with the $L^{2}$-projection variational transfer operator in 3D. Section 2 provides a brief overview of the method, by discussing the mathematical formulation, the discretization, the surface–volume interaction, and the implementation. Section 3 first compares 2D and 3D results to state-of-the-art benchmark results [20]. Next, the method is validated for 3D fracture networks by comparing it to results from conforming mesh simulations, which are used as a reference. For all cases described above, different mesh sizes and mesh convergences are discussed. Finally, we apply the method to realistic fracture networks with hundreds of fractures, demonstrating the capability of the Lagrange multiplier method to model single-phase flow through realistic 3D fracture networks. The presented findings are then summarized in Section 4.

2 Method

In order to accommodate fluid flow through 3D fractured porous rock volumes, the Lagrange multiplier formulation, proposed by Köppel et al. [35], is applied and solved in 3D. The Lagrange multiplier formulation considers fractures as lower-dimensional manifolds (i.e. surface elements). In 3D, this results in surface domains for the fractures and one volume domain for the rock matrix. Accurate transfer of fluid pressure between surface and volume domains is accomplished by using the $L^{2}$-projection variational transfer operator [36] to discretize the Lagrange multipliers. The following subsections discuss the mathematical formulation, followed by a description of the discretization, the $L^{2}$-projection method for surface–volume interaction, and the implementation.

2.1 Mathematical formulation

Following Köppel et al. [35], we formulate a continuous Lagrange multiplier fracture problem.

The matrix domain is designated with $\Omega$ $\subset$ $\mathbb{R}^{n}$, $n=2$ or $3$, and the fracture domain with $\gamma$ $\subset$ $\Omega$ of dimension $n-1$. A normal vector, $\textbf{n}_{\gamma}$, is defined with respect to the fracture surface (Fig. 1). Steady-state fluid flow in the porous-medium matrix, $\Omega$, is governed by

[TABLE]

where $\mathbf{K}$ is the permeability tensor, $p$ is the fluid pressure, and $f$ is the sink/source term.

Flow in the fracture, $\gamma$, is described by

[TABLE]

Above, $\partial\Omega$ and $\partial\gamma$ is the interface boundary between $\Omega$ and $\gamma$. Fluid exchange between $\Omega$ and $\gamma$ is given by $\lambda=\lambda(x),\;x\in\gamma$.

The spaces $V_{\Omega},V_{\gamma},\mathbf{V}$, and $\Lambda$ are defined by:

[TABLE]

with the test functions $q\in V_{\Omega}$, $q_{\gamma}\in V_{\gamma}$, and $\mu\in\Lambda$. The variational formulation is found by multiplying Eqs. (1) and (2) by the test functions, integrating over $\Omega$ and $\gamma$, and using integration by parts on both equations. From that, the variational formulation is given as follows:

Find $(p,p_{\gamma})\in\mathbf{V}$ and $\lambda\in\Lambda$, such that

[TABLE]

and

[TABLE]

Here, Eq. (5) indicates the coupling conditions between the domains $\Omega$ and $\gamma$. The Lagrange multiplier represents the fluid pressure gradient $\lambda=\mathbf{K}\nabla p\cdot\mathbf{n}_{\gamma}$ and ensures the continuity of the fluid pressure and the exchange of the forces between the fracture domain, $\gamma$, and the matrix, $\Omega$, in the direction normal to $\gamma$.

Köppel et al. [35] show that there exists a unique solution to the variational formulation of this problem. Further details on the mathematical proof can be found in their work.

2.2 Discretization

Based on the variational formulation in Eq. (4), the discrete counterpart is formulated in a finite element framework. In order to solve the discrete formulation, three distinct meshes are defined to approximate the matrix $\Omega$, the fracture $\gamma$, and the Lagrange multiplier $\lambda$ in the fracture.

The meshes for the finite element discretization of Eq. (4) are defined to be $\mathcal{M}$ in $\Omega$, $\mathcal{M}_{\gamma}$ in $\gamma$, and $\mathcal{M}_{\lambda}$ in $\lambda$. The respective mesh widths $h_{\mathcal{M}}$, $h_{\mathcal{M},\gamma}$, and $h_{\mathcal{M},\lambda}$ are defined by:

• $h_{\mathcal{M}}:=\max\limits_{1\leq M\leq\mathcal{M}}h_{M}$, where $h_{M}$ = diam $M$,
• $h_{\mathcal{M,\gamma}}:=\max\limits_{1\leq m\leq\mathcal{M}_{\gamma}}h_{m}$, where $h_{m}$ = diam $m$,
• $h_{\mathcal{M,\lambda}}:=\max\limits_{1\leq n\leq\mathcal{M}_{\lambda}}h_{n}$, where $h_{n}$ = diam $n$.

To enhance readability, the respective mesh widths are reduced to $h$, $h_{\gamma}$, and $h_{\lambda}$ for the remainder of this paper.

The shape of elements in each mesh might be of any kind. Hence, the approximation spaces $V_{h,\Omega}$, $V_{h,\gamma}$, and $\Lambda_{h}$ of continuous, piecewise-polynomial functions on $\Omega$, $\gamma$, and $\lambda$ are defined by:

[TABLE]

The definition of the approximation space $\Lambda_{h}$ of the Lagrange multiplier implies that $h_{\lambda}=h_{\gamma}$. To prevent a poorly conditioned system matrix, it is necessary to satisfy $h_{\lambda}\leq\min(h,h_{\gamma})$, which – in combination with the above – results in $h\leq h_{\gamma}$. Generally, by using $q_{\gamma}|_{m}\in\mathbb{P}^{1}(m)$ and $\lambda|_{n}\in\mathbb{P}^{1}(n)$, the Lagrange multiplier $\lambda$ is applied on the fracture mesh. This offers higher flexibility regarding the meshing.

2.3 Volume–surface information transfer

The meshes associated with the spaces $V_{h,\Omega}$ (matrix) and $V_{h,\gamma}$ (fracture network) are generally non-matching (i.e. Fig. 2), which means that the surfaces of the volume elements of the matrix do not necessarily coincide with those of the surface elements of the fracture network.

This leads to mutually non-conforming discretizations which require the use of information transfer techniques such as interpolation or $L^{2}$-projections for handling the coupling terms in Eq. (4). Here the $L^{2}$-projection approach is adopted, as it has been shown to have better approximation properties than interpolation [27].

Ultimately, intersections between the matrix mesh and the fracture mesh have to be found in order to perform quadrature on the coupling term in Eq. (5) up to numerical precision. Then, integration is done on these intersections. For any pair of a volume element $M\in\mathcal{M}$ and a shell element $m\in\mathcal{M}_{\gamma}$, the following procedure is completed:

• •

Computation of the intersection $I=M\cap m$ by using a variant of the Sutherland–Hodgman clipping algorithm [49] as shown in Fig. 3. Here, $M$ is interpreted as a set of half-spaces which are sequentially used to clip $m$.

• •

If $I\neq\emptyset$, $I$ is meshed into the simplicial complex $\mathcal{T}_{I}=\{S\}$, where $S$ is a simplex.

• •

A suitable quadrature rule is mapped to each simplex $S$ which is then mapped to the reference configurations of $M$ and $m$.

The task of detecting pair-wise element intersections is accelerated by employing octree data structures. Additional acceleration can be gained by applying techniques such as spatial-hashing [18]. The applied algorithms are fully automated and no prior knowledge about the relation of the meshes is required. Further details regarding the information transfer procedure can be found in Krause and Zulian [36].

2.4 Implementation

The routines described in this paper are implemented within the open-source software library Utopia [58]. Utopia uses libMesh [34] for the finite element discretization, MOONoLith [57] for the intersection detection, and PETSc [5] with MUMPS [3] for the linear algebra calculations. In the numerical experiments illustrated in Section 3, the size of the algebraic linear system of equations (4) reaches at most one million degrees of freedom, which is solved with the MUMPS direct solver.

3 Numerical results & Discussion

This section discusses various numerical experiments performed to investigate the Lagrange multiplier method in 2D and 3D. First, the implementation is verified in 2D by comparing the results to a benchmark case presented by Flemisch et al. [20]. In a second experiment, the 2D case is extruded in the third dimension. This enables testing the accuracy of the Lagrange multiplier method, combined with the $L^{2}$-projection variational transfer operator (LM–L2) in 3D by comparing it to the 2D benchmark results. Next, a heterogeneous 3D fracture network is built and the LM–L2 method is validated by comparing it with results from a conforming mesh model. Finally, the LM–L2 method is applied to a realistic fracture network with 150 randomly distributed fractures. These final numerical experiments are conducted to investigate different geometrical complexities as well as mesh convergence.

Throughout the numerical experiments, various element shapes are used, following Eqs. (6). The implementation of the LM–L2 method causes no restrictions regarding the element shape, so that the matrix domain, $\Omega$, and the fracture domain, $\gamma$, can therefore be of arbitrary shape. However, for simplicity, the matrix mesh, in all presented numerical experiments, is composed of hexahedral (for 3D) and quadrilateral (for 2D) elements. The respective element spaces for all experiments is second order in $\Omega$ and first order in $\gamma$ and $\lambda$. Nevertheless, the implementation allows consideration of zero-, first-, and second-order elements in all domains. To test convergence and accuracy, various mesh sizes for $\Omega$ and $\gamma$/$\lambda$ are employed throughout this study.

To facilitate comparison with existing studies [20, 35], all physical parameters are normalized throughout this study. All experiments are conducted on a unit domain with an edge length of $1$. Further, the fracture permeability is $\mathbf{K}_{\gamma}=10^{4}\mathbf{I}_{\gamma}$ and the fracture aperture is $a=10^{-4}$, which is incorporated in the applied fracture permeability and the Neumann boundary condition. The permeability in the matrix domain is set to $\mathbf{K}=\mathbf{I}$. These choices ensure that, for the cases studied here, the matrix and fracture network both contribute to a similar extent to the overall fluid flow, i.e. the ratio between a) the average aspect ratio of the fractures and b) the permeability ratio between the matrix and the fractures approximately equal to one [17].

3.1 Benchmark – 2D

First, the implementation of the LM–L2 method is tested by comparing the obtained results to 2D benchmark results presented by Flemisch et al. [20]. Benchmark $1$, which was first introduced by Geiger et al. [24], is chosen as a representative example for this study. Emphasis is placed on the validation of the presented method with reference results from Geiger et al. [24], Flemisch et al. [20], while results of alternative approaches can be found in the listed references. While Köppel et al. [35] used $\mathbb{P}^{0}$ elements for their benchmark comparison, this study employs $\mathbb{P}^{1}$ elements for the Lagrange multiplier mesh and $h_{\gamma}=h_{\lambda}$. The fracture elements consist of line segments and the matrix of quadrilateral elements. Further, a mesh-size convergence study is conducted for the fracture and the matrix mesh.

The fracture network at hand contains six fractures in perpendicular orientation as shown in Fig. 4.

A non-homogeneous Neumann boundary condition on the left domain boundary serves as a fluid source and a non-homogeneous Dirichlet boundary condition on the right domain boundary serves as a fluid sink, which yields fluid flow across the domain from the left to the right boundary. Homogeneous Neumann boundary conditions at the top and bottom boundaries enforce no-flow conditions across those domain boundaries.

The resulting pressure distribution for $h=1/129$, and $h_{\gamma}=h_{\lambda}=1/128$ is depicted in Fig. 5, and shows good agreement with the reference benchmarks [24, 20].

As expected, the central, horizontal fracture serves as a channel of high permeability, facilitating faster fluid flow than the surrounding porous-medium matrix. Matrix regions further removed from the high-permeability fractures (e.g., top and bottom left corner) therefore result in the steepest fluid pressure gradients.

Mesh-size dependency is subsequently investigated in a convergence study. Mesh widths are chosen such that non-conforming meshes for the fractures and the matrix are ensured. The initial realization has the largest mesh width of $h=1/33$ and $h_{\gamma}=h_{\lambda}=1/32$. With each refinement, the mesh resolution is then increased by a factor of two, resulting in $h=1/65$ and $h_{\gamma}=h_{\lambda}=1/64$, $h=1/129$ and $h_{\gamma}=h_{\lambda}=1/128$, and $h=1/257$ and $h_{\gamma}=h_{\lambda}=1/256$.

Figs. 6a and 6b show the results for all realizations and the reference results along the lines AA*′* and BB*′*, yielding good convergence and accuracy for all realizations.

However, for low resolutions of $h_{(\gamma/\lambda)}$, small deviations can be observed along BB*′* between $Arc\ length=0.6$ and $Arc\ length=1.0$. This can be attributed to a lack in resolution of the fracture intersections by the non-conforming mesh configuration. This error increases with decreasing matrix mesh resolution. Köppel et al. [35] suggest that results might be more accurate for $h\geq h_{\gamma,\lambda}$. However, the considered second-order function space in $\Omega$ reduces this effect.

Fig. 6c shows the RMSM in the matrix for $h=1/33$, $h=1/65$, and $h=1/129$, relative to $h=1/257$, which is decreasing approximately linearly with increasing mesh resolution. The error between different results is calculated by the root mean square (RMS) over all elements. To facilitate comparison of non-conforming meshes, results from all resolutions are interpolated on a mesh with the finest resolution. The squared error on a single element in the mesh $\mathcal{M}_{(\gamma)}$ and the RMSM,(γ) error are calculated following:

[TABLE]

where $n$ is a node in $\mathcal{M}$ or $\mathcal{M}_{\gamma}$, and $N=\Sigma n$.

3.2 Benchmark – 3D

The accuracy of the LM–L2 method in 3D is tested by extruding the setup introduced in Section 3.1 in the third dimension, which results in a 3D porous-medium matrix domain and surface domains for the fractures (Fig. 4). This extrusion enables comparison of 3D results and the 2D benchmarks from Flemisch et al. [20]. The extrusion causes the observation lines AA*′* and BB*′* (see Fig. 4) to be located anywhere in the direction of the extrusion. For the present case, their position is chosen at half of the extrusion length. Non-homogeneous Neumann boundary conditions are applied at the fracture and the matrix on the left boundary. A non-homogeneous Dirichlet boundary condition is applied on the right boundary. Homogeneous Neumann boundary conditions are imposed on the remaining boundaries.

Various realizations with different mesh sizes are performed to study convergence. Fig. 7 shows the pressure profile along the lines (a) AA*′* and (b) BB*′* for the reference results in 2D and the 3D results for all realizations. The initial mesh width of $h=1/33$ and $h_{\gamma}=h_{\lambda}=1/32$ is consecutively refined by a factor of two, which results in $h=1/65$ and $h_{\gamma}=h_{\lambda}=1/64$, and $h=1/129$ and $h_{\gamma}=h_{\lambda}=1/128$. The element shapes in this experiment are hexahedrons in $\Omega$ and quadrilaterals in $\gamma/\lambda$. A more detailed convergence study in 3D is presented for the more complex case in Section 3.3. Here, we restrict ourselves to a comparison of the pressure results along the observation lines AA*′* and BB*′* with the reference results (Fig. 7).

While the pressure profile is captured for the mesh widths $h=1/65$ and $h_{\gamma}=h_{\lambda}=1/64$, and $h=1/33$ and $h_{\gamma}=h_{\lambda}=1/32$, they display a slight deviation between $Arc\ length=0.6$ and $Arc\ length=1.0$. As accuracy improves with decreasing mesh size, the realization with $h=1/129$ and $h_{\gamma}=h_{\lambda}=1/128$ only shows minor deviations from the reference results. The more pronounced deviations for these coarser mesh realizations can be explained with a lack of resolution around the fracture intersections in the area of $Arc\ length=0.6$ to $Arc\ length=1.0$. This suggests that the matrix mesh width around the fractures and in particular at the fracture intersections should be smaller. However, this effect in 3D is similar to the one in 2D systems (see Section 3.1).

3.3 Heterogeneous Fracture Network – 3D

Here we generate an artificial heterogeneous fracture network to test the LM–L2 method with more complex geometries. This enables the investigation of accuracy and convergence for fractures with varying orientation, size, and location. Fractures are often assumed to be disc-shaped [17]. Hence, our numerical experiment also approximates fractures as circular surfaces. Fracture tips are generally difficult to represent by non-conforming mesh methods, as they characterize the boundary of a fracture domain. On the tips, the Lagrange multiplier can act in multiple spatial directions, while its primary direction in the center of the fracture is normal to the fracture plane (i.e. flow normal to the fracture plane). This becomes particularly important in 3D, if the fracture edge is not a line (e.g. corners or circular fractures), as this results in a discontinuity of the direction of the Lagrange multiplier. Throughout this experiment, only the matrix-mesh width, $h$, is varied and the results are compared to results from a conforming mesh simulation. This enables investigation of the influence of the matrix mesh width on the accuracy of the results.

The chosen model setup for seven fractures is shown in Fig. 8.

All fractures are located within the matrix domain and are characterized by a random diameter, orientation and location. At the boundaries BC0 and BC1, the applied Dirichlet condition is set to $1$ and [math], respectively. The homogeneous Neumann conditions imposed on the remaining boundaries ensures no-flow.

Convergence is studied with three realizations with a matrix-mesh width of $h=1/33$, $h=1/65$, and $h=1/129$, respectively. The fracture-mesh width for all realizations is constant at $h_{\gamma}=h_{\lambda}=1/129$. The results of these realizations are compared to reference results which are obtained with a conforming mesh using finite element methods in MOOSE [23] with a mesh width of $h=h_{\gamma}=1/129$. For the conforming mesh, the fracture elements are triangles and the matrix elements tetrahedrons. For the LM–L2 method, the element shapes in $\gamma$ are triangles and hexahedrons in $\Omega$.

The observed RMS errors in the matrix and fracture show linear convergence rates with an improved rate for $h=1/129$ (see Fig. 9).

Due to the similar mesh width of this realization to the reference model, this error solely compares the effect of non-conforming meshes. Furthermore, the observed convergence rates are similar for matrix and fractures. Hence, the changing mesh width in the matrix mesh affects the accuracy in the fractures and the matrix similarly. This further suggests that sufficiently high matrix resolution around the fractures is crucial. Generally, the RMS error in the fractures is higher because the entire domain is affected by the non-matching meshes. In contrast, large matrix areas are further away from the fractures and therefore not affected by the non-matching meshes.

Fig. 10 shows the local error at cross sections A, B, and C, which are defined in Fig. 8.

Overall, the error in planes A and B is small, although an increased error is present around the fractures and in particular the fracture tips. The biggest errors can be found on plane C, specifically at the fracture tips of the upper most large fracture. This fracture is cut through its center and has the largest diameter, thereby connecting high-pressure areas on the left with low-pressure areas on the right. These relatively large errors can also be observed for other fractures oriented parallel to the gradient, as they form a shortcut for pressure in the system. This leads to the largest pressure gradients at the fracture tips which are located close to the inflow and outflow of the domain, or which happen to be the farthest apart in the direction of the gradient. This is also where the biggest errors in the fracture domain are found. The difficulties encountered during the generation of the conforming mesh for the reference solution serve as further motivation for the development of non-conforming methods. Particularly, very small elements are necessary to represent the matrix between the fracture intersections with low intersection angle. Specialized meshing tools are required to efficiently mesh high-quality, conforming mixed-dimensional meshes [14, 29, 9].

3.4 Random Fracture Network – 3D

Fracture networks commonly contain hundreds or even thousands of fractures within a rock volume of $100\text{\,}\mathrm{m}$ side length. They are geometrically characterized by several distributions regarding the fracture size, orientation, density or location. It is commonly assumed that the fracture size is distributed following a truncated power law [12]. Fracture orientation, density, and location are site-specific and usually derived from geological information.

These highly heterogeneous geometries are computationally demanding, when generating meshes for numerical simulations. Using the LM–L2 method, the geometry only needs to be represented by the fracture mesh, as the porous-medium matrix is represented by a regular grid. This numerical experiment shows and investigates the application of the LM–L2 method on a realistic fracture network with 150 fractures.

In this experiment, the fracture radius distribution follows a power law, with truncations at 0.1 and 0.4. The 150 fractures are circular, randomly oriented, and distributed in a cubical model domain with a side length of 1. For simplicity, it is further ensured that no fracture intersects the model domain boundary. Fig. 11 shows the fractures embedded in the porous-medium matrix, of which only a single layer is displayed in red, illustrating, how the fracture mesh (green) cuts through the porous-medium matrix mesh (red). The mesh width for this numerical experiment is set to $h=1/33$ in the porous-medium matrix and $h_{\gamma}=h_{\lambda}=1/200$ in the fractures. In accordance with the previous numerical experiment, the element shapes are triangles in $\gamma$ and hexahedrons in $\Omega$. Fluid flow is modeled from left to right by non-homogeneous Dirichlet boundary conditions with values of 1 and 0, respectively. The imposed boundary conditions on the remaining boundaries are homogeneous Neumann boundary conditions.

The resultant dimensionless fluid pressure distribution in the fractures and the porous-medium matrix is depicted in Fig. 12. Additionally, three observation planes A, B, and C are defined to observe the fluid pressure in each plane (Fig. 13). Although the fluid pressure field is clearly influenced by the fracture network, Planes B and C show that the porous-medium matrix still contributes considerably to the overall fluid flow through the entire fractured porous medium. These results demonstrate the presented method’s capability of modeling fluid flow through complex fractured porous media.

4 Conclusion

This study builds on previous work to combine a Lagrange multiplier method with an $L^{2}$-projection variational transfer operator to numerically model single-phase fluid flow problems in 3D fractured porous media. The performance of the method is assessed by comparison with benchmark cases and numerical experiments in 2D and 3D. We also demonstrate the method’s suitability for large-scale, realistic fracture-network realizations in 3D.

Our comparison with benchmark simulations shows good agreement with reference results from the literature. In general, the accuracy of non-conforming methods at fracture intersections is known to depend on the width of the porous-medium matrix mesh. These findings are confirmed in this study by all of our numerical experiments, as we also observe small deviations in fluid pressures at fracture intersections for larger mesh widths. However, the RMS error in the porous-medium matrix mesh decreases linearly with increasing mesh size, as desired.

To study the presented method in 3D, the 2D benchmark case is extruded in the third direction, thereby enabling a comparison of 3D simulations with 2D benchmark results. Although the 3D case provides additional complexity for the transfer operator, results agree well with the benchmark results. Analogous to the 2D results, the 3D results also show the largest deviations from the reference results at fracture intersections.

A 3D fracture network, containing seven fractures, is presented, demonstrating the applicability of circular fractures with random orientations, sizes, and fracture tips inside the model domain. These simulations are compared to reference results obtained with conforming mesh simulations. A convergence study refined the porous-medium matrix mesh while leaving the fracture mesh unaltered. The observed convergence rates show a linear behavior in the porous-medium matrix and the fracture domains, with slightly improved convergence for the highest mesh resolution. The largest errors occur at the location of the steepest fluid pressure gradients, for example, at the tips of fractures that are aligned parallel to the fluid pressure gradient. More specifically, the errors occur at the fracture tips farthest apart from one another and in the direction of the fluid pressure gradient.

Finally, a numerical experiment of a realistic discrete fracture network (DFM), with 150 fractures, demonstrates the capability to model single-phase flow through highly complex fractured porous media. The obtained fluid pressure contours in the 3D domain show the expected behavior for a fractured porous medium, with high-permeability fractures, largely influencing the fluid flow through the numerical model domain. Despite a coarse porous-medium matrix mesh, the presented method achieves detailed representations of fractures embedded in a porous-medium matrix domain.

Generally, the Lagrange Multiplier–L2-projection method shows good agreement with benchmark cases and reference results, while yielding good convergence. In all cases, our results suggest that the matrix mesh of the porous-medium should be locally refined when accuracy is to be improved, particularly around fracture tips and at fracture intersections.

Further research will focus on the development of mesh adaptivity algorithms that allow pressure-gradient- and fracture-location-dependent porous-medium matrix-mesh refinements. Additional extensions of the numerical approach presented here will target transient fluid flow computation for which the parallel-processing algorithm, employed in this study, can be optimized.

Acknowledgment

M.O.S., P.S., A.E., D.V., and S.B.R. thank the Werner Siemens Foundation for their endowment of the Geothermal Energy and Geofluids group at the Institute of Geophysics, ETH Zurich. P.Z., M.G.C.N, and R.K. thank the SCCER-SoE program. We gratefully acknowledge the discussion with Markus Köppel regarding his implementation of the Lagrange multiplier method.

Computer Code Availability

All methods and routines, used for this study, are implemented with the open-source software library Utopia [58]. Utopia’s lead developer is co-author Patrick Zulian at USI Lugano, Switzerland. Co-developers are Alena Kopaničáková, Maria Chiara Giuseppina Nestola, Andreas Fink, Nur Fadel, Victor Magri, Teseo Schneider, and Eric Botter.

The contact address and e-mail of Patrick Zulian are as follows:

• Institute of Computational Science

Università della Svizzera italiana (USI - University of Lugano)

Via Giuseppe Buffi 13

CH-6904 Lugano

• patrick.zulian@usi.ch

Utopia was first available in 2016, the programming language is C++ and it can be accessed through a git repository or a docker container on:

• https://bitbucket.org/zulianp/utopia (18.6 MB),
• https://hub.docker.com/r/utopiadev/utopia.

The software dependencies are as follows:

• PETSc (https://www.mcs.anl.gov/petsc/),

must be compiled with MUMPS enabled

• libMesh for the FE module (https://github.com/libMesh)

There are no hardware requirements given by Utopia. Potential hardware or software requirements of the underlying libraries libMesh and PETSc are not stated here.

ORCID

P. Schädle - https://orcid.org/0000-0002-5485-2028

P. Zulian - https://orcid.org/0000-0002-5822-3288

D. Vogler - https://orcid.org/0000-0002-0974-9240

S. Bhopalam R. - https://orcid.org/0000-0001-7221-695X

M.G.C. Nestola - https://orcid.org/0000-0002-5700-0306

A. Ebigbo https://orcid.org/0000-0003-3972-3786

R.H. Krause - https://orcid.org/0000-0001-5408-5271

M.O. Saar - https://orcid.org/0000-0002-4869-6452