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

(1)

Delft University of Technology

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

method

Kadeethum, T.; Nick, H. M.; Lee, S.; Ballarin, F.

DOI

10.1016/j.advwatres.2020.103620

Publication date

2020

Document Version

Final published version

Published in

Advances in Water Resources

Citation (APA)

Kadeethum, T., Nick, H. M., Lee, S., & Ballarin, F. (2020). Flow in porous media with low dimensional

fractures by employing enriched Galerkin method. Advances in Water Resources, 142, [103620].

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

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)

ContentslistsavailableatScienceDirect

Advances

in

Water

Resources

journalhomepage:www.elsevier.com/locate/advwatres

Flow

in

porous

media

with

low

dimensional

fractures

by

employing

enriched

Galerkin

method

T. Kadeethum

a

_,

_H.M.

_Nick

a,b,∗

_,

_S.

_Lee

c

_,

_F.

_Ballarin

d a Technical University of Denmark, Denmark

b Delft University of Technology, the Netherlands c Florida State University, Florida, USA d mathLab, Mathematics Area, SISSA, Italy

a

r

t

i

c

l

e

i

n

f

o

Keywords:

Fractured porous media Mixed-dimensional Enriched Galerkin Finite element method Heterogeneity Local mass conservative

a

b

s

t

r

a

c

t

ThispaperpresentstheenrichedGalerkindiscretizationformodelingfluidflowinfracturedporousmediausing themixed-dimensionalapproach.Theproposedmethodhasbeentestedagainstpublishedbenchmarks.Since fractureandporousmediadiscontinuitiescansignificantlyinfluencesingle-andmulti-phasefluidflow,the het-erogeneousandanisotropicmatrixpermeabilitysettingisutilizedtoassesstheenrichedGalerkinperformancein handlingthediscontinuitywithinthematrixdomainandbetweenthematrixandfracturedomains.Ourresults illustratethattheenrichedGalerkinmethodhasthesameadvantagesasthediscontinuousGalerkinmethod;for example,itconserveslocalandglobalfluidmass,capturesthepressurediscontinuity,andprovidestheoptimal errorconvergencerate.However,theenrichedGalerkinmethodrequiresmuchfewerdegreesoffreedomthan thediscontinuousGalerkinmethodinitsclassicalform.Thepressuresolutionsproducedbybothmethodsare similarregardlessoftheconductiveornon-conductivefracturesorheterogeneityinmatrixpermeability.This analysisshowsthattheenrichedGalerkinschemereducesthecomputationalcostswhileofferingthesame ac-curacyasthediscontinuousGalerkinsothatitcanbeappliedforlarge-scaleflowproblems.Furthermore,the resultsofatime-dependentproblemforathree-dimensionalgeometryrevealthevalueofcorrectlycapturingthe discontinuitiesasbarriersorhighly-conductivefractures.

1. Introduction

Modelingof fluidflow in fractured porousmedia is essential for awidevarietyof applicationsincludingwaterresourcemanagement (Glaser etal., 2017;Penget al.,2017),geothermal energy(Willems andNick,2019;Salimzadehetal.,2019a;SalimzadehandNick,2019), oil andgas (Wheeler etal., 2019; Kadeethum et al., 2019c; Andri-anovand Nick,2019; Kadeethum et al., 2020c),induced seismicity (Rinaldi andRutqvist, 2019), CO₂ sequestration (Salimzadeh et al., 2018),andbiomedicalengineering(Vinjeetal.,2018;RuizBaieretal., 2019;Kadeethum etal., 2020a).Afractured porousmedium canbe decomposedintothebulkmatrixandfracturedomains,whichare gen-erallyanisotropic,heterogeneous,andhavesubstantiallydiscontinuous materialpropertiesthatcanspanseveralordersofmagnitude(Matthai andNick,2009;Flemischetal.,2018;Jiaetal.,2017;Bisdometal., 2016).These discontinuitiescancriticallyenhanceorhindertheflux withinandbetweenthebulkmatrixandfracturedomains.Accurately

∗_{corresponding}_author_at:_Technical_University_of_Denmark,_Lyngby,_Denmark,_Delft_University_of_Technology,_The_Netherlands.

E-mailaddresses:teekad@dtu.dk(T.Kadeethum),h.m.nick@tudelft.nl(H.M.Nick),lee@math.fsu.edu(S.Lee),francesco.ballarin@sissa.it(F.Ballarin).

capturingtheﬂowbehaviorcontrolledbythesediscontinuitiesin com-plex mediais stillchallenging(NickandMatthai,2011a;De Dreuzy etal.,2013;Flemischetal.,2016;HoteitandFiroozabadi,2008;Zhang etal.,2016).

Therearetwomainapproachestorepresentthefluidflowbetween thematrixandfracturedomain(NickandMatthai,2011b;Flemisch etal.,2018;Juanesetal.,2002).Thefirstmodel,anequi-dimensional model,discretizesthematrixandfracturedomainwithsame dimension-ality(Salinasetal.,2018).Thisapproachisstraightforwardto imple-ment,andnocouplingconditionisrequired.Thisapproachisutilized tomodel,forexample,coupledhydromechanicaloffracturedrocks us-ingthefinite-discreteelementmethod(Lathametal.,2013;2018),and canalsocapturefracturepropagationusinganimmersed-bodymethod (Obeysekaraetal.,2018)orphase-fieldapproach(Santillanetal.,2018; Leeetal., 2016b;2018).Thesecond model,whichwecalla mixed-dimensional modelhereafter,reducesthefracturedomaintoalower dimensionalitybyassumingthefracturethicknessismuchsmaller

com-https://doi.org/10.1016/j.advwatres.2020.103620

Received21December2019;Receivedinrevisedform6May2020;Accepted10May2020 Availableonline31May2020

(3)

Fig. 1. Comparisonofdegreesoffreedomforlinearpolynomialcaseamong(a ) continuousGalerkin(CG),(b )discontinuousGalerkin(DG),and(c )enriched Galerkin(EG)functionspaces.

paredtothesizeofmatrixdomain(Boonetal.,2018;Martinetal., 2005; Berroneetal., 2018).The secondapproach hasseveral bene-ﬁts;forexample,itreducesthedegreesoffreedom(DOF)(Nickand Matthai,2011a)asthefracturedomainisrepresentedastheinterface, whichispartofthematrixdomain,andsubsequently,thisapproachcan improvemeshquality(reducethemeshskewness)(Matthaietal.,2010). Sincethefracturesareinterfaces,onecanusealargermeshsize,which satisﬁesCourant-Friedrichs-Lewy(CFL)conditionmoreeasily(Juanes etal.,2002;NickandMatthai,2011b).

Inthepastdecades,manyapproacheshavebeenproposedtomodel thefracturedporousmediausingthemixed-dimensionalapproach;(1) two-point flux approximationin unstructuredcontrol volume finite-differencetechnique(Karimi-Fardetal.,2004),(2)multi-pointflux ap-proximationusingmixedfiniteelementmethodongeneralquadrilateral andhexahedralgrids(Wheeleretal.,2012),(3)eXtendedfiniteelement combinedwithmixedfiniteelementformulation(D’AngeloandScotti, 2012;PrevostandSukumar,2016;SanbornandPrevost,2011),(4) em-beddeddiscretefracture-matrix(DFM)modelingwithnon-conforming mesh(Hajibeygietal.,2011;Odsaeteretal.,2019),(5)mixed approx-imationsuchasmimeticfinitedifference(FlemischandHelmig,2008; Formaggia etal.,2018),(6)two-fieldformulationusingmixed finite element(MFE)(Martinetal.,2005;Fumagallietal.,2019),and(7) dis-coninuousGalerkin(DG)method(Rivieetal.,2000;Hoteitand Firooz-abadi,2008;Antoniettietal.,2019;Arnoldetal.,2002).

WefocusonthefiniteelementbaseddiscretizationsuchastheDG and MFEmethodsthat ensurethe localmass conservative property. Moreover, they areflexibleenough todiscretize complexsubsurface geometriessuch asintersectionsof fracturesor irregular-shaped ma-trix blocks.Additionally,theaforementioned methodsarecapable of mimickingthefracturepropagationinporoelasticmediausingeither cohesivezonemethodorlinearelasticfracturemechanicsframework (SalimzadehandKhalili,2015;Salimzadehetal., 2019b;Secchiand Schrefler, 2012;SeguraandCarol,2008).However,theMFEmethod requiresanadditionalprimaryvariable(fluidvelocity)whichmay re-quiremorecomputationalresources,especiallyinathree-dimensional

Fig. 2. Comparisonofratio ofdegreesof freedomfor1st,2nd,3rd,4th,and5th poly-nomialdegreecasesontriangularelement amongCG,EG,andDGdiscretizations:(a ) ratioofEGoverCGDOF,(b )ratioofDG overCGDOF,and(c )ratioofDGoverEG DOF.

(4)

Fig. 3. Illustrationof (a ) equi-dimensional and(b ) mixed-dimensionalsettings.Notethat these illustra-tionsarethegraphicalrepresentations;inthe numeri-calmodel𝜕Ω_p,𝜕Ω_q,𝜕Γ_p,or𝜕Γ_qhavetobeimposedon allboundaryfacesofeachdomain.

Fig. 4. IllustrationofEGelements(a )withoutfracture interfaceand(b )withfractureinterface.

domainKadeethumetal.(2019a).Meshadaptivityisalsonot straight-forwardtoimplementLeeandWheeler(2017),anditrequiresthe inver-sionofthepermeabilitytensor,whichmayleadtoanill-posedproblem (ChooandLee,2018).TheDGmethodalsocanbeconsideredasa com-putationallyexpensivemethodasitrequiresalargenumberofDOF(Sun andLiu,2009;Leeetal.,2016a).

Toresolvesomeoftheshortcomingsmentionedabove,wepropose anenrichedGalerkin(EG)discretization(Leeetal., 2016a;Zi etal., 2004;Khoeietal.,2018)tomodelfluidflowinfracturedporousmedia usingthemixed-dimensionalapproach.TheEGmethod,utilizedinthis study,composesoftheCGfunctionspaceaugmentedbya piecewise-constantfunctionatthecenterofeachelement.Thismethodhasthe sameinteriorpenaltytypebilinearformastheDGmethod(SunandLiu, 2009;Leeetal.,2016a).TheEGmethod,however,onlyrequirestohave discontinuousconstantsasillustratedinFig.1,soithasfewerDOFthan theDGmethod.Fig.2presentsthecomparisonoftheDOFratioamong CG,EG,andDGmethods,anditshowsthattheEGmethodrequires halfoftheDOFneededbytheDGmethod(triangularelementswiththe firstpolynomialdegreeapproximation).Notethatthisratiodecreases asthepolynomialdegreeapproximationincrease.TheEGmethodhas beendevelopedtosolvegeneralellipticandparabolicproblemswith dy-namicmeshadaptivity(LeeandWheeler,2017;2018;Leeetal.,2018) andextendedtoaddressthemultiphasefluidflowproblems(Leeand Wheeler,2018).Recently,theEGmethodhasbeenalsoappliedtosolve thenon-linearporoelasticproblem(ChooandLee,2018; Kadeethum etal.,2019b;2020b),andcompareditsperformancewithother two-andthree-fieldformulationmethods(Kadeethumetal.,2019a).Tothe bestofourknowledge,thisisthefirstattempttoapplytheEG discretiza-tioninthemixed-dimensionalsetting.

Therestofthepaperisorganizedasfollows.Themethodology sec-tionincludesmodeldescription,mathematicalequations,andtheir dis-cretizationsfortheEGandDGmethods.Subsequently,theblock struc-ture usedtocomposetheEG functionspaceandthecouplingterms betweenmatrixandfracturedomainsisillustrated.Thenumerical ex-amplessectionpresentsﬁveexamples,andtheconclusionisﬁnally pro-vided.

2. Methodology

2.1. Governingequations

Wefirstbrieflyintroducetheequi-dimensionalmodel,whichisused toderivethemixed-dimensionalmodel.Weareinterestedinsolving steady-stateandtime-dependent singlephase fluidflow in fractured porousmediaonΩ,whichcomposesofmatrixandfracturedomains, Ω_m andΩ_f,respectively. Let Ω⊂ ℝ𝑑 be thedomainof interestind -dimensionalspacewhere𝑑=1,2,or3boundedbyboundary,𝜕Ω.𝜕Ω canbedecomposedtopressureandfluxboundaries,𝜕Ωpand𝜕Ωq, re-spectively.Thetimedomainisdenotedby𝕋 =(0,𝜏],where𝜏 >0isthe finaltime.

Theillustrationoftheequi-dimensionalmodelisshowninFig.3a. Thismodelcomposesoftwomatrixsubdomains,Ω_miwhere𝑖=1,2,and onefracturesubdomain,Ωf.Notethat,forthesakeofsimplicity,this setup isusedtoillustratethegoverningequationwithonlytwo ma-trixsubdomains,butinageneralcase,thedomainmaycomposeofnm subdomains,i.e.𝑖=1,2,…,𝑛_𝑚.Moreover,thedomainmaycontainup tonffractures,wherenmandnfarenumberofmatrixsubdomainand fracture,respectively. Forsimplicity,in thissectionwe willconsider

𝑛𝑓=1.Thefracturesmaynotcutthroughthematrixdomain,whichwe callimmersedfracturesetting.ThistopicwillbediscussedinSection3. Thegoverningsystemwithinitialandboundaryconditionsofthe equi-dimensionalmodelassumingaslightlycompressibleﬂuidforthematrix domainispresentedbelow:

𝑐𝜙𝑚𝑖_𝜕𝑡𝜕(𝑝𝑚𝑖)−∇⋅𝒌_𝜇𝒎𝑖(∇𝑝𝑚𝑖−𝜌𝐠)=𝑔𝑚𝑖inΩ𝑚𝑖×𝕋, for𝑖=1,2, (1)

𝑝𝑚=𝑝𝑚𝐷on𝜕Ω𝑝×𝕋, (2)

−𝒌𝒎𝑖

𝜇 (∇𝑝𝑚𝑖−𝜌𝐠)⋅ 𝒏=𝑞𝑚𝑁on𝜕Ω𝑞×𝕋, (3)

(5)

Fig. 5. (a )geometryusedfortheerrorconvergenceanalysis,(b )illustrationoftheexactsolution,errorconvergenceplotbetweenEGandDGmethodswith polynomialapproximationdegree1and2of(c )matrixand(d )fracturedomains.Notethatfor(c )theslopesofthebestﬁttedlineare2.235(EG),2.022(DG)for thepolynomialapproximationdegree1,and3.137(EG),and3.410(DG)forthepolynomialapproximationdegree2.Theslopesofthebestﬁttedlinefor(d )are 2.289(EG),2.313(DG)forthepolynomialapproximationdegree1,and3.194(EG),and2.948(DG)forthepolynomialapproximationdegree2.

andforthefracturedomainis:

𝑐𝜙𝑓_𝜕𝑡𝜕(𝑝𝑓)−∇⋅𝒌_𝜇𝒇(∇𝑝𝑓−𝜌𝐠)=𝑔𝑓inΩ𝑓×𝕋, (5)

𝑝𝑓=𝑝𝑓𝐷on𝜕Ω𝑝×𝕋, (6)

−𝒌𝒇

𝜇 (∇𝑝𝑓−𝜌𝐠)⋅ 𝒏=𝑞𝑓𝑁on𝜕Ω𝑞×𝕋, (7)

𝑝𝑓=𝑝0𝑓 inΩ𝑓×𝕋 at𝑡=0, (8)

where(· )_irepresentsanindex,cistheﬂuidcompressibility,𝜙mand𝜙f arethematrixandfractureporosity(𝜙fisassumedtobeone),kmandkf

arethematrixandfracturepermeabilitytensor,respectively,𝜇 isﬂuid viscosity,p_mandp_farematrixandfracturepressure,respectively,𝜌 is

ﬂuiddensity,gisthegravitationalvector,gmandgfaresink/sourcefor matrixandfracturedomains,respectively,nisanormalunitvectorto anysurfaces,p_mDandp_fDareprescribedpressureformatrixandfracture domains,respectively,qmN andqfN areprescribedﬂuxformatrixand fracturedomains,respectively,and𝑝_𝑚0

𝑖 and𝑝0𝑓 areprescribedpressure formatrixandfracturedomainsat𝑡=0,respectively.

Toformulatethemixed-dimensionalsettingaspresentedinFig.3b, we integrate along the normal direction to the fracture plane

Martinetal.(2005).AsaresultΩ_fisreducedtoaninterface,Γ.Note thatthegoverningequationofthemixed-dimensionalsettingusedin this studyis proposedbyMartinetal.(2005)inthemixedfinite el-ement formulation, which uses fluid pressure and fluid velocity as the primaryvariables. Themixed-dimensionalsetting has beenused in themixedformulation(Keilegavlenetal.,2017) oradapted to fi-nitevolumediscretization(Glaseretal.,2017;Stefanssonetal.,2018;

(6)

Fig. 6. (a )geometryandboundaryconditions usedforthequarterﬁve-spotproblemand(b ) meshthathasℎ=1.2× 10−1_.

Fig. 7. Pressuresolutionwithℎ=7.2× 10−2_,_p

m,inthematrixofquarterﬁve-spotexample:(a )permeablefracture,(b )impermeablefracturecases,pmplotalong 𝑥=𝑦lineof(c )permeablefracture,and(d )impermeablefracturecases.NotethatwedigitizetheresultsofthereferencesolutionfromAntoniettietal.(2019).Note thattheresultsobtainedwiththeEGandDGmethodsoverlap.

(7)

Fig. 8. Immersed fracture geometry and boundaryconditionsof(a )permeablefracture, (b ) partially permeable fracture, (c ) imper-meablefracturecases,and(d )meshthathas

𝑛𝑒=272.A’isfracturetip.

Glaser et al., 2019) and DG discretization on polytopic grids (Antoniettietal.,2019).Themixed-dimensionstrongformulationand itsboundaryconditionsforthematrixdomainaresimilartothe equi-dimensionalmodel,(1)to(3),butthestrongformulationandits bound-aryconditionsofthefracturedomainare:

𝑐𝑎𝑓_𝜕𝑡𝜕(𝑝𝑓)−∇𝑇⋅ 𝑎𝑓 𝒌𝑻 𝒇 𝜇 ( ∇𝑇𝑝_𝑓−𝜌𝐠)=𝑔_𝑓+−𝒌𝒎 𝜇 ( ∇𝑝_𝑚−𝜌𝐠)inΓ ×𝕋, (9) 𝑝𝑓=𝑝𝑓𝐷on𝜕Γ𝑝×𝕋, (10) − 𝒌𝑻 𝒇 𝜇 (∇𝑝𝑓−𝜌𝐠)⋅ 𝒏=𝑞𝑓𝑁on𝜕Γ𝑞×𝕋, (11) 𝑝𝑓=𝑝0𝑓 inΓ ×𝕋at𝑡=0, (12)

where𝜕Γpand𝜕Γqrepresentpressureandﬂuxboundariesofthe frac-ture domain, respectively, 𝒌𝑻

𝒇 is the tangentialfracture permeability

tensor,∇Tand∇T· arethetangentialgradientanddivergence opera-tors,whicharedeﬁnedas∇𝑇(⋅)=∇(⋅)−𝒏[𝒏⋅ ∇(⋅)]andtr(∇𝑇(⋅)), respec-tively,tr(⋅)istraceoperator,afisafractureaperture,−𝒌_𝜇𝒎

(

∇𝑝_𝑚−𝜌𝐠) representstheﬂuidmasstransferbetweenmatrixandfracturedomains

Martinetal.(2005),⋅isjumpoperator,whichwillbediscussedlater inthediscretizationpart,andpfDandqfNarespeciﬁedpressureandﬂux forthefracturedomain,respectively.

Inthisstudy,if𝑑=3,k_masafulltensorisdeﬁnedas:

𝒌𝒎∶= ⎡ ⎢ ⎢ ⎣ 𝑘𝑥𝑥 𝑚 𝑘𝑥𝑦𝑚 𝑘𝑥𝑧𝑚 𝑘𝑦𝑥𝑚 𝑘𝑦𝑦𝑚 𝑘𝑦𝑧𝑚 𝑘𝑧𝑥 𝑚 𝑘𝑧𝑦𝑚 𝑘𝑧𝑧𝑚 ⎤ ⎥ ⎥ ⎦, (13)

wherealltensorcomponentscharacterizethetransformationofthe com-ponentsofthegradientofﬂuidpressureintothecomponentsofthe ve-locityvector.The𝑘𝑥𝑥

𝑚,𝑘𝑦𝑦𝑚,and𝑘𝑧𝑧𝑚 representthematrixpermeabilityin x-,y-,andz-direction,respectively.kf,ontheotherhand,composesof

twocomponents: 𝒌𝑓∶= [ 𝒌𝑻 𝒇 0 0 𝑘𝑛_𝑓 ] , (14) where𝑘𝑛

𝑓isthenormalfracturepermeability.Notethatwepresenthere ageneralformof𝒌𝑻

𝒇,tobespecific,𝒌𝑻𝒇 isascalarif𝑑=2andtensorif 𝑑=3.Torepresentthefluidmasstransferbetweenmatrixandfracture domains,followingMartinetal.(2005), Antoniettietal.(2019),we definethecouplingconditionsbetweenthematrixandfracturedomain as: ⟨−𝒌_𝜇𝒎(∇𝑝_𝑚−𝜌𝐠)⟩ ⋅ 𝒏= 𝛼𝑓 2 ( 𝑝𝑚1−𝑝𝑚2 ) onΓ, (15) −𝒌𝒎 𝜇 ( ∇𝑝_𝑚−𝜌𝐠)= 2𝛼𝑓 2𝜉 −1 ( ⟨𝑝𝑚⟩ −𝑝𝑓)onΓ, (16) where⟨ · ⟩ isanaverageoperator,whichwillbepresentedlaterinthe discretizationpart,𝛼frepresentsaresistantfactorofthemasstransfer

(8)

Fig. 9. Pressuresolution,pm,inthematrixof

immersedfractureexample using𝑛𝑒=6,698:

(a )permeablefracture,(b )partiallypermeable fracture,and(c )impermeablefracturecases.

betweenthefractureandmatrixdomainsdeﬁnedas:

𝛼𝑓= 2𝑘𝑛_𝑓

𝑎𝑓 , (17)

and𝜉 ∈ (0.5,1.0].Inthispaper,weset𝜉 =1.0forthesakeofsimplicity. Ingeneralcase,𝜉 isusedtorepresentafamilyofthemixed-dimensional model, and more details can be found in Martin et al. (2005),

Antoniettietal.(2019).

2.2. Numericaldiscretization

Inthissection,thediscretizationofthemixed-dimensionalmodelis illustrated.ThedomainΩispartitionedintoneelements,ℎ,whichis thefamilyofelementsT(trianglesin2D,tetrahedronsin3D).Wewill furtherdenotebyeafaceofTasillustratedinFig.4.Wedenoteh_Tasthe diameterofTanddeﬁneh≔ max(hT),whichwemayreferasmeshsize. Letℎdenotesthesetofallfacets,forthematrixdomain,ℎ0theinternal facets,_ℎ𝐷theDirichletboundaryfaces,and_ℎ𝑁denotestheNeumann boundaryfaces.FollowingLeeetal.(2016a)forany𝑒∈0

ℎlet𝑇+and

𝑇− _be_two_neighboring_elements_such_that_𝑒₌_𝜕_𝑇+_∩_𝜕_𝑇−_._Let_𝒏+ _and

𝒏− _be_the_outward_normal_unit _vectors_to_𝜕𝑇+ _and_𝜕𝑇−_,_respectively

(seeFig.4a).ℎ=ℎ0 ∪ ℎ𝐷 ∪ ℎ𝑁,ℎ0 ∩ ℎ𝐷=∅,ℎ0 ∩ ℎ𝑁=∅,and 𝐷

ℎ ∩ℎ𝑁=∅.Wefurtherdeﬁneℎ1∶=ℎ0∪ℎ𝐷.

Thefracturedomainisconformingwith_ℎ,anditisnamedΓh (in-tervalsin2D,trianglesin3Ddomain)hereafteraspresentedinFig.4b.

Wewillfurtherdenotebye_fafaceofΓ_h.LetΛ_hdenotesthesetofall facets,forthefracturedomain,Λ0

ℎ theinternalfacets,Λ𝐷ℎ the Dirich-letboundaryfaces,andΛ𝑁_ℎ denotestheNeumannboundaryfaces,and Λ𝐷_ℎ∩ Λ𝑁_ℎ =∅.LetnbetheoutwardnormalunitvectortoΓ_h,whichis coincidedwith𝒏+ _of_the_𝜕𝑇+_.

Next,wedeﬁnethejumpoperatorofthescalarandvectorvalues as:

𝑋⋅ 𝒏∶=𝑋+𝒏++𝑋−𝒏−and

𝑿⋅ 𝒏∶=𝑿+⋅ 𝒏++𝑿−⋅ 𝒏−, (18) respectively,Moreover,byassumingthatthenormalvectornisoriented from𝑇+_to_𝑇−_,_we_obtain

𝑋∶=𝑋+−𝑋− (19)

FollowingLeeetal.(2016a),Scovazzietal.(2017),theweighted aver-ageisdeﬁnedas:

⟨𝑋⟩𝛿𝑒=𝛿𝑒𝑋++(1−𝛿𝑒)𝑋−, (20) where 𝛿𝑒∶= 𝑘𝑚 − 𝑒 𝑘𝑚+𝑒 +𝑘𝑚−𝑒, (21) and 𝑘𝑚+𝑒 ∶=(𝒏+)𝑇⋅ 𝒌𝒎+⋅ 𝒏+,

(9)

Fig. 10. pmplotalong𝑦=0.75lineofimmersedfractureexample:(a )permeablefracture,(b )partiallypermeablefracture,(c )impermeablefracturecases.Note

thatwedigitizethereferenceresults,Angotetal.2009,fromAngotetal.(2009).AlltheresultsobtainedfromtheEGandDGmethodswithdiﬀerentmeshsizes overlap.

𝑘𝑚−𝑒 ∶=(𝒏−)𝑇⋅ 𝒌𝒎−⋅ 𝒏−, (22)

andaharmonicaverageof𝑘𝑚𝑒+and𝑘𝑚−𝑒 isdeﬁnedas:

𝑘𝑚𝑒∶=

2𝑘_𝑚+_𝑒𝑘_𝑚−_𝑒 (

𝑘𝑚+𝑒 +𝑘𝑚−𝑒

) . (23)

Thearithmeticaverage,𝛿𝑒=0.5,issimplydenotedby⟨ · ⟩.Notethat, ingeneralcase,𝛿e isalsodeﬁnedforthekf;however,forthesakeof

simplicity,weassumethematerialpropertiesofthefracturedomainare homogeneous.Hence,weperformtheseoperatorsinthematrixdomain only.

Remark1. Thenumericaldiscretizationdiscussedin thisstudyonly considersthecaseof aconformingmesh,i.e.thefracturedomain,Γ, elementiscoincidentwithasetoffacesofthematrixdomain,Ω,as illustratedinFig.4b.

Thisstudyfocusesontwofunctionspaces,arisingfromEG,andDG discretizations,respectively.WebeginwithdeﬁningtheCGfunction

spaceforthematrixpressure,pmas CG_𝑘 ℎ ( ℎ)∶={𝜓𝑚∈ℂ0(Ω)∶ 𝜓𝑚||𝑇 ∈ℙ𝑘(𝑇),∀𝑇∈ℎ}, (24) whereCG_𝑘 ℎ (

ℎ)isthespacefortheCGapproximationwithkthdegree polynomialsforthe pm unknown,ℂ0(Ω)denotesthespaceof scalar-valuedpiecewisecontinuouspolynomials,ℙ𝑘(𝑇)isthespaceof polyno-mialsofdegreeatmostkovereachelementT,and𝜓mdenotesageneric functionofCG_𝑘

ℎ (

ℎ).Furthermore,wedeﬁnethefollowingDGfunction spaceforthematrixpressure,pm:

DG_𝑘 ℎ ( ℎ)∶={𝜓𝑚∈𝐿2(Ω)∶𝜓𝑚||𝑇∈ℙ𝑘(𝑇),∀𝑇∈ℎ}, (25) whereDG_𝑘 ℎ (

ℎ)isthespacefortheDGapproximationwithkthdegree polynomialsforthepmspaceandL2(Ω)isthespaceofsquareintegrable functions.Finally,wedeﬁnetheEGfunctionspaceforpmas:

EG_𝑘 ℎ ( ℎ)∶=ℎCG𝑘 ( ℎ)+ℎDG0 ( ℎ), (26) whereEG_𝑘 ℎ (

ℎ)isthespacefortheEGapproximationwithkthdegree polynomialsforthepspaceandDG0

ℎ (

(10)

approx-Fig. 11. Regularfracturenetworkexample:(a ) geometryandboundaryconditions(fractures areshowninred),pressuresolution,pm,inthe

matrixusing𝑛𝑒=2046of(b )permeable

frac-ture,(c )impermeablefracturecases,and(d ) meshthathas𝑛𝑒=184.(Forinterpretationof

thereferencestocolourinthisﬁgurelegend, thereaderisreferredtothewebversionofthis article.)

imationwith0thdegreepolynomials,inotherwords,apiecewise con-stantapproximation.NotethattheEGdiscretizationisexpectedtobe beneficialforanaccurateapproximationofthepmdiscontinuityacross interfaceswherehighpermeabilitycontrast,eitherbetween𝒌+ _𝒎and𝒌− _𝒎 ork_mandk_f,isobserved.Ontheotherhand,thematerialpropertiesof thefracturedomainareassumedtobehomogeneous,whichleadstono discontinuitywithinthefracturedomain.Therefore,theCG discretiza-tionofthep_funknownwillsufficeinthefollowing.Thep_ffunctionspace isdefinedas: CG_𝑘 ℎ ( Γ_ℎ)∶={𝜓_𝑓∈ℂ0(Γ)∶𝜓_𝑓|| |𝑒∈ℙ𝑘(𝑒),∀𝑒∈ Γℎ } , (27) whereCG_𝑘 ℎ (

Γ_ℎ)isthespacefor theCGapproximationwithkth de-greepolynomialsforthepfunknown,ℂ0(Γ)denotesthespaceof scalar-valuedpiecewisecontinuouspolynomials,ℙ𝑘(𝑒)isthespaceof polyno-mialsofdegreeatmostkovereachfacete,and𝜓fdenotesageneric functionof_ℎCG𝑘(Γ_ℎ).

Remark2. Inthisstudy,weonlyfocusonthemixedfunctionspace be-tweenthematrixandfracturedomainsarisingfromeitherEG_𝑘

ℎ ( ℎ)× CG_𝑘 ℎ ( Γ_ℎ) orDG_𝑘 ℎ (

ℎ)×_ℎCG𝑘(Γ_ℎ). The fracturedomain can be dis-cretizedbyeitherEG_𝑘 ℎ ( Γ_ℎ)orDG_𝑘 ℎ (

Γ_ℎ)ifthereareanydiscontinuities insidethefracturemedium.

Thetimedomain,𝕋=(0,𝜏],ispartitionedintoNtopensubintervals suchthat, 0=∶𝑡0<𝑡1<⋯<𝑡𝑁𝑡∶=𝜏. Thelength of the subinterval, Δtn_,_is _deﬁned_as_Δ_𝑡𝑛₌_𝑡𝑛₋_𝑡𝑛−1 _where_n _represents_the_current_time

step.Inthisstudy,implicitﬁrst-ordertimediscretizationisutilizedfor atimedomainof(1)and(5)asshownbelowforboth𝑝𝑛

𝑚,ℎand𝑝𝑛𝑓,ℎ: 𝜕𝑝𝑚,ℎ 𝜕𝑡 ≈ 𝑝𝑛 𝑚,ℎ−𝑝𝑛𝑚,ℎ−1 Δ𝑡𝑛 , and 𝜕𝑝𝑓,ℎ 𝜕𝑡 ≈ 𝑝𝑛 𝑓,ℎ−𝑝𝑛𝑓,ℎ−1 Δ𝑡𝑛 . (28)

WedenotethatthetemporalapproximationofthefunctionΦ(· ,tn₎_by Φn.

Withgiven𝑝𝑛−1

𝑚,ℎ and𝑝𝑛𝑓,ℎ−1,wenowseektheapproximatedsolutions

𝑝𝑛 𝑚,ℎ∈ EG_𝑘 ℎ ( ℎ)and𝑝𝑛𝑓,ℎ∈ CG_𝑘 ℎ ( Γ_ℎ)ofp_m(· ,tn₎_and_p f(· ,tn), respec-tively,satisfying ((𝑝𝑛 𝑚,ℎ,𝑝𝑛𝑓,ℎ;𝑝𝑛𝑚,ℎ−1,𝑝𝑛𝑓,ℎ−1 ) ,(𝜓𝑚,𝜓𝑓)) + (( 𝑝𝑛 𝑚,ℎ,𝑝𝑛𝑓,ℎ ) ,(𝜓𝑚,𝜓𝑓))−(𝜓_𝑚,𝜓_𝑓)=0, ∀𝜓_𝑚∈_ℎEG𝑘(_ℎ)and∀𝜓_𝑓∈_ℎCG𝑘(Γ_ℎ). (29) First,thetemporaldiscretizationpartisdeﬁnedas

((𝑝𝑛 𝑚,ℎ,𝑝𝑛𝑓,ℎ;𝑝𝑛𝑚,ℎ−1,𝑝𝑛𝑓,ℎ−1 ) ,(𝜓𝑚,𝜓𝑓))∶=𝑚𝑚 ( 𝑝𝑛 𝑚,ℎ;𝑝𝑛𝑚,ℎ−1,𝜓𝑚 ) +𝑚_𝑓(𝑝𝑛_𝑓,ℎ;𝑝𝑛−1 𝑓,ℎ,𝜓𝑓 ) , (30) where 𝑚𝑚 ( 𝑝𝑛 𝑚,ℎ;𝑝𝑛𝑚,ℎ−1,𝜓𝑚 ) ∶=∑_{𝑇 ∈} ℎ∫𝑇𝑐𝜙𝑚 𝑝𝑛 𝑚,ℎ−𝑝𝑛𝑚,ℎ−1 Δ𝑡𝑛 𝜓𝑚𝑑𝑉, (31) and 𝑚𝑓 ( 𝑝𝑛 𝑓,ℎ;𝑝𝑛𝑓,ℎ−1,𝜓𝑓 ) ∶=∑_𝑒∈Γ_ℎ∫_𝑒𝑐𝑎𝑓 𝑝𝑛 𝑓,ℎ−𝑝𝑛𝑓,ℎ−1 Δ𝑡𝑛 𝜓𝑓𝑑𝑆. (32)

(11)

Fig. 12. pmplotofregularfracturenetworkexample:(a )along𝑦=0.7lineofpermeablefracture,(b )along𝑥=0.5lineofpermeablefracture,(c )along(𝑥=0.0,𝑦=

0.1)to(𝑥=0.9,𝑦=1.0)lineofimpermeablefracturecases.AlloftheresultsobtainedfromtheEGandDGmethodsofpermeablefractureswithdifferentmeshsizes overlap.Theresultsofimpermeablefracturescasewithdifferentmeshsizes,however,illustratesomedifferences.

Here,∫T· dVand∫e· dSrefertovolumeandsurfaceintegrals, respec-tively. Next,wedeﬁne((𝑝𝑛 𝑚,ℎ,𝑝𝑛𝑓,ℎ ) ,(𝜓𝑚,𝜓𝑓))as: ((𝑝𝑛 𝑚,ℎ,𝑝𝑛𝑓,ℎ ) ,(𝜓𝑚,𝜓𝑓))∶=𝑎 ( 𝑝𝑛 𝑚,ℎ,𝜓𝑚 ) +𝑏 ( 𝑝𝑛 𝑓,ℎ,𝜓𝑚 ) +𝑐(𝑝𝑛_𝑚,ℎ,𝜓_𝑓)+𝑑(𝑝𝑛_𝑓,ℎ,𝜓_𝑓), (33) where 𝑎(𝑝𝑛 𝑚,ℎ,𝜓𝑚 ) ∶= ∑ 𝑇 ∈ℎ∫𝑇 𝒌𝒎 𝜇 ( ∇𝑝𝑛_𝑚,ℎ−𝜌𝐠 ) ⋅ ∇𝜓𝑚𝑑𝑉 − ∑ 𝑒∈1 ℎ⧵Γℎ∫𝑒 ⟨𝒌_𝜇𝒎(∇𝑝𝑛_𝑚,ℎ−𝜌𝐠)⟩_𝛿 𝑒⋅ 𝜓𝑚𝑑𝑆 +𝜃 ∑ 𝑒∈1 ℎ⧵Γℎ∫𝑒 ⟨𝒌𝒎 𝜇 ∇𝜓𝑚⟩𝛿𝑒⋅ 𝑝𝑛𝑚,ℎ𝑑𝑆 + ∑ 𝑒∈1 ℎ⧵Γℎ∫𝑒 𝛽 ℎ𝑒 𝑘𝑚𝑒 𝜇 𝑝𝑛𝑚,ℎ⋅ 𝜓𝑚𝑑𝑆 + ∑ 𝑒∈Γ_ℎ∫𝑒 𝛼𝑓 2 𝑝 𝑛 𝑚,ℎ⋅ 𝜓𝑚𝑑𝑆 + ∑ 𝑒∈Γ_ℎ∫𝑒 2𝛼_𝑓 2𝜉 −1⟨𝑝 𝑛 𝑚,ℎ⟩⟨𝜓𝑚⟩ 𝑑𝑆, (34) 𝑏(𝑝𝑛 𝑓,ℎ,𝜓𝑚 ) ∶=−∑ 𝑒∈Γ_ℎ∫𝑒 2𝛼_𝑓 2𝜉 −1𝑝 𝑛 𝑓,ℎ⟨𝜓𝑚⟩ 𝑑𝑆, (35) 𝑐(𝑝𝑛 𝑚,ℎ,𝜓𝑓 ) ∶=−∑ 𝑒∈Γ_ℎ∫𝑒 2𝛼_𝑓 2𝜉 −1⟨𝑝 𝑛 𝑚,ℎ⟩𝜓𝑓𝑑𝑆, (36)

(12)

Fig. 13. ﬂuxatsurfaceA’comparisonforbothpermeableandimpermeable casesbetweentheEGandDGmethods.AlloftheresultsobtainedfromtheEG andDGmethodswithdiﬀerentmeshsizesoverlap.

and 𝑑(𝑝𝑛 𝑓,ℎ,𝜓𝑓 ) ∶= ∑ 𝑒∈Γ_ℎ∫𝑒 𝑎𝑓 𝒌𝑻 𝒇 𝜇 ( ∇𝑝𝑛_𝑓,ℎ−𝜌𝐠 ) ⋅ ∇𝜓𝑓𝑑𝑆 +∑ 𝑒∈Γ_ℎ∫𝑒 2𝛼_𝑓 2𝜉 −1𝑝 𝑛 𝑓,ℎ𝜓𝑓𝑑𝑆, (37)

Wenotethatthecouplingconditions,(15)and(16),areembeddedin theabovediscretizedequations.Inparticular,theconditions(15)and

(16)arediscretizedas 1 ( 𝑝𝑛 𝑚,ℎ,𝜓𝑚 ) ∶= ∑ 𝑒∈Γ_ℎ∫𝑒 𝛼𝑓 2𝑝 𝑛 𝑚,ℎ⋅ 𝜓𝑚𝑑𝑆, ∀𝜓𝑚∈ℎEG𝑘 ( ℎ), (38) and 2 (( 𝑝𝑛 𝑚,ℎ,𝑝𝑛𝑓,ℎ ) ,(𝜓𝑚,𝜓𝑓))∶= ∑ 𝑒∈Γ_ℎ∫𝑒 2𝛼_𝑓 2𝜉 −1⟨𝑝 𝑛 𝑚,ℎ⟩⟨𝜓𝑚⟩ 𝑑𝑆 − ∑ 𝑒∈Γ_ℎ∫𝑒 2𝛼_𝑓 2𝜉 −1⟨𝑝 𝑛 𝑚,ℎ⟩𝜓𝑓𝑑𝑆− ∑ 𝑒∈Γ_ℎ∫𝑒 2𝛼_𝑓 2𝜉 −1𝑝 𝑛 𝑓,ℎ⟨𝜓𝑚⟩ 𝑑𝑆 + ∑ 𝑒∈Γ_ℎ∫𝑒 2𝛼_𝑓 2𝜉 −1𝑝 𝑛 𝑓,ℎ𝜓𝑓𝑑𝑆, ∀𝜓𝑚∈ℎEG𝑘 ( ℎ)and∀𝜓𝑓∈ℎCG𝑘 ( Γ_ℎ), (39) respectively.Finally,wedeﬁne(𝜓𝑚,𝜓𝑓)as:

(𝜓𝑚,𝜓𝑓)∶=𝓁𝑚(𝜓𝑚)+𝓁𝑓(𝜓𝑓) (40) where 𝓁𝑚(𝜓𝑚)∶= ∑ 𝑇 ∈ℎ∫𝑇 𝑔𝑚𝜓𝑚𝑑𝑉+ ∑ 𝑒∈𝑁 ℎ ∫𝑒𝑞𝑚𝑁𝜓𝑚𝑑𝑆 +𝜃 ∑ 𝑒∈𝐷 ℎ ∫_𝑒𝒌𝜇𝒎∇𝜓𝑚⋅ 𝑝𝑚𝐷𝒏𝑑𝑆+ ∑ 𝑒∈𝐷 ℎ ∫_𝑒ℎ𝛽𝑒 𝑘𝑚𝑒 𝜇 𝜓𝑝⋅ 𝑝𝑚𝐷𝒏𝑑𝑆 (41) and 𝓁𝑓(𝜓𝑓)∶=∑𝑒∈Γ_ℎ∫𝑒𝑎𝑓𝑔𝑓𝜓𝑓𝑑𝑆+∑𝑒𝑓∈Λ_ℎ𝑁∫𝑒𝑞𝑓𝑁𝜓𝑓𝑑𝑆. (42) Here,thechoicesoftheinteriorpenaltymethodisprovidedby𝜃.The discretizationbecomesthesymmetricinteriorpenaltyGalerkinmethod (SIPG),when𝜃 =−1,theincompleteinteriorpenaltyGalerkinmethod

(IIPG), when𝜃 =0, andthenon-symmetricinteriorpenalty Galerkin method(NIPG)when𝜃 =1Riviere(2008).Inthisstudy,weset𝜃 =−1 forthesimplicity.Theinteriorpenaltyparameter, 𝛽,isafunctionof thepolynomialdegree,k,andthecharacteristicmeshsize,he,whichis deﬁnedas:

ℎ𝑒∶=

meas(𝑇+)₊_meas₍_𝑇−₎

2meas(𝑒) , (43)

wheremeas(· )representsameasurementoperator,measuringlength, area,orvolume.Somestudiesfortheoptimalchoiceof𝛽 isprovidedin

Leeetal.(2019),Riviere(2008).

Remark3. TheNeumannboundaryconditionisnaturallyappliedon theboundaryfacesthatbelongtotheNeumannboundarydomainfor boththematrixandfracturedomains,𝑒∈_ℎ𝑁and𝑒_𝑓∈ Λ𝑁_ℎ.The Dirich-letboundarycondition,ontheotherhand,isweaklyenforcedonthe Dirichletboundaryfaces,𝑒∈_ℎ𝐷,forthematrixdomain,butstrongly enforcedontheDirichletbondaryfacesofthefracturedomain,𝑒𝑓∈ Λ𝐷_ℎ.

Let{𝜓_𝑚,𝜋(𝑖1)_}𝑝𝑚,𝜋

𝑖1=1 denotethesetofbasisfunctionsof

𝜋𝑘 ℎ ( ℎ),i.e. 𝜋𝑘 ℎ ( ℎ)=span{𝜓_𝑚(𝑖1)_}𝑝𝑚,𝜋

𝑖1=1 ,havingdenotedby𝑝𝑚,𝜋 thenumberof

DOFforthe𝜋 scalar-valuedspace,where𝜋 canmeaneitherEGorDG.In asimilarway,let{𝜓(𝑖3)

𝑓,CG} 𝑝𝑓,CG

𝑖3=1 bethesetofbasisfunctionsforthespace

CG_𝑘

ℎ (

Γ_ℎ)._𝑝_𝑓_,CG thenumberofDOFfortheCGscalar-valuedspace.

Hence,twomixedfunctionspaces,(1)EG_k×CG_kand(2)DG_k×CG_k wherekrepresentsthedegreeofpolynomialapproximation,are possi-ble,andwillbecomparedinthenumericalexamplesinthenextsection. Thematrixcorrespondingtotheleft-handsideof(29)isassembled com-posingthefollowingblocks:

[ 𝜋𝑘×𝜋𝑘 𝑚𝑚 ] 𝑖𝑖𝑖2 ∶=𝑚_𝑚(𝜓_𝑚(𝑖2),𝜓 𝑚(𝑖1))₊𝑎(𝜓 𝑚(𝑖2),𝜓 𝑚(𝑖1)),𝑖 1 =1,…,_𝑝_𝑚_,𝜋,𝑖2=1,…,𝑝𝑚,𝜋, [ 𝜋𝑘×CG_𝑘 𝑚𝑓 ] 𝑖𝑖𝑖3 ∶=𝑏(𝜓(𝑖3) 𝑓,CG,𝜓𝑚,𝜋 (𝑖1) ) ,𝑖1=1,…,𝑝𝑚,𝜋,𝑖3=1,…,𝑝𝑓,CG, [ CG_𝑘×𝜋_𝑘 𝑓𝑚 ] 𝑖4𝑖2 ∶=𝑐(𝜓_𝑚,𝜋(𝑖2),𝜓(𝑖4) 𝑝𝑓,CG ) ,𝑖4=1,…,𝑝𝑓,CG,𝑖2=1,…,𝑝𝑚,𝜋, [ CG_𝑘×CG_𝑘 𝑓𝑓 ] 𝑖4𝑖3 ∶=𝑚_𝑓(𝜓(𝑖3) 𝑓,CG,𝜓 (𝑖4) 𝑓,CG ) +𝑑(𝜓(𝑖3) 𝑓,CG,𝜓 (𝑖4) 𝑓,CG ) 𝑖4 =1,…,_𝑝 𝑓,CG,𝑖3=1,…,𝑝𝑓,CG. (44)

Inasimilarway,theright-handsideof(29)givesrisetoablockvector ofcomponents [ 𝓁𝜋𝑘 𝑚]𝑖1∶=𝓁𝑚 ( 𝜓𝑚,𝜋(𝑖1)), 𝑖 1=1,…,𝑝𝑚,𝜋, [ 𝓁CG_𝑘 𝑓 ] 𝑖3 ∶=𝓁_𝑓(𝜓(𝑖3) 𝑓,CG ) 𝑖3=1,…,𝑓,CG. (45)

Theresultingblockstructureisthus [ 𝜋𝑘×𝜋_𝑘 𝑚𝑚 𝑚𝑓𝜋𝑘×CG𝑘 CG_𝑘×𝜋𝑘 𝑓𝑚  CG_𝑘×CG_𝑘 𝑓𝑓 ]{ 𝑝𝜋𝑘 𝑚,ℎ 𝑝CG_𝑘 𝑓,ℎ } = { 𝓁𝜋𝑘 𝑚 𝓁CG_𝑘 𝑓 } . (46) where𝑝𝜋𝑘 𝑚,ℎand𝑝 CG_𝑘

𝑓,ℎ collectthedegreesoffreedomformatrixand frac-turepressure,respectively.Finally,weremarkthat(owingto(26))the case𝜋 =EGcanbeequivalentlydecomposedintoa(CG_k×DG₀) ×CG_k mixedfunctionspace,resultingin:

⎡ ⎢ ⎢ ⎢ ⎣ CG_𝑘×CG_𝑘 𝑚𝑚 𝑚𝑚CG𝑘×DG0 𝑚𝑓CG𝑘×CG𝑘 DG₀×CG_𝑘 𝑚𝑚 𝑚𝑚DG0×DG0 𝑚𝑓DG0×CG𝑘 CG_𝑘×CG_𝑘 𝑓𝑚  CG_𝑘×DG0 𝑓𝑚  CG_𝑘×CG_𝑘 𝑓𝑓 ⎤ ⎥ ⎥ ⎥ ⎦ ⎧ ⎪ ⎨ ⎪ ⎩ 𝑝CG_𝑘 𝑚,ℎ 𝑝DG₀ 𝑚,ℎ 𝑝CG_𝑘 𝑓,ℎ ⎫ ⎪ ⎬ ⎪ ⎭ = ⎧ ⎪ ⎨ ⎪ ⎩ 𝓁CG_𝑘 𝑚 𝓁DG₀ 𝑚 𝓁CG_𝑘 𝑓 ⎫ ⎪ ⎬ ⎪ ⎭ . (47)

ThisformulationmakestheEGmethodologyeasilyimplementable inanyexistingDGcodes.MatricesandvectorsarebuiltbyFEniCSform compiler(Alnaesetal.,2015).Theblockstructureissetupusing multi-phenicstoolbox(BallarinandRozza,2019).Randomﬁeldof permeabil-ity(km)ispopulatedusingSciPypackage(Jonesetal.,2001).𝛽,penalty

(13)

Fig. 14. Lowkmcase:kmdistributionof:(a )quarterﬁve-spot(𝑛𝑒=6,568),(b )regularfracturenetwork(𝑛𝑒=26,952)examples,histogramofkm of(c )quarter

ﬁve-spot,and(d )regularfracturenetworkexamples.

Remark4. WenotethatCG,DG,andEGmethodsarebasedonGalerkin method,whichcouldbeextendedtoconsideradaptivemeshesthat con-tainhangingnodes.Inaddition,there arevariousadvanced develop-mentforeachmethodstoenhancetheeﬃciency,includingvariable ap-proximationordertechniques.Especially,forEGmethod,anadaptive enrichment,i.e.,thepiecewise-constantfunctionsonlyaddedtothe el-ementswherethesharpmaterialdiscontinuities(e.g.,betweenmatrix andfracturedomains)areobserved,canbedeveloped.However,inour followingnumericalexamples,wefocusontheclassicalformofeach methodsforthecomparisonbysimulatingtheproposedmixed dimen-sionalapproachformodelingfractures.

3. Numericalexamples

WeillustratethecapabilityoftheEGmethodusingsevennumerical examples.Webeginwithananalysisoftheerrorconvergencerate be-tweentheEGandDGmethodstoverifythedevelopedblockstructurein themixed-dimensionalsetting.WealsoinvestigatetheEGperformance inmodelingthequarterﬁve-spotpatternandhandlingthefracturetip intheimmersedfracturegeometry.Next,wetesttheEGmethodina regularfracturenetworkwithandwithoutaheterogeneityinmatrix permeabilityinput.Lastly,weapplytime-dependentproblemsfortwo three-dimensionalgeometries;theﬁrstonerepresentsthecasewhere

fracturesareorthogonaltotheaxes,andanotherrepresentsgeometry wherefracturesaregivenwitharbitraryorientationswiththeir interac-tionsinathree-dimensionaldomain.

3.1. Errorconvergenceanalysis

Toverifytheimplementationoftheproposedblockstructure uti-lizedtosolvethemixed-dimensionalmodelusingtheEGmethod,we illustratetheerrorconvergencerateoftheEGmethodandcomparethis valuewiththeDGmethod.Theexampleusedinthisanalysisisadapted fromAntoniettietal.(2019).WetakeΩ =[0,1]2,andchoosetheexact solutioninthematrix,Ω,andfracture,Γ ={(𝑥,𝑦)∈ Ω ∶𝑥+𝑦=1},as: ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ 𝑝𝑚=𝑒𝑥+𝑦 inΩ1, 𝑝𝑚= 𝑒 𝑥+𝑦 2 + ⎛ ⎜ ⎜ ⎝ 1 2+ 3𝑎_𝑓 𝑘𝑛 𝑓 √ 2 ⎞ ⎟ ⎟ ⎠ 𝑒 inΩ2, 𝑝𝑓=𝑒 ( 1+√2𝑎𝑓 𝑘𝑛 𝑓 ) inΓ. (48)

Bychoosing𝒌_𝒎=𝑰,(48)satisﬁesthesystemofequations,(1),(9),(15), and(16),presentedinthemethodologysectionwithsink/sourceterms,

(14)

Fig. 15. Highkmcase:kmdistributionof:(a )quarterﬁve-spot(𝑛𝑒=6,568),(b )regularfracturenetwork(𝑛𝑒=26,952)examples,histogramofkmof(c )quarter

ﬁve-spot,and(d )regularfracturenetworkexamples.

g,asfollows: ⎧ ⎪ ⎨ ⎪ ⎩ 𝑔𝑚=−2𝑒𝑥+𝑦 inΩ1, 𝑔𝑚=−𝑒𝑥+𝑦 inΩ2, 𝑔𝑓 =√𝑒 2 inΓ, (49)

All other physical parametersare setto one, andthe homogeneous boundaryconditionsareappliedtoallboundaries.Thegeometryused inthisanalysisandtheillustrationoftheexactsolutionarepresented inFig.5a-b,respectively.

WecalculateL2_norm_of_the_diﬀerence_between_the_exact_solution,_p_,

andapproximatedsolution,ph,andtheresultsarepresentedinFig.5 c-dformatrixandfracturedomains,respectively.Forbothmatrixand fracturedomains,theEGandDGmethodsprovidetheexpected conver-gencerateoftwoandthreeforpolynomialdegreeapproximation,k,of oneandtwo,respectively(Antoniettietal.,2019;Babuska,1973).

3.2. Quarterﬁve-spotexample

ThisnumericalexampleteststheEGmethodperformanceinan in-jection/productionsetting usingﬁve-spotpattern,andweadoptthis examplefromChaveetal.(2018),Antonietti etal.(2019). The ﬁve-spot pattern,whereone injection well is located in the middleand

fourproducersarelocatedateachcornerofthesquare,iscommonly usedinundergroundenergyextractionChenetal.(2006).Duetothe symmetryofthisgeometry,onlyaquarterofthedomain(Ω =[0,1]2) is simulated. Theinjection well is located at (0,0), andthe produc-tion well islocated at (1,1).We placethefracturewith 𝑎𝑓=0.0005 atΓ ={(𝑥,𝑦)∈ Ω ∶𝑥+𝑦=1}.Thegeometry,boundaryconditions,and meshwithℎ=1.2× 10−1_applied_in_this_analysis_are_shown_in_Fig.₆_.

Thefollowingsourcetermisapplied totheentirematrixdomain includingtheinjectionandproductionwells:

𝑔𝑚(𝑥,𝑦)=10.1tan ( 200 ( 0.2−√𝑥2₊_𝑦2)) −10.1tanh ( 200 ( 0.2−√(𝑥−1)2₊₍_𝑦₋₁₎2))_. ₍₅₀₎

Toinvestigatetheeﬀectoffractureconductivity,weperformtwo sim-ulations using diﬀerent fracture conductivity inputs. (i) We choose, 𝒌𝒎=𝑰,𝑘𝑛𝑓=1,and𝒌𝑻𝒇 =100𝑰forthepermeablefracturecase.(ii)We

assumethefractureisimpermeableandset𝒌𝒎=𝑰,𝑘𝑛_𝑓=1× 10−2_,_and

𝒌𝑻

𝒇 =𝑰.Alloftheremainingphysicalparametersaresettoone.

ResultsoftwocasesarepresentedinFig.7a-bforthepressurevalue andFig.7c-dforthepressureprofilealong𝑥=𝑦line.Thepressure pro-fileofthepermeablefracturecaseissmoothlydecreasingfromthe in-jectionwelltowardstheproductionwell.Thepressureprofileofthe

(15)

im-Fig. 16. Thepmsolutionofheterogeneousquarterﬁve-spotexample(lowkm):(a )permeablefracture,(b )impermeablefracturecases,pmplotalong𝑥=𝑦lineof

(c )permeablefracture,and(d )impermeablefracturecases.AlloftheresultsobtainedfromtheEGandDGmethodsapproximatelyoverlap.

permeablefracturecase,ontheotherhand,illustratesajumpacrossthe fractureinterface.Theseﬁndingscomplywiththeresultsofthe previ-ousstudies(Chaveetal.,2018;Antoniettietal.,2019).Ourresults con-vergetothereferencesolutionasthehisreducedfromℎ=1.2× 10−1to

ℎ=7.2× 10−2_._Note_that₍_Antonietti_et_al.,₂₀₁₉₎_perform_these

numer-icalexperimentsusingthesecond-orderDGmethodwithℎ=7.5× 10−2

onpolytopicgrids(Antoniettietal.,2019).Thereisnosigniﬁcant dif-ferencebetweentheEGandDGresultsforbothhvalues(Fig.7c-d).

3.3. Immersedfractureexample

The numerical examples discussedso far contain a fracture that cut through the matrix domain. To testthe EG discretization capa-bilityin the immersedfracturesetting, weadopt this example from

Angotetal.(2009).Sinceweassumethatthefracturetipis substan-tiallysmall,thereisnoﬂuidmasstransferbetweenthematrixand frac-turedomainsacrossthefracturetip,seeA’inFig.8.Hence,thefracture boundary,Λ𝑁_ℎ,thatintersectswith thebulkmatrixmaterialinternal boundary,0

ℎ,isenforcedwithno-ﬂowboundarycondition,𝑞𝑓𝑁=0. Inthisexample,wetakethebulkmatrix,Ω =[0,1]2,andthe frac-turewith𝑎𝑓=0.01,Γ ={(𝑥,𝑦)∈ Ω ∶𝑥=0.5,𝑦≥0.5}.Weperformthree

simulationsusingdiﬀerentfractureconductivityinputs;(i)weassume 𝒌𝒎=𝑰,𝑘𝑛𝑓=1× 102,and𝒌𝑻𝒇=1× 106𝑰forthepermeablefracturecase,

anditsgeometryandboundaryconditionsarepresentedinFig.8a.(ii) Thepartiallypermeablefracturecaseutilizes𝒌𝒎=𝑰,𝑘𝑛𝑓 =1× 102,and 𝒌𝑻

𝒇 =1× 102𝑰.Thiscase geometryandboundaryconditionsare

illus-trated inFig.8b.(iii)Impermeable fracture,geometryandboundary conditionsareshowninFig.8c,andituses𝒌𝒎=𝑰,𝑘𝑛𝑓=1× 10−7,and 𝒌𝑻

𝒇 =1× 10−7𝑰.Allotherphysicalparametersareequaltoone.Example

ofmeshthatcontains𝑛𝑒=272isshowninFig.8d.

InFig.9,thevaluesofpmwith𝑛𝑒=1,741forallcasesarepresented. Thepressureplotalong𝑦=0.75ispresentedinFig.10.Thepermeable andpartiallypermeable fracturecasesillustratethecontinuityof the pressurewhile theimpermeablefracturecase showsthejumpacross thefractureinterface.Ourresultsconverge(𝑛𝑒=272,𝑛𝑒=1,741,and

𝑛𝑒=6,698)tothesolutionsprovidedbyAngotetal.(2009).Moreover, theEGandDGmethodsprovidesimilarresults.Notethatthereference solutionsareperformedontheﬁnitevolumemethodwith65,536 con-trolvolumes.

(16)

Fig. 17. Thepmsolutionofheterogeneousquarterﬁve-spotexample(highkm):(a )permeablefracture,(b )impermeablefracturecases,pmplotalong𝑥=𝑦lineof

(c )permeablefracture,and(d )impermeablefracturecases.AlloftheresultsobtainedfromtheEGandDGmethodsapproximatelyoverlap. 3.4. Regularfracturenetworkexample

Weincreasethecomplexityoftheproblembyincreasingthenumber offracturesasshowninFlemischetal.(2018).Thisexample,however, iscalledregularfracturenetworksinceallfracturesareorthogonalto theaxes(xory).Thegeometryusedinthisexamplewasutilizedby

Geigeretal.(2013)foranalyzingmulti-ratedual-porositymodel.Details formodelgeometryandboundaryconditionsareshown inFig.11a. Weset𝒌𝒎=𝑰forallthematrixdomain,Ω =[0,1]2,and𝑎𝑓=1× 10−4 forallfractures,Γ. Twofractureconductivityinputsareused;(i)we choose𝑘𝑛_𝑓=1× 104 _and_𝒌𝑻

𝒇 =×104𝑰 forthepermeablefracturecase.

(ii)Fortheimpermeablefracturecase,weassume𝑘𝑛

𝑓=1× 10−4,and 𝒌𝑻

𝒇=1× 10−4𝑰.Alloftheremainingphysicalparametersaresettoone.

Exampleofmeshthatcontains𝑛𝑒=184isshowninFig.11d.

Thepmresultsofthepermeableandimpermeablefracturesare pre-sentedinFig.11b-c,respectively.Thepermeablefracturecaseshows thesmoothpmproﬁlewhiletheimpermeablefracturecaseclearly illus-tratesthejumpofpmacrossthefractureinterface.Theseresultsarethe sameasthereferencesolutionsprovidedbyFlemischetal.(2018).

Figs.12a-bpresentthepressureplotsalongthelines𝑦=0.7and𝑥= 0.5ofthepermeablefracturecase.Thepressureplotalong(𝑥=0.0,𝑦=

0.1)to(𝑥=0.9,𝑦=1.0)lineofimpermeablefracturecaseisshownin

Fig.12c.Ourresultsusing𝑛𝑒=184,𝑛𝑒=382,𝑛𝑒=2,046,and𝑛𝑒=26,952 convergetothereferencesolutionsprovidedbyFlemischetal.(2018). The reference solution is simulatedbased on ﬁnite volume method with equi-dimensional setting, and it contains 1,175,056 elements

Flemischetal.(2018).

Besidestheevaluationofpressureresults,wealsoinvestigatethe ﬂuxcalculatedatfaceA’(seeFig.11a)asfollows:

flux=−∑ 𝑒∈A′∫𝑒

𝒌𝑚

𝜇 (∇𝑝𝑚−𝜌𝐠)⋅ 𝒏𝑑𝑆onA′, (51)

AspresentedinFig.13,thedifferencebetweeneachnecaseis insignif-icant,i.e.thedifferentbetween𝑛𝑒=26,952and𝑛𝑒=184casesisless than1%.Furthermore,thereisnodifferencebetweentheEGandDG methods.TheDOFcomparisonbetweenEGandDGmethodsisshown in Table1. TheEGmethod requirestheDOF(inthematrixdomain) approximatelyhalfofthatoftheDGmethod.NotethattheDOFinthe fracturedomainisthesamebecausewediscretizethefracturedomain usingtheCGmethod.

(17)

Fig. 18. Thep_msolutionofheterogeneousregularfracturenetworkexample(lowk_m):(a )permeablefracture,(b )impermeablefracturecases,p_mplotalong(c )

𝑦=0.7and𝑥=0.5linesofpermeablefracturecase,and(d )(𝑥=0.0,𝑦=0.1)to(𝑥=0.9,𝑦=1.0)lineofimpermeablefracturecase. Table 1

Degreesoffreedom(DOF)comparisonbetweenEGandDGmethodsofregular fracturenetworkexample.

ne EG DG

matrix domain fracture domain matrix domain fracture domain

26,952 40,629 13,677 80,856 13,677

2,046 3,123 1,077 6,138 1,077

382 595 213 1,146 213

184 293 109 552 109

3.5. Theheterogeneousinbulkmatrixpermeabilityexample

Thenumericalexamplespresentedsofaronlyconsideran homoge-neousmatrixpermeabilityvalue.Inthissection,weexaminethe capa-bilityoftheEGmethodinhandlingthediscontinuitynotonlybetween thefractureandmatrixdomainsbutalsowithinthematrixdomain(e.g.

NickandMatthai,2011b)byemployingtheheterogeneouspermeability inthebulkmatrix.Weadaptthequarterﬁve-spotasinSection3.2and regularfracturenetworkasinSection3.4.Wechoosetheﬁnestmesh frombothexamplestotesttheEGmethodcapabilitycomparedtothat

oftheDGmethodinhandlingthesharpdiscontinuitybetweenthe max-imumnumberofinterfaces.Thegeometries,boundaryconditions,and inputparametersareutilizedasSections3.2and3.4exceptforthek_m value.

Inthis study,𝒌𝒎=𝑘_𝑚𝑰,km valueisrandomlyprovided valuefor eachcell.Wewilldistinguishinparticulartwodiﬀerentcases,named lowk_mcaseandhighk_mcaseinthefollowing.Thelowk_mcaseis char-acterizedbyalog-normaldistributionwithaverage𝑘̄𝑚=1.0,variance var(𝑘_𝑚)=40,limitedtominimum value𝑘_𝑚min=1.0× 10−2,and

maxi-mumvalue𝑘𝑚max=1.0× 102.Thehighkmcaseuses𝑘̄𝑚=30.0,var(𝑘𝑚)= 90,𝑘_𝑚min=1.0× 10−1,and𝑘𝑚max=1.5× 102.Thisheterogeneousﬁelds

forbothexamplesarepopulatedusingSciPypackage(Jonesetal.,2001) asshowninFigs.14and15forlowk_mandhighk_mcases,respectively.

3.5.1. Quarterﬁve-spotexample

Theresultsforthelowkm casewithpermeable(𝑘𝑛_𝑓=1and𝒌𝑻_𝒇=

100𝑰)andimpermeable(𝑘𝑛

𝑓=1× 10−2and𝒌𝑻𝒇=𝑰)fracturesare

illus-tratedinFig.16.Ingeneral,thepmproﬁlefromthetwocasesaremore dispersethanthehomogeneouskmsetting.TheEGandDGmethods

(18)

Fig. 19. Thepmsolutionofheterogeneousregularfracturenetworkexample(highkm):(a )permeablefracture,(b )impermeablefracturecases,pmplotalong(c ) 𝑦=0.7and𝑥=0.5linesofpermeablefracturecase,and(d )(𝑥=0.0,𝑦=0.1)to(𝑥=0.9,𝑦=1.0)lineofimpermeablefracturecase.

(19)

Fig. 21. Thepmsolutionofheterogeneousandanisotropykminregularfracturenetworkexample:(a )permeablefracture,(b )impermeablefracturecases,thepm

plotalong(c )𝑦=0.7and𝑥=0.5linesofpermeablefracturecase,and(d )(𝑥=0.0,𝑦=0.1)to(𝑥=0.9,𝑦=1.0)lineofimpermeablefracturecase.

regardingpermeableandimpermablefracturesettingsareprovidedas follows:

1. Lowkmwithpermeablefractureresultillustratestheapproximately

smoothpmsolutionbecausethe𝑘̄𝑚=1.0isequalto𝑘𝑛𝑓.Theplot along𝑥=𝑦line,asexpected,showspmgraduallydecreasesfrom theinjectionwelltotheproductionwell.Thisresultcomplies withthatofthehomogeneouskmsetting.

2. Lowkm withimpermeablefractureresultandtheplotalong𝑥= 𝑦lineexhibitalittlejumpofp_m across thefractureinterface. Thisbehaviorisdiﬀerentfromthehomogeneouskmsettingsince 𝑘𝑚min=𝑘𝑛_𝑓,whichleadtothelesspermeabilitycontrastbetween

thefractureandmatrixdomains.

Theresultsforthehighk_m casewithpermeableandimpermeable fracturesarepresentedinFig.17.Using𝑛𝑒=6,568,theEGandDG re-sultsareapproximatelythesame.Theobservationsconcerningthe frac-turepermeabilityarepresentedbelow:

1. Highkmwithpermeablefracturedisplaysajumpacrossthefracture

domainbecausek_m aroundthefractureontheplottinglineis higherthan𝑘𝑛

𝑓;asaresult,thefractureinterfaceactslikeaﬂow barrier.Theplotalong𝑥=𝑦linesupportsthisobservationaspm jumpsacrossthefractureinterface.

2. Highkmwithimpermeablefractureillustratesahugejumpacross

thefractureinterface,ascanbeobservedfromthepmplotalong

𝑥=𝑦line,sincekmismuchhigherthan𝑘𝑛_𝑓.Thepmvariationis lesspronouncedcomparedtothelowkm one,seeFigs.7b

(ho-mogeneous)and16b(heterogeneous),becausetheﬂuidﬂowis blockedbythefractureinterface(sharpmaterialdiscontinuity).

3.5.2. Regularfracturenetworkexample

The pressure results of the low km case between the permeable

(𝑘𝑛

𝑓 =1× 104 and 𝒌𝑻𝒇 =×104𝑰) and impermeable (𝑘𝑛𝑓 =1× 10−4 and 𝒌𝑻

𝒇 =1× 10−4𝑰)fracturesareillustratedinFig.18.Similartothe

ﬁve-spotexample,theEGandDGresultsaresimilarwith𝑛_𝑒=26,952,andpm resultsaremoredispersivethanthehomogeneousk_msettingasshownin

(20)

Fig. 22. Three-dimensionalregularfracturenetwork.Theillustratedsurfaceswithedges(bluein a andredin b andc )indicatethefractures.(a )geometryand initial/boundaryconditionsarepresented.pmisenforcedaszeroontwocornerpointsshowninredandno-ﬂowconditionisappliedtoallboundaries.Initial

conditions(ICs)forpmandpfaresetasone.(b )presentspressuresolutiononplaneC,pm,matrixvelocity,vm,andfracturevelocity,vf,withthepermeablefracture

(casei),(c )presentspmonplaneD,vm,andvfwithimpermeablefracturecases(caseii).Notethatthematrixpressureisonlyshowninthecrosssectionbetween

thefartwoedgesofthemodel.(Forinterpretationofthereferencestocolourinthisﬁgurelegend,thereaderisreferredtothewebversionofthisarticle.)

Fig.19a-b.Thediscussionsregardingpermeableandimpermable frac-turesettingsareprovidedasfollows:

1. Lowkm withpermeablefractureillustratesthefracturedominate

ﬂow regime,even though 𝑘̄𝑚=1.0, kmmin is setto 1.0× 10−2,

whichismuchlessthan𝑘𝑛

𝑓.Thissettingreducestheimpactof thematrixdomain.pmplotalong𝑦=0.7and𝑥=0.5lines sup-portthisobservationbyshowingthehighpressuregradientin thematrixdomain,butpmbecomesmuchlessvariedwhen en-teringthefracturedomain.

2. Lowk_mwithimpermeablefracturealsopresentsthefracture dom-inate ﬂow regime because 𝑘𝑛

𝑓=1× 10−4, which is less than

𝑘𝑚min=1.0× 10−2.Therefore,theﬂowinthematrixisblocked

bythefracturedomain.Theplotalong(𝑥=0.0,𝑦=0.1)to(𝑥= 0.9,𝑦=1.0) line illustrates jumps across the fracture domain, whichissupportingourobservation.

Theresultsforthehighkm casewithpermeableandimpermeable

fracturesarepresentedinFig.19.With𝑛𝑒=26,952,theEGandDG

re-sultsareapproximatelythesame.Theobservationsconcerningthe frac-turepermeabilityarepresentedbelow:

1. High km withpermeable fracture showsthatthematrixdomain

gainsmoremomentumcomparingtothelowk_m case.p_m plot along𝑦=0.7and𝑥=0.5linesillustratesalsoapproximately lin-earreductionalongthematrixdomain,whilepressureisalmost constantinthefracturedomain.

2. High km withimpermeable fractureclearlypresents thefracture

domaindominatetheﬂowbecause𝑘𝑛

𝑓ismuchlessthanthekm.

Hence,theﬂowisblockedbythefractures.Theplotalong(𝑥= 0.0,𝑦=0.1)to(𝑥=0.9,𝑦=1.0)linesupportthisobservationby illustratingmultiplejumpsacrossthefracturedomainandnopm variationinsideeachmatrixblock.

3.6. Theheterogeneousandanisotropyinbulkmatrixpermeabilityexample

Thissectionillustratesthecomparisonbetweenthepmsolutionsof EGandDGmethodsusingtheheterogeneousandanisotropick_m.In

(21)

con-Fig. 23. Comparisonof(a )theaveragepminthewholedomainandblocksAandB,and(b )thepmproﬁlesfort=25andt=75along(𝑥=0.0,𝑦=0.0,𝑧=0.0)to

(𝑥=1.0,𝑦=1.0,𝑧=1.0)line.

trastwiththepreviousexample,weutilizethecoarsestmesh(𝑛𝑒=184) fromthe regularfracturenetworkexample.The matrixpermeability heterogeneity,km,isgeneratedwiththesamespecificationasthelow k_mcaseintheprevioussection.However,inthisstudy,theoff-diagonal termsof (13)arenot zero,andthekm ofeach elementisdefinedas

follows: 𝒌𝒎∶= [ 𝑘𝑚 0.1𝑘𝑚 0.1𝑘_𝑚 𝑘_𝑚 ] , (52)

ThegeneratedheterogeneousﬁeldisshowninFig.20a-bforboth diag-onalandoﬀ diagonalterms.

The pressure results of the low km case between the permeable

(𝑘𝑛

𝑓=1× 104 and 𝒌𝑻𝒇=×104𝑰) and impermeable (𝑘𝑛𝑓=1× 10−4 and 𝒌𝑻

𝒇 =1× 10−4𝑰)fracturesarepresentedinFig.21.TheresultsofEGand

DGmethodsareapproximatelysimilar forbothfracturepermeability settings.These resultsillustratethattheEGmethodcapturesthe dis-continuitiesandallowusingacoarsemeshtomaintainaccuracy.

3.7. Thetime-dependentproblems

Fromthissection,weconsidertime-dependentproblemwherecand

𝜙miarenonzero(1)–(8).Besides,weextendthespatialdomainto three-dimensionalspacetofurtherillustratetheapplicabilityoftheproposed EGmethod.Here,theﬁrstexampleconsidersthegeometrycontaining onlytheorthogonalfracturestotheaxes(x,y,orz),andthesecond ex-ampleassumesthegeometrywitharbitraryorientatednaturalfractures. Notethatwe,here,presentonlytheresultsoftheEGmethod.Theresults oftheDGmethodarecomparabletothoseoftheEGmethod.

3.7.1. Three-dimensionalregularfracturenetworkexample.

Inthisexample,weconsiderathree-dimensionaldomain,whichis ananalogofthetwo-dimensionalcasepresentedinSection3.5.2.This domaincontainsasetofwell-interconnectedfracturesandmeshedwith tetrahedralandtriangularelements(forfractures)with𝑛𝑒=9,544.The fracturegeometryisbasedontheexampleinBerreetal.(2020).The detailsofgeometrywithinitialandboundaryconditionsillustratedin

Fig.22a.

Here,weconsidertwodiﬀerent scenarios:casei)permeable frac-turecasewith𝑘𝑛

𝑓=1× 106and𝒌𝑻𝒇 =1× 106𝑰;andcaseii)impermeable

fracturecasewith𝑘𝑛

𝑓=1× 10−12,and𝒌𝒇𝑻 =1× 10−12𝑰.Forbothcases,

weemployapermeableporousmediumbysetting𝒌𝒎=[1.0,1.0,0.1]𝑰,

andalloftheremainingphysicalparametersaresettooneforthe sim-plicity.Thetemporaldomainisgivenas𝕋=[0,100]whereanuniform timestepsizeΔ𝑡𝑛=1.

InFig.22b-c,thenumerical results ofthepm,vm,andvf forthe

permeable(casei)andimpermeablefracturecases(caseii)at𝑡=100 arepresented.Amerevisualexaminationoftheseresultsalreadyshows thatthefracturepermeabilitycontrolstheflowfield.Forthepermeable fracturecase,thevelocityatthecornerofblockBishigherthanthat ontheoppositecornerinblockAasthefracturesareclosertotheopen cornerpointinblockB.Fortheimpermeablefracturecasethevelocity atthecornerofblockBislowerthanthatontheoppositecornerin blockAsinceblockAislargerthanblockB,anditcansupportflowfor alongertime.

SimilarbehaviorsofthepressurevaluesareobservedinFig.23a, wheretheaveragepressurevaluesofthefulldomainandblockAand B areplottedfor𝕋 =[0,100]. Itis clearthat theaveragepressurein blockBdropsfasterthaninblockAfortheimpermeablecase.Moreover,

Fig.23billustratesthevalueofpmalong(𝑥=0.0,𝑦=0.0,𝑧=0.0)to(𝑥= 1.0,𝑦=1.0,𝑧=1.0)lineat𝑡=25and75.

3.7.2. Algrøynaoutcropexample

The ﬁnal example is a three-dimensional fracture network built basedonanoutcropmapinAlgrøyna,Norway(Fumagallietal.,2019). Themodelhasasizeof850 × 1400 × 600mandcontains52 inter-sectingfractures(themodelisdescribedindetailinBerreetal.,2020). Theﬁniteelementmeshisdiscretizedbytetrahedral(rockmatrix)with

𝑛𝑒=163,575andtriangularelements(fractures)with𝑛𝑒=329,080.See

Fig.24aformoredetails.AsshowninthisﬁgureDirichletboundary conditionsareappliedontwoedgesofthemodeltorepresentan in-jectionandaproductionwell.Allotherboundariesareconsideredas no-ﬂow.Therockmatrixandthefracturesareconsideredpermeable:

𝑘𝑛

𝑓=1× 102,m2𝒌𝒇𝑻 =1× 102𝑰 m2,and𝒌𝒎=[1.0,1.0,0.1]𝑰 m2.Allof

theremainingphysical parametersaresettoone,and𝕋=[0,1× 106_]

secusinganuniformΔtn_of₁_×₁₀5_sec.

Fig.24bandcshowthesimulationresultsincludingthepressure field,thepressureiso-surfaces,andvelocityvectorsintherockmatrix at𝑡=500,000sec.Thepressureprofilesalongalinebetweentwo op-positecornersofthemodelatdifferenttimesareplottedinFig.24d. ThisexampleillustratestheapplicabilityofthepresentedEGmethod

(22)

Fig. 24. Algrøynaoutcropexample:(a )geometryandinitial/boundaryconditions(no-ﬂowconditionisappliedtoallboundaries(exceptstheproductionand injectionwells).Initialconditions(ICs)forpmandpfaresetasone.), b )pressuresolution,pm,matrixvelocity(inblackarrow),vm,andfracturesmesh,(c )pmand

itsiso-surfacesshowningrey,and(d )pressureproﬁlesalongalinefromthetopofinjectionwelltothebottomofproductionwell.

foracomplexthree-dimensionalfracturenetworkwitharbitrary orien-tations.

4. Conclusion

ThisstudypresentstheEGdiscretizationforsolvingasingle-phase ﬂuidﬂowinthefracturedporousmediausingthemixed-dimensional approach.Ourproposedmethodhasbeentested againstseveral

pub-lishedbenchmarksandsubsequentlyassesseditsperformanceinthetest caseswiththeheterogeneousandanisotropicmatrixpermeability.Our resultsillustratethatthepressuresolutionsresultedfromtheEGandDG method,withthesamemeshsize,areapproximatelysimilar. Further-more,theEGmethodenjoysthesamebeneﬁtsastheDGmethod;for instance,preserveslocalandglobalconservationforﬂuxes,canhandle discontinuitywithinandbetweenthesubdomains,andhastheoptimal