Algebraic dynamic multilevel method for embedded discrete fracture model (F-ADM)

Delft University of Technology

Algebraic dynamic multilevel method for embedded discrete fracture model (F-ADM)

HosseiniMehr, Mousa; Cusini, Matteo; Vuik, Cornelis; Hajibeygi, Hadi

DOI

10.1016/j.jcp.2018.06.075

Publication date

2018

Document Version

Final published version

Published in

Journal of Computational Physics

Citation (APA)

HosseiniMehr, M., Cusini, M., Vuik, C., & Hajibeygi, H. (2018). Algebraic dynamic multilevel method for

embedded discrete fracture model (F-ADM). Journal of Computational Physics, 373, 324-345.

https://doi.org/10.1016/j.jcp.2018.06.075

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.

Otherwise as indicated in the copyright section: the publisher

is the copyright holder of this work and the author uses the

Dutch legislation to make this work public.

Contents lists available atScienceDirect

Journal

of

Computational

Physics

www.elsevier.com/locate/jcp

Algebraic

dynamic

multilevel

method

for

embedded

discrete

fracture

model

(F-ADM)

Mousa HosseiniMehr

a

,

b

,

∗

,

Matteo Cusini

b

,

Cornelis Vuik

a

,

Hadi Hajibeygi

b

a_{FacultyofElectricalEngineering,MathematicsandComputerScience,DepartmentofAppliedMathematics,DelftUniversityofTechnology,}

VanMourikBroekmanweg 6,2628XEDelft,theNetherlands

b_{FacultyofCivilEngineeringandGeosciences,DepartmentofGeoscienceandEngineering,DelftUniversityofTechnology,Stevinweg 1,}

2628CV,Delft,theNetherlands

a

r

t

i

c

l

e

i

n

f

o

a

b

s

t

r

a

c

t

Articlehistory:

Received17April2018

Receivedinrevisedform28June2018 Accepted30June2018

Availableonline9July2018

Keywords:

Flowinporousmedia Adaptivemeshrefinement Multiscalebasisfunctions Algebraicmultiscalemethod Multilevelmultiscalemethod Fracturedporousmedia

Scalablephysics-basednonlinearsimulation

Wepresentanalgebraicdynamicmultilevelmethodformultiphaseflowinheterogeneous fractured porous media (F-ADM), where fractures are resolved at fine scale with an embeddeddiscretemodellingapproach.Thisfine-scalediscretesystememploys indepen-dentfine-scalecomputationalgridsforheterogeneousmatrixanddiscretefractures,which resultsinlinearsystemsizesoutofthe scopeoftheclassicalsimulationapproaches.To reducethe computational costs, yet provide accuratesolutions, on thishighly resolved fine-scalemesh, F-ADMimposes independent dynamic multilevel coarse gridsfor both matrixandlower-dimensionaldiscretefractures.Thefully-implicitdiscretesystemisthen mappedintothisadaptivedynamicmultilevelresolution forallunknowns(i.e.,pressure and phase saturation). The dynamic resolution aims for resolving sharp fronts for the transport unknowns, thus constant interpolators are used to map the saturation from coarsetofine gridsbothinmatrix and fractures. However,dueto the globalnature of thepressureunknowns,localmultilevelbasisfunctionsforbothmatrixandfractureswith flexiblematrix-fracturecouplingtreatmentareintroducedforthepressure.Theassembly ofthefullsetsofbasisfunctionsallowsformappingthesolutionsupanddownbetween anyresolutions.Duetoitsadaptive multilevelresolution, F-ADMdevelopsanautomatic integratedframeworktohomogeniseorexplicitlyrepresentafracturenetworkatacoarser levelbyselectionofthemultilevelcoarsenodes ineachsub-domain.Varioustestcases, includingmultiphase flow in 2D and 3D media, are studied, where only a fraction of the fine-scale grids is employed to obtain accurate nonlinear multiphase solutions. F-ADMcastsapromisingapproachforlarge-scalesimulationofmultiphaseflowinfractured media.

©2018ElsevierInc.Allrightsreserved.

1. Introduction

Many geo-engineering applications (e.g.,hydrocarbon and geothermal energy production) are driven by underground fluid flow in large-scale (O(103) [m]) heterogeneous fractured porous media [1]. Even at what is known as Darcy’s

*

Correspondingauthorat:FacultyofElectricalEngineering,MathematicsandComputerScience,DepartmentofAppliedMathematics,DelftUniversityof Technology,VanMourikBroekmanweg 6,2628XEDelft,theNetherlands.

E-mailaddresses:M.Hosseinimehr@tudelft.nl(M. HosseiniMehr),M.Cusini@tudelft.nl(M. Cusini),C.Vuik@tudelft.nl(C. Vuik),H.Hajibeygi@tudelft.nl

(H. Hajibeygi).

https://doi.org/10.1016/j.jcp.2018.06.075

simulations by tacklingdifferent aspects ofthe entirecomplexity map. Multiscalemethods have beendeveloped to effi-cientlysolvetheelliptic(orparabolic)pressureequation,whichhighlyheterogeneouscoefficients,bysolvingthesystemon acoarsegridwhilepreserving thefine-scaleheterogeneities.Ontheotherhand,DLGRtechniquesadaptthegridresolution throughoutthetime-dependentsimulationtoemployahigh-resolutiongridwherenecessary(i.e.,theadvancingsaturation front), andare, therefore,transport-orientedmethods.The AlgebraicDynamicMultilevel(ADM) method[20] has been in-troducedtoaddressthemulti-scalemultilevelcoexistenceofthepressure(ellipticorparabolic) andtransport(hyperbolic) unknowns,atthesametime,withinbothfullyimplicit(FIM)andsequential(implicitandexplicit)simulationframeworks. ADMdevelopsadynamic multilevelsystemforall unknowns,throughanalgebraic formulation,wheretheresolutionsare connectedthroughsetsofbasisfunctions.ADMhasbeenrecentlydevelopedformultiphaseflowsinheterogeneous forma-tions,withcompositionalandcapillaryeffects[21].

ADMextends theapplicability of themultiscale methods tofully-implicit (stable) simulations, allows forcrossing the scalesforallunknownswithamultileveldynamicmesh,doesnotrequirereconstructionofconservativefluxfield,employs thebasis functionswhichare computedonlyatthebeginning ofthesimulation,anddoesnotrely onanysmoothing it-erativeprocedure.Thedevelopment ofsuch adynamicmultilevelschemeforfracturedmedia hasnotyetbeenaddressed, despitetheirhighimportanceinthegeo-scientificcommunityandextensiveliteratureforfine-scaleconsistentdiscrete rep-resentations[22–32].Suchadynamicmultilevelapproachallowsforcapturingexplicitfracturesattheirrelevantresolution whilemaintainingthescalabilityofthesimulationforreal-fieldapplications.

Inthiswork,an ADMmethodforsimulationofmultiphaseflow inheterogeneousfracturedporousmediaisdeveloped (F-ADM).Thisisachievedbydevisingmultilevelmultiscalebasisfunctionsforembeddeddiscretefracturemodelling(EDFM) approach[33–38].Moreprecisely,first,anEDFMfine-scalediscretesystemisobtained[39–41,34,42–45],whichallowsfor independent grids forfracturesandmatrix.This makes itconvenient to treatcomplex fracture geometrieswith multiple intersections[46,47].Giventhisfine-scalediscretesystem,withfullyimplicitflow-transportcouplingtreatment,theF-ADM mapsittoadynamicmultilevelnestedresolution,bothforthematrixandthefractures.Suchadynamicmultilevelgridis obtainedtoemployfine-scaleresolutiononlywhereneeded(e.g.,attheadvancingsaturationfronts).

Similar asin the fine-scale system, the ADM dynamic grid is chosen independently for matrix and fractures, based onafront-trackingcriterionthataimstominimisethecost-accuracytrade-off.Mappingthesolutionsacrossdifferentgrid resolutionsisperformedthroughsequencesofrestrictionandprolongationoperators.Finite-volumerestrictionoperatorsare employedforall unknownsto ensuremassconservationatall levels.Ontheotherhand,differentinterpolationstrategies are consideredforthe twomain unknowns(pressure andsaturation). Specifically,thedevised multilevelmultiscale basis functionsareemployedaspressureinterpolatorswhereas piece-wiseconstantfunctionsareusedtointerpolatesaturation values,asthegridrefinementstrategyavoidscrossingthescalesatthesaturationfrontlocations.

Such a development allows for an automatic framework to explicitly oreffectively representa fracture at anycoarse level,i.e.,throughtheselectionofthelocationofthecoarsenodes.Moreprecisely,ifnocoarsenodeforafracturenetwork inside a dual coarse grid is chosen, the matrix basis functionswill automatically homogenise the effect of the fracture atthe corresponding coarse scale. However, if coarse nodesare selected inside the fracture networkof that dual coarse cell, thefractureswill be explicitlyrepresented atthecorresponding coarser resolution.The accuracy andthe sensitivity oftheF-ADMtotheerrorcriterionandtothepressurebasis functionsarestudiedthroughasetof2Dand3Dtest cases. Thedevised F-ADMstrategyprovides adynamic treatmentofhighly fracturedmedia, andcastsa promisingapproachfor real-fieldapplications.

The rest ofthe paperis organised asfollows. The governingequations, fine-scale discrete systemandthe simulation strategy are briefly presented inSection 2. Then, the F-ADM methodfor flow in porous media withembedded discrete fractures isdescribed inSection 3.Numerical results are presented in Section 4. Finally,the paper is concludedin Sec-tion5.

2. Fine-scale equations and solution strategy 2.1. Governingequations

Massconservationforphase

α

intheabsenceofmass-exchangebetweenphases, capillary,andgravitationaleffects,in porousmediawithnfracdiscreteembeddedfracturesreads

Fig. 1. Theindependentmatrix3Dgridandthefractures2Dgridsareshown.Notethateachdomainhasitsowngridandthatanyfractureorientationcan beconsidered.

∂

∂

t



φ

m

ρ

_αmSm_α



− ∇ ·



ρ

_αm

λ

m

· ∇

pm



=

nfrac



i=1

ρ

α

Q

mfαi

+

ρ

αqmα on



m

⊆ 

n (1)

forthematrixand

∂

∂

t



φ

fi

_ρ

fi αSαfi



− ∇ ·



ρ

fi α

λ

fi

· ∇

pfi



=

ρ

α

Q

αfim

+

nfrac



j=1



ρ

α

Q

fifj α



j=i

+

ρ

αqαfi on



fi

⊆ 

n−1

_∀

_i

_{∈ {}

₁

_{, ...,}

_n frac

}

(2)

for the lower-dimensionalfractures. There existnα phases. Moreover, thesuperscripts m and f in Eqs. (1)–(2) indicate, respectively,therockmatrixandthefractures.Here,

φ

istheporosityofthemedium,

ρ

α , Sα,

λ

are,respectively,thedensity,

saturation,andmobilityofphase

α

.Moreover,

λ

=

krα

μαK holds, wherekr,

μ

and

K are

phaserelativepermeability,viscosity

androckabsolutepermeabilitytensor,respectively.Also,qα isthephasesourceterm(i.e.,wells).Finally,

Q

mfi

α and

Q

αfimare

thephaseflux exchangedbetweenmatrixandthe i-thfracture,whereas

Q

fifj

α representstheinfluxofphase

α

from j-th

fracturetothei-thfracture.Notethatmassbalanceenforces



V

Q

mfi α dV

= −



A_fi

Q

fim α d A and



A_fi

Q

fifj α d A

= −



A_{f j}

Q

fjfi α d A.

Equations (1)–(2),subjecttoproperinitial andboundary conditions,formawell-posed systemfornα unknowns,once the

nα



α=1

Sα

=

1 constraintisemployedtoeliminateoneofthephasesaturationunknowns.Here,thissystemofequationsis solvedfortwophaseflowwithprimaryunknownsofp and S1(fromnowonindicatedasS).

2.2. Fine-scalediscretesystem

Thecoupledsystemofnon-linearequations,i.e.,Eqs.(1)–(2),isdiscretizedinspacewithatwo-point-flux-approximation (TPFA) finite-volume schemein spaceanda backward (implicit) Eulerscheme intime. Independent structured grids are generatedforathree-dimensional(3D)porousmediumand2Dfractures.AnillustrationispresentedinFig.1.Theadvective TPFAfluxbetweencontrolvolumesi and j reads

F_α,i j

=

ρ

α∗ k∗_rα

μ

∗_αTi j

(

pi

−

pj

).

(3) Here, Ti j

=

Ai j di jK H

i j is thetransmissibility betweencells i and j. Ai j is theinterface area betweencells i and j,di j is the

distancebetweenthecellscentresandK_{i j}H istheharmonicaverageofthetwopermeabilities.Thetermsindicatedwiththe superscript

∗

areevaluatedusingaphasepotential upwindscheme.FollowingEDFM[34,46],thefluxesbetweenamatrix celliandafracturecelljaremodelledas

F

mf α,i j

= −

F

f m α,i j

= −

ρ

α∗ k∗_rα

μ

∗_αT mf i j

(

pmi

−

p f j

),

(4)

where the transmissibility Tmf_{i j}

=

K_{i j}HC Ii j. K_{i j}H is the harmonicaverage permeabilitybetweenthe overlapping matrixand

fracture elements,andC Ii j istheconnectivityindexbetweenthem.TheEDFMmodelsthematrix-fractureconnectivityas

C Ii j

=

Amf_{i j}

di j where A

mf

i j istheareafractionoffracturecell j overlappingwithmatrixcelli (seeFig.2),and

d

i jistheaverage

Fig. 2. Illustration of overlapping fracture and matrix cells, where the overlap section forms an irregular hexagon.

Fig. 3. Illustrationoftwointersectingfractureplates(left)anddiscreteelements(right),withtheintersectionlinesshowninred.(Forinterpretationofthe coloursinthefigure(s),thereaderisreferredtothewebversionofthisarticle.)

a consistent embedded formulation [46], however, since we consider highly conductive fractures the difference between pEDMFandEDFMisminimal.

Similarly, the flux exchange between intersecting fracture elements i (belonging to fracture f ) and j (belonging to fracture g)ismodelledas

F

f g α,i j

= −

F

g f α,i j

= −

ρ

α∗ k_rα∗

μ

∗_αT f g i j

(

p f i

−

p g j

).

(5)

Thetransmissibility T_{i j}f g betweenthetwocellsiscomputedbasedonalower dimensionalconnectivityindexformulation. The intersectionof2D fractureplates isa line,asshownin Fig.3.Eachtwo intersecting fracture cellsintersect ata line segment Ii j (the red line in Fig. 3 (right)) with the average distances from the intersection segment of

d

i Ifi j

=

d

g j Ii j. Therefore,T_{i j}f g iscalculatedinaharmonic-averageformulationas

T_{i j}f g

=

K_{i j}H C I_{i I}f i j

×

C I g j Ii j C I_{i I}f i j

+

C I g j Ii j

,

(6) where,e.g.,C I_{i I}f

i j istheconnectivityindexbetweenthe2Dfractureelementi andthe1Dintersectionsegment Ii j.Thus,at eachtime-stepthefollowingsystemofequationsissolved

φ

m

ρ

m αSmα

n+1 i

−

φ

m

ρ

m αSmα

n i



t

+

Nn



j=1 F_α,i j

+ +

nfrac



k=1



N_fk



j=1

F

mfk α,i j



=

ρ

αqm_α,i

,

∀

i

∈ {

1

, ...,

Nm

}

(7)

inthematrixand

φ

fh

_ρ

fh α Sαfh

n+1 i

−

φ

fh

_ρ

fh α Sαfh

n i



t

+

Nn



j=1 Fα,i j

+

Nm



j=1

F

fhm α,i j

+

+

nfrac



k=1



N_fk



j=1

F

fhfk α,i j



=

ρ

αqαfh,i

,

∀

i

∈ {

1

, ...,

Nfh

}

(8) ineach fractureh.Here, Nmand Nfk arethenumberofcellsofthematrixandoffracture k. Nn indicatesthenumberof neighbouringcells(2in1D,4in2D,6in3D).Equations(7)–(8) canbewritteninresidualformas

(

rm_α,i

)

n+1

=

ρ

αqmα,i

−

φ

m

ρ

m αSmα

n+1 i

−

φ

m

ρ

m αSmα

n i



t

−

Nn



j=1 Fα,i j

− −

nfrac



k=1



N_fk



j=1

F

mfk α,i j



∀

i

∈ {

1

, ...,

Nm

}

(9) fortherockmatrix,and

(

rfh α,i

)

n+1

₌

_ρ

αq_αfh_,i

−

φ

fh

_ρ

fh α Sαfh

n+1 i

−

φ

fh

_ρ

fh α Sαfh

n i



t

−

Nn



j=1 F_α,i j

−

Nm



j=1

F

fhm α,i j

−

−

nfrac



k=1



N_fk



j=1

F

fhfk α,i j



,

∀

i

∈ {

1

, ...,

Nfh

}

(10)

forfracture fh.Letusdefine rn

= [(

rm

)

n

,

(

rf1

)

n

...(

rfnfrac

)

n

]

T where

(

rk

)

n istheresidualvector ofmediumk at time-stepn.

Similarly, pn _{and S}n _{indicate thevectorsofpressureandsaturationunknowns(of allmedia).The residualr}n+1 _isa non-linearfunctionoftheprimaryunknowns pn+1 and Sn+1.Thus,ateach time-stepaNewton–Raphsonmethodisemployed tosolvethenon-linearsystemiteratively,i.e.

rν+1

=

rν

+

∂

r

∂

p





ν

δ

pν+1

+

∂

r

∂

S





ν

δ

Sν+1

=

0 (11)

where thesuperscript

ν

istheiteration index. Consequently,ateach Newton’siteration a linear-systemJν

δ

xν+1

= −

rν is

solved.Here,

Jν is

theJacobianmatrixwith

δ

xν+1

_{= [δ}

_p

_,

_δ

_S

_]

T_{.Therefore,assumingtwophases(w ando representingwater}

andoilrespectively),thelinearsystemofequationscanbewrittenas

⎛

⎜

⎜

⎜

⎝

Jmm op Jopmf Jmmos Josmf Jfmop Jopff Josfm Jffos Jmm_wp Jmfwp Jmmws J mf ws Jfm_wp Jff_wp Jfm_ws Jff_ws

⎞

⎟

⎟

⎟

⎠

ν







Jν₀

⎛

⎜

⎜

⎝

δ

pm

δ

pf

δ

Sm_w

δ

Sf_w

⎞

⎟

⎟

⎠

ν+1







δxν₀+1

= −

⎛

⎜

⎜

⎝

rm_o rof rm_w rf_w

⎞

⎟

⎟

⎠

ν

  

rν₀ (12)

Inthiswork,non-linearconvergenceisreachedwhenthefollowingconditionsaresatisfied:

||

r

||

_∞

<

1

∧

||δ

p

||

_∞

||

p

||

_∞

<

2

∧ ||δ

S

||

∞

<

2

.

(13)

Here, 1and 2areuser-definedtolerances.

Forreal-fieldapplications,asmentioned before,thesolutionofthisnonlinearsystemofequationsischallenging.Next, theF-ADMispresentedasapromisingframeworkforresolvingthischallenge.

3. F-ADM method 3.1. Solutionstrategy

At each Newton’s iteration,F-ADM constructsan algebraically reduced systembased on Eq.(12) on a multilevel grid whichisdynamicallydefinedatthebeginningofeachtime-step.AschematicofF-ADMsimulationisshowninFig.4.

As a first step,the F-ADMgrid needs to be described. Sets ofnm andnfi nestedCartesian grids are imposed on the matrixandfracturefine-scalegridcells. Thelevelatwhichthefine-scalecomputational domainisobtainedisreferred to asl

=

0.Let Nl

m bethe numberofgrid cellsin theporousmatrixandNlfi the numberofgridcells infracturei, bothat levell.Thecoarseningratio,

γ

l _{canthereforebedefinedas}

γ

l

=



γ

_ml

,

γ

_{f 1}l

,· · · ,

γ

_nl_frac



=

⎛

⎝

Nlm Nlm−1

,

Nl f1 Nl_f−1 1

,

· · · ,

Nl_f nfrac Nl_f−1 nfrac

⎞

⎠ .

(14)

Fig. 4. Schematic description of F-ADM reservoir simulation.

Notethatthecountofcoarseninglevelsandthecoarseningratiosforthematrixandeveryindividualfractureare indepen-dent,whichleadstoaflexibleframework.Therefore,theF-ADMsolutiongridateachtime-stepisconstructedbycombining grid-cellsatdifferentresolutions.

F-ADMmapsthefine-scalesystemtothedynamicmultilevelgridbyapplyingsequences ofrestriction(R)and prolon-gation(P)operators,whichareconstructedbasedonthevaluescomputedatthebeginningofthesimulation(Fig.4).More precisely,ateachiteration,theF-ADMsystemreads

ˆ

Rl_l−1

. . . ˆ

R0₁J0P

ˆ

10

. . . ˆ

Pll−1







JA D M

δˆ

xl

= − ˆ

Rl_l−1

. . . ˆ

R0₁r0







rl

,

(15) whereR

ˆ

l

l+1 istherestrictionoperatormappingpartofthevectorofsolutionswhichareatresolutionl (

δ

ˆ

xl)toresolution

l

+

1 (δx

ˆ

l+1).Similarly,P

ˆ

ll+1 istheprolongationoperatormappingpartoftheentiresolutionvectorwhichareatlevell

+

1

to levell. Oncethe F-ADM systemis solved, the approximatedsolution at fine-scaleresolution

δ

x ₀ (referencefine-scale solutionisrepresentedas

δ

x0)isgivenby

δ

x0

≈ δ

x 0

= ˆ

P10

. . . ˆ

Pll−1

δ

xl

.

(16)

The F-ADMprolongationoperator Pi_i−1 betweenevery levels i and i

−

1 isobtainedfortheentiredomain, though at each timestep onlya fractionofthedomainneeds togo throughthismap, whichisillustrated by P

ˆ

i

i−1 inEq.(15).The prolongationoperator

P

i

Pi_i−1

=

⎛

⎜

⎜

⎜

⎝

[(

Pp

)

ii−1

]

mm

[(

Pp

)

i i−1

]

mf 0 0

[(

Pp

)

i_i−1

]

_{f m}

[(

Pp

)

i_i−1

]

_{f f} 0 0 0 0

[(

PS

)

i_i−1

]

_mm 0 0 0 0

[(

PS

)

i_i−1

]

_{f f}

⎞

⎟

⎟

⎟

⎠

Ni−1×Ni

.

(17)

Similarly,therestrictionoperatorreads

Ri_i−1

=

⎛

⎜

⎜

⎜

⎝

[

Ri_i−1

]

_m 0 0 0 0

[

Ri_i−1

]

_f 0 0 0 0

[

Ri_i−1

]

_m 0 0 0 0

[

Ri_i−1

]

_f

⎞

⎟

⎟

⎟

⎠

Ni×Ni−1

.

(18)

Afinite-volumerestrictionoperatorisemployedinordertoguaranteemassconservationasconvergenceisreachedon themultilevelF-ADMgrid,i.e.,

Ri_i−1

(

s

,

t

)

=



1 if cell s is inside coarser cell t,

0 otherwise

.

(19)

On theother hand,thepressure andsaturation blocksofthe prolongationoperator,

(

Pp

)

ii−1 and

(

PS

)

ii−1 aredifferent [20] as different interpolation rules are used for each variable. In this work

(

PS

)

i_i−1

=

Ri_i−1

T

, where the superscript T

indicates the transpose operator.The

(

Pp

)

i_i−1 blocks,instead, areconstructed followinga multilevelmultiscaleprocedure

forfracturedmediadevelopedinthisworkanddescribedinthenextsubsection.

3.2. Fracturemultilevelmultiscalebasisfunctions

The sub-block

(

Pp

)

ii−1 intheprolongationoperatorinterpolatesthepressurefromresolutioni to i

−

1.Here,the two-level multiscalestrategy[36] isextended tomultiplelevels.Eventhoughthefine-scalepressureequation isaTPFA-based system, startingfromcoarselevel1,allthecoarse-scalesystemsdevelopMPFA-basedconnectivities,duetothelocalbasis functions. Here,theMPFA-basedpressurematrix Ai_p−1,whichstands asfine-scalesystemforleveli

−

1,isconstructedby themultiscalefinite-elementprocedureas

˜

Ai_p−1

= (

Pi_i−₋1₂

)

T Ai_p−2Pi_i−₋1₂

.

(20)

Then, A

˜

i_p−1 isreducedtoaTPFA-basedsystem,i.e.,

Api−1

=

tp f a

( ˜

Ai_p−1

),

(21)

wherethetp f a

(

•)

extractstheTPFAcomponentsofthematrix(heptadiagonalmatrixfor3DCartesiangrids),byconsidering the transmissibilitybetweencellsi and j equal toentry Ai j oftheoriginal matrix.Thealgebraic multiscaleprocedurefor

the calculation ofbasis functionsis nowperformed onthisTPFA systematlevel i

−

1 [36], inorder toobtain the basis functions atlevel i,i.e.,

(

Pp

)

i_i−1. Notethat the TPFAreduction of theMPFA multilevel systemisdone independently for

eachsub-domainofmatrix,andindividualfractures.

In addition,similar asfor thetwo-level multiscale strategy[36],F-ADM allowsfor full flexibility inthe choiceof the couplingstrategyatallmultilevelresolutions.Ingeneral,thepressureprolongationoperatorofagivenleveli reads

(

Pp

)

i_i−1

=



(

Pm_p

)

i_i−1

(

Ppf

)

ii−1



=



(

Pmm p

)

ii−1

(

P mf p

)

ii−1

(

Ppf m

)

ii−1

(

P f f p

)

ii−1



.

(22)

Here, the blocks

(

Pm

p

)

ii−1 and

(

P

f

p

)

ii−1 contain thebasis functionscorresponding to the matrix andfracture coarse-scale unknowns.Eachone,e.g.,

(

Pm

p

)

ii−1,canhavenonzerovaluesinsidematrixandfracturesub-domains,dependingonthelocal couplingstrategy.Moreprecisely,whilethefully-coupledapproachleadstononzerosub-blocksof

(

Pmfp

)

ii−1 and

(

P

f m p

)

ii−1, the decoupledapproachwouldleadtozerovaluesintheselocalcouplingsub-blocks.Fig.5showsasurface plotofsome matrixandfracturebasisfunctionsattwodifferentcoarseninglevelsfora2Dhomogeneousdomain.Adetaileddescription ofthetwo-levelbasisfunctionentriescanbefoundin[36].

Fig. 5. Exampleofmultilevelcoupledbasisfunctionsfora2Dhomogeneousdomain.Matrix(leftcolumn)andfracture(rightcolumn)basisfunctionsfor coarseninglevels1(toprow)and2(bottomrow)areshown.

3.3. Selectionofthegridresolution

At each time-step, the F-ADM solution grid is constructed by combining grid cells of the previously defined sets of grids both inthe matrix and inthe fractures. The grid resolution attime-step n

+

1 is chosen based on the saturation map at time-step n, thus employing an explicit procedure. A semi-implicit gridding strategy [48] could also be used to avoidmisplacement ofthe grid.Fig. 6showsan example ofa 3D test casewithone fracture plate, forsimplicityofthe illustration. Twocoarse levelsare employed inthe matrixalong withone coarse levelinthe fracture.The grid selection methodisappliedindependentlyinthematrixandinthefractures,therefore,differentcriteriacanbesetforeachmedium. Sincenowell-functionsareemployedindynamicsimulations,fine-scaleresolutioniskeptaroundallwellsatalltime.

4. Numerical results

Numericalresultsarepresentedinthissection.First,theEDFMmodelfortwo-phaseflowisvalidated,thentheaccuracy ofamultilevelmultiscalemethodforsinglephaseflow onhomogeneousandheterogeneoustestcasesisassessed.Finally, theaccuracyandsensitivityoftheproposedF-ADMmethodareinvestigatedforboth2D and3Ddomains.F-ADMresults arecomparedagainstfine-scalesimulations,usedasareference.

Inalltest cases,fluidsandrocksareconsidered tobeincompressibleandquadraticrelative permeabilityfunctionsare used.

4.1. TestCase1:validationofEDFM

To validate the EDFM implementation, a 2D homogeneous reservoir of 9m

×

9m is considered as shown in Fig. 7. A cross-shaped fracture network is present in the middle of the domain. Each fracture segment is 5 [m] long, with a permeability value of 109 timeshigher than that of the rock matrix(Km

=

2.5

·

10−13m2 and Kf

=

1.3

·

10−4m2). The

fracture aperture is 4

·

10−2 _{m. Two incompressible fluid phases are considered. The reservoir has an initial saturation}

S1,init

=

0.1. No flow boundary conditionsare considered at the top andbottom boundaries whereas the left andright

boundarieshave fixedpressures of2.0

·

107 _{Pa and1.0}

_·

₁₀7 _{Pa,respectively. Phase1 is injected fromtheleft boundary} whereasproductionoccursattheright.

Fig. 6. ExampleofF-ADMgridselectionfora3Dtestcasewithonesingleplanarfracture.Twocoarselevelsareusedinmatrix(l=1,2)andonecoarse levelisusedinthefracture(l=1).Thebottomrightfigureshowsanexampleofadynamicgrid.

Fig. 7. Visualization of the 2D domain with a cross-shaped fracture network at the center.

Thereferencesolution(referredtoasDNS)isobtainedbyimposinga225

×

225 gridthatallowstofullyresolvetheflow inside thefracture.EDFMsimulationsareperformedwiththreedifferentmatrixgridresolutionsof15

×

15,25

×

25,and 45

×

45.Thegrid-cellsinsidethefracturesarechosentohavesimilardimensionsasthematrixcells.Thetime-stepsizeis 10−4 dayforall simulations.The pressureandsaturation mapsafter0.0235 [Days]are showninFig. 8.Notethat errors decreasesuponrefinementoftheEDFMsolutiongrid.

Figs. 9aand9b presentthepressureandsaturationerrorsasfunctionsofthe simulationtime.Givenavariable x (i.e.,

pressureorsaturation),theerror,ex,iscalculatedas

e

x

=

||xDNS_||x−_DNSxEDFM_|| ||2

Fig. 8. ComparisonofthefullyresolvedDNS(leftcolumn)using225×225 cellsandEDFMresultswithdifferentgridresolutions(lefttoright:45×45, 25×25,15×15)after0.0235 [Days].Pressureisshownonthetop,whilethesaturationmapisillustratedinthebottom.

Fig. 9. Pressure and saturation errors of EDFM simulations with respect to the reference solution (DNS).

4.2. TestCase2:F-ADMforsingle-phaseflowina2Dhomogeneousfracturedreservoir

Here,theaccuracyoftheF-ADMmethodforincompressiblesingle-phaseflowinfracturedporousmediaisassessed.In thiscase, no gridadaptivity isemployed andF-ADMis equivalent toa multilevelmultiscale method.Thus, theobjective ofthistest caseis toverifytheproposed multilevelinterpolation strategy.A 2Dfractured75m

×

75m reservoir with35 fractures is considered. A 375

×

375 grid is imposed on the rock matrix. Fractures have different lengths butidentical aperturesof5

·

10−3_{m.Thefracturenetworkconsistsof4420 gridcells(foratotal145045 cells).Thematrixpermeability} is10−14m2 whereasfracturespermeabilityis2.08

·

10−6m2.Twoinjectionwellsinbottomleftandtopleftcorners(pinj

=

2

·

107_{Pa)andtwoproductionwellsinbottomrightandtoprightcorners(p}

prod

=

1

·

107Pa)arepresent.

Themultilevelmultiscalemethodisrunwithdifferentnumberofcoarselevelsandcoarseningratios.Here,fully-coupled multiscalebasisfunctionsareconsidered.The errorofthemultiscalepressuresolution(pms)withrespecttothefine-scale pressuremap(pfs)iscalculatedas

e

_p

=

||pfs−pms||2

||pfs||2

.

(24)

Thesettingsofeachmultiscalerun,thepressureerror(ep)andthenumberofdegreesoffreedomarepresentedinTable1.

Remarkthattheerrorislargerwhenahighernumberoflevelsorlargercoarseningratiosareemployed.However,theorder ofmagnitudeofthe multiscaleerroristhe sameforall multilevelmultiscalesettings. Moreover,the coarseninglevelsin

Table 1

Settingsanderrorsofthemultilevelmultiscaleruns.

Num. coarse levels γm Total DOF ep[%]

Case 1 1 5×5 6509 0.11

Case 2 1 25×25 1109 0.22

Case 3 2 5×5 877 0.27

Case 4 1 125×125 893 0.64

Case 5 3 5×5 661 0.98

Fig. 10. Testcase2–comparisonofthe pressuredistributionobtainedbyfine-scalesimulationsand F-ADMsingle-phasesolverwith enforcedstatic multilevelmultiscale,wherethepressureissolvedeverywhereatthecoarsestlevel.Topplotisthefine-scalesolution.Middlerowshows1 levelmultiscale solutionwithγm=25×25 (left)and2 levelmultiscalesolutionwithγm=5×5 (right).Showninthebottomrowis1 levelmultiscalesolutionwith

γm=125×125 (left)and3 levelmultiscalesolutionwithγm=5×5 (right).NotethattheerrornormarepresentedinTable1.

matrixandfracturesare independentofeach other.Therefore,themaximumcoarseninglevelineachdomain(i.e.,matrix andeveryindividualfracture)canbedifferentdependingonitssizeanditsnumberofgridcells.

Thesurface plotsofthepressuredistributionsobtainedwithfine-scalesimulationsandwithdifferentmultilevel multi-scalesettingsarepresentedinFig.10.

Fig. 11. Testcase3–comparisonofthepressure distributionobtainedbyfine-scalesimulations andF-ADMsingle-phasesolverwithenforced static multilevelmultiscaleonheterogeneousdomain,wherethepressureissolvedeverywhereatthecoarsestlevel.Shownonthetoprowarethepermeability map(right)andthefine-scalesolution(left).Middlerowshows1-levelmultiscalesolutionwithγm=9×9 (left)and2 levelmultiscalesolutionwith

γm=3×3 (right).Bottomrowincludes1 levelmultiscalesolutionwithγm=27×27 (left)and3 levelmultiscalesolutionwithγm=3×3 (right).Note thattheerrornormarepresentedinTable2.

4.3. TestCase3:F-ADMforsingle-phaseflowina2Dheterogeneousfracturedreservoir

A heterogeneous 2D domain of 75m

×

75m with a network of 35 fractures is considered. The matrix and fracture grids are 135

×

135 and 1665,respectively, resulting in 19890 total grid cells. The top right plot on Fig. 11 showsthe heterogeneouspermeabilitymap forthistestcasewiththeminimumandmaximumvaluesof1.3

×

10−16_m2 _and9.9

_×

10−13m2,respectively. Wellsarelocated andhave thesameconstraintsasinthe previoustest case. Differentcoarsening levelsandratiosareconsidered.

Thepressureerrorwithrespecttothereference(i.e.,fine-scalesolution)iscalculatedbasedonEq.(24),whichisshown in Table 2 fordifferent numbers of degreesof freedom (DOF). All multiscale runs employ fully-coupled basis functions. As itcanbe seeninTable 2, forthe correspondingcoarseningratios,the multilevelmultiscaleandthe1-levelmultiscale strategiesdelivercomparableresults.Notethatnoiterationshavebeenemployedtoimprovethemultiscaleresults[36].

4.4. TestCase4:multi-phaseflowina2Dhomogeneousfracturedreservoir

The same fractured homogeneous reservoir presented in the previous test case is considered. The rock matrix is discretized with a 135

×

135 fine-scale grid and 1665 grid cells are employed for the fractures (total domain grid

Table 2

Settingsanderrorsofthemultilevelmultiscalerunsforheterogeneoustestcase.

Num. coarse levels γm Total DOF ep[%]

Case 1 1 3×3 2580 12.8

Case 2 1 9×9 780 12.7

Case 3 2 3×3 410 13.6

Case 4 1 27×27 580 11.0

Case 5 3 3×3 210 12.1

Fig. 12. Pressureplotsfor2DF-ADMtestcase.Shownarethefine-scalesolution(a),decoupledF-ADMwiththresholds0.1 (b)and0.8 (c),andcoupled F-ADMwiththresholds0.1 (d)and0.8 (e),respectively.

Fig. 13. Saturationdistributionoffine-scale(a),decoupledF-ADMwiththresholds0.1 (b)and0.8 (c),andcoupledF-ADMwiththresholds0.1 (d)and 0.8 (e),respectively.

cells: 19890). F-ADM employs 2 coarse levels with coarsening ratio of 3 in each direction both for the matrix and the fractures. Both decoupled and fully-coupled multiscale basis functions are considered aspressure interpolators. The saturation difference between neighbouring cells is used as the F-ADM coarsening criterion with four different thresh-old values:



S

= {

0.1,0.3,0.5,0.8

}

.Three equidistant injection wells(pinj

=

2

·

107Pa) and three equidistant production wells (pprod

=

1

·

107Pa) are present at the left and right boundaries, respectively. The total simulation time is set to 1000 days.

Figs.12and13show thepressureandthesaturationmapsattheendofthesimulationforfine-scaleandF-ADMwith 2 differentthresholdvaluesforthecoarseningcriterion.

Fig. 14. Pressureandsaturationerrors:Fig.14aandFig.14bshowF-ADMerrorofpressureoversimulationtime-stepsfordecoupledandcoupledbasis functionsapproaches.Figs.14cand14d areplotsofsaturationerroroversimulationtime-stepsfordecoupledandcoupledapproachrespectively.Theerrors arecalculatedbasedonEq.(25).

Theerrorforpressureandsaturationiscalculatedas

e

x

=

||xms_||x−_fs||x₂fs||2 (25)

where,x iseitherpressureorsaturationandthesubscriptsms andfs indicatemultiscaleandfine-scale,respectively. Figs.14a,14band14c,14d showtheF-ADM(bothusingdecoupledandcoupledbasisfunctions)pressureandsaturation errorsateachtime-stepforthefourvaluesofthecoarseningcriterion.Figs.15aand15b illustratetheamountofactivegrid cells (as apercentageofthe fine-scalenumberofcells)foreach time-steps. Notethatgrid cells aroundwellsare always keptatthefine-scaleresolution.Here,864 gridcellsarekeptatfine-scaleresolutionduetonear-wellrefinement,whichis 4.74 percentoutofthe8.61 percentofactivegridcellsatthefirsttime-step.Figs.15cand15d providetheaveragepressure andsaturation errorsalong withthe averagepercentageofgrid cells employed byF-ADMasfunctionsofthe coarsening criterionthresholdvalue.

F-ADMresultsshowalowsensitivitytothetypeofbasisfunctionsemployed.Notethatcoupledbasisfunctionsaremore computationally expensiveandprovidewithadenserprolongationoperator.However, incaseofusingstaticprolongation operators, asbasisfunctionsofalllevelsare computedonlyonceatthebeginningofsimulation,theextracomputational costsinvolvedincalculationofcoupledbasisfunctionscanbeconsiderednegligible.

4.5. TestCase5:multi-phaseflowina2Dheterogeneousfracturedreservoir

This testcasediffers fromtest case4 onlyinpermeabilityofmatrixasit isheterogeneous.The restofthe input pa-rameters andconfigurations are identical to Test Case 4. The permeability map is identicalto one used in test case 3, withheterogeneitycontrastof7.7

×

103_{.Thewellpositioningfollowsthesameconfiguration(linedrivewith3equidistant} injectors and producers). Moreover, coarsening level,coarsening ratio and coarseningcriteria are the same. As for their

Fig. 15. Activegridcellsandaveragedparameters:Fig.15aandFig.15bdemonstratethepercentageofactivegridcellsforthetwomentionedapproaches respectively.Thedataisobtainedoversimulationtime-stepsandforfourdifferentF-ADMtolerances.Figs.15cand15d showtheaverageofpressureerror, saturationerrorandpercentageofactivegridcellsoverthewholesimulationtimeforbothapproaches.

similarperformancetodecoupledapproach,onlytheresultsoffullycoupledbasisfunctionsarepresentedforpressureand saturation,inFigs.16and17,respectively.

TheF-ADMerroriscalculatedusingthe sameformulationasintheprevious test cases,i.e.Eq. (25), andpresented in Figs.18aand18b,respectively,forpressureandsaturation.Inaddition,Fig.18c showstheamountofactivegridcellsduring thesimulationfordifferentcoarseningtolerencevalues.Finally,Fig.18d providestheaveragepressureandsaturationerrors alongwiththeaveragepercentageofgridcellsasfunctionsofthecoarseningthresholdvalue.

4.6. TestCase6:multi-phaseflowina3Dhomogeneousfracturedreservoir

A3D75m

×

75m

×

30m containing26 fractureplatesofdifferentgeometricalpropertiesasshowninFig.19.Thematrix isdiscretized witha99

×

99

×

27 Cartesiangridforatotal of264,627 cells.Independent2D gridsare imposed oneach fractureplate.Thegridcellsinthefractureshavesizes

ξ

and



η

suchthatmean

(ξ,



η

)

=

2

×

mean

(

xm

,



ym

,



zm

).

Here,



xm

,



ym

,



zm indicate thesizes ofthe matrix grid cells ineach direction. The total numberof grid cells inthe

fracturenetworkis2,646 gridcells.Therockmatrixpermeabilityis10−14_m2 _{whereasthefracturesoneis2.08}

_·

₁₀−6_m2_. Three injectionwells arepresentinthebottom left,middleleft andtopleft boundary, allperforatingthrough thewhole thickness of the reservoir (pinj

=

2

·

107 Pa). Similarly three productionwells are located on the right-hand side ofthe domain(pprod

=

1

·

107 Pa).Thesimulationtimeis1000 days.

F-ADMemploys2 levelsofcoarseningwithcoarseningratiosequalto 3ineach direction.Fig.20shows, forexample, thesaturationdistributionattime-step15 alongwiththeF-ADMgridemployedatthattime-step.Thesameerrormeasures usedfortheprevioustestcaseareconsidered.Figs.21aand21b presentpressureandsaturationerrorsateachtime-step forfour different thresholdvalues ofthe coarsening criterion,i.e. 0.1, 0.3, 0.5 and 0.8. The evolution ofthe number of active grid cells throughoutF-ADM simulations is shownin Fig. 21c.Averagepressure error, saturation error andactive grid cells over the wholesimulation time versus



S thresholdsare given inFig. 21d.In thistest case, 23328 gridcells

Fig. 16. Pressureplotsof2DF-ADMtestcase(heterogeneous).Shownarefine-scale(a)andF-ADMresultswiththresholdsofS=0.1 (b),0.3 (c),0.5 (d) and0.8 (e),respectively.

near the wellsare kept atfine-scale resolution that is 8.73 percent out of10.51 percentof active gridcells atthe first time step. An increase in thresholdcorrelates with a decrease insize ofthe F-ADM systemas lessactive grid cells are used.

5. Conclusions

Analgebraicdynamic multilevelmethodforfully-implicitsimulationsofmultiphaseflowinporousmediawith embed-deddiscretefractures,namelyF-ADM,ispresentedinthiswork.Builtontheembeddeddiscretefracturemodel(EDFM),the fine-scalefullyimplicitsystemwasmappedintoamultileveldynamicgrid,definedindependentlyformatrixandfractures. This wasachievedby development ofthesequences ofprolongation(localbasis functions)operators forfractured media. Theselocalmultilevelbasisfunctionswereintroducedafterselectionofthecoarsenodesonbothmatrixandfracture sub-domains,withflexiblematrix-fracturecoupling.Thefront-detectionstrategyisusedtoemployfine-scalegridsonlywhere

Fig. 17. Saturationdistributionof2DF-ADMtestcase(heterogeneous).Shownarefine-scale(a)andF-ADMsolutionsforthresholdsS=0.1 (b),0.3 (c), 0.5 (d)and0.8 (e),respectively.

needed,whilemultilevelmultiscalegridisusedelsewhere.Theuseofmultiscalebasisfunctionsguaranteestheaccuracyof theglobal(pressure)unknownswherecoarsegridsareimposed.

Numerical results for both 2D and3D test cases were presented to validate the EDFM fully-implicit implementation and investigate its accuracy. Then, the sensitivity of F-ADM to the type of pressure interpolator (i.e., with andwithout matrix-fracture coupling) and the fraction of active grid cells chosen indirectly by different threshold values was stud-ied.Thesingle-phaseflowresultsshowthattheproposed multilevelinterpolationstrategy forembedded discretefracture systems isaccurate compared tothe original 1-level multiscale finite-volumemethod. Multiphase flow results,with dif-ferent amountof active dynamic grids, show that F-ADMis able toprovide accurate results forflow in fractured media by employing only a fraction of the fine-scale grid-cells both in the rock matrix and in the fractures. It is expected that the greaterthe size ofthe domain,thelower the percentageof activegrid cells compared to theglobalnumber of

Fig. 18. F-ADMerrorhistoryforpressure(a)andsaturation(b)calculatedbasedonEq.(25) duringthesimulation.Thepercentageofactivegridcells(c) andtheaverageerrorsoversimulationtime(d)arealsopresented.

Fig. 20. Illustration of the F-ADM mesh at time-step 15 (left) and saturation distribution (right). Only saturation values higher than 0.5 are shown (right).

Fig. 21. F-ADMerrorhistoryforpressure(a)andsaturation(b)calculatedbasedonEq.(25) duringthesimulation.Thepercentageofactivegridcells(c) andtheaverageerrorsoversimulationtime(d)arealsopresented.

fine-scalegrids. Assuch, F-ADMcastsapromising approachforsimulationofmultiphase flowinreal-field fractured me-dia.

Acknowledgements

The authors thank all members of the Delft Advanced Reservoir Simulation (DARSim) research team for the fruitful discussionsduringthedevelopmentoftheF-ADMmethod.

References

[1]B.Berkowitz,Characterizingflowandtransportinfracturedgeologicalmedia:areview,Adv.WaterResour.25(2002)861–884.

[2]T.Y.Hou,X.-H.Wu,Amultiscalefiniteelementmethodforellipticproblemsincompositematerialsandporousmedia,J.Comput.Phys.134(1997) 169–189.

[3]Y.Efendiev,T.Hou,T.Strinopoulos,Multiscalesimulationsofporousmediaflows inflow-basedcoordinatesystem,Comput.Geosci.12 (3)(2008) 257–272.

[4] Y.Efendiev,S.Lee,G.Li,J.Yao,N.Zhang,Hierarchicalmultiscalemodelingforflowsinfracturedmediausinggeneralizedmultiscalefiniteelement method,GEM,Int.J.Geomath.6 (2)(2015)141–162,https://doi.org/10.1007/s13137-015-0075-7.

[5]P.Jenny,S.H.Lee,H.A.Tchelepi,Multi-scalefinite-volumemethodforellipticproblemsinsubsurfaceflowsimulation,J.Comput.Phys.187(2003) 47–67.

[6]P.Jenny,S.H.Lee,H.A.Tchelepi,Adaptivefullyimplicitmulti-scalefinite-volumemethodformulti-phaseflowandtransportinheterogeneousporous media,J.Comput.Phys.217(2006)627–641.

[7]H.Hajibeygi,G.Bonfigli,M.Hesse,P.Jenny,Iterativemultiscalefinite-volumemethod,J.Comput.Phys.227(2008)8604–8621.

[8]I.Lunati,S.H.Lee,Anoperatorformulationofthemultiscalefinite-volumemethodwithcorrectionfunction,MultiscaleModel.Simul.8 (1)(2009) 96–109.

[9]Y.Wang,H.Hajibeygi,H.A.Tchelepi,Monotonemultiscalefinitevolumemethod,Comput.Geosci.(2015)1–16.

[10]O.Moyner,K.-A.Lie,Amultiscalerestriction-smoothedbasismethodforhighcontrastporousmediarepresentedonunstructuredgrids,J.Comput. Phys.304(2016)46–71.

[11]M.Berger,J.Oliger,Adaptivemeshrefinementforhyperbolicpartialdifferentialequations,J.Comput.Phys.53(1984)484–512.

[12]D.Han,C.Yan,L.Peng,Amoreflexibleapproachofdynamiclocalgridrefinementforreservoirmodeling,in:SPE16014,SPESymposiumonReservoir Simulation,1–4February,SanAntonio,Texas,USA,1987.

[13]G.Schmidt,F.Jacobs,Adaptivelocalgridrefinementandmulti-gridinnumericalreservoirsimulation,J.Comput.Phys.77(1988)140–165.

[14]W.Mulder,R.G.Meyling,Numericalsimulationoftwo-phaseflowusinglocallyrefinedgridsinthreespacedimensions,in:SPE21230,in:SPEAdv. Technol.Ser.,vol. 1,1991,pp. 36–41.

[15]M.Biterge,T.Ertekin,Developmentandtestingofastatic/dynamiclocalgrid-refinementtechnique,J.Pet.Technol.44(1992)487–495.

[16]M.Edwards,M.Christie,DynamicallyadaptiveGodunov schemeswithrenormalizationinreservoirsimulation,in:SPESymposiumonReservoir Simu-lation,28February–3March,NewOrleans,Louisiana,1993,SPEpaper25268.

[17]M.Edwards,Ahigher-orderGodunov schemecoupledwithdynamicallocalgridrefinementforflowinaporousmedium,Comput.MethodsAppl. Mech.Eng.131(1996)287–308.

[18]Z.Heinemann,G.Gerken,G.vonHantelmann,Usinglocalgridrefinementinamultiple-applicationreservoirsimulator,in:SPEReservoirSimulation Symposium,15–18November,SanFrancisco,California,USA,1983,SPEpaper12255.

[19]B.Faigle,R.Helmig,I.Aavatsmark,B.Flemisch,EfficientmultiphysicsmodellingwithadaptivegridrefinementusingaMPFAmethod,Comput.Geosci. 18(2014)625–636.

[20]M.Cusini,C.vanKruijsdijk,H.Hajibeygi,Algebraicdynamicmultilevel(ADM)methodforfullyimplicitsimulations ofmultiphaseflowinporous media,J.Comput.Phys.314(2016)60–79.

[21]M.Cusini,B.Fryer,C.vanKruijsdijk,H.Hajibeygi,Algebraicdynamicmultilevelmethodforcompositionalflowinheterogeneousporousmedia,J. Comput.Phys.354(2018)593–612.

[22]P.Dietrich,R.Helmig,M.Sauter,H.Hotzl,J.Kongeter,G.Teutsch,FlowandTransportinFracturedPorousMedia,Springer,2005.

[23]M.Karimi-Fard,L.Durlofsky,K.Aziz,Anefficientdiscrete-fracturemodelapplicableforgeneral-purposereservoirsimulators,SPEJ.(2004)227–236.

[24]R.Ahmed,M.G.Edwards,S.Lamine,B.A.Huisman,M.Pal,Control-volumedistributedmulti-pointfluxapproximationcoupledwithalower-dimensional fracturemodel,J.Comput.Phys.284(2015)462–489.

[25]V.Reichenberger,H.Jakobs, P.Bastian,R.Helmig,Amixed-dimensionalfinitevolumemethodfortwo-phaseflowinfracturedporousmedia,Adv. WaterResour.29(2006)1020–1036.

[26]A.Moinfar,A.Varavei,K.Sepehrnoori,R.T.Johns,Developmentofanefficientembeddeddiscretefracturemodelfor3dcompositionalreservoir simu-lationinfracturedreservoirs,SPEJ.19(2014)289–303.

[27]T.Sandve,I.Berre,J.Nordbotten,Anefficientmulti-pointfluxapproximationmethodfordiscretefracture–matrixsimulations,J.Comput.Phys.231 (2012)3784–3800.

[28] T.Sandve,E.Keilegavlen,J.Nordbotten,Physics-basedpreconditionersforflowinfracturedporousmedia,WaterResour.Res.50(2014)1357–1373,

https://doi.org/10.1002/2012WR013034.

[29] S.Geiger-Boschung,S.Matthäi,J.Niessner,R.Helmig,Black-oilsimulationsforthree-component,three-phaseflowinfracturedporousmedia,SPEJ. 14 (2)(2009)338–354,https://doi.org/10.2118/107485-PA.

[30] A.Fumagalli,S.Zonca,L.Formaggia,Advancesincomputationoflocalproblemsforaflow-basedupscalinginfracturedreservoirs,Math.Comput. Simul.137(2017)299–324,https://doi.org/10.1016/j.matcom.2017.01.007.

[31] A. Sakhaee-Pour,M.Wheeler, Effectiveflowproperties forcellscontainingfracturesofarbitrarygeometry,SPEJ.21 (3)(2016)965–980,https:// doi.org/10.2118/167183-PA.

[32] A.Fumagalli,L.Pasquale,S.Zonca,S.Micheletti,Anupscalingprocedureforfracturedreservoirswithembeddedgrids,WaterResour.Res.52 (8)(2016) 6506–6525,https://doi.org/10.1002/2015WR017729.

[33] J.Natvig,B.Skaflestad,F.Bratvedt,K.Bratvedt,K.-A.Lie,V.Laptev,S.Khataniar,Multiscalemimeticsolversforefficientstreamlinesimulationof fracturedreservoirs,SPEJ.16 (4)(2011)880–888,https://doi.org/10.2118/119132-PA.

[34]H.Hajibeygi,D.Karvounis,P.Jenny,Ahierarchicalfracturemodelfortheiterativemultiscalefinitevolumemethod,J.Comput.Phys.230 (24)(2011) 8729–8743.

[35]T.H.Sandve,I.Berre,E.Keilegavlen,J.M.Nordbotten,Multiscalesimulationofflowandheattransportinfractured geothermalreservoirs:inexact solversandimprovedtransportupscaling,in:Thirty-EighthWorkshoponGeothermalReservoirEngineering,Stanford,California,USA,2013.

[36]M.Tene,M.S.AlKobaisi,H.Hajibeygi,Algebraicmultiscalemethodforflowinheterogeneousporousmediawithembeddeddiscretefractures(F-AMS), J.Comput.Phys.321(2016)819–845.

[37]S.Shah,O.Møyner,M.Tene,K.-A.Lie,H.Hajibeygi,Themultiscalerestrictionsmoothedbasismethodforfracturedporousmedia(F-MSRSB),J.Comput. Phys.318(2016)36–57.

[44]J.H.Norbeck,M.W.McClure,J.W.Lo,R.N.Horne,Anembeddedfracturemodelingframeworkforsimulationofhydraulicfracturingandshear stimula-tion,Comput.Geosci.20 (1)(2016)1–18.

[45]S.Pluimers,HierarchicalFractureModelingApproach,Master’sthesis,TUDelft,DepartmentofGeoscienceandEngineering,2015.

[46]M.Tene,S.B.Bosma,M.S.AlKobaisi,H.Hajibeygi,Projection-basedembeddeddiscretefracturemodel(PEDFM),Adv.WaterResour.105(2017)205–216.

[47]B.Flemisch,I.Berre,W.Boon,A.Fumagalli,N.Schwenck,A.Scotti,I.Stefansson,A.Tatomir,Benchmarksforsingle-phaseflowinfracturedporous media,Adv.WaterResour.111(2018)239–258.

[48] D.vanBatenburg,M.Bosch,P.Boerrigter,A.deZwart,J.Vink,ApplicationofdynamicgriddingtechniquestoIOR/EOR-processes,in:SPEReservoir SimulationSymposium,21–23February,TheWoodlands,Texas,USA,2011,2011,pp. 1–16,https://doi.org/10.2118/141711-MS,SPEpaper141711.