Multiscale Hybrid-Mixed Finite Element Method for Flow Simulation in Fractured Porous Media

The multiscale hybrid-mixed (MHM) method is applied to the numerical approximation of two-dimensional matrix fluid flow in porous media with fractures. The two-dimensional fluid flow in the reservoir and the one-dimensional flow in the discrete fractures are approximated using mixed finite elements. The coupling of the two-dimensional matrix flow with the one-dimensional fracture flow is enforced using the pressure of the one-dimensional flow as a Lagrange multiplier to express the conservation of fluid transfer between the fracture flow and the divergence of the one-dimensional fracture flux. A zero-dimensional pressure (point element) is used to express conservation of mass where fractures intersect. The issuing simulation is then reduced using the MHM method leading to accurate results with a very reduced number of global equations. A general system was developed where fracture geometries and conductivities are specified in an input file and meshes are generated using the public domain mesh generator GMsh. Several test cases illustrate the effectiveness of the proposed approach by comparing the multiscale results with direct simulations.


Introduction
Modeling and simulation of fluid flow in naturally and hydraulically fractured subsurface systems has been a popular research topic in petroleum engineering. According to [Schlumberger (2017)], more than half of the world total oil reserves reside in fractured carbonate reservoirs. Many such reservoirs have their main productivity channel from a network of connected fractures. Oil recovery mechanism in fractured reservoirs, especially shale or tight formations with hydraulically induced fractures, is strongly effected by the accuracy of fracture description. Mathematical models for describing recovery process by numerical studies; see Moinfar [Moinfar (2013)] and references therein for details. It should be noted that the number of degrees of freedom is still beyond the scope of classical computational methods with homogenizing natural fractures [Ţene, Al Kobaisi and Hajibeygi (2016)]. This motivates the development of an efficient multiscale method for heterogeneous fractured porous media. The main idea of multiscale Galerkin methods goes back to the work by Babuška et al. [Babuška and Osborn (1983); Babuška, Caloz and Osborn (1994)] and they construct basis functions by solving local flow problems to incorporate fine-scale effects into coarse-scale equations. Thereafter, this idea has been further developed by many researchers and spawned a series of relevant methods, including multiscale finite element method (MsFEM) [Hou and Wu (1997)], multiscale mixed finite element method (MsMFEM) [Chen and Hou (2002); Aarnes (2004)], heterogeneous multiscale method (HMM) [Weinan and Björn (2005)], generalized multiscale finite element method (GMsFEM) [Efendiev, Galvis and Hou (2013)], and multiscale finite-volume method (MsFVM) Jenny, Lee and Tchelepi (2003). During the past few years, multiscale methods have been extended to simulate fluid flow in fractured media. For example, DFM has been integrated into MsFEM by Zhang et al. [Zhang, Huang, Yao et al. (2017)] and GMsFEM by Akkutlu et al. [Akkutlu, Efendiev and Vasilyeva (2016)]. One of the main contributions of this work is the combination of modeling discrete fracture networks with the Multiscale Hybrid-Mixed (MHM) method. The MHM method was originally developed by Valentin, Harder, and Paredes (see Harder et al. [Harder, Paredes and Valentin (2013); Paredes, Valentin and Versieux (2017); Araya, Paredes, Valentin et al. (2013)] for more details) and is a numerical approximation technique geared towards the approximation of conservation laws that incorporate multiple scales. The MHM method has been extended to mixed finite element approximations in Duran et al. [Duran, Devloo, Gomes et al. (2018)], in which the MHM method is described and compared with other multiscale methods such as the Hybrid Discontinuous Galerkin method [Cockburn (2016)] and the Hybrid High-Order method [Di Pietro, Ern and Lemaire (2016)]. The extension of the MHM method to the numerical simulation of discrete fracture networks is to extend the boundary fluxes over each polygonal domain with fluxes associated with the fracture flow. In this paper, we propose a multiscale hybrid-mixed finite element method for DFM and describe its implementation in the NeoPZ finite element library [Devloo (1997)]. The rest of this paper is organized as follows. The basis of representing fluid flow through the interfaces of the elements using mixed finite elements is a hybridized version of mixed finite element approximations and is presented in Section 2. The fluid flow in the fracture is modeled using a lower dimensional flux is described in Section 3. Finally numerical examples in Section 4 demonstrate the simulation of flow through porous medium where volumetric flow is combined with fluid flow through fractures.
In this section, we describe a mathematical formulation for approximating discrete fracture models coupled with the porous media flow. We propose a hybridized mixed finite element formulation, which substantially decreases the size of the discrete linear system. In this formulation, we use a hierarchal grid of two length scales: (1) a coarse-scale structured or unstructured grid based on which global degrees of unknowns are defined, and (2) a finescale unstructured grid is employed to describe long fractures. We note that such a method can also be combined with continuum models in order to take short and medium fractures into account. Let T d ⊂ R d be a polygonal computational domain. The mixed approximation of the Laplace equation on T d can be formulated as : Here n is the outer unit normal vector on element boundaries.

Hybridized mixed finite element approximations
Mixed finite element approximations have the distinct feature that the numerical results of the hybridized formulation are identical to the original formulation. Moreover, when statically condensing the internal fluxes and pressures onto the pressure on the interfaces, the issuing global matrix is symmetric positive definite (SPD). For simplicity, we restrict our discussion in two-dimensional (d = 2) matrix domain with one-dimensional fractures in this paper. We denote the partition of macro elements as T 2 := {T 2 : T 2 ⊂ Ω} and the interfaces between macro elements as T 1 in this case. A hybrid mixed finite element approximation of a flow through porous medium is described as: The meaning of σ 2 · n and p 1 is illustrated graphically in Fig. 1: the pressure p 1 acts as a Lagrange multiplier to enforce the continuity of σ 2 · n. Eq. (3) represents the Darcy's constitutive law of the fluid flow in porous media. The third term corresponds to the pressure at the interfaces of the macro fluxes that acts as a (weak) Dirichlet condition for the macro domain. Eq. (4) represents the two dimensional conservation law. Eq. (5) establishes that the sum of all normal fluxes at the interface between macro elements has to be weakly zero. Furthermore, the H div (T d ) approximations can be hybridized one side at a time.

Figure 1:
Hybridizing an interface between two H(div)-elements

The multiscale hybrid-mixed method
The multiscale hybrid mixed (MHM) method is a technique developed towards the numerical approximation of partial differential equations whose solution exhibit multiscale features. Within the MHM framework the normal component of the flux over the macro elements is approximated by piecewise continuous functions. The extension of these flux functions in the interior of the macro domains constitute the MHM basis functions.
To illustrate the concept of the extension of boundary fluxes in conjunction with discrete fracture networks, a square macro domain is used with a fracture as illustrated in Fig. 2. The permeability of the matrix is unitary and the permeability of the fracture times fracture width is 100.   Observing the function extensions for this domain with one fracture, the capability of the MHM method for simulating multiscale problems is apparent-the singularity caused by the fracture point is isolated in the interior of the domain. The fine mesh used to capture the details of the flow pattern is confined in the interior of the domain. The size of global system of equations is determined uniquely by the number of boundary fluxes and is, therefore, not affected by the resolution of the interior mesh.

MHM and H(div) approximations
In MHM approximations the computational domain is discretized in two levels: a coarse polygonal mesh which is referred to as macro domains and a fine discretization on micro elements meshing the interior of each macro domain. This two-level discretization is illustrated in Fig. 5. The macroscopic fluxes between the macro elements represent the coarse-grained response of the fluid flow problem. The approximations using fine meshes at the interior of each macro element capture the detailed behaviour of the fluid flow.  (2015)] for example). As such the formulation of the MHM method is identical to the mixed finite element formulation. The difference is that in MHM the flux approximation space is partitioned between macro fluxes associated with the boundary of the macro domaines and internal fluxes and pressures associated with the interior.

Coupling fractures with surrounding porous media
Our work is based on the DFM formulation, which represents fractures as simplified (n − 1)-dimensional entities in an n-dimensional domain. As a result, fractures distributed in 2D space will be discretized in 1D form. Fig. 6 shows a porous rock with one fracture to illustrate the idea of DFM. The matrix domain is represented by Ω m , and the fracture domain by Ω f . The assumption of DFM is that all variables remain constant in the cross direction of the fractures. Therefore, the only difference between the single-porosity model and DFM is the integral computation for the fracture domain. In DFM, the fracture aperture can be taken as a factor in front of the 1D integral to simplify the problem considerably.

Modeling fluid flow in a discrete fracture network
The fluid flow in the fractures is modeled using a mixed approximation. This formulation is based on previous work of the authors published in [Castro, Devloo, Farias et al. (2016)].
The fluid flow in a single fracture can be modeled as: Find ( σ 1 , p 1 ) ∈ H div (T 1 ) × L 2 (T 1 )  (2001)]) such that When fractures intersect, the conservation of mass states that the sum of the normal fluxes at the intersection should sum to zero. The conservation at the intersection is imposed by the introduction of a Lagrange multiplier which has the physical quantity of the pressure at the intersection. The statement then becomes: Eq. (8) expresses the constitutive law of one-dimensional fluid flow in the fracture. p 0 represents the pressure at the fracture intersections. Eq. (10) expresses the conservation of mass at the fracture intersection. The flow in intersecting fractures is illustrated in Fig. 7: the red squares illustrate intersection points where 2, 3, or 4 fractures meet.

Coupling flow in porous media and discrete fracture network
If there is fluid exchange between the two dimensional fluid flow and the fracture, the conservation of mass requires that the fluid flow that leaves the two dimensional domain enters as a source term for the one dimensional conservation law. Reciprocally, the pressure of the fluid in the fracture acts as a pressure boundary condition for the two dimensional matrix fluid flow. The simulation spaces that will interact are listed as follows: • two dimensional pressures p 2 and weight functions ϕ 2 , • one dimensional pressures p 1 and weight, functions ϕ 1 , • one dimensional fluxes σ 1 and weight functions ψ 1 , • zero dimensional pressures p 0 and corresponding weight functions ϕ 0 .
In turn, we can now write down the issuing variational statement as Each equation in the above system corresponds to either a physical conservation law or a constitutive relation: 1. The first equation implements the constitutive law in two dimension indicating that the pressures in the fractures act as pressure boundary conditions.
3. The third equation implements the constitutive law in the fractures where the pressures at the end of the fractures (zero dimensional pressures) act as pressure boundary conditions. 4. The fourth equation implements mass conservation in the fracture where the normal fluxes of the two dimensional elements act as external sources.
5. The fifth equation represents mass conservation at the points where multiple fractures intersect. More than two subdomains can be connected through an interface (hence the generic Σ). This is the case when fluxes along one-dimensional manifolds in the computational domain come together. For example, Fig. 7 shows seven line segments linked together by three point-wise pressures.
The system of Eqs. (11)-(15) looks more complicated than it actually is. For simplicity, we may write these equations in the following saddle-point form: Here we denote t ij the j-th term of equation i (index starting from 1) then the correspondence between the integrals and NeoPZ integrals are • t 11 , t 12 , t 21 , r 2 are computed by a regular two dimensional H(div) element; • t 33 , t 34 , t 43 , r 1 are computed by a one dimensional H(div) element; • t 14 , t 41 are computed by a one dimensional interface element which lies between the two dimensional H(div) element and the one dimensional H(div) element; • r 3 is computed by a one dimensional element associated with a Dirichlet boundary condition; • t 35 , t 53 are computed by an interface element between the one dimensional H(div) elements and point elements at the intersection of fractures.
It is worthy to observe that static condensation can be applied at different levels to reduce the size of the global system of equations: • For each micro element, the internal fluxes and all but one pressure equation can be condensed on the boundary fluxes of the element.
• For each macro element, all internal fluxes and all but one pressure equation can be condensed on the boundary fluxes of the macro domain. The flux in the fractures where they intersect with boundaries of the macro domain can not be condensed.
Hence the global system of equations corresponds to the following variables: • the macro fluxes, • average pressure of each macro domain, • zero dimensional fluxes at the intersections of the fractures with the macro domain boundaries (one equation for each intersection).
This means that, after static condensation, we obtain a much smaller global linear system to solve; see the numerical comparison in the next section. Comparing to the MHM approximations without fracture flow, the global system of equations is incremented by one equation for each intersecting fracture.

Numerical experiments
In this section, we perform a few preliminary numerical experiments in order to validate the proposed numerical method, which shall be denoted as MHDFM for convenience.
Simulations are run on a personal computer. For comparison, we use a commercial software to obtain DFM and full-resolution direct simulation results. We always set the permeabilities of fractures and surrounding porous medium to be K 1 = K f = 10 5 and K 2 = K p = 1, respectively. The proposed algorithm is implemented using the object oriented framework for development of finite element algorithms NeoPZ. The NeoPZ environment implements one, two and three dimensional hp-adaptive finite element approximations. Modules are dedicated to the geometric approximation, the definition of the approximation space and the definition of the differential equation to be approximated. In this work we combine one and two dimensional H(div) approximation spaces through the use of interface elements.

Orthogonal fractures
A homogeneous 2D problem with two orthogonal fractures is considered (see Fig. 8(a)), where a "+"-shaped fracture network located at the center of the quadratic domain. No-flow boundary conditions are applied at the top and bottom, while the pressure values are set to p = 1 and p = 0 at the left and right boundaries. According to the length scales specified in Fig. 8(a), the fractures with aperture a = 0.04 are fully resolved using at least 225 × 225 grid cells (fully resolved direct simulation). For the DFM method, an unstructured grid is employed; see Fig. 8(b). For the MHDFM method, a 13 × 13 coarse grid is used; see But the numerical values are very close to each other. In order to make a quantitative comparison, we give Fig. 9 that represents the pressure along the horizontal center line of the domain. The seemingly horizontal line segment with little pressure loss in the plot matches with the fracture domain from (2, 4.5) to (7, 4.5) due to the fracture permeability is very large compared with the matrix permeability. The pressure obtained by MHDFM shows excellent agreement with DFM and direct simulation results. Note that the role of the vertical fracture is less important than that of the horizontal line one. The results from MHDFM, DFM, and direct simulation are close; see Tab. 1. Compared to the fully resolved fine-scale solution and conventional DFM result, the total flow rate computed by MHDFM is almost identical with a discrepancy less than 2.56 × 10 −3 .

Nonorthogonal fractures
Another homogeneous 2D case with two nonorthogonal intersecting fractures (see Fig. 11) is also tested. The boundary condition and other parameters are the same as in the first case and the only difference is the fracture distribution. The two fractures are respectively located at (3, 8; 6, 2) and (5, 7; 2, 3). The total flow rates computed by DFM and MHDFM are 1.2724 and 1.2661, respectively. The relative difference is 4.95 × 10 −3 .

Disjoint fracture networks
We summarize the size of global linear systems in different methods in Tab. 3. Except the first example, it is difficult for us to obtain computational meshes for the direct simulation  Figure 11: A DFM example with two nonorthogonal fractures method in order to resolve the fracture fully. We can see that we are able to reduce the size of global systems to solve by using the proposed multiscale hybrid-mixed method. Moreover, we make a few comments on the method: • In MHDFM, we use H(div) approximation for the flux and pressure variables. This requires more degree of freedom on each element, but gives better accuracy; • For convenience and load balance, we use uniform finer meshes in the MHDFM simulation. This can be replaced by locally refined meshes to reduce computational cost; • In this section, we focus on validation of the method and did not pay attention to the actual computational efficiency or parallelization of MHDFM. These aspects are the potential advantages of MHDFM and deserve further studies in the future.

Conclusions
In this paper, a discrete fracture model is approximated using mixed finite elements. The fluid flow in the fracture is modeled using lower dimensional elements applied to a Laplace-Beltrami operator. Fracture intersections are modeled using a Lagrange multiplier enforcing local conservation. The pressure in the fracture acts as a boundary condition for the two dimensional flow. The Multiscale Hybrid Method is applied to separate the local features of the fracture-reservoir coupling from the global features of fluid flow. The resulting numerical model leads to a very reduced global system of equations. In turn we are able to reduce computational cost to an acceptable level. Results of the proposed numerical model are compared with the standard DFM simulation and the finescale simulations using finite volume techniques. In all cases the difference in flow rate are less than 1%.