Adaptive Dynamic Multilevel Simulation of Fractured Geothermal Reservoirs

Delft University of Technology

Adaptive Dynamic Multilevel Simulation of Fractured Geothermal Reservoirs

HosseiniMehr, Mousa; Vuik, Cornelis; Hajibeygi, Hadi

DOI

10.1016/j.jcpx.2020.100061

Publication date

2020

Document Version

Final published version

Published in

Journal of Computational Physics: X

Citation (APA)

HosseiniMehr, M., Vuik, C., & Hajibeygi, H. (2020). Adaptive Dynamic Multilevel Simulation of Fractured

Geothermal Reservoirs. Journal of Computational Physics: X, 7, 1-24. [100061].

https://doi.org/10.1016/j.jcpx.2020.100061

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.

Contents lists available atScienceDirect

Journal

of

Computational

Physics:

X

www.elsevier.com/locate/jcpx

Adaptive

dynamic

multilevel

simulation

of

fractured

geothermal

reservoirs

Mousa HosseiniMehr

a

,

b

,

∗

,

Cornelis Vuik

a

,

Hadi Hajibeygi

b

a_{FacultyofElectricalEngineering,MathematicsandComputerScience,DepartmentofAppliedMathematics,DelftUniversityofTechnology,} VanMourikBroekmanweg6,2628XEDelft,Netherlands

b_{FacultyofCivilEngineeringandGeosciences,DepartmentofGeoscienceandEngineering,DelftUniversityofTechnology,Stevinweg1,} 2628CV,Delft,Netherlands

a

r

t

i

c

l

e

i

n

f

o

a

b

s

t

r

a

c

t

Articlehistory:

Received17January2020

Receivedinrevisedform11May2020 Accepted12May2020

Availableonline22May2020 Keywords:

Flowinporousmedia Adaptivemeshrefinement Multiscalesimulation Fracturedporousmedia Mass-heattransfercoupling Geothermalenergy

Analgebraicdynamicmultilevel(ADM)methodforfully-coupledsimulationofflowand heattransportinheterogeneousfracturedgeothermalreservoirsispresented.Fracturesare modeledexplicitlyusingtheprojection-basedembeddeddiscretemethod(pEDFM),which accuratelyrepresentsfractureswithgeneric conductivityvalues,frombarriersto highly-conductivemanifolds.Afullyimplicitschemeisusedtoobtainthecoupleddiscretesystem includingmassandenergybalanceequationswithtwomainunknowns(i.e.,pressureand temperature) atfine-scale level. The ADM methodis then developed to map the fine-scale discretesystemtoadynamicmultilevelcoarsegrid, independentlyfor matrixand fractures.To obtaintheADMmap,multilevelmultiscale coarsegridsare constructedfor matrix as well as for each fractureat allcoarsening levels.On this hierarchicalnested grids, multilevel multiscale basisfunctions (forflow and heat)are solvedlocally atthe beginning of the simulation. They are used during the ADM simulation to allow for accurate multilevel systems in presence of parameter heterogeneity. The resolution of ADMsimulationsisdefineddynamicallybasedonthesolutiongradient(i.e.fronttracking technique) usingauser-definedthreshold. The ADMmappingoccurs algebraicallyusing theso-calledADMprolongationandrestriction operators,forallunknowns.Avariety of 2Dand3Dfractured testcaseswithhomogeneousandheterogeneouspermeabilitymaps are studied. It is shown that ADM is able to model the coupled mass-heat transport accurately by employing only a fraction of fine-scale grid cells. Therefore, it promises an efficient approach for simulation of large and real-field scale fractured geothermal reservoirs.Allsoftwaredevelopmentsofthispaper ispubliclyavailableathttps://gitlab. com/DARSim2simulator.

©2020TheAuthors.PublishedbyElsevierInc.ThisisanopenaccessarticleundertheCC BYlicense(http://creativecommons.org/licenses/by/4.0/).

1. Introduction

Geothermalenergyresourcesarebecomingmoreimportantinthetransitiontowardsrenewableenergy.Duetolow car-bonfootprints[1–4] andhighpotentialofsustainability,thedemandforgeothermalenergyhasbeenrisingandisexpected toincreasefurtherwithinthenextdecades.Accuratescalablemodelingoffluidandheattransportingeothermalreservoirs

*

Correspondingauthorat:FacultyofElectricalEngineering,MathematicsandComputerScience,DepartmentofAppliedMathematics,DelftUniversityof Technology,VanMourikBroekmanweg6,2628XEDelft,Netherlands.

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

https://doi.org/10.1016/j.jcpx.2020.100061

2590-0552/©2020TheAuthors.PublishedbyElsevierInc.ThisisanopenaccessarticleundertheCCBYlicense (http://creativecommons.org/licenses/by/4.0/).

isofparamountimportanceinordertofulfilbothscientificandsocietalexpectationsonsuccessfulfielddevelopmentplans [5,6].

While maintainingfield scalability,simulation modelsshould accurately capturehighcontrastsinboth massandheat transport propertiesinsidetheporousmedia.Moreover,alargenumberofreservoirscontaincomplexgeologicalnetworks offracturesandfaultswithabroadrangeofconductivities.Suchextremecontrastshavesignificantimpactonflowandheat patterns.Therefore,highfidelityrepresentationofthephysicalphenomenawithintheheterogeneousfracturedreservoirsis crucial[7].However,whileachievingthistask,thecomputationalsciencecommunityfacesavarietyofcomplexchallenges, asthemodelsdemandhigh-resolutiongridstobeimposedontheentirereservoir(inordersofkilometres)[8].Non-linear behavior ofthesystemduetostrongmass-heatcouplingresultsinpoorstabilityandconvergence.Incaseofmulti-phase flow(e.g.,high-enthalpysystems)theseissuesbecomemoresevere[9].Thegeo-mechanicalprocesses(i.e.,elasticandplastic deformation)[10–12],reactivetransport(e.g.,geo-chemicalinteractionbetweenthesubstances)[13–15] andcompositional changes are other notable challenges that the system might face. Expectedly, the presence of fractures and faults [16–

18] with wide heterogeneity contrasts compared to the rest of the formation significantly increases the computational complexity.Thesechallengesintroducehighdemandsfordevelopingadvancedsimulationmethodsthatareabletoprovide efficiency (i.e., applicable to real field-scale problems) while maintaining accuracy atthe desired level. Toaddress these challenges,variousadvancedsimulationschemeshavebeendeveloped.

Discrete Fracture Models (DFM)aim to represent fracturesas conforming lower dimensional domains by discretizing them atthematrixcellinterfaces[19–21]. Ontheother hand,theEmbeddedDiscreteFracture Model(EDFM)makes use ofnon-conforminggridstodiscretizefracturesaslowerdimensionalobjectsfullyindependentfromthematrixgrid struc-ture[22–26]. Conservativeflux transfertermsare introducedbetweenoverlappingmatrixandfracture elements.EDFMis applicable to provide acceptable solutions forhighly conductive fractures. Even forhigh-conductive fractures, it leads to very small, yet unphysical, flux leakage. More importantly, EDFMcannot accurately represent flow barriers (e.g., sealing fracturesandimpermeablefaults).Inordertoresolvethislimitation,projection-basedEDFM(pEDFM)wasintroduced[27] where consistent connectivitiesbetweenmatrix andfractureswere developed. Therefore,pEDFM can be appliedto frac-tured porousmediawithgenericrangeofconductivitycontrastbetweenfracture andthehosting rock.EvenwithpEDFM independentgrid,thesizeofthereservoiranddensityofthefracturesmakeanyfine-scalesimulationapproachimpractical. Multiscale Finite-Volume(MsFV) [28–33] andMultiscale Finite-Element(MsFE)[34–36] methodshavebeendevised to accomplishefficientsimulationsinpresenceofheterogeneity,bysolvingcoarse-scalesystems,whileavoidingupscaling,yet preservingfine-scaleheterogeneities[37,38].DynamicLocalGridRefinement(DLGR)[39–47] isanothermethod,which em-ploys fine-scaleresolutioninpartsofthedomainwhereneeded(i.e.,moving fronts).Algebraic dynamicmultilevel(ADM) methodhasbeendevelopedtoaddressthemultilevelmultiscalechallengesonequationswithelliptic,parabolicand hyper-bolic systemsofPDEs [48]. Byextending thegrid refinementstrategy andtakingtheadvantages ofmultilevel multiscale basis functions, the ADMmethod isable to algebraicallymap the fine-scale discrete systems(either fullyimplicit or se-quential)to adynamic multilevel systemforall unknowns. Oncethe systemissolvedatADMresolution, itisprolonged back tofine-scale resolution.Incaseoffullyimplicit(FIM)formulation, asthesystemissolvedatADMresolutionforall unknowns, reconstruction of conservative flux field is not required. This method has been further developed for multi-phaseflowinheterogeneousformationswithcompositionalandcapillaryeffects[49],andfracturedreservoirsusingEDFM scheme[52,50].Recentdevelopmentsalsoincludemultiscalefiniteelementmethodforflow-heat transfermodeling using EDFM approach[37,38].However, such methodhasnot yet beendevelopedforcoupled mass-heattransport infractured porousmediawithpEDFMmethod.

In this article,an Algebraic Dynamic Multilevel (ADM) method forsimulation of non-isothermal single-phase flow in fractured porousmediaispresented.Duetothecomplexnon-linearityofthesystem(strongmass-heatcoupling),twosets ofequations(massandenergyconservation)aresolvedusingtheFIMapproach.Intheenergybalanceequation,conduction terms aredescribed forboth fluidandrock.Projection-based EDFM(pEDFM)is employedin ordertoexplicitlyrepresent fracturesandtoprovideindependentgriddingofmatrixandfractures.Here,theapplicabilityofthepEDFMimplementation [27] hasbeenextendedtoa fullygeneric3D geometrywhereit allowsforincludingfractures(orflow barriers)withany orientation ontheCartesiancoordinatesystem.Thisiscrucialforpracticalapplications.Whilemaintainingthegeometrical flexibility ofEDFM,the matrix-matrixandfracture-matrixconnectivitiesareadjusted toaccountforprojection offracture plates on the matrix transmissibilities. This allows for consistent modeling of fractureswith all conductivity (from high to low) compared withthat of the matrix.The fine-scale discrete systemincluding both matrix andfracture network is obtainedfortwo mainunknowns, namely,pressureandtemperature. Wenote thatlocalthermalequilibriumisassumed, meaning that the temperatures influid and rockare considered to be equal.This assumption canbe recognized as safe for the majority (or perhaps all) of the geo-science applications [51] due to large contact area between the liquid and solidphases.Oncethefine-scalesystemisdiscretizedusingpEDFM,itismappedintoadynamic multilevelresolutiongrid (i.e.,ADMgrid).TheADMgridisconstructedbasedonahierarchicallynestedmultilevelmultiscalegridresolution,where locally computed multiscale basis functionsare employed at each coarsening level [34,28,30]. Thesebasis functions are computedonlyatthebeginning ofthe simulationasa pre-processingstep.Thisprovides highercomputationalefficiency, whilecapturingthefine-scaleheterogeneouspropertiesaccurately.

In order toconstruct themultilevel multiscaleoperators for thegeothermal fractured system, different typesofbasis functionsare computedforpressure andtemperatureunknowns. Forpressure,asaglobalvariable,multiscalebasis func-tions are computed. However, temperature resembles a local variable for which constant basis functions are used. This

distinctionincomputationofbasisfunctionsispreservedon ahierarchicalpatternforallthecoarseninglevelsacrossthe domain.Tocapturetheeffectofpermeableandimpermeablefractures, allbasisfunctionsarefullycoupled,meaning that extrabasisfunctionsarecomputedforfracturesaswell.Theconstructionofallbasisfunctionsaredonefullyalgebraically. Inthisprocess, themultilevel prolongationoperatoris assembledso thatthe partition ofunityholds [53,54]. The multi-levelfinite-volume restrictionoperatorensures massconservationinside thedomain.Ateach time-step,the ADMmap is obtainedbyusingafronttrackingtechniquetodetecttheregionwherethesolutiongradientishigherthanauser-defined threshold.Thesesub-domainsofsharpgradientareresolvedatfine-scalelevel.Moreover,regionsnearthesourcetermsare alwayskeptatfine-scaleresolutiontomaintainaccuracy.Therestofthedomainisresolvedathierarchicallynestedcoarser resolution.ThisADMmapisabletoadaptitselftothemovingfrontdynamicallyandefficiently,whilecapturingthephysics ofinterest inbothflow andheattransport withan accuracysetby the user-definedthreshold.The Newtonlinearization schemeisemployedtofindnonlinearsolutionsiteratively.OncetheADMsolutionisfound,itisprolongedbacktothefine scaleusingprolongationoperators.

Varioushomogeneousandheterogeneoustestcasesareconsideredandtheirnumericalresultsarepresentedto demon-strate the accuracy andapplicability of the devised method.Furthermore, its sensitivitytowards fractures withdifferent rangesofconductivityisstudied.

Thisarticleisorganizedasfollows.ThegoverningequationsandFIMapproacharedescribedinSection1.Thefine-scale discretization using pEDFM is covered in Section 2. The ADM method as the simulation strategy forfractured geother-malreservoirs is presentedinSection 3.The test casesandtheir numericalresultsare showninSection 4. Thepaperis concludedinSection5.

2. Fine-scaleequationsandsolutionstrategy 2.1. Governingequations

Inthissection,thetwomaingoverningequationsaredescribed. Here,thetemperaturesoffluidandrockare assumed tobeidentical.

2.1.1. Massbalance

Massbalanceequationfornon-isothermalsingle-phasefluidflowinporousmediawithnfracdiscreteembeddedfractures reads

∂

∂t



φ

m

ρ

f



− ∇ ·



ρ

f 1

μ

f Km

· ∇

pm



=

ρ

fqmw

+

nfrac



i=1

ρ

∗_f

Q

mfi

_,

_on

_

m

⊆ 

n

,

(1)

fortherockmatrix(m)and

∂

∂t



φ

fi

_ρ

f



− ∇ ·



ρ

f 1

μ

f Kfi

· ∇

_pfi



=

ρ

fqfiw

+

ρ

∗f

Q

fim

+

nfrac



j=1



ρ

∗_f

Q

fifj



j=i

,

on



fi

⊆ 

n−1

_∀

_i

_{∈ {}

₁

_{, ...,n}

frac

},

(2) forthelowerdimensionalfracture( fi).Here,themainunknownisthepressurep.Moreover,

φ

andK aretheporosityand permeabilityofeachmedium.K isatensortoaccountforagenericanisotropiccase.Superscriptsm, fi andw correspond tomatrix,fracturei andwell,respectively.Subscripts f andr denotefluidandrock.

ρ

and

μ

arethedensityandviscosity ofthefluidrespectively.Inaddition,qmw andqfiw _{arethesourceterms(i.e.,wells)onmatrixm andfracture fi.Moreover,}

Q

mfi _and

Q

fim _{arethefluxexchangebetweenmatrixm andoverlappingfracture f}

i correspondingtothegridcells where overlapoccurs.

Q

fifj _{isthefluxexchangefrom j-thfracture tothei-thfractureontheintersectingelements.Thismeans}

that thementionedflux exchangetermsare non-zeroonlywherematrix-fracture overlaporfracture-fracture intersection exists.Duetomassconservation,

˝

V

Q

mfi_dV

= −

˜

A_fi

Q

fim_{d A,and}

˜

A_fi

Q

fifj_{d A}

= −

˜

A_{f j}

Q

fjfi_{d A hold.}

ThePeacemanwellmodelisusedtoobtainthewellsourcetermsformatrix qmw

=

P I

· 

∗

_{· (}

_pw

₋

_pm

₎

V

(3) andfractures qfiw

=

P I

· 

∗

_{· (}

_pw

₋

_pfi

₎



A

.

(4)

Here, P I isthewellproductivityindexand



∗ istheeffectivemobility(



=

K

/

μ

)betweenthewellandthepenetrated gridcellinthemedium.



V and



A arethecontrolvolumeandcontrolareausedinthediscretesystemformatrixm and

fracture fi respectively. The flux exchange terms

Q

mfi,

Q

fim (matrix-fracture connectivities) and

Q

fifj (fracture-fracture connectivities)arewrittenas:

Q

mfi

=

_{C I}mfi

· 

∗

· (

_pfi

−

_pm

₎

₍₅₎

Q

mfi

=

_{C I}fim

· 

∗

· (

_pm

−

_pfi

₎

₍₆₎

Q

fifj

=

_{C I}fifj

· 

∗

· (

_pfj

−

_pfi

_),

₍₇₎

where C I denotesthe connectivity indexbetween each two non-neighboring elements. As an example, the connectivity indexbetweeni-thmatrixelementand j-thfractureelementiscalculatedasC Ii j

=

Ai j

di j,where Ai j isthearea fractionof

fracturecell j overlappingwithmatrixcelli and

d

i j

istheaveragedistancebetweenthesecells.

2.1.2. Energybalance

Assuminglocalequilibrium,energybalanceontheentiredomainreads

∂

∂t



(

ρ

U)m_{e f f}



− ∇ ·



ρ

fHmf 1

μ

f Km

· ∇

pm



− ∇ ·





m_{e f f}

· ∇

Tm



=

=

ρ

fHmfqmw

+

nfrac



i=1

ρ

∗_fH∗_f

Q

mfi

+

nfrac



i=1

R

mfi

_,

_on

_

m

⊆ 

n

,

(8) intherockmatrix(m)and

∂

∂t



(

ρ

U)fi e f f



− ∇ ·



ρ

fHffi 1

μ

f Kfi

· ∇

_pfi



− ∇ ·





fi e f f

· ∇

T fi



=

=

ρ

fHffiqfiw

+

ρ

∗fH∗f

Q

fim

+

nfrac



j=1



ρ

∗_fH∗_f

Q

fifj



j=i

+

R

fim

+

nfrac



j=1



R

fifj



j=i

,

on



fi

⊆ 

n−1

_∀

_i

_{∈ {}

₁

_{, ...,n}

frac

},

(9) inthelowerdimensionalfracture( fi).Here,besidepressurep,thesecondmainunknownisT asthetemperatureinboth fluidandsolidrock.

(

ρ

U

)

e f f istheeffectivepropertydefinedas

(

ρ

U)e f f

= φρ

fUf

+ (

1

− φ)ρ

rUr

,

(10)

where Uf andUr denotethespecific internalenergyinfluid androckrespectively.Distinctly,

(

ρ

U

)

m_{e f f}

= φ

m

ρ

fUf

+ (

1

−

φ

m

)

ρ

rUr and

(

ρ

U

)

e f ffi

= φ

fi

ρ

fUf

+ (

1

− φ

fi

)

ρ

rUr.Moreover, Hf isthespecificfluidenthalpy.Thethreementionedterms can be expressed asnon-linear functionsof pressure andtemperature.



e f f is theeffective thermal conductivityof the saturatedrockdefinedas



e f f

= φ 

f

+ (

1

− φ)

r

.

(11)

Here,



f and



rarethethermalconductivitiesinfluidandrock,respectively.Thesubscripts f andr indicatefluidand solidrock.Notethat,



m_{e f f}

= φ

m



f

+(

1

−φ

m

)



rand



_{e f f}fi

= φ

fi



f

+(

1

−φ

fi

)



r.Lastly,

R

mfi and

R

fimaretheconductive heat flux exchange between matrixm and the overlapping fracture fi.

R

fifj _{denotes the conductive heat flux exchange}

from j-th fracture tothe i-th fracture wherethe intersection occurs. Similar to massflux exchange, the conductive flux exchange termsarenon-zeroonlyfortheexistingmatrix-fractureoverlapsorfracture-fractureintersections.

˝

V

R

mfi_dV

=

−

˜

A_fi

R

fim_{d A,and}

˜

A_fi

R

fifj_{d A}

= −

˜

A_{f j}

Q

fjfi_{d A holdaswelltohonortheconservationofenergy.}

Toobtaintheconductiveheatfluxexchanges,i.e.,

R

mfi_,

R

fim_{(matrix-fractureconnectivities)and}

R

fifj _{(fracture-fracture}

connectivities),theembeddeddiscreteschemeisused,i.e.,

R

mfi

=

_{C I}mfi

· 

∗ e f f

· (

Tfi

−

Tm

)

(12)

R

mfi

=

_{C I}fim

· 

∗ e f f

· (

T m

₋

_Tfi

₎

₍₁₃₎

R

fifj

=

_{C I}fifj

· 

∗ e f f

· (

Tfj

−

Tfi

),

(14)

where



∗_{e f f} isobtainedasharmonically-averagedpropertybetweenthetwo non-neighboringelements. C I is identicalto theconnectivityindexusedinEq.(5).

Fig. 1. Theindependentmatrix3Dgridandthefractures2Dgridsareshown.Notethateachdomainhasitsowngridandthatanyfractureorientationcan beconsidered.

2.1.3. Correlations

Thefollowingcorrelationsareusedtocomputecertainfluidandrockproperties.

Fluid viscosity: Viscosityasafunctionoftemperatureiscalculatedviathefollowingcorrelation[55],

μ

f

(T)

=

2

.

414

×

10−5

×

10 247.8

T−140

_.

Fluid density: Thedensityoffluidisafunctionofpressureandtemperature,definedas[51]

ρ

f

(P

,

T)

=

ρ

f s

(T)

[1

+

cw

(T) (

P

−

Ps

)

]

,

wherePS

=

1bar.cw

(

T

)

and

ρ

f s

(

T

)

arecomputedfromthefollowingempiricalcorrelations[8,56] cw

(T

)

=



0

.

0839 T2

+

652

.

73 T

−

203714



×

10−12

ρ

f s

(T)

= −

0

.

0032 T2

+

1

.

7508 T

+

757

.

5

.

Fluid Entalphy: Fluidenthalpyisobtainedvia[51] Hf

(

P,T

)

=

uws

+

Cp f

(T

−

Ts

)

+

P

ρ

f

,

whereuws

=

420000J

/

kg.

2.2. Fine-scalediscretesystemusingpEDFM

The coupled system of non-linear equations described above, i.e., Eqs. (1)-(9) with two main unknowns (i.e., p, T )

isdiscretizedspatiallyusingtwo-point-flux-approximation(TPFA)finite-volumeschemeandtemporallyusingthebackward (implicit)Eulerscheme[57].Thediscretizationresultsinsetofstructuredgridsforathree-dimensional(3D)porousmedium with2D fractureplanesincludedthataregeneratedindependently.Avisualexampleofthegriddingstructureisshownin Fig.1.Theobtainednon-linearFIMsystemissolvedusingtheNewton-Raphsoniterativemethod[57].

2.2.1. Massandheatfluxes

Themassfluxexchangebetweeneachtwoneighboringcontrolvolumesi and j (insideonemedium)usingTPFAscheme canbewrittenas Fi j

=

ρ

∗_f

μ

∗Ti j

(p

i

−

pj

),

(15) where Ti j

=

Ai j di jK H

i j denotes thetransmissibilitybetweengridcells i and j. Ai j istheinterface areabetweenthetwogrid cells,di j isthedistancebetweentheircellcentersandKi jH istheharmonicaverageofthepermeabilitiessetattheinterface ofgrid cells i and j. The propertieswithsuperscript

∗

are obtainedusing theupwindscheme. Notethat density

ρ

and viscosity

μ

arefunctionsofpressureandtemperature.

Themassfluxbetweeneachtwo non-neighboringgridcells i and j (eithermatrix-fractureorfracture-fracture connec-tivities)isobtainedvia EDFMformulation. Thisflux betweenamatrix(denotedasm)celli andafracture(denoted as f )

Fig. 2. Figureonthe left:Anillustration ofa2Dfractureplateoverlappingamatrix cell.Inthisexample,the overlappingsectionformsanirregular hexagon.Figureontheright:Avisualizationoftwointersectingfractureplates(left)anddiscreteelementsofdifferentintersectingfractures(right),with theintersectionlinesshowninred.

F

i jmf

= −

F

f m i j

= −

ρ

∗_f

μ

∗T mf i j

(p

mi

−

p f j

),

(16)

where the transmissibility Tmf_{i j}

=

K_{i j}HC Ii j. K_{i j}H is theharmonically averaged permeabilitybetweenthe overlapping matrix andfractureelements,andC Ii j istheconnectivityindexbetweenthem.TheEDFMmodelsthematrix-fractureconnectivity indexasC Ii j

=

Amf_{i j}

di j whereA

mf

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

d

i j

istheaveragedistancebetweenthesecells[24].

Themassfluxexchangebetweentwointersectingfractureelementsi (onfracture f )and j (onfractureg)isgivenas

F

f g i j

= −

F

g f i j

= −

ρ

∗_f

μ

∗T f g i j

(p

f i

−

p g j

),

(17)

where T_{i j}f g isthetransmissibility betweenthetwo cellsobtainedvia alowerdimensionalconnectivityindexformulation. Mindthattheintersectionof2Dfractureplatesisalineandfor1Dfractureline-segmentsitisapoint.Fig.2(right)shows a visualization exampleofintersectionbetweentwo non-neighboring2D fracture elements.The intersectionformsa line segment Ii j (showninredcolor)withtheaveragedistancesfromtheintersectionsegment of

d

_{i I}f_{i j}

=

d

g_{j Ii j}.Thus, T_{i j}f g is obtainedusinganharmonic-averageformulation,i.e.,

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

,

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

i j representstheconnectivityindexbetweenthe2Dfractureelementi andthe1Dintersectionlinesegment

Ii j,asshowninFig.2.

Theconvectiveheatfluxexchangebetweentheneighboringcontrolvolumesi and j reads

˘

Fi j

=

ρ

∗_fH∗_f

μ

∗ Ti j

(p

i

−

pj

).

(19)

Here, H∗_f istheenthalpyoffluiddeterminedattheinterfacebetweengridcellsi and j.OnecanconcludeF

˘

i j

=

H∗_f Fi j. Similarly, the convective heat flux exchange betweennon-neighboring elements can be obtainedvia multiplication of their massfluxexchange(

F

i j)by theeffectiveenthalpy(H∗f)determinedattheintersectionoftwooverlappingelements. Namely,

F

˘

i j

=

H∗_f

F

i j.

Theconductiveheatfluxexchangebetweeneachtwoneighboringcellsi and j iswrittenas

Gi j

= T

cond,i j

(T

i

−

Tj

),

(20)

where,

T

cond,i j

=

Ai j

di j

(

e f f

)

H

i j is thetransmissibilitybetweengridcells i and j. Ai j,di j and

(

De f f

)

Hi j aretheinterface area, thedistancefromthecellcentersandtheharmonicaverageoftheeffectiveconductivitysetattheinterfacebetweengrid cells i and j respectively.

Theconductiveheatfluxexchangebetweennon-neighboringelementsi and j isobtainedas

G

mf i j

= −

G

f m i j

= −T

mf cond,i j

(T

m i

−

T f j

),

(21)

formatrix-fractureconnectivitiesand

G

f g i j

= −

G

f g i j

= −T

f g cond,i j

(T

f i

−

T g j

),

(22)

forfracture-fractureconnectivities

Here,

T

cond,i j

= (e f f

)

i jH C Ii j,where

(

e f f

)

i jH is theharmonic average ofthe effective conductivitybetween the over-lapping matrix and fracture elements. A similar procedure as in Eqs (18) will be used to obtain the fracture-fracture connectivityindices.

Ateachtime-stepthefine-scalediscretemassbalanceequationreads



φ

m

ρ

f



n+1 i

−



φ

m

ρ

f



n i

t

+

Nn



j=1 Fi j

+

nfrac



k=1



N_fk



j=1

F

mfk i j



=



ρ

fqmw



i

,

∀

i

∈ {

1

, ...,

Nm

}

(23)

forelementi inmatrixm and



φ

fh

_ρ

_f



n+1 i

−



φ

fh

_ρ

_f



n i

t

+

Nn



j=1 Fi j

+

Nm



j=1

F

fhm i j

+

nfrac



k=1



N_fk



j=1

F

fhfk i j



=



ρ

fqfhw



i

,

∀

i

∈ {

1

, ...,

Nfh

}

(24) forelementi infracture fh.NmandNfk denotethenumberofelementsinmatrixm andfracture fkrespectively.Nnisthe

numberofneighboringgridcells.

Similarly,thefine-scalediscreteformoftheenergybalanceequationateachtime-stepiswrittenas:



(

ρ

U

)

m_{e f f}



n+1 i

−



(

ρ

U)m_{e f f}



n i

t

+

Nn



j=1

˘

Fi j

+

nfrac



k=1



N_fk



j=1

˘

F

mfk i j



+

+

Nn



j=1

˘

Gi j

+

nfrac



k=1



N_fk



j=1

˘

G

mfk i j



=



ρ

fHmfqmw



i

,

∀

i

∈ {

1

, ...,

Nm

}

(25) forelementi inmatrixm and



(

ρ

U

)

fh e f f



n+1 i

−



(

ρ

U)fh e f f



n i

t

+

Nn



j=1

˘

Fi j

+

Nm



j=1

˘

F

fhm i j

+

nfrac



k=1



N_fk



j=1

˘

F

fhfk i j



+

+

Nn



j=1

˘

Gi j

+

Nm



j=1

˘

G

fhm i j

+

nfrac



k=1



N_fk



j=1

˘

G

fhfk i j



=



ρ

fH_ffhqfhw



i

,

∀

i

∈ {

1

, ...,

Nfh

}

(26) for element i in fracture fh. The residual form of equations (23)-(26) can be written as subtracting the left-hand side fromtheright-hand sideineach equation.Theresidualofthemassbalanceequationforeachelementi in matrixm (i.e.,

(

rm M B

)

n+1 i )reads

(r

m_{M B}

)

n_i+1

=



ρ

fqmw



i

−



φ

m

ρ

f



n+1 i

−



φ

m

ρ

f



n i

t

−

−

Nn



j=1 Fi j

−

Nm



j=1

F

fhm i j

−

nfrac



k=1



N_fk



j=1

F

fhfk i j



=

0

,

∀

i

∈ {

1

, ...,

Nm

}.

(27)

The residualsforthe other threeequations(i.e.,the massbalancein fractures

(

rfh

M B

)

n+1

i ,the energybalancein matrix

(

rm E B

)

n+1

i andtheenergybalanceinfractures

(

r fh

E B

)

n+1

i )isobtainedsimilartoEq.(27).Thefullresidualvector attime-step

n isobtainedbyaverticalconcatenationoftheresidualspermedium.Namely,

rn

= [(

rm_{M B}

)

n

, (r

f1 M B

)

n

, (r

f2 M B

)

n

, ..., (r

f_nfrac E B

)

n

, (r

mM B

)

n

, (r

f1 E B

)

n

, (r

f2 E B

)

n

, ..., (r

f_nfrac E B

)

n

]

T

.

(28)

As theresidualattime-step n

+

1 (rn+1) isa non-linearfunction ofpressureandtemperatureunknownsattime-step

n

+

1 (i.e.,pn+1 _andTn+1_{),theNewton-Raphsoniterationmethodisappliedtosolvethesystemateachtime-step.Inother} words,

rν+1

=

rν

+

∂r

∂

p





ν

δp

ν+1

+

∂r

∂

T





ν

δT

ν+1

=

0

.

(29)

Here,thesuperscript

ν

denotestheiterationindex.AteachNewtoniteration,thelinearizedsystemofequationiswritten asJν

δ

xν+1

= −

rν .Inthissystem,Jν istheJacobianmatrixincludingallderivatives.Eachblock Jeα containsthederivatives ofequation e withrespecttounknown

α

,i.e., Jeα

= ∂

re

/∂

α

(containingthelistof, ∂r∂M Bp

|

ν ,

∂rM B ∂T

|

ν , ∂rE B ∂p

|

ν and ∂rE B ∂T

|

ν )and

δ

xν+1_{isthefullvectorofupdatestotheunknowns,precisely,}

_δ

_xν+1

_{= [δ}

_p

_,

_δ

_T

_]

T_.Clearly,

_δ

_p

_{= [δ}

_pm

_,

_δ

_pf1

_,

_...,

_δ

_pfnfrac

]

T _and

δ

T

= [δ

Tm

_,

_δ

_Tf1

_,

_...,

_δ

_Tfnfrac

]

T_{.SubscriptsM B andE B refertoequationsofmassbalance,andenergybalance.Thefulllinear} systemcanbeillustratedas:

⎛

⎜

⎜

⎜

⎜

⎜

⎜

⎝

⎛

⎝

JmmM Bp J mf M Bp Jfm_{M B} p J ff M Bp

⎞

⎠

Jmm M BT J mf M BT Jfm_{M B} T J ff M BT



⎛

⎝

JmmE Bp J mf E Bp Jfm_{E B} p J ff E Bp

⎞

⎠

Jmm_{E B} T J mf E BT Jfm_{E B} T J ff E BT



⎞

⎟

⎟

⎟

⎟

⎟

⎟

⎠

ν







Jν₀

⎛

⎜

⎜

⎜

⎜

⎝

δ

pm

δ

pf

δT

m

δT

f

⎞

⎟

⎟

⎟

⎟

⎠

ν+1







δxν₀+1

= −

⎛

⎜

⎜

⎜

⎜

⎝

rm M B rf_{M B} rm_{E B} rf_{E B}

⎞

⎟

⎟

⎟

⎟

⎠

ν

  

rν₀ (30)

Inordertoachieveaconvergedsolution,thefollowingconditionshavetobeassured:

||

r(M B)

||

2

||

r₍0_{M B)}

||

2

<



(M B)

∧

||

r(E B)

||

2

||

r0_{(E B)}

||

2

<



(E B)

∧

||δ

p

||

2

||

p

||

2

<



(p)

∧

||δ

T

||

2

||

T

||

2

<



(T) (31)

Here,eachthreshold(



x)isauser-definedtolerancesetinitiallyasinputatthebeginningofthesimulation.

To improvethe convergence, atthe beginning ofeach time-step, the mass balanceand energybalance equations are solvedasdecoupledsystemsinasequentialmanner.Inotherwords,first,decoupledlinearsystemofmassbalanceequation (32) issolvedandavectorofpressureupdates(i.e.,

δ

p)isobtained.Afterupdating thepressure,alldependentproperties andderivativesare re-calculated.Similarly,inthe nextstep,the decoupledlinearsystemofenergybalanceequation (33) is solved andavector of temperatureupdatesisacquired. Consequently, thetemperatureunknownsare updated andall dependent parameters and derivatives are recalculated. After thesesteps, the fully coupled linear system(30) issolved iteratively.

⎛

⎝

JmmM Bp J mf M Bp Jfm_{M B} p J ff M Bp

⎞

⎠

ν







JMBpν₀

δp

m

δp

f



ν+1







δpν₀+1

= −

rm_{M B} rf_{M B}



ν

  

rM Bν0 (32)

Jmm_{E B} T J mf E BT Jfm_{E B} T J ff E BT



ν







J_EBTν 0

δ

Tm

δ

Tf



ν+1







δT ν₀+1

= −

rm_{E B} rf_{E B}



ν

  

rE Bν0 (33)

TheextensiontopEDFMisexplainedinsection2.2.2.

Obtaining thesolutionofthislinearizedsystemofequationsiscomputationallychallengingduetothestrongcoupling ofmassandheattransfer,especiallyforrealfield-scaleapplications.Insection 3thesimulationframeworkofthisworkis presentedwithwhichresolvingthecomputationalchallengeisaimedusingtheADMmethod.

2.2.2. pEDFMimplementation

TocorrectforEDFMlimitationonfractureswithgenericconductivityduetoparalleltransmissibilities,thematrix-matrix, fracture-matrix andfracture-fracture connectivitiesare modified intheoverlapping regions. The mentionedmodifications eliminate the parallel transmissibilities, such that pEDFM is applicable to anyconductivity contrast between the matrix andfractures.Initially,allconnectivitiesbetweentwoneighboringmatrixcellsthataredisconnectedduetotheoverlapping fracturesaredetected.Duetogeometricalalgorithmdevisedduringthedevelopmentofthismethod,acontinuousprojection path(visibleinFig.3assolidlinesinlightbluecolor)isautomaticallyobtainedontheinterfacesassuchitdisconnectsthe neighboringconnectionslettingthefluxoccuronlyononeroute(i.e.,throughmatrix-fracture-matrix).

Let fracture element f overlap matrix grid cell i (as shownin Fig.3) withan area fraction Ai f. A set ofprojections is definedon the interface betweenthe overlappedmatrix grid cell i and its neighboring grid cells (in orange)that are affectedbythecrossing(i.e., j andk).Pleasenotethatinthe3Ddimensionalcase,therewillbethreeprojections.Foreach dimension(i.e.,x,y andz)theprojectionareafractionsareobtainedvia:

Fig. 3. pEDFMillustrationfora2Dmatrixwithstructuredgridanda1Doverlappingfracture.Theoverlappedmatrixcellsarehighlightedinyellowcolor. WithclassicalEDFMmethod,thesecellsareconnectedtoeachoverlappingfractureelementviatheconnectivityindexinthematrix-fracturefluxexchange. pEDFMintroducesextranon-neighboringconnectionsbetweenfractureelementsandmatrixcellshighlightedinorange.

Ai f⊥xe

=

Ai f

×

cos

(

γ

),

xe

∈ {

x,y,z

},

(34)

where

γ

istheanglebetweenthefractureelementandtheinterfaceconnectingthematrixgridcelli andtheneighboring gridcellinthecorrespondingdimension(showninFig.3onthezoomed-insectionandhighlightedinredcoloras Ai f⊥x and Ai f⊥y).Newtransmissibilitiesaredefinedtoconnectthefractureelement f toeachnon-neighboringmatrixgridcells (i.e., j andk inthe2DexampleshowninFig.3):

Tief

=

Ai f⊥xe

d

ief



ief

,

xe

∈ {

x,y,z

},

(35)

where,

d

i

_ef is theaverage distance betweenthe fracture element f andmatrix grid cell ie.



ief is the effectivefluid

mobilitybetweenthesecells. Asaresultofthenewtransmissibilities,theconnectivitybetweenthematrixgridcelli and

itscorrespondingneighboringcellsismodified: Tiie

=

Aiie

−

Ai f⊥xe

x

e



ief

, X

∈ {

x,y,z

}.

(36)

Forsimplicityoftheimplementation,themodifiedtransmissibilities areobtainedbymultiplicationofcoefficient

α

asa fractionoftheprojectedcross-sectionarea, andthecross-section areaofthecorresponding interface.Pleasenote thatfor alloverlappingfractureelements(exceptfortheboundariesofthefractures),theprojectionwillcoverthewholeinterface. Therefore,

α

is1

.

0 formostofthecases,resultinginzerotransmissibilitybetweenthematrixgridcells(i.e.,Tiie

=

0),thus

removingtheparalleltransmissibilities[27].

3. ADMmethodforflowandheattransfer 3.1. Solutionstrategy

Atthebeginningofevery Newtoniteration,thefine-scalefully-implicitlinearizedsystemofequationsfromEq.(30) is used to algebraically assemble a reducedlinear system ona dynamic multilevel grid resolutionwhich is definedat the beginningofeachtime-step.Atthefirststage(i.e.,obtainingthedynamicmultilevelgrid),setsofNml andNlfi hierarchically

nested grids are imposed on matrix m and fracture fi. Please note that l refers to the coarsening level at which the correspondingCartesiangridisspecified.Here,thefine-scalelevelisreferredtoasl

=

0,meaningN0_mandN0_f

i gridcellsare

definedatfine-scalelevelformatrixandfracturerespectively. Itisassumedthatthefine-scaleresolution isasufficiently accuraterepresentationofthedomaintocapturethephysicalphenomena.Thecoarseningratio(i.e.,theratiobetweenthe resolutionsoflevell andl

−

1 iswrittenas

γ

ml

=

Nl−1 m Nml and

γ

l_f i

=

Nl_f−1 i Nl_f i (37) formatrixm andfracture fi.The coarseningratios

γ

l areconstant onthehierarchicalsetofcoarseninglevels l,butboth propertiesaredefinedindependentlyformatrixandeachfracture, providingpracticalflexibility.Combinationofgridcells at different resolutions results in construction of the ADM grid. To map from the fine-scale resolution to this dynamic multilevelresolution, ADMemploys aseriesof restriction(R) andprolongation(P)operators, which areconstructed and computedonlyatthebeginningofthesimulation(seetheofflinestageofFig.4).TheADMsystemiswrittenas

Fig. 4. Schematic description of ADM strategy.

ˆ

R_ll−1

. . . ˆ

R0₁J0P

ˆ

10

. . . ˆ

Pll−1







JA D M

δ

x

ˆ

A D M

= − ˆ

Rl_l−1

. . . ˆ

R01r0







rA D M

.

(38)

Inthissystem,R

ˆ

l_l−1 istheso-calledrestrictionoperatorwhichmapsthesectionsofthesolutionvector(

δ

x

ˆ

A D M)atlevel

l

−

1 tolevell.Additionally,P

ˆ

l

l−1 istheprolongationoperatormappingpartsofthesolutionvectorfromlevell tolevell

−

1 [48,50].AftersolvingtheADMsystem,thesolutionacquiredatADMlevel(

δ

x

ˆ

A D M)isprolongedbacktofine-scaleresolution toobtainanapproximatedsolution(i.e.,

δ

x ₀),namely

Notethat

δ

x0 representsthereferencefine-scalesolution.TheADMRestriction(R

ˆ

l_l−1)andprolongation(P

ˆ

l_l−1)operators areobtainedfromthestaticmultilevelmultiscalerestriction(R_ll−1)andprolongation(Pl_l−1)operatorsrespectively.Thestatic prolongationoperatorPl_l−1 isconstructedbyassemblingthelocallycomputedbasisfunctions(seesection3.2)atlevell and

iswrittenas Pl_l−1

=

⎛

⎜

⎜

⎜

⎜

⎜

⎝

[(

Pp

)

l_l−1

]

_mm

[(

Pp

)

l_l−1

]

_mf 0 0

[(

Pp

)

l_l−1

]

_{f m}

[(

Pp

)

l_l−1

]

_{f f} 0 0 0 0

[(

PT

)

l_l−1

]

_mm

[(

PT

)

l_l−1

]

_mf 0 0

[(

PT

)

l_l−1

]

_{f m}

[(

PT

)

l_l−1

]

_{f f}

⎞

⎟

⎟

⎟

⎟

⎟

⎠

Nl−1×Nl

.

(40) NotethatPl

l−1 consistsoftwomaindiagonalblocksforthetwounknowns(i.e.,pressure p andtemperatureT ).Within each main block, the off-diagonalterms refer to the couplingterms betweenmatrix andfractures. The staticrestriction operatorRl_l−1 reads Rl_l−1

=

⎛

⎜

⎜

⎜

⎜

⎝

[

Rl_l−1

]

_m 0 0 0 0

[

Rl_l−1

]

_f 0 0 0 0

[

Rl_l−1

]

_m 0 0 0 0

[

Rl_l−1

]

_f

⎞

⎟

⎟

⎟

⎟

⎠

Nl×Nl−1

.

(41)

Toassuremassconservation,afinite-volumerestrictionoperatorisused.Therefore,Rl_l−1 elementscontaineither1 or 0 where

Rl_l−1

(u,

v)

=



1 if cell v is inside coarser cell u,

0 otherwise

.

(42)

Forthetestcasesthatarepresentedinthisarticle,differentinterpolatorsareemployedforeachofthevariables.The pro-longationoperatorforthepressure(i.e.,

(

Pp

)

ll−1)usesmultilevelmultiscalebasisfunctionsasinterpolator(seesection 3.2), whereasthe temperatureprolongationoperator containsconstant basisfunctions,namely

(

PT

)

l_l−1

= [

Rl_l−1

]

T

where super-script T refers to transpose operator. Therefore, the matrix-fracture coupling blocks of prolongation operator related to temperature(i.e.,

[(

PT

)

l_l−1

]

_mf and

[(

PT

)

l_l−1

]

_{f m})arezero.

3.2. Multilevelmultiscalebasisfunctions

The pressure prolongation operator

(

Pp

)

l_l−1 used in this work employs multiscale basis functionsat each coarsening level.Itcontainsdiagonalsub-blocks

[(

Pp

)

l_l−1

]

_mm,

[(

Pp

)

l_l−1

]

_{f f} andoff-diagonalsub-blocks

[(

Pp

)

l_l−1

]

_mf and

[(

Pp

)

l_l−1

]

_{f m}.The ADMmethodallowsfordifferentmatrix-fracturecouplingstrategiesforcalculationofbasisfunctions.Incaseofdecoupled basis functions, where the fracturesand matrixdo not influence the localsolutions of each other, theoff-diagonal sub-blocks consist of onlyzero elements. However, by choosing the fully-coupled approach, each medium will influence the localsolution(i.e.,thebasisfunctions)oftheoverlappingdomain,hencetheoff-diagonalblockscontainnon-zeroelements aswell[54,50].

The prolongation operator is constructed fullyalgebraically [54] and it includes the basis functionsat all coarsening levels. To compute the basis functions, first the fine-scale pressure system of equationsis assembled. Consequently, the coarse grids are imposed on the fine-scale computational grid cells formatrix and fractures. The approach to construct thehierarchicallynestedprimalanddualcoarsegridsisflexibleinthiswork.Oneoptionistoimposethecoarsegridsby determiningidenticalprimal coarsegrids ateachcoarseninglevel,namelyeachprimal coarsegrid atlevell has identical sizeof

γ

l

₌

_γ

l

x

×

γ

yl

×

γ

zl.Inthiscasethecoarsenodesarenotlocatedonthecornersoredgesoftheglobaldomain.Fig.5 showsanexampleoftheprimalandcoarsegridsofthisoptionfora2D domainof75

×

75 gridincluding3 fractures.The otheroptionistoincludethecoarsenodesonthecornersandedgesoftheglobaldomain,resultinginprimalcoarsegrids withdifferentsizes.Fig.6illustratesthecoarsegridconstructionofthisapproach.Pleasenotethatthelevelatwhicheach ofthemedia(eithermatrixorfractures)iscoarsenedcanbedefinedindependently.

Once the coarse grids are constructed, the basis functions are calculated algebraically for each coarsening level. The multiscaleprocedureconvertsthe two-pointfluxapproximation (TPFA)connectivitiesofthefine-scalesystem(i.e., A0_)to multi-pointfluxapproximation(MPFA)schemeatcoarse-scalesystemsforeachcoarseninglevell (i.e., Al).Therefore,prior toconstructionofcoarse-scalesystematlevell,thesystematlevell

−

1 isreducedtotheTPFAscheme,namely

Fig. 5. Coarse gridstructurewithidenticalprimalcoarsegridsateachcoarseninglevelfora2D 75×75 matrixand3 fractures.Thecoarseningratio inmatrixisγl

m,x=5 andγml,y=5 inx andy directionrespectively.Thenumberofcoarsegridatlevell=2 isNlm,x=_(γ75l m,x)2 =

3 inx directionand Nl

m,y=_(γ75l

m,y)2=3 iny direction.

Fig. 6. Coarsegridstructureincludingprimalcoarsenodesonthecornersandedgesfora2D76×76 matrixand3 fractures.Thecoarseningratioinmatrix isγl

m,x=5 andγml,y=5.Thenumberofcoarsegridatlevell=2 isNlm,x=_(γ76l−1

m,x)2+1=4 andN

l m,y=_(γ76l−1

m,y)2+1=4 inx andy directionrespectively.

Here, the T P F A_O perator removes the MPFA components from the entries of the coarse-scale system matrix. This procedureisperformedindependentlyoneachmedium[50].Fig.7showscoupledbasisfunctionsforamatrixcoarsenode (left)andtwofracturecoarsenodes(right).

3.3. Selectionofthegridresolution

The fine-scale grid resolution is mapped into the ADM grid resolution using the ADM operators. This dynamic grid resolutionisobtainedbyauser-definedgridselectioncriterionasafronttrackingtechnique.Ateachtime-stepn,theADM grid is selected basedon thesolution attime-step n

−

1 (explicit scheme) usinga thresholdset asthe input parameter forthe gridselectioncriterion. Thecondition forwhichthecriterion isconsidered, comparesthedifference ofthevalues ofa mainunknown (inthiscase, thetemperature)betweenneighboringcellsandthegrid cellsina coarsenode (spatial gradientoftheunknown)[48].Let



I

l and



J

l bethesetoftwoneighboringcoarsegridcells I and J atcoarseninglevell. The indicesoffine-scalegridcells belongingtothecoarsegridcells



_lI and



_lJ areindicatedbyi and j respectively. The temperaturedifference

 ¯

TI J isobtainedas

 ¯

TI J

=

max



|

Ti

−

Tj

|



Fig. 7. Exampleofcoupledbasisfunctionsfora2Dhomogeneousdomain.Thefigureontheleftshowsthebasisfunctionofamatrixcoarsenodeandthe figureontherightillustratesthebasisfunctionsoftwocoarsenodesonafracture.

Fig. 8. ExampleofanADMgridfora3Ddomainwithonefractureplateataspecifictime-stepduringthesimulation.Twocoarseninglevelsareemployed forthematrixandonecoarseninglevelisusedinthefracture.Here,themovingfrontofthefluidiscapturedinfine-scaleresolution.

Thecoarsegridblock I isrefinedfromcoarseninglevell toresolutionl

−

1 ifthecondition

 ¯

TI N

>

tol (45)

holds.Here,N referstoall coarsegridcellsneighboringcoarsegridcell I (atcoarseninglevell).Fig.8showsanexample of a 3D domain including one fracture with two coarsening levels over matrix and one coarsening level over fracture. Pleasenotetheindependencyofcoarseningstrategyandgridselectionprocedureformatrixandfracture.Itcanberealized thatthe movingfrontoftheflooding fluidwillbecapturedwithhigherresolutionasitcorresponds tosharpergradients. Moreover,the gridcells aroundthewellsare always keptatfine-scaleresolutionto guaranteeaccuratecaptureofsource termfluxes.

4. Numericalresults

Numericalresultsare presentedinthissection. First,thepEDFMmodel isvalidated. Thereafter,ADMresultsare com-paredagainstfine-scalesimulations,usedasareference.

4.1. Testcase1:validationofpEDFM

Thistest casedemonstratesthe validationofpEDFM implementation.Forthispurpose,a2D homogeneous domainof 3 m

×

3 m is considered (see Fig. 9). Twofractures intersect in themiddle of the domain forming a cross-shape. Both fracturesare1

.

5 [m]longwithaperture of8

×

10−3 _{[m].Intermsoffracturesconductivity,twoscenarios aretakeninto} account. In the first scenario,fractures are af

=

108 times more permeable than the rock matrix (Km

=

10−14 m2 and

Fig. 9. Visualization of test case 1 geometry, a 2D domain with a cross-shaped fracture network at the center.

Table 1

Inputparametersoftestcase1forfluidandrockproperties.

Property Value

Rock thermal conductivity (r) 4[W/m K]

Fluid thermal conductivity (f) negligible

Rock density (ρr) 2750[kg/m3]

Fluid specific heat (Cpf) 4200[J/kg K]

Rock specific heat (Cpr) 790[J/kg K]

Matrix porosity (φ) 0.2 Matrix permeability 10−14_[m2_] High-perm fractures permeability 10−6_[m2_] Low-perm fractures permeability 10−22_[m2_] Fractures aperture 8×10−3_[m]

Kf

=

10−6 m2) andinthesecondscenario,fracturesare108 timeslesspermeablethantherockmatrix(Km

=

10−14 m2 and Kf

=

10−22 m2). The initial pressure and temperature of the reservoir are p0

=

1

.

5

×

107

[

pa

]

and T0

=

400

[

K

]

respectively.Coldwaterwithpressurepinj

=

2

×

107

[

pa

]

andtemperatureTinj

=

400

[

K

]

isinjectedfromtheleftboundary andhotwaterisproducedfromtherightboundarywithpressure pprod

=

107

[

pa

]

.No-flowboundaryconditionsappliesto topandbottomboundaries.Table1showstheinputparametersofthistestcase.

Toobtainthereferencesolution(denotedasDNS),a375

×

375 gridresolutionisimposedonthedomain.Thisallowsfor fullyresolvingtheflowinsidethefracture.Inthiscasefracturesaredefinedaschannelsalongthe middlerowandcolumn ofthediscretizeddomain.EDFMandpEDFMsimulationsarerunwithfourdifferentmatrixgriddingresolutionof75

×

75, 45

×

45 and25

×

25.Atallruns,thefracturesapertureisidenticalandsettobeaf

=

8

×

10−3

[

m

]

whichisapproximately equal tothematrixgridcell sizeofDNSdiscretization.The sizeofgridcells insidethefractureshavesimilardimensions as matrixgrid cells. Incases withhighly conductive fractures, thesimulation isrun for 6

[

hours

]

andinthe caseswith low-permeabilityfractures,thesimulationisrunfor12

[

hours

]

.Theresultsareprovidedat50 isochronalintervals.

In thefirstscenario,fracturesareconsidered tobe highlyconductive (seeTable 1formoredetails).The resultsofthe simulationsasplotsofpressureandtemperaturesolutionsareshownintheFig.10and11.

TheplotsinFig.12presentthepressureandtemperatureerrorsasfunctionsofthesimulationtime.Theerrorofvariable

x (i.e.,pressureortemperature),denotedasex,iscalculatedvia

ex

=

||

xDNS

−

xEDFM/pEDFM

||

2

||

xDNS

||

2

.

(46)

Atthecenterofthedomainwherethefracturesintersect,adifferenceinthetemperaturedistributionbetweentheDNS, EDFMandpEDFMmethods canbe observed.Please note thatthe differencearisesfromthediscretization approaches. In DNSmethod,fracturesareactuallychannelswithaperture ofone gridcell.Asa result,thosematrixgridcellsareflooded withthecoldfluid.ThisisnotthecaseforEDFMandpEDFM.

Inthesecondscenario,bothfracturesaresettobe108 _{timeslesspermeablethanthematrixrock.Figs.13and14show} thesimulationresultsofthisscenarioaspressureandtemperaturesolutionsrespectively.

Fig.15showsthepressureandtemperatureerrorsversusthesimulationtime forthesecond scenario.AstheEDFMis incapableofcapturingthelowpermeablefracturesandtheresultsarenotaccurateandrepresentative,theerrorsareonly shownforpEDFM.

Fig. 10. Testcase1:High-permeablefractures:comparisonofpressuresolutionbetweenthefullyresolvedDNS(singleplotatleftcolumn)using375×375 gridcells,pEDFMresults(toprow)andEDFMresults(bottomrow).BothpEDFMandEDFMarerunwithdifferentgridresolutions(lefttoright:75×75, 45×45 and25×25)at25th timeinterval(i.e.,t=3hours).

Fig. 11. Testcase1:High-permeablefractures:comparisonoftemperaturesolutionbetweenthefullyresolvedDNS,pEDFMandEDFM.Similartoprevious figure,thesingleplotatleftcolumnshowsDNSsolution,toprowshowsthepEDFMresultsandbottomrowshowstheEDFMresults.

Fig. 12. Testcase1:High-permeablefractures:PressureandTemperatureerrorsofpEDFMandEDFMsimulationswithrespecttothereferencesolution (DNS).

4.2. Testcase2:ADMwitha2Dhomogeneousfracturedreservoir

Inthistestcase,theaccuracyoftheADMmethodformass-heattransportinahomogeneousdomaincontainingfractures withmixedconductivitiesisstudied.Forthispurpose,a2Dfractured100m

×

100m domainwith30 fracturesisconsidered. Thelength ofeachfractureisdifferentbutthesizeoftheir apertureisidenticalandsetto af

=

5

·

10−3 m. A136

×

136

Fig. 13. Testcase1:Low-permeablefractures:comparisonofpressuresolutionbetweenthefullyresolvedDNS(singleplotatleftcolumn)using375×375 gridcells,pEDFMresults(toprow)andEDFMresults(bottomrow).BothpEDFMandEDFMarerunwithdifferentgridresolutions(lefttoright:75×75, 45×45 and25×25)at25th timeinterval(i.e.,t=3hours).

Fig. 14. Testcase1:Low-permeablefractures:comparisonoftemperaturesolutionbetweenthefullyresolvedDNS,pEDFMandEDFM.Similartoprevious figure,thesingleplotatleftcolumnshowsDNSsolution,toprowshowsthepEDFMresultsandbottomrowshowstheEDFMresults.

Fig. 15. Test case 1: Low-permeable fractures: Pressure and Temperature errors of pEDFM compared to reference solution (DNS).

grid isimposed on therockmatrix andthefracture networkconsistsof 315 gridcells (intotal 18811 cells).The matrix permeabilityisKm

=

10−14m2andthepermeabilityofthefracturenetworkhastherangeofKfmin

=

10−20m

2_andK fmax

=

10−7 m2.Twoinjectionwellsarelocatedatthebottomleftandtopleftcornerswithinjectionpressureofpinj

=

2

·

107Pa. Additionally,therearetwoproductionwellsatthebottomrightandthetoprightcornerswithpressureofpprod

=

1

·

107Pa. Table2demonstratestheinputparametersofthistestcase.

The coarseninglevel for matrix islm

=

2 andfor fracturesis lf

=

1. Over the entire computational grids, the coars-ening ratio is

γ

=

3

×

3. The ADM simulations are run for three different thresholds of coarsening criterion, namely,

Table 2

Inputparametersoftestcase2forfluidandrockproperties.

Property Value

Rock thermal conductivity (r) 4[W/m K]

Fluid thermal conductivity (f) negligible

Rock density (ρr) 2750[kg/m3]

Fluid specific heat (Cpf) 4200[J/kg K]

Rock specific heat (Cpr) 790[J/kg K]

Matrix porosity (φ) 0.2 Matrix permeability 10−14_[m2_] fractures permeability (min) 10−20_[m2_] fractures permeability (max) 10−7_[m2_] Fractures aperture 5×10−3_[m]

Fig. 16. Testcase2:pressureplots.Fig.16aillustratesthefracturenetworkofthistestcase.Thefigureontoprightisthefine-scalesolution(16b),andthe figuresonthebottomrowshowtheADMsolutionwiththresholdsofT=5 (16c),T=10 (16d)andT=20 (16e)respectively.



T

=

5

,

10

,

20

[

◦K

]

andthe ADMresultsare comparedto resultoffine-scale simulationasreferencesolution. TheADM errorforpressureandtemperatureisobtainedvia

ex

=

||

xFS

−

xADM

||

2

||

xFS

||

2

(47) where, x iseitherpressure ortemperatureandthe subscriptsADM and FS denote ADMandfine-scale,respectively. The simulationisrunfor1000

[

days

]

andtheresultsareprintedat100 isochronalintervals.

Figs.16and17show thepressureandtemperaturesolutionoffine-scaleandADMsimulations.Theerrors(e) andthe percentageofactivegridcells( AGC )arementionedundertheADMplots.

Fig. 18 provides more details of the error analysis for test case 2. On the top rowpressure andtemperature errors (Fig.18aand18b)versussimulationtime-stepsareshown.Fig.18dplotsthepercentageofactivegridcellsusedduringthe simulation.Lastly,averagederrorsandaverageactivegridcellsareplottedovertheentiresimulationforallthreethresholds. Pleasenotethattheregionsaroundthewellsarealwayskeptatfine-scaleresolutionandasaresult(dependingonthe sizeofthedomain)anoticeableportionofactivegridcellsareimposed.Inthistestcase,approximately5 percentoftotal gridcellsarekeptatfine-scaleduetonear-wellrefinement.

Pleasenotethatinthistestcase,intersectionbetweenpermeablefracturesandimpermeablefracturesoccur.Bydefault, theimpermeablefracturesareconsideredas“non-sealing”flowbarriers,meaningthatthecrossingofthehighlypermeable

Fig. 17. Testcase2:temperatureplots.Fig.17bontoprightshowsthefine-scalesolution.TheADMsolutionsareshownonthefiguresofbottomrowwith thresholdsofT=5 (17c),T=10 (17d)andT=20 (17e).

Fig. 18. Testcase2:erroranalysis.Fig.18aand18b showADMerrorofpressureandtemperatureoversimulationtime-steps.Theerrorsarecalculated basedonEq.(47).Fig.18cshowsthepercentageofactivegridcellsusedduringsimulationversusthetime-steps.Lastly,Fig.18ddemonstratestheaveraged valuesofthementionedparametersovertheentiresimulationtimeforthreethresholds.