Projection-based Embedded Discrete Fracture Model (pEDFM)

This work presents a new discrete fracture model, namely the Projection-based Embedded Discrete Frac- ture Model (pEDFM). Similar to the existing EDFM approach, pEDFM constructs independent grids for the matrix and fracture domains, and delivers strictly conservative velocity ﬁelds. However, as a signiﬁcant step forward, it is able to accurately model the effect of fractures with general conductivity contrasts relative to the matrix, including impermeable ﬂow barriers. This is achieved by automatically adjusting the matrix transmissibility ﬁeld, in accordance to the conductivity of neighboring fracture networks, alongside the introduction of additional matrix-fracture connections. The performance of pEDFM is investigated extensively for two- and three-dimensional scenarios involving single-phase as well as multiphase ﬂows. These numerical experiments are targeted at determining the sensitivity of the model towards the frac- ture position within the matrix control volume, grid resolution and the conductivity contrast towards the matrix. The pEDFM signiﬁcantly outperforms the original EDFM and produces results comparable to those obtained when using DFM on unstructured grids, therefore proving to be a ﬂexible model for ﬁeld-scale simulation of ﬂow in naturally fractured reservoirs.


Introduction
Accurate and efficient simulation of flow through subsurface formations is essential for effective engineering operations (including production, storage optimization and safety assessments). Alongside their intrinsic heterogeneous properties, the target geological formations often contain complex networks of naturallyformed or artificially-induced fractures, with a wide range of conductivity properties. Given their significant impact on flow patterns, the accurate representation of these lower-dimensional structural features is paramount for the quality of the simulation results ( Berkowitz, 2002 ).
Discrete Fracture Models (DFM) reduce the dimensionality of the problem by constraining the fractures, as well as any inhibiting flow barriers, to lie at the interfaces between matrix rock cells ( Ahmed et al., 2015 ;Karimi-Fard et al., 2004 ;Reichenberger et al., 2006 ). Then, local grid refinements are applied, where a higher level of detail is necessary, leading to a discrete representation of the flow equations on, sometimes complex, unstructured grids ( Karimi-Fard and Durlofsky, 2016 ;Matthäi et al., 2007 ;Sahimi et al., 2010 ;Schwenck et al., 2015 ;Tatomir et al., 2011 ). Although the DFM approach has been extended to include complex fluids and rock physics -e.g. compositional displacements ( Moortgat et al., 2016 ;Moortgat and Firoozabadi, 2013 ) and geomechanical effects ( Garipov et al., 2016 ) -its reliance on complex computational grids may raise important challenges in realfield applications.
This has led to the emergence of models which make use of non-conforming grids w.r.t. fracture-matrix connections, such as eXtended Finite Element Methods (XFEM, see Flemisch et al., 2016 ) and Embedded Discrete Fracture Models (EDFM, introduced in Lee et al., 20 0 0 ; Li and Lee, 20 08 ). This paper focuses on the latter, which are especially appealing due to their intrinsic ability to deliver mass-conservative flux fields. To this end, the lowerdimensional structural features with relatively small lengths (i.e. fully contained in a single fine-scale matrix cell) are first homogenized, by altering the effective permeability of their support rock ( Pluimers, 2015) . Then, the remaining fracture networks are discretized on separate numerical grids, defined independently from that of the matrix ( Deb and Jenny, 2016 ;Karvounis and Jenny, 2016 ). A comprehensive comparison between DFM and EDFM, along with other fracture models, is performed by Flemisch et al. (2017) .
The EDFM has been applied to reservoirs containing highlyconductive fractures with complex geometrical configurations, while considering compositional fluid physics ( Moinfar et al., 2014 )  and plastic and elastic deformation ( Norbeck et al., 2016 ). It has seen used as an upscaling technique  and was successfully paired with multiscale methods for efficient flow simulation ( Hajibeygi et al., 2011 ;Shah et al., 2016 ;Ţ ene et al., 2016a ;2016b ). However, the experiments presented in this paper show that, in its current formulation, the model is not suitable in cases when the fracture permeability lies below that of the matrix. In addition, even when fractures coincide with the interfaces of matrix cells, the existing EDFM formulation still allows for independent flow leakage (i.e. disregarding the properties of the fracture placed between neighbouring matrix cells). To resolve these important limitations, this work proposes a new formulation for embedded fracture approaches, namely, the projection-based EDFM (pEDFM). The pEDFM is shown to successfully accommodate to lower-dimensional structural features with a wide range of permeability contrasts towards the matrix. This includes highly conductive fractures and flow barriers with small apertures, relative to the reservoir scale, which allows their representation as 2D plates. For the remainder of the paper, they will be referred to, simply, as fractures , regardless of their conductive properties.
The devised pEDFM formulation retains the geometric flexibility of the classic EDFM procedure. More specifically, once the fracture and matrix grids are independently defined, and the crossmedia communication points are identified, pEDFM adjusts the matrix-matrix and fracture-matrix transmissibilities in the vicinity of fracture networks. This ensures that the conductivity of the fracture networks, which can be several orders of magnitude below or above that of the matrix, are automatically taken into account when constructing the flow patterns. Finally, when fractures are explicitly placed at the interfaces of matrix cells, pEDFM automatically provides identical results to DFM.
The paper is structured as follows. First, the governing equations are presented and discretized according to the pEDFM approach. Then, a series of numerical experiments are presented, targeted at validating pEDFM by comparing it to DFM (i.e. using unstructured grids with fractures being confined at the interfaces) and EDFM. The sensitivity of pEDFM with respect to fracture position and orientation, grid resolution and the conductivity contrast towards the matrix is studied extensively. Finally, the results are summarized, conclusions are drawn and possible directions for future work are discussed.

pEDFM formulation
In order to accommodate fractures with a wide range of conductivity contrasts towards the matrix, pEDFM extends the classic EDFM discretization of the governing flow equations by automatically scaling the matrix-matrix connections in the vicinity of fracture networks. At the same time, additional fracture-matrix connections are added in order to keep the system of equations wellposed in all possible scenarios. This is explained in detail in the following subsections.

Governing equations
The mass-conservation equations for isothermal Darcy flow in fractured media, without compositional effects, can be written as for the matrix (superscript m ) and for the fracture (superscript f ) spatial domains. Here, φ is the rock porosity, p the pressure, while s i , λ i and ρ i are the phase saturation, mobility and density, respectively. The q mf and q fm stand for the cross-media connections, while Q m and Q f are source terms, e.g. due to perforating wells, capillary and gravity effects, in the matrix and fracture domain, respectively.

Discretization
In order to solve the coupled system of Eqs. (1) and (2) , independent grids are generated for the rock and fracture domains (See Fig. 1 ). This approach alleviates complexities related to grid generation, since, unlike in DFM, fractures do not need to be confined to the interfaces between matrix grid cells.
The advection term from Eqs. (1) and (2) is defined for each (matrix-matrix and fracture-fracture) grid interface, following the two-point-flux approximation (TPFA) finite volume discretization of the flux F ij between each pair of neighbouring cells i and j as (3) d ij is the distance between the cell centers, λ i j is the effective fluid mobility at the interface between i and j (absolute permeabilities are harmonically averaged, fluid properties are upwinded ( Chen, 2007 )).
The fracture-matrix coupling terms are modelled similar to Hajibeygi et al. (2011) ;Li and Lee (2008) , i.e., for matrix cell i (with volume V i ) connected to a fracture cell f (of area A f ), where T i f = CI i f λ i f = T f i is the cross-media transmissibility. In addition, λ if is the effective fluid mobility at the interface between matrix and fracture (just as before, absolute permeabilities are harmonically averaged, fluid properties are upwinded), while the CI if is the conductivity index, defined as where S if is the surface area of the connection (to be further specified below) and d if is the average distance between the points contained in the rock control volume V i and the fracture surface A f ( Hajibeygi et al., 2011 ;Li and Lee, 2008 ), i.e., where d if stands for the distance between finite volume dv i and fracture plate. Appendix A gives an analytical method for its computation on 2D structured grids.

EDFM
Consider the fractured medium from Fig. 2 , which is discretized on a structured grid. Let A if be the area of intersection between fracture cell f and matrix volume i (highlighted in yellow in Fig. 2 ). The classical EDFM formulation ( Hajibeygi et al., 2011 ;Li and Lee, 2008 ) defines the transmissibility as is the effective cross-media mobility. The transmissibility of the where A ij , A ik are the areas and λ ij , λ ik the effective mobilities of the corresponding matrix interfaces.

pEDFM
The paper proposes an extension to the EDFM formulation, by modifying the matrix-matrix and fracture-matrix in the vicinity of fractures. This enables the development of a general embedded discrete fracture modeling approach (pEDFM), applicable in cases with any conductivity contrast between fractures and matrix. To this end, first a set of matrix-matrix interfaces is selected, such that they define a continuous projection path of each fracture network on the matrix domain (highlighted in red on the right side of Fig. 2 ). While, devising a generic algorithm for the construction of these paths lies outside the scope of this paper, it is important to ensure their continuity for each fracture network.
Consider fracture cell f intersecting matrix volume i on an ndimensional structured grid over a surface area, A if . Let A i f ⊥ x e be its corresponding projections on the path, along each dimension, e = 1 , . . . , n (depicted in red on the left side of Fig. 2 ). Also, let i e be the matrix control volumes which reside on the opposite side of the interfaces affected by the fracture cell projections (highlighted in orange in Fig. 2 ). Then, the following transmissibilities are defined where A ii e are the areas of the matrix interfaces hosting the fracture cell projections and λ if , λ i e f , λ ii e are effective fluid mobilities between the corresponding cells. Notice that the projected areas, A i f ⊥ x e , are eliminated from the matrix-matrix transmissibilities and, instead, make the object of stand-alone connections between the fracture and the non-neighbouring (i.e. not directly intersected) matrix cells i e . Also, the matrix-matrix connectivity T ii e will be eventually zero if the fracture elements (belonging to one or multiple fractures) cross through the entire matrix cell i . Finally, note that, for fractures that are explicitly confined to lie along the interfaces between matrix cells, the pEDFM formulation, as given in Eq. (10) , naturally reduces to the DFM approach on unstructured grids, while the EDFM does not.
Given the above TPFA finite-volume discretization of the advection and source terms from Eqs. (1) and (2) , after applying backward Euler time integration, the coupled system is linearized with the Newton-Raphson scheme and solved iteratively.

Numerical results
This section presents the results of numerical experiments of single-and two-phase incompressible flow through two-and three-dimensional fractured media. Their aim is to validate pEDFM, whose formulation was presented in the previous section, and study its sensitivity to fracture position, grid resolution and fracture-matrix conductivity contrast, respectively. The reference solution for these studies is obtained on a fully resolved grid, i.e. where the size of each cell is equal to the fracture aperture. This allows the following model error measurement, where N coarse is the number of grid cells used by pEDFM and p f ine is the corresponding fully-resolved pressure, interpolated to the coarse scale, if necessary. Some of the experiments were repeated for the classic EDFM, as well as unstructured DFM, for comparison purposes.
For simplicity, but without loss of generality, the flow in these experiments is driven by Dirichlet boundary conditions, instead of injection and production wells, while capillary and gravity effects are neglected. Finally, the simulations were performed using the DARSim 1 inhouse simulator, using a sequentially implicit strategy for the multiphase flow cases.

pEDFM validation
In order to validate pEDFM as a fine-scale model suitable to accommodate fractures with a wide range of permeabilities, a 2D homogeneous domain ( k m = 1 ) is considered, having a + −shaped fracture network, located in the middle. In order to drive the incompressible single-phase flow, Dirichlet boundary conditions with non-dimensional pressure values of p = 1 and p = 0 are imposed on the left and right boundaries of the domain, respectively, while the top and bottom sides are subject to no-flow conditions.
As shown in Fig. 3 , the study is first conducted in a scenario where the fractures are 8 orders of magnitude more conductive than the matrix. The reference solution, in this case, is computed on a 1001 × 1001 structured cartesian grid. From the bottom of Fig. 3 , it is clear that both EDFM and pEDFM, on a coarser 11 × 11 domain, can reproduce the behavior of the flow as dictated by the highly-conductive embedded fracture network.
As shown in Fig. 4 , the same experiment was rerun for the case where the fracture permeability lies 8 orders of magnitude below that of the host matrix. The results expose the limitations of EDFM, where the impermeable fractures are simply by-passed by the flow through the (unaltered) matrix, resulting in a pressure field corresponding to a reservoir with homogeneous (non-fractured) permeability. On the other hand, through its new formulation, pEDFM is able to reproduce the effect of the inhibiting flow barrier (see bottom of Fig. 4 ), confirming its applicability to this case.
These experiments confirm that pEDFM is a suitable extension of EDFM to a wider range of geological scenarios, being able to reproduce the correct flow behaviour in the presence of both high and low permeable fractures, embedded in the porous matrix.

Sensitivity to the fracture position within the grid cell
Given that pEDFM typically operates on much coarser grids than the fully resolved case, it is of interest to elicit its sensitivity to the fracture position within the host grid cell. To this end, the + −shaped fracture configuration is considered; the reference solution is computed on a 3 7 × 3 7 (i.e., 2187 × 2187) cell grid, while pEDFM grid operates at 10 × 10 resolution. From Fig. 3 , it appears that in the case when the fracture network is highly conductive, the horizontal fracture is the one that dictates the flow. Consequently, successive simulations are conducted for both EDFM and pEDFM, while moving the horizontal fracture from top to bottom, as shown in Fig. 5 . Their accuracy is measured using Eq. (11) .
The results show that EDFM is more accurate when fractures are placed at the cell center, rather than when they are close to the interface. However, once the fracture coincides with the interface, EDFM connects it to both matrix cells (each, with a CI calculated using S i f = A i f in Eq. (6) , instead of 2 A if as was the case in Eq. (8) ), thus explaining the abrupt dip in error. In contrast, the pEDFM error attains its peak when fractures are placed at the cell centers and does not exhibit any jumps over the interface. The error of both methods lies within similar bounds (still pEDFM is more accurate) showing that they applicable to the case when fractures are highly conductive. The consistent aspect of pEDFM is that, its results for the case when fractures coincide with the matrix interfaces, its results are identical to the DFM method, while -as explained before-this is not the case for EDFM.
When the network is nearly impermeable, the location of the vertical fracture is critical to the flow ( Fig. 4 ). As such, for the purposes of the current experiment, it will be shifted from left to right, as shown in Fig. 6 . The resulting error plot shows a dramatic increase for EDFM, when compared to Fig. 6 , due to its inability to handle fractures with conductivities that lie below that of the matrix. pEDFM, on the other hand shows a similar behavior and error  range as was observed in the case with highly conductive fractures, i.e., it retains its high accuracy.
These results show a promising trend for pEDFM, which is able to maintain reasonable representation accuracy of the effect of the embedded fractures. The slight increase in error for fractures placed near the matrix cell centers may be mitigated by employing moderate local grid refinements.

Sensitivity to the grid resolution
Another important factor in assessing the quality of an embedded fracture model is its order of accuracy with respect to the grid resolution. A series of nested matrix grids for the + −shaped fracture test case of Figs. 3 and 4 was constructed. The number of cells over each axis is gradually increased using N x = N y = 3 i for-  Fig. 7. Grid resolution sensitivity of pEDFM, EDFM and DFM on the case with highly conductive fractures. The sequence of plots on the top shows pressure error maps for the three methods, when N x = N y = 3 6 (or h = 0 . 0015 ) holds for pEDFM and EDFM, and DFM employs comparable total number of elements. Fig. 8. Grid resolution sensitivity of pEDFM, EDFM and DFM on the case with nearly impermeable fractures. The sequence of plots on the top shows pressure error maps for the three methods, when N x = N y = 3 6 (or h = 0 . 0015 ) holds for pEDFM and EDFM, and DFM employs comparable total number of elements. mula, where i = 2 , 3 , . . . , up to a reference grid resolution, where i = 7 . At the same time, the fracture grid is also refined accordingly such that its step size approximately matches the one in the matrix, h = x = y . The measure of accuracy for this case is similar to Eq. (11) , where, this time, no interpolation is necessary, since the cell centroids are inherited from one level to another in the nested grid hierarchy.
For a better comparison, alongside pEDFM and EDFM, the same sequence of geological scenarios was simulated using DFM on a 2D unstructured grid ( Karimi-Fard et al., 2004 ), where the number of triangles was tweaked to match N = N x × N y as closely as possible and without imposing any preferential grid refinement around the fractures.
The results of this study, in the case when the fractures are highly conductive, are depicted in Fig. 7 . It follows that all three methods experience a linear decay in error with increasing grid resolution. The three error snapshots, which were taken when N x = N y = 3 6 (or h = 0 . 0015 ), show that the pressure mismatch is mainly concentrated around the tips of the horizontal fracture, which represent the network's inflow and outflow points, respectively. For EDFM, the error decays radially for points further away from these fracture tips. For pEDFM, the contour curves are slightly skewed, depending on the choice between upper and lower matrix interfaces for the fracture projection (both are equally probable since the horizontal fracture crosses the grid cell centroids). Finally, for DFM the error distribution shows some heterogeneity, which is a consequence of using unstructured grids in a medium which, except for the neighborhood of the fractures, is homogeneous.
The scenario when the fracture network is considered almost impermeable can not be properly handled by EDFM, regardless of which grid resolution is used ( Fig. 8 ). This serious limitation is, once again, successfully overcome by using pEDFM, which, similar to DFM, maintains its linear scalability with grid refinement on this challenging test case. The error snapshots depict that, this time, the pressure is inaccurate around the tips, as well as the body, of the vertical barrier. This can be explained by the fact that an embedded model on a coarse grid can have difficulty in placing the sharp discontinuity in the pressure field at exactly the right location. Still, the pressure mismatch decays with increasing grid resolution, suggesting that local grid refinements around highly contrasting fractures can benefit pEDFM, in a similar manner to DFM.
To conclude, pEDFM shows a similar convergence behavior, in terms of grid resolution, to the widely used DFM approach. This confirms that, in order to diminish the model representation error, moderate local grid refinements can be applied near fractures.

Sensitivity to the fracture-matrix conductivity contrast
This last sensitivity study is aimed at determining the response of pEDFM while changing the conductivity contrast between the + −shaped fracture network ( k f = 10 −8 , . . . , 10 8 ) and the matrix ( k m = 1 ). To this end, a coarse grid resolution of N x = N y = 3 5 was used and the resulting pressure was compared to that from the reference case, where N x = N y = 3 7 , using Eq. (11) .
The results are depicted in Fig. 9 and are in line with the conclusions from previous sections. Namely, for fracture logpermeabilities on the positive side of the spectrum, the results of EDFM and pEDFM are in agreement. As the permeability contrast passes 5 orders of magnitude, the pressure error plateaus, since, at beyond this stage, the fractures are the main drivers of the flow. However, for fracture permeabilities close to or below that of the matrix, the error of EDFM increases dramatically. pEDFM, on the Fig. 9. Sensitivity of pEDFM to the fracture-matrix conductivity contrast on the + −shaped fracture network case with a grid resolution of 3 5 × 3 5 . The sequence of plots on the top show pressure error maps for pEDFM, when k f = 10 −5 , 1 and 10 5 , respectively.   other hand is able to cope with these scenarios, due to its formulation, its behavior showing an approximately symmetric trend, when compared to that of the positive side of the axis.
The snapshots on the top of Fig. 9 , taken for lower, similar and higher fracture permeabilities w.r.t. the matrix, show the error in the pressure produced by pEDFM. It is clear that the model inaccuracy is concentrated around the tips of fractures which actively influence the flow. Also note that there is a small error even in the case when k f = k m , since the pEDFM discretization ( Section 2 ) is slightly different than that of a homogeneous reservoir. Of course, this result is only presented here for academic purposes -in realistic scenarios, when the contrast is not high enough, such fractures can be homogenized into the matrix field.
This concludes the sensitivity studies conducted in this paper. The test cases presented in the following subsections are meant to test the applicability of pEDFM to more complex 2D and 3D fractured media.

Time-lapse 2D multiphase results
This set of experiments is aimed at determining the performance of pEDFM in multiphase flow scenarios on 2D porous media with increasingly complex fracture geometries and heterogeneities.

Homogeneous matrix
pEDFM is first applied in an incompressible 2-phase flow scenario through a 2D homogeneous domain which is crossed by a set of fractures with heterogeneous properties, as shown in Fig. 10 . The boundary conditions are similar to those used for previous experiments, namely Dirichlet with non-dimensional values of p = 1 and p = 0 on the left and right edges, respectively, while the top and bottom sides are subject to no-flow conditions.
The low permeable fractures inhibit the flow, leaving only two available paths: through the middle of the domain and along the bottom boundary. As can be seen in the time-lapse saturation maps presented in Fig. 10 , the front, indeed, respects these embedded obstacles. The injected fluid is mostly directed through the permeable X-shaped network and surpasses the vertical barrier, in the lower right part of the domain, on its way to the production boundary.
This result reinforces the conclusion that the conservative pressure field obtained using pEDFM leads to transport solutions which honour a wide range of matrix-fracture conductivity contrasts.

Heterogeneous matrix
The following experiment compares the behaviour of EDFM and pEDFM for simulating 2-phase incompressible flow through a 2D porous medium with heterogeneous (i.e. patchy) matrix permeability ( Fig. 11 ). The interplay between the large-(matrix-matrix) and small-scale (fracture-matrix) conductivity contrasts raises additional numerical challenges ( Hamzehpour et al., 2016 ) and is a stepping stone in assessing the model's applicability to realistic cases.
The embedded fracture map used for this test case (top of Fig. 11 ) is based on the Brazil I outcrop from Bertotti and Bisdom ;Bisdom et al. (2016) . The conductivities of the fractures forming the North-West to South-East diagonal streak, were set to 10 −8 , thus creating an impermeable flow barrier across the domain (noticeable in dark blue on the top-right of Fig. 11 ). For the rest of the fractures, permeabilities were randomly drawn from a loguniform distribution supported on the interval [10 −8 , 10 8 ] . Finally, similar to previous experiments, fixed pressure boundary condi-tions p = 1 and p = 0 are set on the left and right edges, respectively, while the top and bottom sides are subject to no-flow conditions.
The middle row of plots from Fig. 11 show the pressure field and time-lapse saturation results obtained using EDFM. Note that the injected fluid is allowed to bypass the diagonal flow barrier, towards the production boundary. This, once again shows the limitation of EDFM, which is only able to capture the effect of fractures with permeabilities higher than the matrix. However, by disregarding flow barriers, EDFM delivers an overly optimistic and nonrealistic production forecast.
In contrast to EDFM, the pressure field obtained using pEDFM shows sharp discontinuities (bottom-left of Fig. 11 ). The accompanying saturation plots confirm that the injected phase is confined by the diagonal barrier and forced to flow through the bottom of the domain, thus significantly delaying its breakthrough towards the production boundary.
These results confirm that pEDFM outperforms to EDFM, due to its applicability in cases with complex and dense fracture geometries and in the presence of matrix heterogeneities.

Comparison between 3D pEDFM and unstructured DFM
Finally, a test case on a 3D domain containing 3 layers of fractures, stacked along the Z axis ( Fig. 12 ) is conducted. The top layer is a vertically extruded version of the 2D fracture map from Fig. 10 . The second layer consists of a single fracture network whose intersecting plates have highly heterogeneous properties. Finally, the third layer is represented by 3 large individual plates, with a cluster of small parallel fractures packed in between.
In this scenario, the incompressible single-phase flow is driven from the left boundary, where the pressure is set to the nondimensional value of p = 1 , towards the right, where p = 0 , while all the other boundaries of the domain are subject to no-flow conditions. No other source terms are present and gravity and capillary effects are neglected. The results of pEDFM, on a matrix grid with N x = N y = N z = 100 and a total of 23381 fracture cells, are compared to those obtained using DFM on an unstructured grid ( Fig. 13 ), where the number of tetrahedra (matrix) and triangles (fractures) were chosen to approximately match the degrees of freedom on the structured grid.
The two pressure fields are plotted in Fig. 14 using iso surfaces for equidistant values, and are in good agreement, for decisionmaking purposes.
This last numerical experiment shows that pEDFM has good potential for field-scale application.

Conclusions and future work
In this paper, a novel Projection-based Embedded Discrete Fracture Model (pEDFM) was devised, for flow simulation through fractured porous media. It inherits the grid flexibility of the classic EDFM approach. However, unlike its predecessor, its formulation allows it to capture the effect of fracture conductivities ranging from highly permeable networks to inhibiting flow barriers.
The new model was validated on 2D and 3D test cases, while studying its sensitivity towards fracture position within a matrix cell, grid resolution and the cross-media conductivity contrast. The results show that pEDFM is scalable and able to handle dense and complex fracture maps with heterogeneous properties in single-, as well as multiphase flow scenarios. Finally, its results on structured grids were found comparable to those obtained using the DFM approach on unstructured, fracture-conforming meshes.
In conclusion, pEDFM is a flexible model, its simple formulation recommending it for implementation in next-generation simulators for fluid flow through fractured porous media.   A15. Analytical calculation of the average distance between a fracture and a matrix cell on 2D structured grids. Four right triangles are constructed by intersecting the cell's edges with the fracture line. Note that triangles 3 and 4 overlap with triangles 1 or 2 (a,b). When the fracture coincides with the cell diagonal, triangles 3 and 4 have zero area (c). If the fracture is aligned with one of the axes, two rectangles are formed instead (e). Finally, the same procedure is followed when the fracture lies outside the cell (d).
Possible directions for future work include the extension of pEDFM to unstructured grids. Furthermore, systematic benchmarking studies (including CPU time comparisons) between pEDFM and alternative mass-conservative DFM approaches on unstructured grids can provide valuable insights for its application and performance in real-field cases.
connections from Eq. (10) . In addition, this procedure can be directly applied to 3D structured grids where fractures are extruded along the Z axis, while a generalization for fracture plates with any orientation is the subject of future research.