Projection-based Embedded Discrete Fracture Model (pEDFM)

Accepted Manuscript

Projection-based Embedded Discrete Fracture Model (pEDFM) Matei T¸ene, Sebastian B.M. Bosma, Mohammed Saad Al Kobaisi, Hadi Hajibeygi PII: DOI: Reference:

S0309-1708(17)30099-4 10.1016/j.advwatres.2017.05.009 ADWR 2849

To appear in:

Advances in Water Resources

Received date: Revised date: Accepted date:

30 January 2017 12 May 2017 13 May 2017

Please cite this article as: Matei T¸ene, Sebastian B.M. Bosma, Mohammed Saad Al Kobaisi, Hadi Hajibeygi, Projection-based Embedded Discrete Fracture Model (pEDFM), Advances in Water Resources (2017), doi: 10.1016/j.advwatres.2017.05.009

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Highlights • The Projection-based Embedded Discrete Fracture Model (pEDFM) is

CR IP T

developed. • pEDFM is able to accurately represent fractures of any conductivity contrasts (with respect to the matrix rock), in a strictly mass-conservative formulation.

• It automatically reduces to Discrete Fracture Model (DFM) formulation

AN US

when fractures are at the interfaces.

• pEDFM is tested on 2D and 3D test cases with increasing degree of heterogeneities and fracture density.

• Systematic studies are conducted to elicit the sensitivity of pEDFM towards fracture-matrix conductivity contrast and grid resolution, while

AC

CE

PT

ED

M

comparing with the classic EDFM and unstructured DFM methods.

1

ACCEPTED MANUSCRIPT

Projection-based Embedded Discrete Fracture Model (pEDFM)

a Delft

CR IP T

Matei T ¸ enea,∗, Sebastian B.M. Bosmaa , Mohammed Saad Al Kobaisib , Hadi Hajibeygia University of Technology, P.O. Box 5048, 2600 GA Delft, the Netherlands. Petroleum Institute, P.O. Box 2533, Abu Dhabi, United Arab Emirates.

b The

Abstract

AN US

This work presents a new discrete fracture model, namely the Projection-based

Embedded Discrete Fracture Model (pEDFM). Similar to the existing EDFM approach, pEDFM constructs independent grids for the matrix and fracture domains, and delivers strictly conservative velocity fields. However, as a significant step forward, it is able to accurately model the effect of fractures with

M

general conductivity contrasts relative to the matrix, including impermeable flow barriers. This is achieved by automatically adjusting the matrix transmissibility field, in accordance to the conductivity of neighboring fracture networks,

ED

alongside the introduction of additional matrix-fracture connections. The performance of pEDFM is investigated extensively for two- and three-dimensional

PT

scenarios involving single-phase as well as multiphase flows. These numerical experiments are targeted at determining the sensitivity of the model towards the fracture position within the matrix control volume, grid resolution and the

CE

conductivity contrast towards the matrix. The pEDFM significantly outperforms the original EDFM and produces results comparable to those obtained when using DFM on unstructured grids, therefore proving to be a flexible model

AC

for field-scale simulation of flow in naturally fractured reservoirs. Keywords: fracture modelling, embedded discrete fracture model, fractured ∗ Corresponding

author Email addresses: [email protected] (Matei T ¸ ene), [email protected] (Sebastian B.M. Bosma), [email protected] (Mohammed Saad Al Kobaisi), [email protected] (Hadi Hajibeygi)

Preprint submitted to Advances in Water Resources

May 15, 2017

ACCEPTED MANUSCRIPT

reservoirs, heterogeneous geological properties, flow in porous media

1

1. Introduction Accurate and efficient simulation of flow through subsurface formations is

3

essential for effective engineering operations (including production, storage op-

4

timization and safety assessments).

5

properties, the target geological formations often contain complex networks of

6

naturally-formed or artificially-induced fractures, with a wide range of conduc-

7

tivity properties. Given their significant impact on flow patterns, the accurate

8

representation of these lower-dimensional structural features is paramount for

9

the quality of the simulation results [1].

CR IP T

2

AN US

Alongside their intrinsic heterogeneous

Discrete Fracture Models (DFM) reduce the dimensionality of the problem

11

by constraining the fractures, as well as any inhibiting flow barriers, to lie at

12

the interfaces between matrix rock cells [2, 3, 4]. Then, local grid refinements

13

are applied, where a higher level of detail is necessary, leading to a discrete

14

representation of the flow equations on, sometimes complex, unstructured grids

15

[5, 6, 7, 8, 9]. Although the DFM approach has been extended to include

16

complex fluids and rock physics – e.g. compositional displacements [10, 11] and

17

geomechanical effects [12] – its reliance on complex computational grids may

18

raise important challenges in real-field applications.

PT

ED

M

10

This has led to the emergence of models which make use of non-conforming

20

grids w.r.t. fracture-matrix connections, such as eXtended Finite Element Meth-

21

ods (XFEM) [13] and Embedded Discrete Fracture Models (EDFM) [14, 15].

22

This paper focuses on the latter, which are especially appealing due to their

23

intrinsic ability to deliver mass-conservative flux fields. To this end, the lower-

AC

CE

19

24

dimensional structural features with relatively small lengths (i.e. fully contained

25

in a single fine-scale matrix cell) are first homogenized, by altering the effective

26

permeability of their support rock [16]. Then, the remaining fracture networks

27

are discretized on separate numerical grids, defined independently from that of

28

the matrix [17, 18]. A comprehensive comparison between DFM and EDFM,

3

ACCEPTED MANUSCRIPT

29

along with other fracture models, is performed in [19]. The EDFM has been applied to reservoirs containing highly-conductive frac-

31

tures with complex geometrical configurations, while considering compositional

32

fluid physics [20] and plastic and elastic deformation [21]. It has seen used as an

33

upscaling technique [22, 23] and was successfully paired with multiscale methods

34

for efficient flow simulation [24, 25, 26, 27]. However, the experiments presented

35

in this paper show that, in its current formulation, the model is not suitable in

36

cases when the fracture permeability lies below that of the matrix. In addition,

37

even when fractures coincide with the interfaces of matrix cells, the existing

38

EDFM formulation still allows for independent flow leakage (i.e. disregarding

39

the properties of the fracture placed between neighbouring matrix cells). To

40

resolve these important limitations, this work proposes a new formulation for

41

embedded fracture approaches, namely, the projection-based EDFM (pEDFM).

42

The pEDFM is shown to successfully accommodate to lower-dimensional struc-

43

tural features with a wide range of permeability contrasts towards the matrix.

44

This includes highly conductive fractures and flow barriers with small apertures,

45

relative to the reservoir scale, which allows their representation as 2D plates.

46

For the remainder of the paper, they will be referred to, simply, as fractures,

47

regardless of their conductive properties.

ED

M

AN US

CR IP T

30

The devised pEDFM formulation retains the geometric flexibility of the clas-

49

sic EDFM procedure. More specifically, once the fracture and matrix grids are

50

independently defined, and the cross-media communication points are identi-

51

fied, pEDFM adjusts the matrix-matrix and fracture-matrix transmissibilities

CE

PT

48

in the vicinity of fracture networks. This ensures that the conductivity of the

53

fracture networks, which can be several orders of magnitude below or above

54

that of the matrix, are automatically taken into account when constructing the

AC

52

55

flow patterns. Finally, when fractures are explicitly placed at the interfaces of

56

matrix cells, pEDFM automatically provides identical results to DFM.

57

The paper is structured as follows. First, the governing equations are pre-

58

sented and discretized according to the pEDFM approach. Then, a series of

59

numerical experiments are presented, targeted at validating pEDFM by com4

ACCEPTED MANUSCRIPT

paring it to DFM (i.e. using unstructured grids with fractures being confined at

61

the interfaces) and EDFM. The sensitivity of pEDFM with respect to fracture

62

position and orientation, grid resolution and the conductivity contrast towards

63

the matrix is studied extensively. Finally, the results are summarized, conclu-

64

sions are drawn and possible directions for future work are discussed.

65

2. pEDFM formulation

CR IP T

60

In order to accommodate fractures with a wide range of conductivity con-

67

trasts towards the matrix, pEDFM extends the classic EDFM discretization of

68

the governing flow equations by automatically scaling the matrix-matrix con-

69

nections in the vicinity of fracture networks. At the same time, additional

70

fracture-matrix connections are added in order to keep the system of equations

71

well-posed in all possible scenarios. This is explained in detail in the following

72

subsections.

73

2.1. Governing equations

77

ED

dia, without compositional effects, can be written as  m ∂(φρi si ) − ∇ · (ρi λi · ∇p) = Qm + [ρi q]mf ∂t

PT

76

M

The mass-conservation equations for isothermal Darcy flow in fractured me-

74

75

AN US

66

for the matrix (superscript 

(1)

on Ωf ⊂ Rn−1

(2)

) and

f ∂(φρi si ) − ∇ · (ρi λi · ∇p) = Qf + [ρi q]f m ∂t

CE

78

m

on Ωm ⊂ Rn

for the fracture (superscript f ) spatial domains. Here, φ is the rock porosity, p

80

the pressure, while si , λi and ρi are the phase saturation, mobility and density,

AC

79

81

respectively. The q mf and q f m stand for the cross-media connections, while Qm

82

and Qf are source terms, e.g. due to perforating wells, capillary and gravity

83

effects, in the matrix and fracture domain, respectively.

5

ACCEPTED MANUSCRIPT

84

2.2. Discretization In order to solve the coupled system of Eqs. (1)-(2), independent grids

86

are generated for the rock and fracture domains (See Fig. 1). This approach

87

alleviates complexities related to grid generation, since, unlike in DFM, fractures

88

do not need to be confined to the interfaces between matrix grid cells.

AN US

CR IP T

85

(a) Matrix grid

(b) Fracture grid

Figure 1: In pEDFM, independent grids are defined separately for the matrix and fracture

M

domains.

The advection term from Eqs. (1)-(2) is defined for each (matrix-matrix

90

and fracture-fracture) grid interface, following the two-point-flux approxima-

91

tion (TPFA) finite volume discretization of the flux Fij between each pair of

92

neighbouring cells i and j as

PT

ED

89

93

Aij ¯ dij λij

Fij = Tij (pi − pj ).

(3)

is the transmissibility, Aij is the interfacial area, dij is

Here, Tij =

95

¯ ij is the effective fluid mobility at the the distance between the cell centers, λ

96

interface between i and j (absolute permeabilities are harmonically averaged,

AC

CE 94

97

98

99

100

fluid properties are upwinded [28]). The fracture-matrix coupling terms are modelled similar to [15, 24], i.e., for

matrix cell i (with volume Vi ) connected to a fracture cell f (of area Af ), Z mf Fif = qif dV = Tif (pf − pi ) (4) Vi

6

ACCEPTED MANUSCRIPT

101

and Ff i =

102

Z

Af

qff im dA = Tf i (pi − pf ),

(5)

where Tif = CIif λif = Tf i is the cross-media transmissibility. In addition, λif

104

is the effective fluid mobility at the interface between matrix and fracture (just

105

as before, absolute permeabilities are harmonically averaged, fluid properties

106

are upwinded), while the CIif is the conductivity index, defined as CIif =

107

Sif , hdiif

CR IP T

103

(6)

where Sif is the surface area of the connection (to be further specified below)

109

and hdiif is the average distance between the points contained in the rock control

110

111

AN US

108

volume Vi and the fracture surface Af [15, 24], i.e., Z 1 dif dvi , hdiif = Vi Ωi

(7)

where dif stands for the distance between finite volume dvi and fracture plate.

113

Appendix A gives an analytical method for its computation on 2D structured

114

grids.

115

2.2.1. EDFM

M

112

Consider the fractured medium from Fig. 2, which is discretized on a struc-

117

tured grid. Let Aif be the area of intersection between fracture cell f and

118

matrix volume i (highlighted in yellow in Fig. 2). The classical EDFM formu-

119

lation [15, 24] defines the transmissibility as

PT

ED

116

Tif =

CE

120

2Aif λif hdiif

(8)

where, in this case, Sif = 2Aif for computing CI in Eq. (6) and λif is the

122

effective cross-media mobility. The transmissibility of the matrix-matrix con-

123

nections in the neighborhood of the fracture (between control volumes i and j,

124

k, respectively) are left unmodified from their TPFA finite volumes form, i.e

AC

121

125

Tij =

Aij λij ∆x

Tik =

Aik λik . ∆y

(9)

126

where Aij , Aik are the areas and λij , λik the effective mobilities of the corre-

127

sponding matrix interfaces. 7

ACCEPTED MANUSCRIPT

128

2.2.2. pEDFM

CR IP T

k A if ⊥ X A if ⊥ Y

A if

i

AN US

j

f

Figure 2: Illustration of pEDFM on a 2D structured grid. The matrix cells highlighted in yellow are connected directly to the fracture, as defined in the classic EDFM. The cells highlighted in orange take part in the additional non-neighbouring connections between fracture

M

and matrix grid cells, as required by pEDFM.

The paper proposes an extension to the EDFM formulation, by modifying

130

the matrix-matrix and fracture-matrix in the vicinity of fractures. This enables

131

the development of a general embedded discrete fracture modeling approach

132

(pEDFM), applicable in cases with any conductivity contrast between fractures

133

and matrix. To this end, first a set of matrix-matrix interfaces is selected, such

134

that they define a continuous projection path of each fracture network on the

135

matrix domain (highlighted in red on the right side of Fig. 2). While, devising

CE

PT

ED

129

a generic algorithm for the construction of these paths lies outside the scope of

137

this paper, it is important to ensure their continuity for each fracture network.

138

Consider fracture cell f intersecting matrix volume i on an n-dimensional

139

structured grid over a surface area, Aif . Let Aif ⊥xe be its corresponding projec-

140

tions on the path, along each dimension, e = 1, . . . , n (depicted in red on the left

141

side of Fig. 2). Also, let ie be the matrix control volumes which reside on the op-

142

posite side of the interfaces affected by the fracture cell projections (highlighted

AC

136

8

ACCEPTED MANUSCRIPT

143

in orange in Fig. 2). Then, the following transmissibilities are defined Tif =

144

Aif λif , hdiif

T ie f =

Aif ⊥xe λi f hdiie f e

and

Tiie =

Aiie − Aif ⊥xe λiie , ∆xe (10)

where Aiie are the areas of the matrix interfaces hosting the fracture cell pro-

146

jections and λif , λie f , λiie are effective fluid mobilities between the correspond-

147

ing cells. Notice that the projected areas, Aif ⊥xe , are eliminated from the

148

matrix-matrix transmissibilities and, instead, make the object of stand-alone

149

connections between the fracture and the non-neighbouring (i.e. not directly

150

intersected) matrix cells ie . Also, the matrix-matrix connectivity Tiie will be

151

eventually zero if the fracture elements (belonging to one or multiple fractures)

152

cross through the entire matrix cell i.

AN US

CR IP T

145

Finally, note that, for fractures that are explicitly confined to lie along the

154

interfaces between matrix cells, the pEDFM formulation, as given in Eq. (10),

155

naturally reduces to the DFM approach on unstructured grids, while the EDFM

156

does not.

M

153

Given the above TPFA finite-volume discretization of the advection and

158

source terms from Eqs. (1)-(2), after applying backward Euler time integration,

159

the coupled system is linearized with the Newton-Raphson scheme and solved

160

iteratively.

161

3. Numerical results

PT

ED

157

This section presents the results of numerical experiments of single- and

163

two-phase incompressible flow through two- and three-dimensional fractured

164

media. Their aim is to validate pEDFM, whose formulation was presented in

165

the previous section, and study its sensitivity to fracture position, grid resolution

AC

CE

162

166

and fracture-matrix conductivity contrast, respectively. The reference solution

167

for these studies is obtained on a fully resolved grid, i.e. where the size of each

168

cell is equal to the fracture aperture. This allows the following model error

169

measurement,

170

kk1 = Ncoarse

PN 1

|p0f ine − pcoarse | Ncoarse 9

(11)

ACCEPTED MANUSCRIPT

where Ncoarse is the number of grid cells used by pEDFM and p0f ine is the

172

corresponding fully-resolved pressure, interpolated to the coarse scale, if neces-

173

sary. Some of the experiments were repeated for the classic EDFM, as well as

174

unstructured DFM, for comparison purposes.

CR IP T

171

175

For simplicity, but without loss of generality, the flow in these experiments

176

is driven by Dirichlet boundary conditions, instead of injection and production

177

wells, while capillary and gravity effects are neglected.

Finally, the simulations were performed using the DARSim 1 in-house sim-

178

ulator, using a sequentially implicit strategy for the multiphase flow cases.

180

3.1. pEDFM validation

AN US

179

8 7 6 5

M

4

AC

ED

1 0

(b) log10 (k) on the 11 × 11 pEDFM grid

PT

(a) log10 (k) on the 1001 × 1001 fully

0.5

1

1

0.5

0.5

0 0

2

resolved grid

CE

1

3

0.5

1 0.5 1 0

(c) Fully resolved pressure

1 0 0

0.5 0.5

1

0

(d) EDFM pressure

1 0.5

0 0

0.5

1

0

(e) pEDFM pressure

Figure 3: Fully-resolved (c) EDFM (d) and pEDFM (e) pressure solutions in a homogeneous reservoir containing a +-shaped highly-conductive fracture network (top).

10

ACCEPTED MANUSCRIPT

In order to validate pEDFM as a fine-scale model suitable to accommo-

182

date fractures with a wide range of permeabilities, a 2D homogeneous domain

183

(k m = 1) is considered, having a +-shaped fracture network, located in the mid-

184

dle. In order to drive the incompressible single-phase flow, Dirichlet boundary

185

conditions with non-dimensional pressure values of p = 1 and p = 0 are imposed

186

on the left and right boundaries of the domain, respectively, while the top and

187

bottom sides are subject to no-flow conditions.

CR IP T

181

As shown in Fig. 3, the study is first conducted in a scenario where the

189

fractures are 8 orders of magnitude more conductive than the matrix. The

190

reference solution, in this case, is computed on a 1001×1001 structured cartesian

191

grid. From the bottom of Fig. 3, it is clear that both EDFM and pEDFM, on a

192

coarser 11 × 11 domain, can reproduce the behavior of the flow as dictated by

193

the highly-conductive embedded fracture network.

AN US

188

As shown in Fig. 4, the same experiment was rerun for the case where the

195

fracture permeability lies 8 orders of magnitude below that of the host matrix.

196

The results expose the limitations of EDFM, where the impermeable fractures

197

are simply by-passed by the flow through the (unaltered) matrix, resulting in

198

a pressure field corresponding to a reservoir with homogeneous (non-fractured)

199

permeability. On the other hand, through its new formulation, pEDFM is able

200

to reproduce the effect of the inhibiting flow barrier (see bottom of Fig. 4),

201

confirming its applicability to this case.

PT

ED

M

194

These experiments confirm that pEDFM is a suitable extension of EDFM to

203

a wider range of geological scenarios, being able to reproduce the correct flow

CE

202

behaviour in the presence of both high and low permeable fractures, embedded

205

in the porous matrix.

AC

204

206

3.2. Sensitivity to the fracture position within the grid cell

207

Given that pEDFM typically operates on much coarser grids than the fully

208

resolved case, it is of interest to elicit its sensitivity to the fracture position

209

within the host grid cell. To this end, the +-shaped fracture configuration is

210

considered; the reference solution is computed on a 37 × 37 (i.e., 2187 × 2187) 11

ACCEPTED MANUSCRIPT

0 -1 -2 -3

CR IP T

-4 -5 -6 -7 -8

(b) log10 (k) on the 11 × 11 pEDFM

resolved grid

grid

1

1

0.5

0.5

AN US

(a) log10 (k) on the 1001 × 1001 fully

1

0.5

1 0 0

0.5 0.5

0.5

0

(c) Fully resolved pressure

0.5

0 0

0.5

1 0

M

1

1

1

0 0

(d) EDFM pressure

0.5 1

0

(e) pEDFM pressure

ED

Figure 4: Fully-resolved (c) EDFM (d) and pEDFM (e) pressure solutions in a homogeneous reservoir containing a +-shaped nearly impermeable flow barrier (top).

cell grid, while pEDFM grid operates at 10 × 10 resolution.

PT

211

From Fig. 3, it appears that in the case when the fracture network is highly

213

conductive, the horizontal fracture is the one that dictates the flow. Conse-

214

quently, successive simulations are conducted for both EDFM and pEDFM,

215

while moving the horizontal fracture from top to bottom, as shown in Fig. 5.

216

Their accuracy is measured using Eq. (11).

AC

CE

212

217

The results show that EDFM is more accurate when fractures are placed

218

at the cell center, rather than when they are close to the interface. However,

219

once the fracture coincides with the interface, EDFM connects it to both matrix

220

cells (each, with a CI calculated using Sif = Aif in Eq. (6), instead of 2Aif as

221

was the case in Eq. (8)), thus explaining the abrupt dip in error. In contrast,

12

ACCEPTED MANUSCRIPT

8

4

1

CR IP T

0

10 -2

0

AN US

pEDFM EDFM

center

interface

center

Position

M

Figure 5: Sensitivity of pEDFM to the position of highly conductive fractures, embedded within the matrix grid cells. To this end, the horizontal fracture of the +-shaped network is successively moved from top to bottom over a 2 grid cell window (top), while monitoring the

ED

pressure mismatch towards the corresponding fully resolved simulation (bottom).

the pEDFM error attains its peak when fractures are placed at the cell centers

223

and does not exhibit any jumps over the interface. The error of both methods

224

lies within similar bounds (still pEDFM is more accurate) showing that they

225

applicable to the case when fractures are highly conductive. The consistent

CE

PT

222

aspect of pEDFM is that, its results for the case when fractures coincide with

227

the matrix interfaces, its results are identical to the DFM method, while –as

228

explained before– this is not the case for EDFM.

AC

226

229

When the network is nearly impermeable, the location of the vertical fracture

230

is critical to the flow (Fig. 4). As such, for the purposes of the current exper-

231

iment, it will be shifted from left to right, as shown in Fig. 6. The resulting

232

error plot shows a dramatic increase for EDFM, when compared to Fig. 6, due

13

ACCEPTED MANUSCRIPT

0

-4

5

CR IP T

-8

10 -2

pEDFM EDFM

4

AN US

3

2

1

0

center

interface

center

Position

M

Figure 6: Sensitivity of pEDFM to the position of nearly impermeable fractures, embedded within the matrix grid cells. The vertical fracture of the +-shaped network is successively moved from left to right over a 2 grid cell window (top), while monitoring the pressure mis-

ED

match towards the corresponding fully resolved simulation (bottom).

to its inability to handle fractures with conductivities that lie below that of the

234

matrix. pEDFM, on the other hand shows a similar behavior and error range

235

as was observed in the case with highly conductive fractures, i.e., it retains its

236

high accuracy.

CE

PT

233

These results show a promising trend for pEDFM, which is able to maintain

238

reasonable representation accuracy of the effect of the embedded fractures. The

239

slight increase in error for fractures placed near the matrix cell centers may be

240

mitigated by employing moderate local grid refinements.

241

3.3. Sensitivity to the grid resolution

AC

237

242

Another important factor in assessing the quality of an embedded fracture

243

model is its order of accuracy with respect to the grid resolution. A series of 14

ACCEPTED MANUSCRIPT

nested matrix grids for the +-shaped fracture test case of Figs. 3 and 4 was

245

constructed. The number of cells over each axis is gradually increased using

246

Nx = Ny = 3i formula, where i = 2, 3, . . . , up to a reference grid resolution,

247

where i = 7. At the same time, the fracture grid is also refined accordingly such

248

that its step size approximately matches the one in the matrix, h = ∆x = ∆y.

249

The measure of accuracy for this case is similar to Eq. 11, where, this time, no

250

interpolation is necessary, since the cell centroids are inherited from one level

251

to another in the nested grid hierarchy.

CR IP T

244

For a better comparison, alongside pEDFM and EDFM, the same sequence

253

of geological scenarios was simulated using DFM on a 2D unstructured grid [2],

254

where the number of triangles was tweaked to match N = Nx × Ny as closely

255

as possible and without imposing any preferential grid refinement around the

256

fractures.

AN US

252

The results of this study, in the case when the fractures are highly conduc-

258

tive, are depicted in Fig. 7. It follows that all three methods experience a linear

259

decay in error with increasing grid resolution. The three error snapshots, which

260

were taken when Nx = Ny = 36 (or h = 0.0015), show that the pressure mis-

261

match is mainly concentrated around the tips of the horizontal fracture, which

262

represent the network’s inflow and outflow points, respectively. For EDFM,

263

the error decays radially for points further away from these fracture tips. For

264

pEDFM, the contour curves are slightly skewed, depending on the choice be-

265

tween upper and lower matrix interfaces for the fracture projection (both are

266

equally probable since the horizontal fracture crosses the grid cell centroids).

CE

PT

ED

M

257

Finally, for DFM the error distribution shows some heterogeneity, which is a

268

consequence of using unstructured grids in a medium which, except for the

269

neighborhood of the fractures, is homogeneous.

AC

267

270

The scenario when the fracture network is considered almost impermeable

271

can not be properly handled by EDFM, regardless of which grid resolution is

272

used (Fig. 8). This serious limitation is, once again, successfully overcome by

273

using pEDFM, which, similar to DFM, maintains its linear scalability with grid

274

refinement on this challenging test case. The error snapshots depict that, this 15

(a) Error pEDFM

(b) Error EDFM

10 -1

10 -2

ED

10 -4

M

10 -3

pEDFM EDFM DFM

(c) Error DFM

AN US

10 -2

CR IP T

ACCEPTED MANUSCRIPT

Figure 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

PT

three methods, when Nx = Ny = 36 (or h = 0.0015) holds for pEDFM and EDFM, and DFM employs comparable total number of elements.

time, the pressure is inaccurate around the tips, as well as the body, of the

276

vertical barrier. This can be explained by the fact that an embedded model

277

on a coarse grid can have difficulty in placing the sharp discontinuity in the

AC

CE

275

278

pressure field at exactly the right location. Still, the pressure mismatch decays

279

with increasing grid resolution, suggesting that local grid refinements around

280

highly contrasting fractures can benefit pEDFM, in a similar manner to DFM.

281

To conclude, pEDFM shows a similar convergence behavior, in terms of grid

282

resolution, to the widely used DFM approach. This confirms that, in order to

16

ACCEPTED MANUSCRIPT

10 -2 9 8 7 5 4 3 2 1

(a) Error pEDFM

(b) Error EDFM

10 -2

10 -1

10 -2

ED

10 -4

M

10 -3 pEDFM EDFM DFM

(c) Error DFM

AN US

10 -1

CR IP T

6

Figure 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

PT

three methods, when Nx = Ny = 36 (or h = 0.0015) holds for pEDFM and EDFM, and DFM employs comparable total number of elements.

diminish the model representation error, moderate local grid refinements can be

284

applied near fractures.

285

3.4. Sensitivity to the fracture-matrix conductivity contrast

AC

CE

283

286

This last sensitivity study is aimed at determining the response of pEDFM

287

while changing the conductivity contrast between the +-shaped fracture network

288

(k f = 10−8 , . . . , 108 ) and the matrix (k m = 1). To this end, a coarse grid

289

resolution of Nx = Ny = 35 was used and the resulting pressure was compared

17

ACCEPTED MANUSCRIPT

290

to that from the reference case, where Nx = Ny = 37 , using Eq. 11. The results are depicted in Fig. 9 and are in line with the conclusions from

292

previous sections. Namely, for fracture log-permeabilities on the positive side

293

of the spectrum, the results of EDFM and pEDFM are in agreement. As the

294

permeability contrast passes 5 orders of magnitude, the pressure error plateaus,

295

since, at beyond this stage, the fractures are the main drivers of the flow. How-

296

ever, for fracture permeabilities close to or below that of the matrix, the error of

297

EDFM increases dramatically. pEDFM, on the other hand is able to cope with

298

these scenarios, due to its formulation, its behavior showing an approximately

299

symmetric trend, when compared to that of the positive side of the axis.

AN US

CR IP T

291

The snapshots on the top of Fig. 9 , taken for lower, similar and higher frac-

301

ture permeabilities w.r.t. the matrix, show the error in the pressure produced

302

by pEDFM. It is clear that the model inaccuracy is concentrated around the

303

tips of fractures which actively influence the flow. Also note that there is a

304

small error even in the case when k f = k m , since the pEDFM discretization

305

(Section 2) is slightly different than that of a homogeneous reservoir. Of course,

306

this result is only presented here for academic purposes – in realistic scenarios,

307

when the contrast is not high enough, such fractures can be homogenized into

308

the matrix field.

ED

M

300

This concludes the sensitivity studies conducted in this paper. The test

310

cases presented in the following subsections are meant to test the applicability

311

of pEDFM to more complex 2D and 3D fractured media.

312

3.5. Time-lapse 2D multiphase results

CE

PT

309

This set of experiments is aimed at determining the performance of pEDFM

314

in multiphase flow scenarios on 2D porous media with increasingly complex

AC

313

315

fracture geometries and heterogeneities.

316

3.5.1. Homogeneous matrix

317

pEDFM is first applied in an incompressible 2-phase flow scenario through

318

a 2D homogeneous domain which is crossed by a set of fractures with hetero-

18

ACCEPTED MANUSCRIPT

10 -3

10 -3

10 -3

9

9

8

8

8

7

7

7

6

6

6

5

5

5

4

4

3

3

2

2

1

1

(a) Error kf = 10−5

CR IP T

9

4 3 2 1

(b) Error kf = 1

(c) Error kf = 105

AN US

10 -1

10 -2

10 -4

-7

-6

-5

-4

-3

-2

-1

0

1

2

3

4

5

6

7

8

ED

-8

pEDFM EDFM

M

10 -3

Figure 9: Sensitivity of pEDFM to the fracture-matrix conductivity contrast on the +-shaped fracture network case with a grid resolution of 35 × 35 . The sequence of plots on the top show

319

PT

pressure error maps for pEDFM, when kf = 10−5 , 1 and 105 , respectively.

geneous properties, as shown in Fig. 10. The boundary conditions are similar to those used for previous experiments, namely Dirichlet with non-dimensional

321

values of p = 1 and p = 0 on the left and right edges, respectively, while the top

322

and bottom sides are subject to no-flow conditions.

AC

CE 320

323

The low permeable fractures inhibit the flow, leaving only two available

324

paths: through the middle of the domain and along the bottom boundary. As

325

can be seen in the time-lapse saturation maps presented in Fig. 10, the front,

326

indeed, respects these embedded obstacles. The injected fluid is mostly directed

327

through the permeable X-shaped network and surpasses the vertical barrier, in 19

CR IP T

ACCEPTED MANUSCRIPT

(d) 0.2 PVI

(e) 0.3 PVI

(f) 0.4 PVI

(g) 0.5 PVI

M

(c) 0.1 PVI

(b) Pressure

AN US

(a) log10 (k)

Figure 10: Fracture permeabilities (top-left), pressure field (top-right) and time-lapse sat-

ED

uration results (bottom) on a 2D test case with homogeneous matrix conductivity, under incompressible 2-phase flow conditions.

328

the lower right part of the domain, on its way to the production boundary. This result reinforces the conclusion that the conservative pressure field ob-

330

tained using pEDFM leads to transport solutions which honour a wide range of

331

matrix-fracture conductivity contrasts.

CE

PT

329

332

3.5.2. Heterogeneous matrix

AC

333

The following experiment compares the behaviour of EDFM and pEDFM

334

for simulating 2-phase incompressible flow through a 2D porous medium with

335

heterogeneous (i.e. patchy) matrix permeability (Fig. 11). The interplay be-

336

tween the large- (matrix-matrix) and small-scale (fracture-matrix) conductivity

337

contrasts raises additional numerical challenges [29] and is a stepping stone in

338

assessing the model’s applicability to realistic cases. 20

CR IP T

ACCEPTED MANUSCRIPT

(b) log10 (kf )

AN US

(a) log10 (km )

M

EDFM

(c) Pressure

(d) 0.1 PVI

(e) 0.3 PVI

(f) 0.5 PVI

CE

PT

ED

pEDFM

(g) Pressure

(h) 0.1 PVI

(i) 0.3 PVI

(j) 0.5 PVI

Figure 11: Heterogeneous matrix and fracture permeability maps (top), pressure and time-

AC

lapse saturation results for EDFM (middle) and pEDFM (bottom) on a 2D densely fractured test case, under incompressible 2-phase flow conditions.

339

The embedded fracture map used for this test case (top of Fig. 11) is based

340

on the Brazil I outcrop from [30, 31]. The conductivities of the fractures forming

341

the North-West to South-East diagonal streak, were set to 10−8 , thus creating

21

ACCEPTED MANUSCRIPT

an impermeable flow barrier across the domain (noticeable in dark blue on the

343

top-right of Fig. 11). For the rest of the fractures, permeabilities were randomly

344

drawn from a log-uniform distribution supported on the interval [10−8 , 108 ].

345

Finally, similar to previous experiments, fixed pressure boundary conditions

346

p = 1 and p = 0 are set on the left and right edges, respectively, while the top

347

and bottom sides are subject to no-flow conditions.

CR IP T

342

The middle row of plots from Fig. 11 show the pressure field and time-lapse

349

saturation results obtained using EDFM. Note that the injected fluid is allowed

350

to bypass the diagonal flow barrier, towards the production boundary. This,

351

once again shows the limitation of EDFM, which is only able to capture the

352

effect of fractures with permeabilities higher than the matrix. However, by

353

disregarding flow barriers, EDFM delivers an overly optimistic and nonrealistic

354

production forecast.

AN US

348

In contrast to EDFM, the pressure field obtained using pEDFM shows sharp

356

discontinuities (bottom-left of Fig. 11). The accompanying saturation plots con-

357

firm that the injected phase is confined by the diagonal barrier and forced to flow

358

through the bottom of the domain, thus significantly delaying its breakthrough

359

towards the production boundary.

ED

M

355

These results confirm that pEDFM outperforms to EDFM, due to its appli-

361

cability in cases with complex and dense fracture geometries and in the presence

362

of matrix heterogeneities.

363

3.6. Comparison between 3D pEDFM and unstructured DFM

PT

360

Finally, a test case on a 3D domain containing 3 layers of fractures, stacked

365

along the Z axis (Fig. 12) is conducted. The top layer is a vertically extruded

366

version of the 2D fracture map from Fig. 10. The second layer consists of a single

AC

CE 364

367

fracture network whose intersecting plates have highly heterogeneous properties.

368

Finally, the third layer is represented by 3 large individual plates, with a cluster

369

of small parallel fractures packed in between.

370

In this scenario, the incompressible single-phase flow is driven from the left

371

boundary, where the pressure is set to the non-dimensional value of p = 1, 22

ACCEPTED MANUSCRIPT

8 6 4 2

CR IP T

0 -2 -4 -6 -8

(a) log10 (kf ) layer 1

(b) log10 (kf ) layer 2

(c) log10 (kf ) layer 3

AN US

Figure 12: 3D porous medium containing 3 layers of fractures with heterogeneous properties.

towards the right, where p = 0, while all the other boundaries of the domain are

373

subject to no-flow conditions. No other source terms are present and gravity

374

and capillary effects are neglected. The results of pEDFM, on a matrix grid

375

with Nx = Ny = Nz = 100 and a total of 23381 fracture cells, are compared to

376

those obtained using DFM on an unstructured grid (Fig. 13), where the number

377

of tetrahedra (matrix) and triangles (fractures) were chosen to approximately

378

match the degrees of freedom on the structured grid.

AC

CE

PT

ED

M

372

Figure 13: Fracture-conforming unstructured grid constructed by DFM for the 3D test case.

379

380

The two pressure fields are plotted in Fig. 14 using iso surfaces for equidistant values, and are in good agreement, for decision-making purposes. 23

AN US

CR IP T

ACCEPTED MANUSCRIPT

(a) DFM pressure

(b) pEDFM pressure

Figure 14: Comparison between the pressures obtained using pEDFM and unstructured DFM (using similar grid resolutions) for a 3D incompressible single-phase test case with 3 layers of heterogeneous fractures.

M

This last numerical experiment shows that pEDFM has good potential for

381

field-scale application.

383

4. Conclusions and future work

ED

382

In this paper, a novel Projection-based Embedded Discrete Fracture Model

385

(pEDFM) was devised, for flow simulation through fractured porous media. It

386

inherits the grid flexibility of the classic EDFM approach. However, unlike its

387

predecessor, its formulation allows it to capture the effect of fracture conduc-

388

tivities ranging from highly permeable networks to inhibiting flow barriers.

CE

PT

384

The new model was validated on 2D and 3D test cases, while studying its

AC

389

390

sensitivity towards fracture position within a matrix cell, grid resolution and

391

the cross-media conductivity contrast. The results show that pEDFM is scal-

392

able and able to handle dense and complex fracture maps with heterogeneous

393

properties in single-, as well as multiphase flow scenarios. Finally, its results

394

on structured grids were found comparable to those obtained using the DFM

24

ACCEPTED MANUSCRIPT

395

approach on unstructured, fracture-conforming meshes. In conclusion, pEDFM is a flexible model, its simple formulation recommend-

397

ing it for implementation in next-generation simulators for fluid flow through

398

fractured porous media.

CR IP T

396

Possible directions for future work include the extension of pEDFM to un-

399

structured grids.

Furthermore, systematic benchmarking studies (including

401

CPU time comparisons) between pEDFM and alternative mass-conservative

402

DFM approaches on unstructured grids can provide valuable insights for its

403

application and performance in real-field cases.

404

Acknowledgements

AN US

400

The financial support of PI/ADNOC is acknowledged. Also, many thanks

406

are due towards the members of the DARSim (Delft Advanced Reservoir Simu-

407

lation) research group of TU Delft, for the useful discussions during the devel-

408

opment of the pEDFM modeling approach.

409

Appendix A. Average distance calculation

ED

M

405

The computation of the average distance between a matrix control volume

411

and a fracture surface, which appears in Eqs. (6) and (10), may involve numer-

412

ical integration for arbitrarily shaped cells. For 2D structured grids, however,

413

analytical formulas were given in [24] for a few specific fracture orientations.

414

This subsection, an adaptation of Chapter 2.3.2 from [16], provides a general

CE

PT

410

415

procedure to handle fracture lines with arbitrary orientation. The interfaces of each cell intersected by a fracture are extended until they

417

intersect the fracture line, resulting in four right triangles with surfaces A1 to

AC

416

418

A4 , as shown in Fig. A.15. Then, given the average distance between each

419

triangle and its hypotenuse, hdi1 to hdi4 , as (see [24]),

420

Lxi · Lyi hdii = p , 3 Lxi 2 + Lyi 2 25

(A.1)

ACCEPTED MANUSCRIPT

4

2

2

4

2

1

3 (a)

2

1

3

(b)

1

(c)

2

1

AN US

4

3

CR IP T

1

(d)

(e)

Figure A.15: 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

M

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).

where Lxi and Lyi are the lengths of the axis-aligned sides of triangle i, the

422

average distance between grid cell i and fracture line f is obtained,

ED

421

hdiif =

PT

423

A1 hdi1 + A2 hdi2 − A3 hdi3 − A4 hdi4 . A1 + A2 − A3 − A4

(A.2)

Note that no modification is required to the formula in the case when frac-

425

tures lie outside the cell, i.e. for the non-neighbouring connections from Eq. 10.

CE

424

In addition, this procedure can be directly applied to 3D structured grids where

427

fractures are extruded along the Z axis, while a generalization for fracture plates

428

with any orientation is the subject of future research.

429

References

AC

426

430

431

[1] B. Berkowitz, Characterizing flow and transport in fractured geological media: A review, Adv. Water Resour. 25 (8–12) (2002) 861–884. 26

ACCEPTED MANUSCRIPT

432

[2] M. Karimi-Fard, L. J. Durlofsky, K. Aziz, An efficient discrete fracture

433

model applicable for general purpose reservoir simulators, SPE J. 9 (2)

434

(2004) 227–236. [3] V. Reichenberger, H. Jakobs, P. Bastian, R. Helmig, A mixed-dimensional

436

finite volume method for two-phase flow in fractured porous media, Adv.

437

Water Resour. 29 (7) (2006) 1020–1036.

CR IP T

435

[4] R. Ahmed, M. G. Edwards, S. Lamine, B. A. Huisman, M. Pal, Three-

439

dimensional control-volume distributed multi-point flux approximation

440

coupled with a lower-dimensional surface fracture model, J. Comput. Phys.

441

303 (2015) 470–497.

AN US

438

[5] S. K. Matth¨ ai, A. A. Mezentsev, M. Belayneh, Finite element node-centered

443

finite-volume two-phase-flow experiments with fractured rock represented

444

by unstructured hybrid-element meshes, SPE Reserv. Eval. Eng. 10 (2007)

445

740–756.

M

442

[6] M. Sahimi, R. Darvishi, M. Haghighi, M. R. Rasaei, Upscaled unstructured

447

computational grids for efficient simulation of flow in fractured porous me-

448

dia, Transport in porous media 83 (1) (2010) 195–218.

ED

446

[7] N. Schwenck, B. Flemisch, R. Helmig, B. Wohlmuth, Dimensionally re-

450

duced flow models in fractured porous media: crossings and boundaries,

451

Computational Geosciences 19 (1) (2015) 1219–1230.

PT

449

[8] A. Tatomir, A. Szymkiewicz, H. Class, R. Helmig, Modeling two phase flow

453

in large scale fractured porous media with an extended multiple interacting

454

continua method, Computer Modeling in Engineering and Sciences 77 (2)

AC

CE

452

455

(2011) 81.

456

[9] M. Karimi-Fard, L. J. Durlofsky, A general gridding, discretization, and

457

coarsening methodology for modeling flow in porous formations with dis-

458

crete geological features, Advances in Water Resources 96 (6) (2016) 354–

459

372. 27

ACCEPTED MANUSCRIPT

460

[10] J. Moortgat, A. Firoozabadi, Three-phase compositional modeling with

461

capillarity in heterogeneous and fractured media, SPE J 18 (2013) 1150 –

462

1168. [11] J. Moortgat, M. Amooie, M. Soltanian, Implicit finite volume and dis-

464

continuous galerkin methods for multicomponent flow in unstructured 3d

465

fractured porous media, Advances in Water Resources 96 (2016) 389–404.

CR IP T

463

[12] T. T. Garipov, M. Karimi-Fard, H. A. Tchelepi, Discrete fracture model for

467

coupled flow and geomechanics, Computational Geosciences 20 (1) (2016)

468

149–160.

AN US

466

469

[13] B. Flemisch, A. Fumagalli, A. Scotti, A review of the xfem-based approx-

470

imation of flow in fractured porous media, in: Advances in Discretization

471

Methods, Springer, 2016, pp. 47–76.

[14] S. H. Lee, C. L. Jensen, M. F. Lough, Efficient finite-difference model for

473

flow in a reservoir with multiple length-scale fractures, SPE J. 3 (2000)

474

268–275.

M

472

[15] L. Li, S. H. Lee, Efficient field-scale simulation of black oil in naturally

476

fractured reservoir through discrete fracture networks and homogenized

477

media, SPE Reserv. Eval. Eng. 11 (2008) 750–758.

PT

478

ED

475

[16] S. Pluimers, Hierarchical fracture modeling, Msc thesis, Delft University of Technology, The Netherlands (2015).

CE

479

[17] D. Karvounis, P. Jenny, Adaptive hierarchical fracture model for enhanced

481

geothermal systems, Multiscale Modeling & Simulation 14 (1) (2016) 207–

482

231.

AC

480

483

[18] R. Deb, P. Jenny, Numerical modeling of flow-mechanics coupling in a frac-

484

tured reservoir with porous matrix, Proceedings of the 41st Workshop on

485

Geothermal Reservoir Engineering Stanford University, Stanford, Califor-

486

nia, February 22-24, SGP-TR-209. (2016) 1–9.

28

ACCEPTED MANUSCRIPT

487

[19] B. Flemisch, I. Berre, W. Boon, A. Fumagalli, N. Schwenck, A. Scotti,

488

I. Stefansson, A. Tatomir, Benchmarks for single-phase flow in fractured

489

porous media, arXiv:1701.01496 (2017). [20] A. Moinfar, A. Varavei, K. Sepehrnoori, R. T. Johns, Development of an

491

efficient embedded discrete fracture model for 3d compositional reservoir

492

simulation in fractured reservoirs, SPE J. 19 (2014) 289–303.

CR IP T

490

[21] J. H. Norbeck, M. W. McClure, J. W. Lo, R. N. Horne, An embedded

494

fracture modeling framework for simulation of hydraulic fracturing and

495

shear stimulation, Computational Geosciences 20 (1) (2016) 1–18.

AN US

493

496

[22] A. Fumagalli, L. Pasquale, S. Zonca, S. Micheletti, An upscaling procedure

497

for fractured reservoirs with embedded grids, Water Resources Research

498

52 (8) (2016) 6506–6525.

[23] A. Fumagalli, S. Zonca, L. Formaggia, Advances in computation of local

500

problems for a flow-based upscaling in fractured reservoirs, Mathematics

501

and Computers in Simulation 137 (2017) 299–324.

M

499

[24] H. Hajibeygi, D. Karvounis, P. Jenny, A hierarchical fracture model for

503

the iterative multiscale finite volume method, J. Comput. Phys. 230 (2011)

504

8729–8743.

PT

ED

502

[25] M. T ¸ ene, M. S. Al Kobaisi, H. Hajibeygi, Algebraic multiscale method

506

for flow in heterogeneous porous media with embedded discrete fractures

507

(F-AMS), J. Comput. Phys. 321 (2016) 819–845.

CE

505

[26] S. Shah, O. Moyner, M. T ¸ ene, K.-A. Lie, H. Hajibeygi, The multiscale

509

restriction smoothed basis method for fractured porous media (F-MsRSB),

510

J. Comput. Phys. 318 (2016) 36–57.

AC

508

511

[27] M. T ¸ ene, M. Al Kobaisi, H. Hajibeygi, Multiscale projection-based Embed-

512

ded Discrete Fracture Modeling approach (F-AMS-pEDFM), in: ECMOR

513

XV-15th European Conference on the Mathematics of Oil Recovery, 2016.

29

ACCEPTED MANUSCRIPT

514

[28] Z. Chen, Reservoir simulation: mathematical techniques in oil recovery, SIAM, 2007.

515

[29] H. Hamzehpour, M. Asgari, M. Sahimi, Acoustic wave propagation in het-

517

erogeneous two-dimensional fractured porous media, Physical Review E

518

93 (6) (2016) 063305.

519

[30] G.

Bertotti,

K.

(NE

Bisdom,

Brazil),

Fracture

CR IP T

516

patterns

520

Fm.

521

be07fe95-417c-44e9-8c6a-d13f186abfbb.

in

the

Jandeira

http://data.4tu.nl/repository/uuid:

[31] K. Bisdom, G. Bertotti, H. M. Nick, The impact of different aperture distri-

523

bution models and critical stress criteria on equivalent permeability in frac-

524

tured rocks, Journal of Geophysical Research: Solid Earth 121 (5) (2016)

525

4045–4063.

AC

CE

PT

ED

M

AN US

522

30