Projection-based Embedded Discrete Fracture Model (pEDFM)

(1)

Projection-based Embedded Discrete Fracture Model (pEDFM)

Tene, Matei; Bosma, Sebastian; Al Kobaisi, Mohammed S.; Hajibeygi, Hadi

DOI

10.1016/j.advwatres.2017.05.009

Publication date

2017

Document Version

Final published version

Published in

Advances in Water Resources

Citation (APA)

Tene, M., Bosma, S., Al Kobaisi, M. S., & Hajibeygi, H. (2017). Projection-based Embedded Discrete

Fracture Model (pEDFM). Advances in Water Resources, 105, 205-216.

https://doi.org/10.1016/j.advwatres.2017.05.009

Important note

To cite this publication, please use the final published version (if applicable).

Please check the document version above.

Copyright

Other than for strictly personal use, it is not permitted to download, forward or distribute the text or part of it, without the consent of the author(s) and/or copyright holder(s), unless the work is under an open content license such as Creative Commons. Takedown policy

Please contact us and provide details if you believe this document breaches copyrights. We will remove access to the work immediately and investigate your claim.

This work is downloaded from Delft University of Technology.

(2)

a

r

t

i

c

l

e

i

n

f

o

Article history:

Received 30 January 2017 Revised 12 May 2017 Accepted 13 May 2017 Available online 15 May 2017

Keywords:

Fracture modelling

Embedded discrete fracture model Fractured reservoirs

Heterogeneous geological properties Flow in porous media

a

b

s

t

r

a

c

t

Thisworkpresentsanewdiscretefracturemodel,namelytheProjection-basedEmbeddedDiscrete Frac-tureModel(pEDFM).SimilartotheexistingEDFMapproach,pEDFMconstructsindependentgridsforthe matrixandfracturedomains,anddeliversstrictlyconservativevelocityfields.However,asasignificant stepforward,itisabletoaccuratelymodeltheeffectoffractureswithgeneralconductivitycontrasts rel-ativetothematrix,includingimpermeableflowbarriers.Thisisachievedbyautomaticallyadjustingthe matrixtransmissibilityfield, inaccordancetotheconductivityofneighboringfracturenetworks, along-sidetheintroductionofadditionalmatrix-fractureconnections.TheperformanceofpEDFMisinvestigated extensivelyfortwo-andthree-dimensionalscenariosinvolvingsingle-phaseaswellasmultiphaseflows. Thesenumericalexperimentsaretargetedatdeterminingthesensitivityofthemodeltowardsthe frac-turepositionwithinthematrixcontrolvolume,gridresolutionandtheconductivitycontrasttowardsthe matrix.ThepEDFMsignificantlyoutperformstheoriginalEDFMandproducesresultscomparabletothose obtainedwhenusingDFMonunstructuredgrids,thereforeprovingtobeaflexiblemodelforfield-scale simulationofflowinnaturallyfracturedreservoirs.

1. Introduction

Accurate and efficient simulation of flow through subsurface formations is essential for effective engineering operations (in-cluding production,storageoptimizationandsafety assessments). Alongside their intrinsic heterogeneous properties,the target ge-ological formations oftencontain complex networksof naturally-formed 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 structuralfeatures is paramountforthe qualityof thesimulation results(Berkowitz,2002).

Discrete Fracture Models (DFM) reduce the dimensionality of the problemby constrainingthefractures, aswell asany inhibit-ing flow barriers, to lie at the interfaces between matrix rock cells(Ahmedetal.,2015; Karimi-Fardetal.,2004;Reichenberger et al., 2006). Then, local grid refinements are applied, where a higher level of detail is necessary, leading to a discrete repre-sentationof theflow equations on,sometimes complex, unstruc-turedgrids(Karimi-FardandDurlofsky,2016;Matthäietal.,2007;

∗ _{Corresponding author.}

E-mail addresses: m.tene@tudelft.nl (M. ¸T ene), s.b.m.bosma@student.tudelft.nl

(S.B.M. Bosma), malkobaisi@pi.ac.ae (M.S. Al Kobaisi), h.hajibeygi@tudelft.nl (H. Ha- jibeygi).

Sahimi et al., 2010; Schwenck etal., 2015; Tatomir et al., 2011).

Although the DFMapproach has been extended to include com-plex ﬂuids and rock physics – e.g. compositional displacements

(Moortgat et al., 2016; Moortgat and Firoozabadi, 2013) and

ge-omechanicaleffects (Garipov et al., 2016) – its reliance on com-plex computational grids mayraise important challengesin real-ﬁeldapplications.

This hasled to the emergence of models which make use of non-conforming grids w.r.t. fracture-matrix connections, such as eXtendedFiniteElementMethods(XFEM,seeFlemischetal.,2016) andEmbeddedDiscreteFractureModels(EDFM,introduced inLee et al., 2000; Li and Lee, 2008). This paper focuses on the lat-ter, which are especially appealing due to their intrinsic ability to deliver mass-conservative flux fields. To this end, the lower-dimensional structural features withrelatively small lengths (i.e. fullycontainedinasinglefine-scalematrixcell)arefirst homog-enized,byalteringtheeffectivepermeabilityoftheirsupportrock

(Pluimers, 2015). Then, the remaining fracture networksare

dis-cretizedon separate numericalgrids, deﬁnedindependentlyfrom that of the matrix (Deb and Jenny, 2016; Karvounis and Jenny, 2016). A comprehensive comparison between DFM and EDFM, alongwithother fracturemodels, isperformedbyFlemischetal. (2017).

The EDFM has been applied to reservoirs containing highly-conductive fractures with complex geometrical conﬁgurations, whileconsideringcompositionalﬂuidphysics(Moinfaretal.,2014)

http://dx.doi.org/10.1016/j.advwatres.2017.05.009

(3)

and plastic and elastic deformation (Norbeck et al., 2016). It hasseen used asan upscaling technique (Fumagalli et al., 2016;

2017) and was successfully paired with multiscale methods for eﬃcientﬂowsimulation (Hajibeygi etal.,2011; Shahetal., 2016;

¸Teneetal.,2016a;2016b).However, theexperimentspresentedin

thispapershow that,inits currentformulation,themodelisnot suitableincaseswhenthefracturepermeabilityliesbelowthatof thematrix.Inaddition,evenwhenfracturescoincidewiththe in-terfacesofmatrixcells,theexistingEDFMformulationstillallows for independent ﬂow leakage (i.e. disregarding the properties of thefractureplacedbetweenneighbouringmatrixcells).Toresolve theseimportantlimitations,thisworkproposesanewformulation for embedded fracture approaches, namely, the projection-based EDFM(pEDFM).ThepEDFMisshowntosuccessfullyaccommodate tolower-dimensionalstructuralfeatureswithawiderangeof per-meabilitycontrasts towards thematrix. Thisincludes highly con-ductivefracturesandﬂowbarrierswithsmallapertures,relativeto thereservoirscale,whichallowstheirrepresentationas2Dplates. Fortheremainderofthepaper,theywillbereferredto,simply,as fractures,regardlessoftheirconductiveproperties.

The devised pEDFMformulationretainsthegeometric flexibil-ityoftheclassicEDFMprocedure.Morespecifically,oncethe frac-tureand matrix grids are independently defined, andthe cross-media communication points are identified, pEDFM adjusts the matrix-matrixand fracture-matrix transmissibilities in the vicin-ityoffracturenetworks. Thisensuresthat theconductivityofthe fracturenetworks, which can be severalorders ofmagnitude be-lowor abovethat ofthe matrix,are automaticallytakeninto ac-countwhenconstructingtheflowpatterns.Finally,whenfractures areexplicitlyplacedattheinterfacesofmatrixcells,pEDFM auto-maticallyprovidesidenticalresultstoDFM.

The paper is structured as follows. First, the governing equa-tions are presented anddiscretized according to the pEDFM ap-proach.Then,aseriesofnumericalexperimentsarepresented, tar-getedatvalidatingpEDFMbycomparingitto DFM(i.e.using un-structured grids with fractures being conﬁned at the interfaces) andEDFM.ThesensitivityofpEDFMwithrespecttofracture posi-tionandorientation,gridresolutionandtheconductivitycontrast towards thematrix is studied extensively. Finally, theresults are summarized,conclusionsaredrawnandpossibledirectionsfor fu-tureworkarediscussed.

2. pEDFMformulation

In orderto accommodate fractureswitha wide rangeof con-ductivitycontraststowards thematrix,pEDFMextendstheclassic EDFMdiscretizationofthegoverningﬂowequationsby automati-callyscalingthematrix-matrixconnections inthevicinityof frac-turenetworks. At the same time, additional fracture-matrix con-nectionsareaddedinordertokeepthesystemofequations well-posedin all possible scenarios. This isexplained in detail in the followingsubsections.

2.1.Governingequations

The mass-conservation equationsforisothermal Darcy ﬂow in fracturedmedia, withoutcompositional effects,can bewritten as

_∂

₍

_φρ

isi

)

∂

t −

∇

·

(

ρ

i

λ

i·

∇

p

)

m =Qm_{+ [}

_ρ

iq]m f on

m⊂ Rn (1) forthematrix(superscriptm₎_and

_∂

₍

_φρ

isi

)

∂

t −

∇

·

(

ρ

i

λ

i·

∇

p

)

f =Qf₊_[

_ρ

iq]f m on

f⊂ Rn−1 (2)

forthefracture(superscriptf₎_spatial_domains._Here,

φ

_is_the_rock

porosity,pthepressure,whilesi,

λ

i and

ρ

i arethe phase

satura-tion,mobilityanddensity,respectively.Theqmf _and_qfm _stand_for

the cross-media connections, while Qm _and _Qf _are _source _terms,

e.g. due to perforating wells,capillary andgravity effects, in the matrixandfracturedomain,respectively.

2.2. Discretization

Inordertosolve thecoupledsystemofEqs.(1)and(2), inde-pendentgridsaregeneratedfortherockandfracturedomains(See

Fig.1).Thisapproachalleviatescomplexitiesrelatedtogrid gener-ation,since,unlikeinDFM,fracturesdonotneedtobeconﬁnedto theinterfacesbetweenmatrixgridcells.

The advection term fromEqs. (1) and (2)is defined for each (matrix-matrix andfracture-fracture) gridinterface, following the two-point-flux approximation (TPFA) finite volume discretization ofthe flux Fij betweeneach pairof neighbouringcellsi and j as

Fi j =Ti j

(

pi− pj

)

. (3)

Here, Ti j= Ai j

di j

¯

λ

i j is thetransmissibility, Aij is theinterfacial area,

dijisthedistancebetweenthecellcenters,

λ

¯i j istheeffectiveﬂuid

mobilityattheinterface betweeni andj(absolute permeabilities are harmonicallyaveraged, ﬂuid properties are upwinded (Chen, 2007)).

The fracture-matrix coupling terms are modelled similar to

Hajibeygietal.(2011);LiandLee(2008),i.e.,formatrixcelli(with

volumeVi)connectedtoafracturecellf(ofareaAf),

Fi f= Vi qm f_{i f} dV=Ti f

(

pf− pi

)

(4) and Ff i= Af qf m_{f i} dA=Tf i

(

pi− pf

)

, (5)

whereTi f=CIi f

λ

i f=Tf iisthecross-mediatransmissibility.In

ad-dition,

λ

if is theeffective ﬂuid mobilityatthe interface between

matrixandfracture(justasbefore,absolutepermeabilitiesare har-monicallyaveraged, ﬂuid propertiesareupwinded), whiletheCIif

istheconductivityindex,deﬁnedas

CIi f= Si f

d

i f,

(6) whereSif isthesurfaceareaoftheconnection(tobefurther

spec-iﬁedbelow) and

d

if isthe averagedistance betweenthe points

containedintherockcontrolvolumeV_iandthefracturesurfaceA_f

(Hajibeygietal.,2011;LiandLee,2008),i.e.,

d

i f= 1 Vi i di f d

v

i, (7)

where dif standsfor the distance betweenﬁnite volume dvi and

fractureplate.AppendixAgivesananalyticalmethodforits com-putationon2Dstructuredgrids.

2.2.1. EDFM

ConsiderthefracturedmediumfromFig.2,whichisdiscretized on a structured grid.Let Aif be the area of intersectionbetween

fracturecellfandmatrixvolumei(highlightedinyellowinFig.2). TheclassicalEDFMformulation(Hajibeygi etal.,2011;LiandLee, 2008)deﬁnesthetransmissibilityas

Ti f=

2Ai f

d

i f

λ

i f

(8) where,inthiscase, S_{i f}=2A_{i f} forcomputingCI inEq.(6)and

λ

if

(4)

Fig. 1. In pEDFM, independent grids are deﬁned separately for the matrix and fracture domains.

Fig. 2. Illustration of pEDFM on a 2D structured grid. The matrix cells highlighted in yellow are connected directly to the fracture, as deﬁned in the classic EDFM. The cells highlighted in orange take part in the additional non-neighbouring connections between fracture and matrix grid cells, as required by pEDFM. (For inter- pretation of the references to color in this ﬁgure, the reader is referred to the web version of the article.)

matrix-matrixconnectionsintheneighborhoodofthefracture (be-tweencontrolvolumesiandj,k,respectively)areleftunmodiﬁed fromtheirTPFAﬁnitevolumesform,i.e

Ti j= Ai j

x

λ

i j Tik= Aik

y

λ

ik. (9)

whereAij, Aik are theareasand

λ

ij,

λ

ik theeffectivemobilities of

thecorrespondingmatrixinterfaces.

2.2.2. pEDFM

The paperproposesan extension totheEDFMformulation,by modifyingthematrix-matrixandfracture-matrixinthevicinityof fractures. This enables the development of a general embedded discrete fracture modeling approach (pEDFM), applicable incases with any conductivitycontrast between fractures and matrix.To this end, ﬁrst a set of matrix-matrix interfaces is selected, such thattheydeﬁneacontinuousprojectionpathofeachfracture net-work onthe matrixdomain (highlightedinred onthe rightside ofFig.2).While,devisingagenericalgorithmfortheconstruction ofthesepathsliesoutsidethescope ofthispaper,itisimportant toensuretheircontinuityforeachfracturenetwork.

Consider fracture cell fintersecting matrix volume i on an n -dimensionalstructuredgridoverasurfacearea,Aif.LetAi f⊥xebeits

corresponding projectionsonthepath,alongeach dimension,e= 1,...,n(depictedinredontheleftsideofFig.2).Also,letiebethe

matrix control volumeswhichreside on the opposite side ofthe interfaces affectedby the fracture cell projections(highlightedin

orangeinFig.2).Then,thefollowingtransmissibilitiesaredeﬁned

Ti f= Ai f

d

i f

λ

i f ,Tief= Ai f⊥xe

d

ief

λ

iefandTiie= Aiie− Ai f⊥xe

xe

λ

iie, (10) whereAiie are theareasofthematrixinterfaceshostingthe

frac-ture cell projections and

λ

_if,

λ

_i_e_f_,

λ

_ii_e are effective ﬂuid mobili-tiesbetweenthecorrespondingcells.Noticethattheprojected ar-eas,Ai f⊥xe, are eliminated from the matrix-matrix

transmissibili-tiesand, instead,makethe objectofstand-alone connections be-tween the fracture andthe non-neighbouring(i.e. not directly in-tersected)matrixcellsie.Also, thematrix-matrixconnectivityTiie

willbeeventually zeroifthefractureelements(belongingtoone ormultiplefractures)crossthroughtheentirematrixcelli.

Finally,notethat,forfracturesthatareexplicitlyconﬁnedtolie alongtheinterfacesbetweenmatrixcells,thepEDFMformulation, asgiveninEq.(10),naturallyreducestotheDFMapproachon un-structuredgrids,whiletheEDFMdoesnot.

GiventheaboveTPFAﬁnite-volumediscretizationofthe advec-tionandsource termsfromEqs.(1) and(2), afterapplying back-wardEulertimeintegration,thecoupledsystemislinearizedwith theNewton–Raphsonschemeandsolvediteratively.

3. Numericalresults

This section presents the results of numerical experiments of single- and two-phase incompressible ﬂow through two- and three-dimensionalfracturedmedia.TheiraimistovalidatepEDFM, 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 solutionfor thesestudies is obtainedon a fully resolved grid,i.e. wherethesize ofeachcell isequalto thefracture aperture.This allowsthefollowingmodelerrormeasurement,

1 Ncoarse = N 1

|

p f ine− pcoarse

|

Ncoarse (11) whereNcoarseisthenumberofgridcellsusedbypEDFMandp _{f ine}

is the corresponding fully-resolved pressure, interpolated to the coarsescale,ifnecessary.Some oftheexperimentswererepeated fortheclassicEDFM,aswellasunstructuredDFM,forcomparison purposes.

Forsimplicity,butwithoutlossofgenerality,theﬂowinthese experimentsisdrivenbyDirichletboundaryconditions,insteadof injectionandproductionwells,whilecapillaryandgravity effects areneglected.

(5)

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

Finally,thesimulationswereperformedusingtheDARSim1 in-housesimulator,usingasequentiallyimplicitstrategyforthe mul-tiphaseﬂowcases.

3.1.pEDFMvalidation

In order to validate pEDFM as a ﬁne-scale model suitable to accommodatefractureswithawiderangeofpermeabilities, a2D homogeneousdomain(km₌₁₎ _is_considered,_having_a₊_−shaped

fracturenetwork, locatedin themiddle.In orderto drivethe in-compressiblesingle-phaseﬂow,Dirichletboundaryconditionswith non-dimensionalpressure valuesof p=1and p=0are imposed ontheleftandrightboundariesofthedomain,respectively,while thetopandbottomsidesaresubjecttono-ﬂowconditions.

As shown in Fig.3, the studyis ﬁrst conductedin a scenario where the fractures are 8 orders of magnitude more conductive thanthematrix.Thereferencesolution,inthiscase, iscomputed ona 1001 × 1001 structured cartesian grid. From the bottom of

Fig.3,itisclearthatbothEDFMandpEDFM,onacoarser11× 11 domain,canreproducethebehavioroftheﬂowasdictatedbythe highly-conductiveembeddedfracturenetwork.

AsshowninFig.4,thesameexperimentwasrerunforthecase wherethefracture permeabilitylies8ordersofmagnitudebelow thatofthehostmatrix.TheresultsexposethelimitationsofEDFM, wheretheimpermeablefracturesaresimplyby-passedbytheflow throughthe(unaltered)matrix,resultinginapressurefield corre-spondingtoareservoirwithhomogeneous(non-fractured) perme-ability.Ontheotherhand,throughitsnewformulation,pEDFMis abletoreproducetheeffectoftheinhibitingflowbarrier(see bot-tomofFig.4),confirmingitsapplicabilitytothiscase.

Theseexperimentsconﬁrm thatpEDFMisasuitable extension of EDFMto a wider range of geological scenarios, being able to reproducethecorrectﬂowbehaviourinthepresenceofbothhigh andlowpermeablefractures,embeddedintheporousmatrix.

3.2. Sensitivitytothefracturepositionwithinthegridcell

Given that pEDFM typically operates on much coarser grids than thefully resolved case, it isof interest toelicit its sensitiv-ityto thefracture positionwithin the hostgrid cell.Tothis end, the ₊_−shaped fracture conﬁgurationis considered; thereference solution is computed on a 37 _{× 3}7 _(i.e.,₂₁₈₇ _{× 2187)} _cell _grid, whilepEDFMgridoperatesat10× 10resolution.

FromFig.3,itappears thatinthecasewhenthefracture net-work ishighlyconductive, thehorizontalfracture isthe one that dictates the ﬂow. Consequently, successive simulations are con-ducted for both EDFM and pEDFM,while moving the horizontal fracturefromtop tobottom,asshowninFig.5.Theiraccuracyis measuredusingEq.(11).

The results show that EDFMis more accurate when fractures areplacedatthecellcenter,ratherthanwhentheyareclosetothe interface.However, oncethefracturecoincides withtheinterface, EDFMconnectsitto both matrixcells(each,withaCI calculated usingS_{i f}₌A_{i f} inEq.(6),insteadof2A_ifaswasthecaseinEq.(8)), thusexplainingtheabruptdipinerror.Incontrast,thepEDFM er-ror attains its peak when fracturesare placed atthe cell centers and doesnot exhibit anyjumps over the interface. The error of bothmethodslieswithin similarbounds(stillpEDFMismore ac-curate) showing that they applicable to the case when fractures arehighlyconductive. Theconsistent aspectofpEDFMisthat, its resultsforthecasewhenfracturescoincidewiththematrix inter-faces, its results are identical tothe DFM method,while –as ex-plainedbefore– thisisnotthecaseforEDFM.

When the networkis nearly impermeable, thelocation of the vertical fracture is critical to the ﬂow (Fig. 4). As such, for the purposesofthecurrentexperiment, itwillbeshiftedfromleftto right,asshowninFig.6.Theresultingerrorplotshowsadramatic increaseforEDFM,whencomparedtoFig.6,duetoitsinabilityto handlefractureswithconductivitiesthatliebelowthatofthe ma-trix.pEDFM,ontheotherhandshowsasimilarbehavioranderror

(6)

Fig. 4. Fully-resolved (c) EDFM (d) and pEDFM (e) pressure solutions in a homogeneous reservoir containing a + −shaped nearly impermeable ﬂow barrier (top).

Fig. 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 pressure mismatch towards the corresponding fully resolved simulation (bottom).

rangeaswasobservedinthecasewithhighlyconductivefractures, i.e.,itretainsitshighaccuracy.

TheseresultsshowapromisingtrendforpEDFM,whichisable to maintain reasonable representation accuracy of the effect of the embeddedfractures. The slightincrease inerrorfor fractures placed nearthematrixcellcenters maybemitigated by employ-ingmoderatelocalgridreﬁnements.

3.3.Sensitivitytothegridresolution

Anotherimportantfactorinassessingthequalityofan embed-dedfracturemodelisitsorderofaccuracywithrespecttothegrid resolution.Aseriesofnestedmatrixgridsforthe₊_−shaped

frac-ture test caseof Figs. 3 and 4 was constructed. The number of

(7)

for-Fig. 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 mismatch towards the corresponding fully resolved simulation (bottom).

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.

(8)

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,wherei=2,3,...,up toareferencegrid resolution,where i=7.Atthesametime,thefracturegridisalsoreﬁnedaccordingly such thatits stepsizeapproximatelymatchestheone inthe ma-trix,h=

x=

y.Themeasureofaccuracyforthiscaseissimilar to Eq. (11), where,this time, no interpolation is necessary, since the cell centroids are inheritedfrom one level to anotherin the nestedgridhierarchy.

Forabettercomparison,alongsidepEDFMandEDFM,thesame sequenceofgeologicalscenarioswassimulatedusingDFMona2D unstructuredgrid(Karimi-Fardetal.,2004),wherethenumberof triangleswastweakedtomatchN=Nx× Nyascloselyaspossible

andwithoutimposinganypreferentialgridreﬁnementaroundthe fractures.

The results of this study, in the case when the fractures are highlyconductive, are depictedin Fig.7. It followsthat all three methods experience a linear decay in error withincreasing grid resolution. The three error snapshots, which were taken when Nx=Ny=36 (or h=0.0015),show thatthepressure mismatchis

mainly concentrated around the tips of the horizontal fracture, whichrepresentthe network’sinﬂowandoutﬂowpoints, respec-tively.ForEDFM, theerrordecaysradially forpointsfurtheraway fromthesefracturetips.ForpEDFM,thecontourcurvesareslightly skewed, depending on thechoice betweenupperand lower ma-trix interfacesfor the fracture projection (both are equally prob-able since thehorizontal fracturecrosses thegrid cell centroids). Finally,forDFM theerrordistribution showssome heterogeneity, which isaconsequence ofusingunstructuredgrids ina medium which, except for the neighborhood of the fractures, is homoge-neous.

The scenario when the fracture network is considered almost impermeable cannotbe properlyhandledbyEDFM,regardlessof which grid resolution is used (Fig. 8). This serious limitation is,

onceagain,successfullyovercomebyusingpEDFM,which,similar toDFM,maintainsitslinearscalabilitywithgridrefinementonthis challengingtest case. The error snapshots depict that, thistime, the pressure is inaccurate around the tips, as well as the body, ofthe vertical barrier. This can be explained by the fact that an embeddedmodel on a coarse grid can have difficulty in placing the sharp discontinuity inthe pressure field at exactly the right location. Still,the pressure mismatchdecayswithincreasing grid resolution, suggesting that local grid refinements around highly contrasting fractures can benefit pEDFM, in a similar manner toDFM.

Toconclude, pEDFMshows asimilar convergence behavior, in termsof grid resolution, to thewidely usedDFM approach. This conﬁrmsthat,inordertodiminishthemodelrepresentationerror, moderatelocalgridreﬁnementscanbeappliednearfractures.

3.4.Sensitivitytothefracture-matrixconductivitycontrast

This last sensitivity study is aimed at determining the re-sponse of pEDFM while changing the conductivity contrast be-tween the ₊ _−shaped fracture network (kf₌₁₀−8_,_{. . .}_,₁₀8₎ _and thematrix(km₌₁_)._To _this_end, _a_coarse _grid_resolution _of_N

x=

Ny=35wasusedandtheresultingpressurewascomparedtothat

fromthereferencecase,whereNx=Ny=37,usingEq.(11).

The results are depicted in Fig. 9 and are in line with the conclusions from previous sections. Namely, for fracture log-permeabilitiesonthepositive sideofthespectrum,theresultsof EDFMandpEDFMare inagreement.As thepermeabilitycontrast passes5ordersofmagnitude,thepressureerrorplateaus,since,at beyondthisstage, thefractures arethe main driversof theﬂow. However,forfracture permeabilitiesclosetoorbelowthat ofthe matrix,the errorof EDFMincreases dramatically. pEDFM,on the

(9)

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.

Fig. 10. Fracture permeabilities (top-left), pressure ﬁeld (top-right) and time-lapse saturation results (bottom) on a 2D test case with homogeneous matrix conductivity, under incompressible 2-phase ﬂow conditions.

(10)

Fig. 11. Heterogeneous matrix and fracture permeability maps (top), pressure and time-lapse saturation results for EDFM (middle) and pEDFM (bottom) on a 2D densely fractured test case, under incompressible 2-phase ﬂow conditions.

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

other hand is able to cope with these scenarios, due to its for-mulation,itsbehaviorshowinganapproximatelysymmetrictrend, whencomparedtothatofthepositivesideoftheaxis.

ThesnapshotsonthetopofFig.9,takenforlower,similarand higherfracture permeabilities w.r.t.thematrix,show theerrorin thepressure producedby pEDFM.It isclearthatthe model inac-curacyis concentratedaroundthetipsoffractureswhichactively inﬂuencetheﬂow.Alsonotethatthereisasmallerroreveninthe casewhen kf ₌_km_,_since _the _pEDFM_{discretization}₍_Section ₂₎_is

slightlydifferentthanthatofahomogeneousreservoir.Ofcourse, thisresultisonlypresentedhereforacademicpurposes– in realis-ticscenarios,whenthecontrastisnothighenough,suchfractures canbehomogenizedintothematrixﬁeld.

Thisconcludes thesensitivitystudies conductedin thispaper. Thetestcasespresentedinthefollowingsubsectionsaremeantto testtheapplicability ofpEDFMtomorecomplex2D and3D frac-turedmedia.

(11)

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

3.5.Time-lapse2Dmultiphaseresults

This set of experiments is aimed at determining the perfor-manceofpEDFMinmultiphaseﬂowscenarioson2Dporousmedia withincreasinglycomplexfracturegeometriesandheterogeneities.

3.5.1. Homogeneousmatrix

pEDFM isfirstapplied inan incompressible2-phase flow sce-nariothrough a 2D homogeneous domain which is crossed by a setoffractureswithheterogeneousproperties,asshowninFig.10. Theboundaryconditionsaresimilartothoseusedforprevious ex-periments,namelyDirichletwithnon-dimensionalvaluesofp=1 andp=0onthe left andrightedges,respectively, whilethetop andbottomsidesaresubjecttono-flowconditions.

The lowpermeablefracturesinhibittheﬂow, leavingonlytwo availablepaths:through themiddleofthe domainandalongthe bottom boundary. As can be seen in the time-lapse saturation mapspresented in Fig. 10, the front, indeed, respects these em-beddedobstacles.Theinjectedﬂuidismostlydirectedthroughthe permeableX-shapednetworkandsurpassestheverticalbarrier,in thelower rightpartofthe domain,onits waytothe production boundary.

Thisresultreinforcestheconclusionthattheconservative pres-sureﬁeldobtainedusingpEDFMleadstotransportsolutionswhich honourawiderangeofmatrix-fractureconductivitycontrasts.

3.5.2. Heterogeneousmatrix

ThefollowingexperimentcomparesthebehaviourofEDFMand pEDFMforsimulating 2-phaseincompressibleﬂow through a 2D porous medium with heterogeneous (i.e. patchy) matrix perme-ability(Fig. 11). The interplaybetweenthelarge- (matrix-matrix) andsmall-scale (fracture-matrix)conductivitycontrastsraises ad-ditionalnumericalchallenges (Hamzehpour etal., 2016) and isa stepping stone in assessing the model’s applicability to realistic cases.

The embedded fracture map used for this test case (top of

Fig.11)isbasedontheBrazilIoutcropfromBertottiandBisdom;

Bisdom et al.(2016). The conductivities of the fracturesforming

the North-Westto South-East diagonal streak, were set to 10−8_, thuscreatinganimpermeableﬂowbarrieracrossthedomain (no-ticeablein dark blue onthe top-right of Fig.11). Forthe rest of the fractures, permeabilities were randomly drawn from a log-uniformdistributionsupported ontheinterval[10−8,108_]_._Finally, similar to previous experiments, ﬁxed pressure boundary

condi-tions p=1 and p=0 are set on the left and right edges, re-spectively,while thetop andbottom sidesaresubjectto no-ﬂow conditions.

The middle rowofplots from Fig.11 show the pressure field andtime-lapsesaturationresults obtainedusingEDFM.Note that the injected fluid is allowed to bypass the diagonal flow barrier, towardstheproductionboundary.This,onceagainshowsthe limi-tationofEDFM,whichisonlyabletocapturetheeffectoffractures withpermeabilitieshigherthanthematrix.However,by disregard-ingflowbarriers,EDFMdeliversan overlyoptimistic and nonreal-isticproductionforecast.

Incontrast toEDFM, thepressure field obtainedusingpEDFM showssharpdiscontinuities(bottom-leftofFig.11).The accompa-nyingsaturationplotsconfirm thattheinjected phaseisconfined by thediagonalbarrierandforcedtoflowthrough thebottomof the domain, thus significantly delaying its breakthrough towards theproductionboundary.

TheseresultsconﬁrmthatpEDFMoutperformstoEDFM,dueto itsapplicability incaseswithcomplexanddensefracture geome-triesandinthepresenceofmatrixheterogeneities.

3.6. Comparisonbetween3DpEDFMandunstructuredDFM

Finally,atestcaseona3D domaincontaining3layersof frac-tures,stackedalongtheZaxis(Fig.12)isconducted.Thetoplayer isaverticallyextrudedversionofthe2DfracturemapfromFig.10. The second layer consistsof a singlefracture network whose in-tersectingplateshavehighlyheterogeneousproperties.Finally,the thirdlayerisrepresentedby3largeindividualplates,witha clus-terofsmallparallelfracturespackedinbetween.

Inthisscenario,theincompressiblesingle-phaseﬂowisdriven from the left boundary, where the pressure is set to the non-dimensionalvalueofp=1,towardstheright, wherep=0,while alltheotherboundariesofthedomainaresubjecttono-ﬂow con-ditions.No other source termsarepresentandgravity and capil-laryeffects areneglected.The resultsofpEDFM,on amatrixgrid with Nx=Ny=Nz=100 and a total of 23381 fracture cells, are

compared to those obtainedusing DFM on an unstructured grid (Fig. 13), where the numberof tetrahedra (matrix)and triangles (fractures) were chosen to approximately match the degrees of freedomonthestructuredgrid.

ThetwopressureﬁeldsareplottedinFig.14usingisosurfaces for equidistant values, and are in good agreement, for decision-makingpurposes.

ThislastnumericalexperimentshowsthatpEDFMhasgood po-tentialforﬁeld-scaleapplication.

4. Conclusionsandfuturework

Inthispaper,anovelProjection-basedEmbeddedDiscrete Frac-tureModel(pEDFM)wasdevised,forflowsimulationthrough frac-tured 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 fromhighlypermeablenetworkstoinhibitingflowbarriers.

The newmodelwasvalidatedon2D and3D test cases,while studying its sensitivity towards fracture position within a matrix cell,gridresolutionandthecross-mediaconductivitycontrast.The resultsshowthatpEDFMisscalableandabletohandledenseand complexfracturemapswithheterogeneouspropertiesinsingle-,as wellasmultiphaseﬂowscenarios.Finally,itsresultsonstructured gridswerefoundcomparabletothoseobtainedusingtheDFM ap-proachonunstructured,fracture-conformingmeshes.

Inconclusion,pEDFMisaflexiblemodel,itssimpleformulation recommendingitforimplementationinnext-generationsimulators forfluidflowthroughfracturedporousmedia.

(12)

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

Fig. 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 pEDFMtounstructuredgrids.Furthermore,systematic benchmark-ing studies (including CPU time comparisons) between pEDFM and alternative mass-conservative DFM approaches on unstruc-tured grids can provide valuable insights for its application and performanceinreal-ﬁeldcases.

Acknowledgements

TheﬁnancialsupportofPI/ADNOCisacknowledged.Also,many thanks are due towards the members of the DARSim (Delft Ad-vanced Reservoir Simulation) research group of TU Delft, forthe usefuldiscussions duringthe developmentofthe pEDFM model-ingapproach.

AppendixA. Averagedistancecalculation

Thecomputationoftheaveragedistancebetweenamatrix con-trolvolumeanda fracturesurface, whichappearsin Eqs.(6)and

(10),mayinvolvenumericalintegrationforarbitrarilyshapedcells.

For2Dstructuredgrids,however,analyticalformulasweregivenin

Hajibeygietal.(2011)forafewspeciﬁcfractureorientations.This

subsection,an adaptation of Chapter 2.3.2 from Pluimers (2015), providesageneralproceduretohandlefracturelineswitharbitrary orientation.

The interfaces of each cell intersected by a fracture are ex-tendeduntiltheyintersectthefractureline,resultinginfourright triangleswithsurfacesA1 toA4,asshowninFig.A.15.Then,given the average distance between each triangle and its hypotenuse,

d

1to

d

4,as(seeHajibeygietal.(2011)),

d

i=

Lxi· Lyi

3

Lxi2+Lyi2

, (A.1)

whereLxiandLyiarethelengthsoftheaxis-alignedsidesof

trian-glei,theaveragedistancebetweengridcelliandfracturelinefis obtained,

d

i f=

A1

d

1+A2

d

2− A3

d

3− A4

d

4 A1+A2− A3− A4 .

(A.2) Notethatnomodiﬁcationisrequiredtotheformulainthecase when fractureslie outsidethe cell, i.e. for the non-neighbouring

(13)

connections fromEq.(10). In addition,this procedure can be di-rectlyappliedto3Dstructured gridswherefracturesareextruded alongtheZaxis,whileageneralizationforfractureplateswithany orientationisthesubjectoffutureresearch.

References

Ahmed, R. , Edwards, M.G. , Lamine, S. , Huisman, B.A. , Pal, M. , 2015. Three-dimensional control-volume distributed multi-point ﬂux approximation coupled with a lower-dimensional surface fracture model. J. Comput. Phys. 303, 470–497 . Berkowitz, B. , 2002. Characterizing ﬂow and transport in fractured geological media:

a review. Adv. Water Res. 25 (8–12), 861–884 .

Bertotti, G., Bisdom, K., Fracture patterns in the Jandeira Fm. (NE Brazil). http://data. 4tu.nl/repository/uuid:be07fe95- 417c- 44e9- 8c6a- d13f186abfbb .

Bisdom, K. , Bertotti, G. , Nick, H.M. , 2016. The impact of different aperture distribu- tion models and critical stress criteria on equivalent permeability in fractured rocks. J. Geophys. Res. 121 (5), 4045–4063 .

Chen, Z. , 2007. Reservoir simulation: mathematical techniques in oil recovery. SIAM . Deb, R. , Jenny, P. , 2016. Numerical modeling of ﬂow-mechanics coupling in a fractured reservoir with porous matrix. Proc. 41st Workshop Geotherm. Reservoir Eng., Stanford, California, February 22–24, SGP-TR-209. 1–9 .

Flemisch, B., Berre, I., Boon, W., Fumagalli, A., Schwenck, N., Scotti, A., Stefansson, I., Tatomir, A., 2017. Benchmarks for single-phase ﬂow in fractured porous media. ArXiv:1701.01496.

Flemisch, B. , Fumagalli, A . , Scotti, A . , 2016. A review of the XFEM-based approximation of ﬂow in fractured porous media. In: Advances in Discretization Methods. Springer, pp. 47–76 .

Fumagalli, A. , Pasquale, L. , Zonca, S. , Micheletti, S. , 2016. An upscaling procedure for fractured reservoirs with embedded grids. Water Resour. Res. 52 (8), 6506–6525 .

Fumagalli, A. , Zonca, S. , Formaggia, L. , 2017. Advances in computation of local prob- lems for a ﬂow-based upscaling in fractured reservoirs. Math. Comput. Simul. 137, 299–324 .

Garipov, T.T. , Karimi-Fard, M. , Tchelepi, H.A. ,2016. Discrete fracture model for coupled ﬂow and geomechanics. Comput. Geosci. 20 (1), 149–160 .

Hajibeygi, H. , Karvounis, D. , Jenny, P. , 2011. A hierarchical fracture model for the iterative multiscale ﬁnite volume method. J. Comput. Phys. 230, 8729–8743 . Hamzehpour, H. , Asgari, M. , Sahimi, M. , 2016. Acoustic wave propagation in hetero-

geneous two-dimensional fractured porous media. Phys. Rev. E 93 (6), 063305 . Karimi-Fard, M. , Durlofsky, L.J. , 2016. A general gridding, discretization, and coars-

ening methodology for modeling ﬂow in porous formations with discrete geological features. Adv Water Resour. 96 (6), 354–372 .

Karimi-Fard, M. , Durlofsky, L.J. , Aziz, K. , 2004. An eﬃcient discrete fracture model applicable for general purpose reservoir simulators. SPE J. 9 (2), 227–236 .

Karvounis, D. , Jenny, P. , 2016. Adaptive hierarchical fracture model for enhanced geothermal systems. Multiscale Model. Simul. 14 (1), 207–231 .

Lee, S.H. , Jensen, C.L. , Lough, M.F. , 20 0 0. Efficient finite-difference model for flow in a reservoir with multiple length-scale fractures. SPE J. 3, 268–275 .

Li, L. , Lee, S.H. , 2008. Eﬃcient ﬁeld-scale simulation of black oil in naturally fractured reservoir through discrete fracture networks and homogenized media. SPE Reservoir Eval. Eng. 11, 750–758 .

Matthäi, S.K. , Mezentsev, A .A . , Belayneh, M. , 2007. Finite element node-centered finite-volume two-phase-flow experiments with fractured rock represented by unstructured hybrid-element meshes. SPE Reservoir Eval. Eng. 10, 740–756 . Moinfar, A . , Varavei, A . , Sepehrnoori, K. , Johns, R.T. , 2014. Development of an effi-

cient embedded discrete fracture model for 3d compositional reservoir simulation in fractured reservoirs. SPE J. 19, 289–303 .

Moortgat, J. , Amooie, M. , Soltanian, M. , 2016. Implicit ﬁnite volume and discontin- uous galerkin methods for multicomponent ﬂow in unstructured 3d fractured porous media. Adv. Water Resour. 96, 389–404 .

Moortgat, J. , Firoozabadi, A. , 2013. Three-phase compositional modeling with capil- larity in heterogeneous and fractured media. SPE J. 18, 1150–1168 .

Norbeck, J.H. , McClure, M.W. , Lo, J.W. , Horne, R.N. , 2016. An embedded fracture modeling framework for simulation of hydraulic fracturing and shear stimula- tion. Comput. Geosci. 20 (1), 1–18 .

Pluimers, S. , 2015. Hierarchical Fracture Modeling. Delft University of Technology, The Netherlands Msc thesis .

Reichenberger, V. , Jakobs, H. , Bastian, P. , Helmig, R. ,2006. A mixed-dimensional ﬁ- nite volume method for two-phase ﬂow in fractured porous media. Adv. Water Resour. 29 (7), 1020–1036 .

Sahimi, M. , Darvishi, R. , Haghighi, M. , Rasaei, M.R. , 2010. Upscaled unstructured computational grids for eﬃcient simulation of ﬂow in fractured porous media. Transp. Porous Media 83 (1), 195–218 .

Schwenck, N. , Flemisch, B. , Helmig, R. , Wohlmuth, B. , 2015. Dimensionally reduced ﬂow models in fractured porous media: crossings and boundaries. Comput. Geosci. 19 (1), 1219–1230 .

Shah, S. , Moyner, O. , ¸T ene, M. , Lie, K.-A. , Hajibeygi, H. , 2016. The multiscale restric- tion smoothed basis method for fractured porous media (F-MsRSB). J. Comput. Phys. 318, 36–57 .

Tatomir, A . , Szymkiewicz, A . , Class, H. , Helmig, R. , 2011. Modeling two phase ﬂow in large scale fractured porous media with an extended multiple interacting con- tinua method. Comput. Model. Eng. Sci. 77 (2), 81 .

¸T ene, M. , Al Kobaisi, M. , Hajibeygi, H. , 2016. Multiscale projection-based Embedded Discrete Fracture Modeling approach (F-AMS-pEDFM). In: ECMOR XV-15th Eu- ropean Conference on the Mathematics of Oil Recovery .

¸T ene, M. , Al Kobaisi, M.S. , Hajibeygi, H. , 2016. Algebraic multiscale method for ﬂow in heterogeneous porous media with embedded discrete fractures (F-AMS). J. Comput. Phys. 321, 819–845 .

Projection-based Embedded Discrete Fracture Model (pEDFM)