Flow in porous media with low dimensional fractures by employing enriched Galerkin method

This paper presents the enriched Galerkin discretization for modeling ﬂuid ﬂow in fractured porous media using the mixed-dimensional approach. The proposed method has been tested against published benchmarks. Since fracture and porous media discontinuities can signiﬁcantly inﬂuence single- and multi-phase ﬂuid ﬂow, the het- erogeneous and anisotropic matrix permeability setting is utilized to assess the enriched Galerkin performance in handling the discontinuity within the matrix domain and between the matrix and fracture domains. Our results illustrate that the enriched Galerkin method has the same advantages as the discontinuous Galerkin method; for example, it conserves local and global ﬂuid mass, captures the pressure discontinuity, and provides the optimal error convergence rate. However, the enriched Galerkin method requires much fewer degrees of freedom than the discontinuous Galerkin method in its classical form. The pressure solutions produced by both methods are similar regardless of the conductive or non-conductive fractures or heterogeneity in matrix permeability. This analysis shows that the enriched Galerkin scheme reduces the computational costs while oﬀering the same ac- curacy as the discontinuous Galerkin so that it can be applied for large-scale ﬂow problems. Furthermore, the results of a time-dependent problem for a three-dimensional geometry reveal the value of correctly capturing the discontinuities as barriers or highly-conductive fractures.

There are two main approaches to represent the fluid flow between the matrix and fracture domain ( Nick and Matthai, 2011b;Flemisch et al., 2018;Juanes et al., 2002 ). The first model, an equi-dimensional model, discretizes the matrix and fracture domain with same dimensionality . This approach is straightforward to implement, and no coupling condition is required. This approach is utilized to model, for example, coupled hydromechanical of fractured rocks using the finite-discrete element method ( Latham et al., 2013;, and can also capture fracture propagation using an immersed-body method pared to the size of matrix domain Martin et al., 2005;Berrone et al., 2018 ). The second approach has several benefits; for example, it reduces the degrees of freedom (DOF) ( Nick and Matthai, 2011a ) as the fracture domain is represented as the interface, which is part of the matrix domain, and subsequently, this approach can improve mesh quality (reduce the mesh skewness) ( Matthai et al., 2010 ). Since the fractures are interfaces, one can use a larger mesh size, which satisfies Courant-Friedrichs-Lewy (CFL) condition more easily ( Juanes et al., 2002;Nick and Matthai, 2011b ).
We focus on the finite element based discretization such as the DG and MFE methods that ensure the local mass conservative property. Moreover, they are flexible enough to discretize complex subsurface geometries such as intersections of fractures or irregular-shaped matrix blocks. Additionally, the aforementioned methods are capable of mimicking the fracture propagation in poroelastic media using either cohesive zone method or linear elastic fracture mechanics framework ( Salimzadeh and Khalili, 2015;Salimzadeh et al., 2019b;Secchi and Schrefler, 2012;Segura and Carol, 2008 ). However, the MFE method requires an additional primary variable (fluid velocity) which may require more computational resources, especially in a three-dimensional

Fig. 2.
Comparison of ratio of degrees of freedom for 1 st , 2 nd , 3 rd , 4 th , and 5 th polynomial degree cases on triangular element among CG, EG, and DG discretizations: ( a ) ratio of EG over CG DOF, ( b ) ratio of DG over CG DOF, and ( c ) ratio of DG over EG DOF. Fig. 3. Illustration of ( a ) equi-dimensional and ( b ) mixed-dimensional settings. Note that these illustrations are the graphical representations; in the numerical model Ω p , Ω q , Γ p , or Γ q have to be imposed on all boundary faces of each domain. domain Kadeethum et al. (2019a) . Mesh adaptivity is also not straightforward to implement Lee and Wheeler (2017) , and it requires the inversion of the permeability tensor, which may lead to an ill-posed problem ( Choo and Lee, 2018 ). The DG method also can be considered as a computationally expensive method as it requires a large number of DOF ( Sun and Liu, 2009;Lee et al., 2016a ).
To resolve some of the shortcomings mentioned above, we propose an enriched Galerkin (EG) discretization ( Lee et al., 2016a;Zi et al., 2004;Khoei et al., 2018 ) to model fluid flow in fractured porous media using the mixed-dimensional approach. The EG method, utilized in this study, composes of the CG function space augmented by a piecewiseconstant function at the center of each element. This method has the same interior penalty type bilinear form as the DG method ( Sun and Liu, 2009;Lee et al., 2016a ). The EG method, however, only requires to have discontinuous constants as illustrated in Fig. 1 , so it has fewer DOF than the DG method. Fig. 2 presents the comparison of the DOF ratio among CG, EG, and DG methods, and it shows that the EG method requires half of the DOF needed by the DG method (triangular elements with the first polynomial degree approximation). Note that this ratio decreases as the polynomial degree approximation increase. The EG method has been developed to solve general elliptic and parabolic problems with dynamic mesh adaptivity  and extended to address the multiphase fluid flow problems . Recently, the EG method has been also applied to solve the non-linear poroelastic problem ( Choo and Lee, 2018;Kadeethum et al., 2019b;2020b ), and compared its performance with other twoand three-field formulation methods ( Kadeethum et al., 2019a ). To the best of our knowledge, this is the first attempt to apply the EG discretization in the mixed-dimensional setting.
The rest of the paper is organized as follows. The methodology section includes model description, mathematical equations, and their discretizations for the EG and DG methods. Subsequently, the block structure used to compose the EG function space and the coupling terms between matrix and fracture domains is illustrated. The numerical examples section presents five examples, and the conclusion is finally provided.

Governing equations
We first briefly introduce the equi-dimensional model, which is used to derive the mixed-dimensional model. We are interested in solving steady-state and time-dependent single phase fluid flow in fractured porous media on Ω, which composes of matrix and fracture domains, Ω m and Ω f , respectively. Let Ω ⊂ ℝ be the domain of interest in ddimensional space where = 1 , 2, or 3 bounded by boundary, Ω. Ω can be decomposed to pressure and flux boundaries, Ω p and Ω q , respectively. The time domain is denoted by = ( 0 , ] , where > 0 is the final time.
The illustration of the equi-dimensional model is shown in Fig. 3 a. This model composes of two matrix subdomains, Ω mi where = 1 , 2 , and one fracture subdomain, Ω f . Note that, for the sake of simplicity, this setup is used to illustrate the governing equation with only two matrix subdomains, but in a general case, the domain may compose of n m subdomains, i.e. = 1 , 2 , … , . Moreover, the domain may contain up to n f fractures, where n m and n f are number of matrix subdomain and fracture, respectively. For simplicity, in this section we will consider = 1 . The fractures may not cut through the matrix domain, which we call immersed fracture setting. This topic will be discussed in Section 3 . The governing system with initial and boundary conditions of the equidimensional model assuming a slightly compressible fluid for the matrix domain is presented below: and for the fracture domain is: where ( · ) i represents an index, c is the fluid compressibility, m and f are the matrix and fracture porosity ( f is assumed to be one), k m and k f are the matrix and fracture permeability tensor, respectively, is fluid viscosity, p m and p f are matrix and fracture pressure, respectively, is fluid density, g is the gravitational vector, g m and g f are sink/source for matrix and fracture domains, respectively, n is a normal unit vector to any surfaces, p mD and p fD are prescribed pressure for matrix and fracture domains, respectively, q mN and q fN are prescribed flux for matrix and fracture domains, respectively, and 0 and 0 are prescribed pressure for matrix and fracture domains at = 0 , respectively. To formulate the mixed-dimensional setting as presented in Fig. 3 b, we integrate along the normal direction to the fracture plane Martin et al. (2005) . As a result Ω f is reduced to an interface, Γ. Note that the governing equation of the mixed-dimensional setting used in this study is proposed by Martin et al. (2005) in the mixed finite element formulation, which uses fluid pressure and fluid velocity as the primary variables. The mixed-dimensional setting has been used in the mixed formulation ( Keilegavlen et al., 2017 ) 8. Immersed fracture geometry and boundary conditions of ( a ) permeable fracture, ( b ) partially permeable fracture, ( c ) impermeable fracture cases, and ( d ) mesh that has = 272 . A' is fracture tip. Glaser et al., 2019 ) and DG discretization on polytopic grids ( Antonietti et al., 2019 ). The mixed-dimension strong formulation and its boundary conditions for the matrix domain are similar to the equidimensional model, (1) to (3) , but the strong formulation and its boundary conditions of the fracture domain are: where Γ p and Γ q represent pressure and flux boundaries of the fracture domain, respectively, is the tangential fracture permeability tensor, ∇ T and ∇ T · are the tangential gradient and divergence operators, which are defined as represents the fluid mass transfer between matrix and fracture domains Martin et al. (2005) , ⋅ is jump operator, which will be discussed later in the discretization part, and p fD and q fN are specified pressure and flux for the fracture domain, respectively.
In this study, if = 3 , k m as a full tensor is defined as: where all tensor components characterize the transformation of the components of the gradient of fluid pressure into the components of the velocity vector. The , , and represent the matrix permeability in x-, y-, and z-direction, respectively. k f , on the other hand, composes of two components: where is the normal fracture permeability. Note that we present here a general form of , to be specific, is a scalar if = 2 and tensor if = 3 . To represent the fluid mass transfer between matrix and fracture domains, following Martin et al. (2005) , Antonietti et al. (2019) , we define the coupling conditions between the matrix and fracture domain as: where ⟨ · ⟩ is an average operator, which will be presented later in the discretization part, f represents a resistant factor of the mass transfer between the fracture and matrix domains defined as: and ∈ (0.5, 1.0]. In this paper, we set = 1 . 0 for the sake of simplicity. In general case, is used to represent a family of the mixed-dimensional model, and more details can be found in Martin et al. (2005) , Antonietti et al. (2019) .

Numerical discretization
In this section, the discretization of the mixed-dimensional model is illustrated. The domain Ω is partitioned into n e elements,  ℎ , which is the family of elements T (triangles in 2D, tetrahedrons in 3D). We will further denote by e a face of T as illustrated in Fig. 4 . We denote h T as the diameter of T and define h ≔ max ( h T ), which we may refer as mesh size. Let  ℎ denotes the set of all facets, for the matrix domain,  0 ℎ the internal facets,  ℎ the Dirichlet boundary faces, and  ℎ denotes the Neumann boundary faces. Following Lee et al. (2016a) for any ∈  0 ℎ let + and − be two neighboring elements such that = + ∩ − . Let + and − be the outward normal unit vectors to + and − , respectively (see Fig. 4 Fig. 4 b.
We will further denote by e f a face of Γ h . Let Λ h denotes the set of all facets, for the fracture domain, Λ 0 ℎ the internal facets, Λ ℎ the Dirichlet boundary faces, and Λ ℎ denotes the Neumann boundary faces, and Λ ℎ ∩ Λ ℎ = ∅. Let n be the outward normal unit vector to Γ h , which is coincided with + of the + . Next, we define the jump operator of the scalar and vector values as: respectively, Moreover, by assuming that the normal vector n is oriented from + to − , we obtain Scovazzi et al. (2017) , the weighted average is defined as: where and and a harmonic average of + and − is defined as: The arithmetic average, = 0 . 5 , is simply denoted by ⟨ · ⟩. Note that, in general case, e is also defined for the k f ; however, for the sake of simplicity, we assume the material properties of the fracture domain are homogeneous. Hence, we perform these operators in the matrix domain only.
Remark 1. The numerical discretization discussed in this study only considers the case of a conforming mesh, i.e. the fracture domain, Γ, element is coincident with a set of faces of the matrix domain, Ω, as illustrated in Fig. 4 b.
This study focuses on two function spaces, arising from EG, and DG discretizations, respectively. We begin with defining the CG function space for the matrix pressure, p m as is the space for the CG approximation with k th degree polynomials for the p m unknown, ℂ 0 (Ω) denotes the space of scalarvalued piecewise continuous polynomials, ℙ ( ) is the space of polynomials of degree at most k over each element T , and m denotes a generic . Furthermore, we define the following DG function space for the matrix pressure, p m : where is the space for the DG approximation with k th degree polynomials for the p m space and L 2 ( Ω) is the space of square integrable functions. Finally, we define the EG function space for p m as: is the space for the EG approximation with k th degree polynomials for the p space and  is the space for the DG approx- Fig. 11. Regular fracture network example: ( a ) geometry and boundary conditions (fractures are shown in red), pressure solution, p m , in the matrix using = 2046 of ( b ) permeable fracture, ( c ) impermeable fracture cases, and ( d ) mesh that has = 184 . (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) imation with 0 th degree polynomials, in other words, a piecewise constant approximation. Note that the EG discretization is expected to be beneficial for an accurate approximation of the p m discontinuity across interfaces where high permeability contrast, either between + and − or k m and k f , is observed. On the other hand, the material properties of the fracture domain are assumed to be homogeneous, which leads to no discontinuity within the fracture domain. Therefore, the CG discretization of the p f unknown will suffice in the following. The p f function space is defined as: is the space for the CG approximation with k th degree polynomials for the p f unknown, ℂ 0 (Γ) denotes the space of scalarvalued piecewise continuous polynomials, ℙ ( ) is the space of polynomials of degree at most k over each facet e , and f denotes a generic Remark 2. In this study, we only focus on the mixed function space between the matrix and fracture domains arising from either  . The fracture domain can be dis- if there are any discontinuities inside the fracture medium.
The length of the subinterval, Δt n , is defined as Δ = − −1 where n represents the current time step. In this study, implicit first-order time discretization is utilized for a time domain of (1) and (5) as shown below for both ,ℎ and ,ℎ : We denote that the temporal approximation of the function Φ( · , t n ) by Φ n . With given −1 ,ℎ and −1 ,ℎ , we now seek the approximated solutions First, the temporal discretization part is defined as where ( Here, ∫ T · dV and ∫ e · dS refer to volume and surface integrals, respectively.
where ( ,ℎ , We note that the coupling conditions, (15) and (16) , are embedded in the above discretized equations. In particular, the conditions (15) and (16) are discretized as and respectively. Finally, we define  ( , ) as: and Here, the choices of the interior penalty method is provided by . The discretization becomes the symmetric interior penalty Galerkin method (SIPG), when = −1 , the incomplete interior penalty Galerkin method (IIPG), when = 0 , and the non-symmetric interior penalty Galerkin method (NIPG) when = 1 Riviere (2008) . In this study, we set = −1 for the simplicity. The interior penalty parameter, , is a function of the polynomial degree, k , and the characteristic mesh size, h e , which is defined as: where meas( · ) represents a measurement operator, measuring length, area, or volume. Some studies for the optimal choice of is provided in Lee et al. (2019) , Riviere (2008) .
Remark 3. The Neumann boundary condition is naturally applied on the boundary faces that belong to the Neumann boundary domain for both the matrix and fracture domains, ∈  ℎ and ∈ Λ ℎ . The Dirichlet boundary condition, on the other hand, is weakly enforced on the Dirichlet boundary faces, ∈  ℎ , for the matrix domain, but strongly enforced on the Dirichlet bondary faces of the fracture domain, be the set of basis functions for the space .  , CG the number of DOF for the CG scalar-valued space. Hence, two mixed function spaces, (1) EG k × CG k and (2) DG k × CG k where k represents the degree of polynomial approximation, are possible, and will be compared in the numerical examples in the next section. The matrix corresponding to the left-hand side of (29) is assembled composing the following blocks: ) , , CG ) , CG , , CG , In a similar way, the right-hand side of (29) gives rise to a block vector of components The resulting block structure is thus where ,ℎ and CG ,ℎ collect the degrees of freedom for matrix and fracture pressure, respectively. Finally, we remark that (owing to (26) ) the case = EG can be equivalently decomposed into a (CG k × DG 0 ) × CG k mixed function space, resulting in: This formulation makes the EG methodology easily implementable in any existing DG codes. Matrices and vectors are built by FEniCS form compiler ( Alnaes et al., 2015 ). The block structure is setup using multiphenics toolbox ( Ballarin and Rozza, 2019 ). Random field of permeability ( k m ) is populated using SciPy package ( Jones et al., 2001 ). , penalty parameter, is set at 1.1 and 1.0 for DG and EG methods, respectively. Remark 4. We note that CG, DG, and EG methods are based on Galerkin method, which could be extended to consider adaptive meshes that contain hanging nodes. In addition, there are various advanced development for each methods to enhance the efficiency, including variable approximation order techniques. Especially, for EG method, an adaptive enrichment, i.e., the piecewise-constant functions only added to the elements where the sharp material discontinuities (e.g., between matrix and fracture domains) are observed, can be developed. However, in our following numerical examples, we focus on the classical form of each methods for the comparison by simulating the proposed mixed dimensional approach for modeling fractures.

Numerical examples
We illustrate the capability of the EG method using seven numerical examples. We begin with an analysis of the error convergence rate between the EG and DG methods to verify the developed block structure in the mixed-dimensional setting. We also investigate the EG performance in modeling the quarter five-spot pattern and handling the fracture tip in the immersed fracture geometry. Next, we test the EG method in a regular fracture network with and without a heterogeneity in matrix permeability input. Lastly, we apply time-dependent problems for two three-dimensional geometries; the first one represents the case where fractures are orthogonal to the axes, and another represents geometry where fractures are given with arbitrary orientations with their interactions in a three-dimensional domain.

Error convergence analysis
To verify the implementation of the proposed block structure utilized to solve the mixed-dimensional model using the EG method, we illustrate the error convergence rate of the EG method and compare this value with the DG method. The example used in this analysis is adapted from Antonietti et al. (2019) . We take Ω = [ 0 , 1 ] 2 , and choose the exact solution in the matrix, Ω, and fracture, Γ = {( , ) ∈ Ω ∶ + = 1} , as: By choosing = , (48) satisfies the system of equations, (1), (9), (15) , and (16) , presented in the methodology section with sink/source terms, g , as follows: All other physical parameters are set to one, and the homogeneous boundary conditions are applied to all boundaries. The geometry used in this analysis and the illustration of the exact solution are presented in Fig. 5 a-b, respectively. We calculate L 2 norm of the difference between the exact solution, p , and approximated solution, p h , and the results are presented in Fig. 5 cd for matrix and fracture domains, respectively. For both matrix and fracture domains, the EG and DG methods provide the expected convergence rate of two and three for polynomial degree approximation, k , of one and two, respectively ( Antonietti et al., 2019;Babuska, 1973 ).

Quarter five-spot example
This numerical example tests the EG method performance in an injection/production setting using five-spot pattern, and we adopt this example from Chave et al. (2018) , Antonietti et al. (2019) . The fivespot pattern, where one injection well is located in the middle and four producers are located at each corner of the square, is commonly used in underground energy extraction Chen et al. (2006) . Due to the symmetry of this geometry, only a quarter of the domain ( Ω = [ 0 , 1 ] 2 ) is simulated. The injection well is located at (0,0), and the production well is located at (1,1). We place the fracture with = 0 . 0005 at Γ = {( , ) ∈ Ω ∶ + = 1} . The geometry, boundary conditions, and mesh with ℎ = 1 . 2 × 10 −1 applied in this analysis are shown in Fig. 6 .
The following source term is applied to the entire matrix domain including the injection and production wells: To investigate the effect of fracture conductivity, we perform two simulations using different fracture conductivity inputs. (i) We choose, = , = 1 , and = 100 for the permeable fracture case. (ii) We assume the fracture is impermeable and set = , = 1 × 10 −2 , and = . All of the remaining physical parameters are set to one. Results of two cases are presented in Fig. 7 a-b for the pressure value and Fig. 7 c-d for the pressure profile along = line. The pressure profile of the permeable fracture case is smoothly decreasing from the injection well towards the production well. The pressure profile of the im- permeable fracture case, on the other hand, illustrates a jump across the fracture interface. These findings comply with the results of the previous studies ( Chave et al., 2018;Antonietti et al., 2019 ). Our results converge to the reference solution as the h is reduced from ℎ = 1 . 2 × 10 −1 to ℎ = 7 . 2 × 10 −2 . Note that ( Antonietti et al., 2019 ) perform these numerical experiments using the second-order DG method with ℎ = 7 . 5 × 10 −2 on polytopic grids ( Antonietti et al., 2019 ). There is no significant difference between the EG and DG results for both h values ( Fig. 7 c-d).

Immersed fracture example
The numerical examples discussed so far contain a fracture that cut through the matrix domain. To test the EG discretization capability in the immersed fracture setting, we adopt this example from Angot et al. (2009) . Since we assume that the fracture tip is substantially small, there is no fluid mass transfer between the matrix and fracture domains across the fracture tip, see A' in Fig. 8 . Hence, the fracture boundary, Λ ℎ , that intersects with the bulk matrix material internal boundary,  0 ℎ , is enforced with no-flow boundary condition, = 0 . In this example, we take the bulk matrix, Ω = [ 0 , 1 ] 2 , and the fracture with = 0 . 01 , Γ = {( , ) ∈ Ω ∶ = 0 . 5 , ≥ 0 . 5} . We perform three simulations using different fracture conductivity inputs; (i) we assume = , = 1 × 10 2 , and = 1 × 10 6 for the permeable fracture case, and its geometry and boundary conditions are presented in Fig. 8 a. (ii) The partially permeable fracture case utilizes = , = 1 × 10 2 , and = 1 × 10 2 . This case geometry and boundary conditions are illustrated in Fig. 8 b. (iii) Impermeable fracture, geometry and boundary conditions are shown in Fig. 8 c, and it uses = , = 1 × 10 −7 , and = 1 × 10 −7 . All other physical parameters are equal to one. Example of mesh that contains = 272 is shown in Fig. 8 d. In Fig. 9 , the values of p m with = 1 , 741 for all cases are presented. The pressure plot along = 0 . 75 is presented in Fig. 10 . The permeable and partially permeable fracture cases illustrate the continuity of the pressure while the impermeable fracture case shows the jump across the fracture interface. Our results converge ( = 272 , = 1 , 741 , and = 6 , 698 ) to the solutions provided by Angot et al. (2009) . Moreover, the EG and DG methods provide similar results. Note that the reference solutions are performed on the finite volume method with 65,536 control volumes.

Regular fracture network example
We increase the complexity of the problem by increasing the number of fractures as shown in Flemisch et al. (2018) . This example, however, is called regular fracture network since all fractures are orthogonal to the axes ( x or y ). The geometry used in this example was utilized by Geiger et al. (2013) for analyzing multi-rate dual-porosity model. Details for model geometry and boundary conditions are shown in Fig. 11 a. We set = for all the matrix domain, Ω = [ 0 , 1 ] 2 , and = 1 × 10 −4 for all fractures, Γ. Two fracture conductivity inputs are used; (i) we choose = 1 × 10 4 and = ×10 4 for the permeable fracture case.
(ii) For the impermeable fracture case, we assume = 1 × 10 −4 , and = 1 × 10 −4 . All of the remaining physical parameters are set to one. Example of mesh that contains = 184 is shown in Fig. 11 d.
The p m results of the permeable and impermeable fractures are presented in Fig. 11 b-c, respectively. The permeable fracture case shows the smooth p m profile while the impermeable fracture case clearly illustrates the jump of p m across the fracture interface. These results are the same as the reference solutions provided by Flemisch et al. (2018) .
Besides the evaluation of pressure results, we also investigate the flux calculated at face A' (see Fig. 11 a) as follows: As presented in Fig. 13 , the difference between each n e case is insignificant, i.e. the different between = 26 , 952 and = 184 cases is less than 1%. Furthermore, there is no difference between the EG and DG methods. The DOF comparison between EG and DG methods is shown in Table 1 . The EG method requires the DOF (in the matrix domain) approximately half of that of the DG method. Note that the DOF in the fracture domain is the same because we discretize the fracture domain using the CG method.

The heterogeneous in bulk matrix permeability example
The numerical examples presented so far only consider an homogeneous matrix permeability value. In this section, we examine the capability of the EG method in handling the discontinuity not only between the fracture and matrix domains but also within the matrix domain (e.g. Nick and Matthai, 2011b ) by employing the heterogeneous permeability in the bulk matrix. We adapt the quarter five-spot as in Section 3.2 and regular fracture network as in Section 3.4 . We choose the finest mesh from both examples to test the EG method capability compared to that of the DG method in handling the sharp discontinuity between the maximum number of interfaces. The geometries, boundary conditions, and input parameters are utilized as Sections 3.2 and 3.4 except for the k m value.
In this study, = , k m value is randomly provided value for each cell. We will distinguish in particular two different cases, named low k m case and high k m case in the following. The low k m case is characterized by a log-normal distribution with average ̄ = 1 . 0 , variance var ( ) = 40 , limited to minimum value min = 1 . 0 × 10 −2 , and maximum value max = 1 . 0 × 10 2 . The high k m case uses ̄ = 30 . 0 , var ( ) = 90 , min = 1 . 0 × 10 −1 , and max = 1 . 5 × 10 2 . This heterogeneous fields for both examples are populated using SciPy package ( Jones et al., 2001 ) as shown in Figs. 14 and 15 for low k m and high k m cases, respectively.

Quarter five-spot example
The results for the low k m case with permeable ( = 1 and = 100 ) and impermeable ( = 1 × 10 −2 and = ) fractures are illustrated in Fig. 16 . In general, the p m profile from the two cases are more disperse than the homogeneous k m setting. The EG and DG methods provide similar results, with = 6 , 568 and ℎ = 7 . 2 × 10 −2 . The discussions  regarding permeable and impermable fracture settings are provided as follows: 1. Low k m with permeable fracture result illustrates the approximately smooth p m solution because the ̄ = 1 . 0 is equal to . The plot along = line, as expected, shows p m gradually decreases from the injection well to the production well. This result complies with that of the homogeneous k m setting. 2. Low k m with impermeable fracture result and the plot along = line exhibit a little jump of p m across the fracture interface. This behavior is different from the homogeneous k m setting since min = , which lead to the less permeability contrast between the fracture and matrix domains.
The results for the high k m case with permeable and impermeable fractures are presented in Fig. 17 . Using = 6 , 568 , the EG and DG results are approximately the same. The observations concerning the fracture permeability are presented below: 1. High k m with permeable fracture displays a jump across the fracture domain because k m around the fracture on the plotting line is higher than ; as a result, the fracture interface acts like a flow barrier. The plot along = line supports this observation as p m jumps across the fracture interface. 2. High k m with impermeable fracture illustrates a huge jump across the fracture interface, as can be observed from the p m plot along = line, since k m is much higher than . The p m variation is less pronounced compared to the low k m one, see Figs. 7 b (homogeneous) and 16 b (heterogeneous), because the fluid flow is blocked by the fracture interface (sharp material discontinuity).

Regular fracture network example
The pressure results of the low k m case between the permeable ( = 1 × 10 4 and = ×10 4 ) and impermeable ( = 1 × 10 −4 and = 1 × 10 −4 ) fractures are illustrated in Fig. 18 . Similar to the fivespot example, the EG and DG results are similar with = 26 , 952 , and p m results are more dispersive than the homogeneous k m setting as shown in 1. Low k m with permeable fracture illustrates the fracture dominate flow regime, even though ̄ = 1 . 0 , k m min is set to 1 . 0 × 10 −2 , which is much less than . This setting reduces the impact of the matrix domain. p m plot along = 0 . 7 and = 0 . 5 lines support this observation by showing the high pressure gradient in the matrix domain, but p m becomes much less varied when entering the fracture domain. 2. Low k m with impermeable fracture also presents the fracture dominate flow regime because = 1 × 10 −4 , which is less than min = 1 . 0 × 10 −2 . Therefore, the flow in the matrix is blocked by the fracture domain. The plot along ( = 0 . 0 , = 0 . 1) to ( = 0 . 9 , = 1 . 0) line illustrates jumps across the fracture domain, which is supporting our observation.
The results for the high k m case with permeable and impermeable fractures are presented in Fig. 19 . With = 26 , 952 , the EG and DG re-sults are approximately the same. The observations concerning the fracture permeability are presented below: 1. High k m with permeable fracture shows that the matrix domain gains more momentum comparing to the low k m case. p m plot along = 0 . 7 and = 0 . 5 lines illustrates also approximately linear reduction along the matrix domain, while pressure is almost constant in the fracture domain. 2. High k m with impermeable fracture clearly presents the fracture domain dominate the flow because is much less than the k m . Hence, the flow is blocked by the fractures. The plot along ( = 0 . 0 , = 0 . 1) to ( = 0 . 9 , = 1 . 0) line support this observation by illustrating multiple jumps across the fracture domain and no p m variation inside each matrix block.

The heterogeneous and anisotropy in bulk matrix permeability example
This section illustrates the comparison between the p m solutions of EG and DG methods using the heterogeneous and anisotropic k m . In con- trast with the previous example, we utilize the coarsest mesh ( = 184 ) from the regular fracture network example. The matrix permeability heterogeneity, k m , is generated with the same specification as the low k m case in the previous section. However, in this study, the off-diagonal terms of (13) are not zero, and the k m of each element is defined as follows: The generated heterogeneous field is shown in Fig. 20 a-b for both diagonal and off diagonal terms. The pressure results of the low k m case between the permeable ( = 1 × 10 4 and = ×10 4 ) and impermeable ( = 1 × 10 −4 and = 1 × 10 −4 ) fractures are presented in Fig. 21 . The results of EG and DG methods are approximately similar for both fracture permeability settings. These results illustrate that the EG method captures the discontinuities and allow using a coarse mesh to maintain accuracy.

The time-dependent problems
From this section, we consider time-dependent problem where c and mi are nonzero (1) -(8) . Besides, we extend the spatial domain to threedimensional space to further illustrate the applicability of the proposed EG method. Here, the first example considers the geometry containing only the orthogonal fractures to the axes ( x, y , or z ), and the second example assumes the geometry with arbitrary orientated natural fractures. Note that we, here, present only the results of the EG method. The results of the DG method are comparable to those of the EG method.
In this example, we consider a three-dimensional domain, which is an analog of the two-dimensional case presented in Section 3.5.2 . This domain contains a set of well-interconnected fractures and meshed with tetrahedral and triangular elements (for fractures) with = 9 , 544 . The fracture geometry is based on the example in Berre et al. (2020) . The details of geometry with initial and boundary conditions illustrated in Fig. 22 a. Here, we consider two different scenarios: case i) permeable fracture case with = 1 × 10 6 and = 1 × 10 6 ; and case ii) impermeable fracture case with = 1 × 10 −12 , and = 1 × 10 −12 . For both cases, we employ a permeable porous medium by setting = [ 1 . 0 , 1 . 0 , 0 . 1 ] , and all of the remaining physical parameters are set to one for the simplicity. The temporal domain is given as = [0 , 100] where an uniform time step size Δ = 1 .
In Fig. 22 b-c, the numerical results of the p m , v m , and v f for the permeable (case i) and impermeable fracture cases (case ii) at = 100 are presented. A mere visual examination of these results already shows that the fracture permeability controls the flow field. For the permeable fracture case, the velocity at the corner of block B is higher than that on the opposite corner in block A as the fractures are closer to the open corner point in block B. For the impermeable fracture case the velocity at the corner of block B is lower than that on the opposite corner in block A since block A is larger than block B, and it can support flow for a longer time.
Similar behaviors of the pressure values are observed in Fig. 23

Algrøyna outcrop example
The final example is a three-dimensional fracture network built based on an outcrop map in Algrøyna, Norway ( Fumagalli et al., 2019 ). The model has a size of 850 × 1400 × 600 m and contains 52 intersecting fractures (the model is described in detail in Berre et al., 2020 ). The finite element mesh is discretized by tetrahedral (rock matrix) with = 163 , 575 and triangular elements (fractures) with = 329 , 080 . See Fig. 24 a for more details. As shown in this figure Dirichlet boundary conditions are applied on two edges of the model to represent an injection and a production well. All other boundaries are considered as no-flow. The rock matrix and the fractures are considered permeable: = 1 × 10 2 , m 2 = 1 × 10 2 m 2 , and = [ 1 . 0 , 1 . 0 , 0 . 1 ] m 2 . All of the remaining physical parameters are set to one, and = [0 , 1 × 10 6 ] sec using an uniform Δt n of 1 × 10 5 sec. Fig. 24 b and c show the simulation results including the pressure field, the pressure iso-surfaces, and velocity vectors in the rock matrix at = 500 , 000 sec. The pressure profiles along a line between two opposite corners of the model at different times are plotted in Fig. 24 d. This example illustrates the applicability of the presented EG method for a complex three-dimensional fracture network with arbitrary orientations.

Conclusion
This study presents the EG discretization for solving a single-phase fluid flow in the fractured porous media using the mixed-dimensional approach. Our proposed method has been tested against several pub-lished benchmarks and subsequently assessed its performance in the test cases with the heterogeneous and anisotropic matrix permeability. Our results illustrate that the pressure solutions resulted from the EG and DG method, with the same mesh size, are approximately similar. Furthermore, the EG method enjoys the same benefits as the DG method; for instance, preserves local and global conservation for fluxes, can handle discontinuity within and between the subdomains, and has the optimal error convergence rate. However, it has much fewer degrees of freedom compared to that of the DG method in its classical form. We note that this comparison can vary based on advanced developments of each method, e.g., a hybridized discontinuous Galerkin method or variable approximation orders. Besides, the results of the time-dependent problem for a three-dimensional geometry highlight the importance of correctly capturing the discontinuities with conductivity values, from barriers to highly-conductive fractures, present in geological media. This work can be extended to multiphysics problems, e.g., poroelastic and transport phenomena, and general form of the mixed-dimensional abstraction, i.e., coupled between d and − dimensionality, where d and n are any integers and − ≥ 0 .

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.