A benchmark study of the multiscale and homogenization methods for fully implicit multiphase flow simulations

(1)

A benchmark study of the multiscale and homogenization methods for fully implicit

multiphase flow simulations

Hajibeygi, Hadi; Olivares, Manuela Bastidas; HosseiniMehr, Mousa; Pop, Sorin; Wheeler, Mary

DOI

10.1016/j.advwatres.2020.103674

Publication date

2020

Document Version

Final published version

Published in

Advances in Water Resources

Citation (APA)

Hajibeygi, H., Olivares, M. B., HosseiniMehr, M., Pop, S., & Wheeler, M. (2020). A benchmark study of the

multiscale and homogenization methods for fully implicit multiphase flow simulations. Advances in Water

Resources, 143, [103674]. https://doi.org/10.1016/j.advwatres.2020.103674

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

A

benchmark

study

of

the

multiscale

and

homogenization

methods

for

fully

implicit

multiphase

ﬂow

simulations

Hadi

Hajibeygi

a,∗

_,

_Manuela

_Bastidas

_Olivares

b

_,

_Mousa

_HosseiniMehr

a

_,

_Sorin

_Pop

b

_,

Mary

Wheeler

c

a Faculty of Civil Engineering and Geosciences, Delft University of Technology, P.O. Box 5048, 2600 GA Delft, the Netherlands b Faculty of Sciences, Hasselt University, Diepenbeek, Belgium

c Institute for Computational Engineering and Sciences, The University of Texas at Austin, 201 East 24th Street, ACE 5.324, Campus Mail C0200, Austin, TX 78712 US

a

r

t

i

c

l

e

i

n

f

o

Keywords: Multiscale Homogenization

Algebraic dynamic multilevel Adaptive mesh reﬁnement Flow in porous media Fully implicit simulation

a

b

s

t

r

a

c

t

Accuratesimulationofmultiphaseflowinsubsurfaceformationsischallenging,astheformationsspanlarge lengthscales(km)withhigh-resolutionheterogeneousproperties.Todealwiththischallenge,differentmultiscale methodshavebeendeveloped.Suchmethodsconstructcoarse-scalesystems,basedonagivenhigh-resolution fine-scalesystem.Furthermore,theyareamenabletoparallelcomputingandallowfora-posteriorierrorcontrol. Themultiscalemethodsdifferfromeachotherinthewaythetransitionbetweenthedifferentscalesismade. Multiscale(finiteelementandfinitevolume)methodscomputelocalbasisfunctionstomapthesolutions(e.g. pressure)betweencoarseandfinescales.Instead,homogenizationmethodssolvelocalperiodicproblemsto determineeffectivemodelsandparameters(e.g.permeability)atacoarserscale.Itisyetunknownhowthesetwo methodscomparewitheachother,especiallywhenappliedtocomplexgeologicalformations,withnoclearscale separationinthepropertyfields.Thispaperdevelopsthefirstcomparisonbenchmarkstudyofthesetwomethods andextendstheirapplicabilitytofullyimplicitsimulationsusingthealgebraicdynamicmultilevel(ADM)method. Ateachtimestep,onthegivenfine-scalemeshandbasedonanerroranalysis,thefullyimplicitsystemissolvedon adynamicmultilevelgrid.Theentriesofthissystemareobtainedbyusingmultiscalelocalbasisfunctions (ADM-MS),and,respectively,byhomogenizationoverlocaldomains(ADM-HO).Bothsetsoflocalbasisfunctions (ADM-MS)andlocaleffectiveparameters(ADM-HO)arecomputedatthebeginningofthesimulation,withnofurther updatesduringthemultiphaseflowsimulation.Thetwomethodsareextendedandimplementedinthesame open-sourceDARSim2simulator(https://gitlab.com/darsim2simulator),toprovidefairqualitycomparisons.The resultsrevealinsightfulunderstandingofthetwoapproaches,andqualitativelybenchmarktheirperformance.It isre-emphasizedthatthetestcasesconsideredhereincludepermeabilityfieldswithnoclearscaleseparation.The developmentofthispapershedsnewlightsonadvancedmultiscalemethodsforsimulationofcoupledprocesses inporousmedia.

1. Introduction

Geologicalformationsspanlarge(km)lengthscales,having hetero-geneousproperties characterizedathighresolutions(cmandbelow). Asfortheuncertaintywithintheintegratedfielddata,typically, sev-eralequiprobable realizationsoftheproperty fieldsaregeneratedto studyandsimulatethefluidflowandtransport.Classicalsimulation ap-proachesaretooexpensiveforsuchstudies.Therefore,advanced sim-ulationmethodsarerequiredtoallowforanaccuraterepresentation of theheterogeneousproperties.Atthesame time,they should

pro-∗_{Corresponding}_author.

E-mail addresses: H.Hajibeygi@tudelft.nl (H. Hajibeygi), manuela.bastidas@uhasselt.be (M.B. Olivares), S.Hosseinimehr@tudelft.nl (M. HosseiniMehr), sorin.pop@uhasselt.be(S.Pop),mfw@ices.utexas.edu(M.Wheeler).

videaneﬃcientsimulationframeworktostudymultiplerealizations

Jansenetal.(2009);Wachspress(1966).

Modelorderreductiontechniqueshavebeendevelopedtoprovidea meaningfulapproximatesimulationframework.Suchtechniqueshave tobefastenoughtobeappliedtolarge-scalecomputationaldomains. Inthissense, anyadvancedmethodof thistypecanbe seenasﬁeld applicable onlyifitallowsforreducingtheerror belowanydesired thresholdvalueHajibeygietal.(2012).

Here we only consider numerical model order reduction tech-niques,among whichmultiscale EfendievandHou (2009); Hou and

Wu (1997) andhomogenization Weinan (2011)methods stand very

promising.

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

Received10October2019;Receivedinrevisedform17June2020;Accepted2July2020 Availableonline3July2020

(3)

These approaches are diﬀerent in the sense that the multi-scale method deals with crossing the solution, e.g., the pressure, across the scales Aarnes and Hou (2002); Jenny et al. (2003);

Hajibeygi et al. (2008); Chung et al. (2015), whereas in the latter eﬀective, lower-resolution parameters andfunctions like the perme-ability or the transmissibility, arederived Weinan andYue (2004);

Abdulleetal.(2012);Weinanetal.(2007);Lietal.(2020);Singhand Wheeler(2018);Vasilyevaetal.(2020).Moreover,whilethemultiscale basisfunctionshavebeenexpressedinapurelyalgebraicformulation

Wangetal.(2014),thesamedoesnotholdforthehomogenization ap-proach.Speciallytheintegrationofhomogenizedparameterswithinthe fullyimplicitframeworkinanalgebraicmannerhasnotyetbeen devel-opedsofar.Thepresentworkisaﬁrststepinthisdirection.

Atthesametime,thetwomethodshavemanysimilarities.Bothﬁnd theirmappingstrategyvialocalsolutionsoftheoriginalgoverning equa-tionswithlocalboundaryconditions.Multiscalebasisfunctionsoften employreduced-dimensionalboundaryconditionsTeneetal.(2015);

MøynerandLie(2016), whilehomogenization schemesimpose peri-odic boundaryconditionsonlocal problems,andconsiderlocal rep-resentative micro-structures even in the case of non-periodic prop-erties Allaire (1992); Abdulle and Weinan (2003); Arbogast and Xiao(2013);Bastidas etal.(2019);Brownetal.(2013).Both meth-ods are eﬀective for global equations within the fully coupled sys-tem of local-global unknowns, e.g., the global pressure and the lo-calsaturation.Bothhavebeenextendedtononlinearandgeologically complexmodelsAmanbeketal.(2019a);HosseiniMehretal.(2018);

Singhetal.(2019a).Recentdevelopmentsofthesetwoclassesof ap-proacheshaveintroducedafully-implicitdynamic multilevel simula-tion framework(ADM) in which heterogeneousdetailed geo-models aremappedintoadaptivedynamiccoarsermeshCusinietal.(2018);

Faigleetal.(2014);Klemetsdaletal.(2020);Carciopoloetal.(2020). TheADMmethoddevelopsafully-implicitdiscretesystemfor cou-pled flow and transport system of equations, in which each equa-tion can be represented at a different resolution than the defined fine-scale one. More importantly, the procedure can be done fully algebraic, with the dynamic mesh resolution defined based on a front-tracking strategy. In contrast to the rich existing literature of Adaptive MeshRefinement(AMR) methodsPauet al.(2009,2012);

BergerandOliger(1984);SchmidtandJacobs(1988);Edwards(1996);

Sammon(2003);KlemetsdalandLie(2020),ADMcanbedeﬁnedasan adaptivemeshcoarseningstrategywhichisconvenientlyapplicablefor heterogeneousandnonlinearcoupledsystemsCusinietal.(2016).

Irrespectiveof the choiceof thedynamic meshstrategy, it is al-waysachallenge toconstruct adaptive multiscaleentriesof the im-plicitsystems.TheADMmethodsofarhasincludedmultiscalebasis functionsCusini etal. (2016). Inaddition,homogenization methods havealsobeendevelopedformultiphasesimulationsondynamicgrids

Amanbeketal.(2019a);Cusinietal.(2019).Inthiscontext,twoaspects canbeofinterest:thestudyofthehomogenization-basedcoarsersystem entriesandthedevelopmentofabenchmarkstudyofthequalityofthe twoapproachesofADM-multiscale(ADM-MS)andADM-homogenized (ADM-HO)forcoupledimplicitmultiphaseﬂowscenarios.

ThispaperdevelopssuchaunifiedframeworkinwhichtheADM methodisextendedtoaccountforbothmultiscaleandhomogenization schemesformultiphaseflow simulations.Thisdevelopment makesit possibletoallowfordifferentcoarse-scaleentriesfordynamic simula-tions,andimportantlytobenchmarkthetwoclassesofmultiscaleand homogenizationstrategies.Noteworthyis that,oncetheeffective pa-rametersarecomputed,allotherhomogenizationproceduresare imple-mentedalgebraically.Thisisdonebyintroducingconstantunitylocal basis,withthesupportofprimal(non-overlapping)coarse-scale parti-tions.ThemultiscaleADMisimplemented fullyalgebraicsincelocal basisfunctionsarealsosolvedalgebraicallyovertheoverlapping(dual) coarsegriddomainsZhouandTchelepi(2012).Theoutcomeofthis de-velopmentismadeavailabletothepublicviaanopen-sourceDARSim2 simulator,https://gitlab.com/darsim2simulator.

Numericaltestcasesareconsideredforthechallenging,highly het-erogeneousSPE10ChristieandBlunt(2001).Thenumberofactivegrid cells,pressureandsaturationerrors,andthesolutionmapsareall re-portedindetail.Thedevelopmentofthispapershedsnewlightsinthe applicationofmultiscaleandhomogenizationapproachesinadvanced next-generationenvironmentsforﬁeld-relevantsimulationscenarios.

Thepaperisstructuredasfollows.Next,inSection2,the mathe-matical modelis statedbrieﬂy.Section3presentsthecomputational framework for both multiscale and homogenization ADM methods.

Section4presentsthetestcases,andconclusionsaredrawninSection5. TheAppendixgivesmoredetailsonthemultiscaleandthe homogeniza-tionapproaches.

2. Governingequations

Weconsiderﬂowoftwoimmiscibleandincompressiblephasesof

𝛼 and𝛽 throughaheterogeneousporousmedium.AttheDarcyscale, massbalanceforthephasei∈{𝛼,𝛽}reads

𝜕

𝜕𝑡(𝜙𝜌𝑖𝑆𝑖)−∇⋅ (𝜌𝑖𝜆𝑖⋅ (∇𝑝−𝜌𝑖𝑔∇𝑧))=𝜌𝑖𝑞𝑖. (1)

Here,𝜙 istheporosityofthemedium,𝜌i[kg/m3]andSiarethedensity andsaturationofthephasei,respectively.Thephasemobilitytensor

𝜆iisequalto𝐾𝐾_𝑟𝑖∕𝜇𝑖,whereK[m2]istherockabsolutepermeability tensor,and𝐾_𝑟𝑖=𝐾_𝑟𝑖(𝑆_𝑖)isthesaturation-dependentrelative permeabil-ityofphasei.Moreover,𝜇 [Pa.s]isthephaseviscosity.Fortheease ofpresentation,thetwophasepressuresareassumedequal,𝑝=𝑝_𝛼=𝑝_𝛽 [Pa](seee.g.AzizandSettari(2002)).However,theextensionto mod-elsinvolvingasaturationdependentcapillarypressureisalsopossible. Inaddition,g[m/s2_]_is_the_{gravitational}_acceleration_which_acts_on_∇_z

direction,andq[1/s]isthephasesourceterm.

Hereitisassumedthatthetwoﬂuidsareoccupyingcompletelythe porespace,andnootherﬂuidphaseispresent.Thisgivestheconstraint

𝑆𝛼+𝑆_𝛽=1,whichreducesthenumberofunknownsintheabove equa-tionstotwo:S_𝛼(inshortfromhereon,S)andp.Finally,themodelis completedbyinitialconditionsforthesaturation,andwithboundary conditions.Wedonotspecifythemexplicitlysincenoneofthemplaya roleinthemultiscalestrategy.

The fully-implicit coupled simulation approach Aziz and Set-tari (2002)estimates alltheparametersatnexttimestep (𝑛+1). As such,thesemi-discretenonlinearresidualforthephasei∈{𝛼,𝛽}reads

𝑟𝑛+1 𝑖 =[𝜌𝑖𝑞𝑖]𝑛+1− ( 𝜙𝜌𝑖𝑆𝑖)𝑛+1−(𝜙𝜌𝑖𝑆𝑖)𝑛 Δ𝑡 +∇⋅ ( 𝜌𝑖𝜆𝑖⋅(∇𝑝−𝜌𝑖𝑔∇𝑧))𝑛+1. (2) For ﬁnding thesolution pair (𝑝𝑛+1_,_𝑆𝑛+1₎ _one _needs _to_employ _a

linearization scheme. Here we restrict thediscussion totheNewton scheme,whichis2nd-orderconvergentbutrequiresastartingpointthat iscloseenoughtothesolution.Inotherwords,thetimestepmaybe sub-jecttorestrictionsalsodependingonthemeshsize.Alternatively,one mayconsiderapproacheslikethemodiﬁedPicardCeliaetal.(1990)or theL-SchemeRaduetal.(2017),whicharelessdemandingfromthe computationalpointof view,ormorerobust w.r.t.thestartingpoint and meshresolution, butconverge slowerthan theNewton scheme

Bastidasetal.(2019).Appliedto(2),theNewtonlinearizationreads

𝑟𝑛+1_≈_𝑟𝜈₊𝜕𝑟

𝜕𝑝|𝜈𝛿𝑝𝜈+1+𝜕𝜕𝑆𝑟|𝜈𝛿𝑆𝜈+1, (3)

whichcanbeexpressedalgebraicallyas𝐉𝜈_𝛿𝐱𝜈+1₌₋_𝐫𝜈_,_i.e., ⎡ ⎢ ⎢ ⎣ 𝜕𝑟𝛼 𝜕𝑝 𝜕𝑟𝜕𝑆𝛼 𝜕𝑟𝛽 𝜕𝑝 𝜕𝑟𝛽 𝜕𝑆 ⎤ ⎥ ⎥ ⎦ ⏟⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏟ 𝐉 𝜈_[ 𝛿𝑝 𝛿𝑆 ] ⏟⏟⏟ 𝛿𝐱 𝜈+1 =− [ 𝑟𝛼 𝑟𝛽 ] ⏟⏟⏟ 𝐫 𝜈 . (4)

Ineachtimestep,thelinearEq.(4)issolvediteratively(innerloop) severaltimesuntilnonlinearconvergence(outerloop)isreached.The

(4)

Fig.1. ExampleofanADMgrid(4thfromthetop),obtainedbycombining ﬁne-scale(top)andcoarserresolutionsoflevel1(2ndfromthetop)andlevel 2(3rdfromthetop).Alsoshownisthesaturationproﬁlecorrespondingtothe ADMgrid(bottom).

overallcomputationalcomplexityofthesimulationdependshighlyon thecomplexityofthesolutionofthislinearsystem.Advancedmultiscale andhomogenizationmethodsaimatsolvingthislinearsystemona dy-namicmultilevelmesh.Notethat,asshownbeforeCusinietal.(2018), theoveralleﬃciencyofanyadvancedmethodshouldincludenotonly thespeedupofsolvingthelinearEq.(4)butalsothecountoftheNewton (outer)loops.Next,theADMmethodbasedonmultiscaleand homoge-nizationformulationsispresented.

3. Dynamicmultilevelsimulationbasedonmultiscaleand homogenizationmethods

3.1. ADMFrameworkformulation

Thefully-implicitlinearsystem(4)istooexpensivetobesolvedfor realfieldscenarios.Amultileveldynamicmesh,asshowninFig.1,is generatedwithintheADMframework,basedonanerrorestimate strat-egy.Theerrorestimateisdevelopedbasedonafronttrackingcriterion, whichappliesfine-scalegridsonlyatsub-regionswithsharpgradients. Thefine-scalesystemisthenalgebraicallyreducedintothismultilevel grid,throughsequencesofrestrictionandprolongationoperators.To obtaintheADMgrid,first,setsof𝑁𝑙=𝑁_𝑥𝑙×𝑁_𝑦𝑙hierarchicallynested coarsegridsareimposedonthefinemesh.Here,lindicatesthe

coars-eninglevel.Moreover,𝛾l_is_the_coarsening_ratio_which_is_deﬁned_as 𝛾𝑙₌₍_𝛾𝑙 𝑥,𝛾𝑦𝑙)= ( 𝑁𝑙−1 𝑥 𝑁𝑙 𝑥 , 𝑁 𝑙−1 𝑦 𝑁𝑙 𝑦 ) , (5)

fortwo-dimensional(2D)domains.TheADMgridisconstructedby as-semblingacombinationofcellsatdiﬀerentresolutionswithinthe com-putationaldomain.Byusingthesequenceofrestriction(R)and prolon-gation(P)operators,onecanexpresstheADMsystemas

̂𝐑𝑙−1 𝑙 … ̂𝐑01𝐉𝟎̂𝐏 1 0… ̂𝐏𝑙𝑙−1 ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟ 𝐉ADM 𝛿 ̂𝑥ADM=−̂𝐑𝑙−1_𝑙 … ̂𝐑01𝑟0 ⏟⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏟ ̂𝐫ADM . (6)

Here, ̂𝐑𝑙−1_𝑙 is therestrictionoperatorwhichmapsthepartsof the solutionvectorthatareatlevel(𝑙−1)tolevell.Similarly,the prolonga-tionoperator̂𝐏𝑙_𝑙−1mapsthepartsofthesolutionvectorthatareatlevel

ltolevel𝑙−1.OncetheADMsystem(6)issolved,theapproximated ﬁne-scalesolution𝛿𝑥′

0canbeacquiredbyprolongingtheADMsolution

𝛿 ̂𝑥ADM,i.e.

𝛿𝑥0≈𝛿𝑥′0=̂𝐏 1

0… ̂𝐏𝑙𝑙−1𝛿𝑥ADM. (7)

TheADMRestriction ̂𝐑𝑙−1_𝑙 andprolongation ̂𝐏𝑙_𝑙−1operatorsare as-sembledusingthestaticmultilevelmultiscalerestriction𝐑𝑙−1

𝑙 and pro-longation𝐏𝑙

𝑙−1operators,respectively.Theyareconstructedonlyatthe beginningofthesimulationandarekeptunchangedthroughoutthe en-tiresimulation.

Thestaticprolongationoperator𝐏𝑙

𝑙−1isconstructedasanassembly ofthelocallycomputedbasisfunctionsateachcoarseninglevelland reads 𝐏𝑙 𝑙−1= ( (𝑃_𝑝)𝑙_𝑙−1 0 0 (𝑃_𝑆)𝑙_𝑙−1 ) 𝑁𝑙−1×𝑁_𝑙 . (8)

Here,(𝑃_𝑝)𝑙_𝑙−1and(𝑃_𝑆)𝑙_𝑙−1 arethetwomaindiagonalblocks corre-spondingtomainunknowns(i.e.,pressurepandsaturationS).Inthe caseofusingthehomogenizationscheme,i.e.ADM-HO,aswillbe de-scribedinSection3.3),constantbasisfunctionsforpressureareused. However,forthemultiscale-basedADM,i.e.ADM-MS,aswillbe de-scribedinSection3.2,locally-computedbasisfunctionsareused.Note thatthesaturationprolongationoperatorforbothapproachesis con-stantunityfunctionatallcoarseninglevels,whichrepresentsthe con-servativeﬁnite-volumeintegration.

Thestaticrestrictionoperator𝐑𝑙−1 𝑙 reads 𝐑𝑙−1 𝑙 = ( (𝑅)𝑙−1_𝑙 0 0 (𝑅)𝑙−1_𝑙 ) 𝑁𝑙×𝑁𝑙−1 . (9)

Inthiswork,aﬁnite-volumerestrictionoperatorisusedtoguarantee localmassconservation,i.e.

𝑅𝑙−1 𝑙 (𝑖,𝑗)=

{

1 ifcelliisinsidecoarsercellj,

0 otherwise. (10)

3.2. ADMUsingmultiscale(ADM-MS)

IntheADM-MS method,theprolongationoperatorforpressureis foundbasedonmultiscalebasisfunctions.Theselocalbasisfunctions arecomputedalgebraicallyWangetal.(2014),basedonthesteady-state pressureequation.Inthisstudy,theincompressibleﬂowequation (ellip-ticpressureequation)isusedtoconstructthemultiscalebasisfunctions

Teneetal.(2015).AnexampleofabasisfunctionisshowninFig.2. Thecoarse grid constructionandcomputationof multiscalebasis functionsareexplainedinmoredetailsinAppendixA.

3.3. ADMUsinghomogenization(ADM-HO)

Homogenizationisanothermethodthatcanbeappliedtoproblems involvingmultiplescales.Inthismethod,oneusesthemathematical modelsat micro(ﬁne) scale(1) toderive eﬀectiveupscaled models

(5)

Fig.2. Anexampleofabasisfunctionbelongingtothemiddlecoarsenodeof aheterogeneous2Ddomain.

Fig.3. SketchofthecoarsepartitionofΩwhenusingtwocoarseninglevels.

andparametersinwhichtherapidlyoscillating charactersticsare av-eragedout.Indoingso,theupscaledmodelmayhaveadiﬀerent struc-turethantheonesattheﬁnescale.WerefertoAmazianeetal.(2017);

Bourgeatetal.(1996);Hornung(1997);vanDuijnetal.(2007)for the-oreticaldetails.

ThegoalofthisworkistobuildaunifiedADMplatformwherethe multiscaleandhomogenizationmethodscanbecompared.Therefore, here,thehomogenization methodis usedonly toconstructeffective propertiesatthedynamicmultilevelmesh.Inthissetup,the homoge-nizedpropertiesofADM-HOatmultilevelmesharefoundasinADM-MS bysolvinglocalflow(pressure)equationsbasedonanincompressible (elliptic)equation.

Moreprecisely,oneassumesthatascaleseparationholdsand dou-blesthespatialvariableintoafastandaslowone.Themethodrelieson thehomogenizationansatz,meaningthatallquantitiesin(1)canbe ex-pandedregularlyintermsofascaleseparationparameter.Suchideasare employedinBastidasetal.(2019);AbdulleandNonnenmacher(2009);

Amanbeketal.(2019b);Amazianeetal.(1991);Singhetal.(2019b);

Szymkiewiczetal.(2011);Henningetal.(2015,2013)todevelop ef-fectivenumericalsimulationschemesevenincaseofnon-periodic me-dia.Moredetailsaboutthehomogenizationprocedurecanbefoundin

AppendixB.

Inthepresentcontext,foragivenﬁne-scaleeﬀectivepermeability

Kandforeachcoarseninglevell,aneﬀectivepermeabilitytensorKl_is computedlocallyinapre-processingstep.FirstthedomainΩ⊂ ℝ2 _is

dividedintocoarsecellsΩlthatcorrespondtoapartitionofthedomain ΩasshowninFig.3.

Fig.4. Exampleofthelocalsolutions𝜔1 _(top_right,_for_x_-direction)_and_𝜔2 (bottomright,fory-direction)foracoarsecellinsidea2Ddomain.The hetero-geneouspermeabilityﬁeldisalsoshownfortheentiredomain(left).

ForeachcoarsecellΩ_l atlevell,thecomponentsoftheeﬀective permeabilitytensorarecalculateas

𝐊𝑙 𝑖,𝑗|||| Ω_𝑙 = ∫Ω_𝑙 ( 𝐾(𝐞𝑗+∇𝜔𝑗))⋅ 𝐞𝑖𝑑y. (11) for𝑖,𝑗=1,2.Here𝜔j_are_the_periodic_solutions_of_the_pressure_equation onlocaldomains(knownasmicro-cellequationinHOliterature),i.e. −∇⋅(𝐾(∇_𝑦𝜔𝑗+𝐞_𝑗))=0,forally∈ Ω_𝑙. (12) Weremarkthat{𝐞_𝑗}2

𝑗=1isthecanonicalbasisofdimension2andKisthe abovementionedpermeabilitytensor.Toguaranteetheuniquenessof thesolution𝜔j_one_assumes_that_its_average_value_over_the_local_coarse cellΩ_lis0.

Todeterminethevalueoftheeﬀectivepermeabilitytensorateach coarsecellΩl,twolocal(micro-cell)problems(12)aresolvedforeach spatialdirectionin2D.Fig.4providesanillustrationoftheselocal so-lutionsforacoarseelement.

Notethatthelocalproblems(12)capturetherapidlyoscillating char-acteristicswithinacoarse element,completelydecoupledfromother coarse elements.Thehomogenizedparameters,like multiscalebases, arecomputedatthebeginningofthesimulation.Fig.5illustratesthe calculationoftheeﬀectivepermeabilityatdiﬀerentlevels.

Thehomogenizedparametersareusedtoconstructthecoarse sys-tementries.Moreprecisely,thehomogenizedvalueinacoarsecellis distributedequallytothefinecellsconstructingit.Thenthefine-scale Jacobianandresidualarecomputedwiththefine-scalesaturationfield. ThissystemisthenmappedtotheADMresolutionbysetting prolonga-tionoperatorsin(6)tounity.Thisisaconvenientprocedure,developed inthiswork,tointegratethenumericalhomogenizationmethodwith anexistingadvancedsimulator.

NoticethatbasedonthefeaturesofthepermeabilitytensorK,the resulting effectiveparameter Kmaydependonthemacro-scale loca-tion andthesizeof thecoarse-scalepartition. Nevertheless, onecan showthatinpractice,theadaptiverefinementofthemeshisan impor-tantaspectthatimprovethecalculationoftheeffectiveparameters(see

Bastidasetal.(2019)andFig.5).

TheADMprocedureissketchedinFig.6.Moredetailsabouttherole ofthehomogenizationintheoﬄinestageandthecompletealgorithm ofADM-MSandADM-HOcanbefoundinAlgorithm1.

4. Simulationresults

Tobenchmarkthehomogenizationandmultiscalebasedsolutions forthedynamicmeshonheterogeneousmedia,twoheterogeneous non-periodicpermeabilityfieldsfromthetopandbottomlayersoftheSPE 10thComparativeSolutionProjectChristieandBlunt(2001)are con-sidered.Forbothtestcases,thecomputationaldomainentails216×54 gridcellsatfine-scalewithΔ𝑥=Δ𝑦=1[m].Ano-flowconditionis im-posedonallboundaries.Initiallythereservoircontainsonly the2nd phase(e.g.oil),i.e.𝑆=0.The1stphase(e.g.water)isinjectedfroman

(6)

Fig.5.Exampleoffourdiﬀerentlevelsofhomogenizedpermeabilityvalues:ﬁnescale(bottomright),coarselevel1(bottomleft),coarselevel2(topright)and coarselevel3(topleft).

Algorithm 1 The ADM algorithm using multiscale basis functions (ADM-MS)orhomogenization(ADM-HO).

Startofthesimulation;

Readtheinputﬁlesandscanthekeywords;

Givenaﬁnescalepermeability(𝐾)andthenumberofcoarseninglevels (𝐿):

ifMultiscalethen for𝑙=0to𝐿do

ComputethemultiscalebasisfunctionsΦ𝑙_𝑀𝑆;

end else

Homogenizationischosen.

for𝑙=0to𝐿do

Computethehomogenized𝐊𝑙_;

ComputetheconstantbasisfunctionsΦ𝑙_{𝐶𝑜𝑛𝑠𝑡};

end end

fortimestep𝑡𝑛_do

SelectADMgridresolution;

BuildADMprolongationandrestrictionoperators; Take𝑖𝑡𝑒𝑟=1anduseinitialpressureandsaturation;

whileerror≥tolerance&notconvergeddo

Assembleﬁnescalesystem; SolvetheADMsystem;

Prolongsolutionbacktoﬁnescale; Updateproperties;

iferror≤tolerancethen

Converged. else Notconverged. end Nextiteration𝑖=𝑖+1 end

Nexttimestep(𝑡=𝑡+Δ𝑡)

end

Table1

Inputparametersofﬂuidandrockproperties.

Property value

Porosity ( 𝜙) 0.2

Water density ( 𝜌w ) 1000 [Kg/m 3 ]

Oil density ( 𝜌o ) 1000 [Kg/m 3 ]

Water viscosity ( 𝜇w ) 10 −3 [Pa · s]

Oil viscosity ( 𝜇o ) 10 −3 [Pa · s]

Initial pressure ( p 0 ) 10 7 [Pa] Connate water saturation ( S wc ) 0 [-]

Residual oil saturation ( S or ) 0 [-]

Injection pressure ( p inj ) 2 × 10 7 [Pa]

Production pressure ( p prod ) 0 [Pa]

injectionwell,whilethereservoirﬂuidisproducedfromtheproduction well.Thelocationsoftheinjectionandproductionwellsarespeciﬁed ineachtestcase.

Table1showstheinputparametersoftheﬂuidandrockproperties usedinalltestcases.Notethatthedensityandviscosityratiosare as-sumedtobe1.Sincecompressibilityandgravitationalforcesareboth neglected,thedensityvalueshavenoinﬂuenceontheresults.

ThenumericalresultsprovidedbytheADM-MSandADM-HO meth-odsarecomparedtothoseobtainedfromsimulationatﬁnescale (ref-erence).BothADMmethodsemploythecoarseningratioof3×3with twocoarseninglevels.Thisissetaccordingtothesizeofthedomain.

4.1. Testcase1:SPE10toplayer

In thistest case,one injection welland oneproduction well are placed inthebottomleftcornerandtop rightcornerofthedomain, respectively.Thesimulationtimeis𝑡=1000[days]andtheresultsare reportedon100equidistanttimeintervals.Thepermeability distribu-tionoftheSPE10toplayerisshowninFig.7.

Fig.8showsthehomogenizedversionofthepermeabilityattwo dif-ferentlevels.Wehighlightthatthehomogenizedpermeabilityatboth coarselevelspreservesthestructureoftheoriginalﬁne-scale permeabil-ity.Thehighandlowpermeablezonesremainclearlydetectable.

(7)

Fig.6. SchematicdescriptionofADMreservoirsimulation.

Fig.7. Fine-scalepermeability(Log10scale)fromtoplayeroftheSPE10dataset.

Fig.8.HomogenizedpermeabilityofthetoplayeroftheSPE10withcoarsening ratio3.

Fig.9. Saturationproﬁlesat2000days.Thethresholdvalueforthefront track-ingcriterionisΔ𝑆=0.3.

Thesaturationandpressureﬁeldsattheﬁnaltimestepareshown inFig.9andFig.10,respectively.

Fromtheseresults,itisunderstoodthatADM-HOonacoarsecell containinghighandlowpermeablefinecellscanleadtoahigherflow leakage,ascomparedtofine-scaleandADM-MSapproaches.Thiseffect canbeseeninFig.9.Fig.11,illustratingtheadaptivemeshat2000days afterinjection.Noticethattherefinementofthepermeabilityismost dominantatthesaturationfront,duetothechosenmeshrefinement criterion.Forthisfigure, thecoarseningthresholdvalueisΔ𝑆=0.3, i.e.,acellissuccessivelycoarsenedifΔSislowerthan0.3.

TheerrorhistorymapsforbothADM-MSandADM-HOareshown in Fig.12.Therelative errors,presented inFig.12 andFig.14, are expressedintermsoftheL2_norm_over_the_entire_medium,_calculated

withrespecttotheﬁne-scalesolutionas

𝐸𝑟𝑟𝑜𝑟(𝑆)=‖𝑆ref−𝑆ADM‖2 ‖𝑆ref‖2 (13) 𝐸𝑟𝑟𝑜𝑟(𝑃)=‖𝑃ref−𝑃ADM‖2 ‖𝑃ref‖2 . (14)

(8)

Fig.10. Pressureproﬁlesat2000days.Thethresholdvalueforthefront track-ingcriterionisΔ𝑆=0.3.

Fig.11. AdaptivemeshandhomogenizedpermeabilityfortheSPE10toplayer testcase.ThethresholdvalueforthefronttrackingcriterionisΔ𝑆=0.3.

Theresultsindicatethatthehomogenization-basedsimulationshave higher errors compared withthe multiscale-based simulations.They bothhavesimilaraverageusageofactivegridcells,withADM-MS hav-ingslightlyfewergrid cells.Thisis shown inFig.13.Note thatthe gridcellsaroundwellsarekeptattheﬁne-scaleresolutionpermanently. Furthermore,fortightererrortolerancevalues,thequalityofboth ap-proachesbecomecomparable.

Fig.14providestheaveragepressureandsaturationerrorstogether withtheaveragepercentageofactivegridcellsduringthewhole simu-lationtimeasfunctionsofthecoarseningcriterionthreshold.

4.2. Testcase2:SPE10bottomlayer

InthesecondtestcasethepermeabilitydistributionoftheSPE10 bottomlayer,presentedinFig.15,is considered.Thelocationofthe injectionandproductionwellsarethetopleftandthebottomright cor-ners,respectively.Thesimulationtimeis20[days].Allothersimulation parametersremainunchanged.

Fig.16showsthehomogenizedpermeabilityvaluesattwodiﬀerent levels.Duetothemanyhighcontrastchannels,moreactivecellsare employedcomparedwiththeSPEtoplayer,asshowninFig.17.

Thesaturationandpressuremapsattheﬁnaltimestepareshownin

Fig.18andFig.19,respectively.

Similartotheprevioustestcases,Fig.20comparestheerrorbetween thetwoADMapproaches.Moreover,inFig.21,thepercentageofactive gridcellspereachtime-stepisshown.

Fig.12. ComparisonofthesaturationandpressureerrorusingADM-MSand ADM-HOand3diﬀerentvaluesforthefronttrackingcriterion.

Fig.13. ComparisonoftheactivegridcellsusingADM-MSandADM-HOand 3diﬀerentvaluesforthefronttrackingcriterion.

Fig.22illustratestheaveragevalues oftheerrorsin thepressure andthesaturation,andthepercentageoftheactivegridcellsforeach coarseningcriterionthreshold.

Theresults indicateanoticeabledifferenceintheerrorsof ADM-MSandADM-HO.ThepressureerrorinADM-HOissignificantlyhigher sinceADM-HOuseshomogenizedeffectiveparameters.Thisaspectcan be improvedbyemployingfirstordercorrections.However,suchan

(9)

Fig.14. Averageerrorsforthepressureandsaturationandaverageactivegrid cellsforeachstrategy(ADM-MSandADM-HO).

Fig.15. Fine-scalepermeability(Log10scale)frombottomlayeroftheSPE10 testcase.

Fig.16. HomogenizedpermeabilityoftheSPE10bottomlayerwithcoarsening ratio3.

Fig.17. ReﬁnementofthepermeabilityofthebottomlayeroftheSPE10using ADM-HOafter20days.Thethresholdvalueforthefronttrackingcriterionis Δ𝑆=0.3.

Fig.18. Saturationproﬁlesat20days.Thethresholdvalueforthefront track-ingcriterionisΔ𝑆=0.3.

Fig.19. Pressureproﬁlesat20days.Thethresholdvalueforthefronttracking criterionisΔ𝑆=0.3.

approachwoulddeviatefromtheADMframework,andrequiresmore computationaleffort,thereforeitisnotadoptedhere.ADM-MSinstead employsmultiscalebasisfunctions.Duetothemoreaccuratepressure calculations, theADM-MS saturationerror isalso lowerthanthatof ADM-HO.Thedifferenceinthepercentageofactivegridcellsusedin thetwoapproachesislessnoticeablethanthedifferenceintheerrors. However,theADM-HO usesmoreactivegridcells,especiallyinthis SPE10bottomlayertestcase.

(10)

Fig.20. ComparisonofthesaturationandpressureerrorusingADM-MSand ADM-HOand3diﬀerentvaluesforthefronttrackingcriterion.

Fig.21. ComparisonoftheactivegridcellsusingADM-MSandADM-HOand 3diﬀerentvaluesforthefronttrackingcriterion.

Fig.22. Averageerrorsforthepressureandsaturationandaverageactivegrid cellsforbothapproaches(ADM-MSandADM-HO).

5. Conclusion

Homogenizationandmultiscalemethodshavebeendevelopedand evolvedduringthepastdecadeaspromisingadvancedsimulation ap-proachesforlarge-scaleheterogeneoussystems.Inthiswork,thetwo methodswereinvestigated,extendedintoaunifiedfully-implicit frame-work,andbenchmarkedforsimulationofmultiphaseflowinporous me-dia.Itwasshownthatthetwomethodsallowtheconstructionofcoarser levelsystems,andbothrelyonlocalsolutionstofindtheir correspond-ingmaps.Whilehomogenizationmethodsdelivereffectivemodelsand parameters, multiscalemethods findaninterpolation ofthe solution (pressure)across scales.Thisisthemaindifferencebetween thetwo approaches.

Forhighlyheterogeneoustestcases,itwasshownthatthetwo ap-proachesprovideaccuratesolutions.Withthedevelopedmultiscale nu-mericalstrategies,theADM-MSsolutionsaremoreaccuratewhen com-paredtoADM-HO.Theuse ofaconstanteffectiveparameterinstead oflocalmultiscalebasisfunctionsresultsinrelativelyhighererrors.In addition,usingconstantunityprolongationoperatoralongwiththe ef-fectivecoarse-scaleparametersallowsforstraightforward implementa-tionoftheADM-HOmethodfordomainswithnon-periodic permeabil-ityfields.Thestudyofthis papershedsnewlighton theapplication ofmultiscaleandhomogenizationmethodsforreal-fieldsimulationof multiphaseflowinporousmedia.Notethatthecomputationalcostsof thetwoapproacheswerecomparable,astheyappliedalmostthesame activecellsduringthesimulation.Ongoingstudyincludesbenchmark studiesof ADM-HOandADM-MSfor3D fracturedporousmedia,on compilablesimulationplatform,whichallowsscientificCPU compari-sonstudy.

DeclarationofCompetingInterest

Theauthorsdeclarethattheyhavenoknowncompetingﬁnancial interestsorpersonalrelationshipsthatcouldhaveappearedtoinﬂuence theworkreportedinthispaper.

CRediTauthorshipcontributionstatement

HadiHajibeygi:Conceptualization,Supervision,Writing-original draft.ManuelaBastidasOlivares:Methodology,Software,Writing -originaldraft.MousaHosseiniMehr:Methodology,Software,Writing -originaldraft.SorinPop:Conceptualization,Writing-review& edit-ing,Supervision.MaryWheeler:Conceptualization,Writing-review& editing.

(11)

Acknowledgements

Hadi Hajibeygi was sponsored through the Dutch Science Foun-dation (NWO) grant 17509, under Innovational Research Incentives Scheme Vidi (project “ADMIRE”). Sorin Pop and Manuela Bastidas acknowledge the Research Foundation-Flanders (FWO) through the Odysseusprogram(projectGOG1316N).AllauthorsacknowledgeT.U. Delft DARSim group members for the fruitful discussions, specially MatteoCusiniandJeroenRijntjes,fortheirhelpregardingthe DAR-Sim2simulator.DARSim2open-sourcesimulatorcanbe accessedvia

https://gitlab.com/darsim2simulatorlink.

AppendixA. ADMbasedonmultiscalemethod

IntheADM-MSapproach,multiscaleﬁnitevolumemethod(MSFV)

Jennyetal.(2003);CortinovisandJenny(2014)isusedtocompute lo-calbasisfunctionsatmultiplecoarseninglevels.Thecomputationof ba-sisfunctionsΦisdonebysolvingtheincompressibleﬂuidﬂowequation (ellipticpartofthemassbalance)CortinovisandJenny(2017)which reads

−∇⋅(𝜆 ⋅ ∇Φ)=0. (15)

Theincompressiblebasisfunctionsarefoundtobethemosteﬃcient ones,comparedwiththecompressibleandmorecomplexformulations

Teneetal.(2015).Thefirststepistoimposecoarsegridsontopofthe finemesh,forthecoarselevel1.Here,tosimplifythevisualization,a 2D15×15discretedomainisconsidered(seeFig.23).Byconnecting thecentersofthecoarsecells,thedual-coarsegridisobtained.Thedual gridmakesanoverlappingpartitioningofthefine-scaledomain,with3 categoriesofinterior(white),edge(green),andvertex(blue)cells.The coarseningratiointheillustratedexampleofFig.23is5×5.

TheEq.(15)issolvedateachdualcoarsegridhandforeachcoarse node(vertex)k,i.e.,−∇⋅ (𝜆 ⋅ ∇Φℎ_𝑘)=0.Inordertosolvethislocal sys-tem,Dirichletboundaryconditionsof1.0(forthecorrespondingcoarse node)and0.0(fortheotherthreecoarse nodes)areimposed.These Dirichletvaluesallowtosolvethebasisfunctionsontheedges,ifa re-duceddimensional(1D)ellipticproblemisconsidered.Thesolutionat theedgeandvertexcellsarethenimposedasDirichletboundary con-ditionforthefull2Dproblem.Thesolutionofthiswell-posedsystemis thebasisfunctionofthecorrespondingcoarsenodeatthe correspond-ingdualcoarsegrid.Fig.24showsaschematicofthementioneddual coarsegridhandanexampleofabasisfunctionbelongingtothebottom leftcoarsenode(Φℎ₁).

Fig.25showsall thefourbasisfunctionsforthementioned dual coarsegridh.

Thecombinationofthebasisfunctionsatallthedualcoarse grid cellssurroundingthecorrespondingcoarsenodeformsthebasis func-tionbelongingtothatcoarsenode.Fig.26illustratesanexampleofa

Fig.23. Constructionofthecoarseanddual-coarsegridsontheﬁne-scale dis-cretedomain.Finecellsarepartitionedw.r.t.thedualcoarsemeshas:interior, edgeandvertexcells.

Fig.24. Illustrationofadualcoarsegridandabasisfunctionbelongingtothe bottomleftcoarsenode.Asitcanbeseen,thevalueofthebottomleftcoarse nodeissettobe1.0,whiletheotherthreevertexcellsaresetto0.

Fig.25. Allthefourbasisfunctionsbelongingtothedualcoarsegridh.Shown beloweachplotistheDirichletvalueatthecornerofeachdualcoarsecellfor theplottedbasisfunction.

Fig.26. Anexampleofabasisfunctionbelongingtothebottomleftcoarse nodeofaheterogeneousdomainwith27×27gridcells.Thecoarseningratio hereis9×9.

basisfunctionbelongingtothebottomleftcoarsenodeofanexample heterogeneous2Ddomai.

Toobtainthebasisfunctionsathighercoarseninglevels,the hierar-chicallynestedcoarsegridisconstructedonthesamedomain.Thesame procedureisfollowedtocomputethebasisfunctionsathigher coarsen-inglevels.Fig.27showsthecoarsegridconstructionat2consequent coarseninglevels,fora2Ddomainwith75×75ﬁnecells.

Notethat,accordingtothevastmultiscaleliterature,construction ofbasisfunctionscanbedonepurelyalgebraic,oncethewire-basket decompositionoftheﬁnecellsintoVertex,Edge,Face,andInterioris known Teneetal. (2016). Apartitioning methodshould be applied for complex mesh MøynerandLie (2016); Parramoreet al.(2016);

Shah et al. (2016); Gulbransen et al. (2010); Bosma et al. (2017);

Mehrdoost(2019);MehrdoostandBahrainian(2016).

AppendixB. ADMbasedonhomogenizationtheory

ThemainideaoftheADM-HOistouseahomogenizedversionof thepermeabilityKinsteadofvolume-averagedpermeabilitiesat

(12)

diﬀer-Fig.27. Coarsegridconstructionat2consequentcoarseninglevelsfora2D domainwith75ﬁnecells.Thecoarseningratioof5×5ischosen.Thegridsizes ofcoarselevel1and2are15×15and3×3,respectively.

entcoarselevels.Indoingso,twoassumptionsarecommonlymade:the permeabilityisperiodicateachofthecoarseninglevels,andthescales arewellseparated.WerefertoAllaire(1992)fortherigorous mathemat-icalsupportofthisapproach.Althoughthetestcasesconsideredheredo notsatisfythetwoassumptionsstatedbefore,thehomogenizationidea canstillbeconsideredfordevelopingmultiscalesimulationtools,and inthissensewerefertoBastidasetal.(2019);Amanbeketal.(2019a);

Singhetal.(2019a);Amanbeketal.(2019b);Singhetal.(2019b). Ateachcoarseninglevell,wecallmicro-scalecell(i.e.,localcoarse cell)theregionΩ_lwhereintheparameterschangerapidly.ForeachΩ_l, thecharacteristiclengthis𝓁,whereListhecharacteristiclengthforthe macro-scaledomainΩ.Thefactor𝜀∶= 𝓁_𝐿 reflectsthescaleseparation. Toidentifythefastchangesintheparameterswedoublethevariables anddefinethefastvariabley∶= x_𝜀.Inthenon-dimensionalsetting,Ω canbewrittenasthefiniteunionofthelocalcellsΩl.WeletΩ =∪Ω_𝑙 forsomesetofindices.

Forcalculatingthehomogenizedpermeabilityweconsideran auxil-iaryellipticproblem:Auxiliaryproblem.Givenafine-scale permeabil-ityK,findafunctionu𝜖thatsatisfies

∇⋅ (𝐾⋅ ∇𝑢𝜖)=0 inΩ,

𝑢𝜖₌₀ _on_𝜕Ω. ₍₁₆₎

Heretheboundaryconditionsarespeciﬁedforcompleteness.

Weemploythehomogenizationansatz,meaningthattheunknownu𝜖

canbewrittenas

𝑢𝜖₍_𝑥₎₌_̂𝑢

0+𝜖 ̂𝑢1+𝜖2̂𝑢2+…, (17)

whereeacĥ𝑢𝑖(x,y)isperiodic.

Duetothedoublingofvariables,thegradientanddivergence oper-atorsbecome

∇=∇_𝑥+1

𝜖∇𝑦 and div=div𝑥+1_𝜖div𝑦.

Inserting(17)andthetwoscaleoperatorsaboveintheauxiliaryproblem

(16),onegets ( div_𝑥+1 𝜖div𝑦 )(( ∇_𝑥+1 𝜖∇𝑦 )( ̂𝑢0+𝜖 ̂𝑢1+(𝜖2) )) =0, ̂𝑢0+𝜖 ̂𝑢1+(𝜖2)=0.

Collectingthetermswithfactor𝜖−2_we_obtain_for_each_domain_Ω

l ∇_𝑦⋅(𝐾∇̂𝑢0

)

=0forally∈ Ω_𝑙, (18)

Clearly,any𝑢0=𝑢0(x)which doesnotdependonthefastvariabley

isasolutionoftheproblem(18).Thus,onecanprovethatallsolutions dependonlyonx.Thefunctionu₀isinfactthemacro-scaleapproximation

oftheoriginalunknownu𝜖.

Ontheotherhand,forthe𝜖−1_terms_one_has

∇_𝑦⋅(𝐾(∇_𝑥̂𝑢0+∇𝑦̂𝑢1

))

=0forally∈𝑌.

Onecandetermine ̂𝑢1asafunctionof ̂𝑢0andeliminateitfromthe

system.Tothisend, ̂𝑢1 iswrittenasalinearcombinationoffunctions

𝜔j_,_and_with𝜕̂𝑢0 𝜕𝑥𝑗 ascoeﬃcients ̂𝑢1(x,y)= dim ∑ 𝑗=1 𝜕̂𝑢0(x) 𝜕𝑥𝑗 𝜔 𝑗₍_x_,_y₎_.

Thefunctions𝜔j_are_solutions_of_the_micro-cell_(local_domain) prob-lemsdeﬁnedoneachΩ_l:

−∇_𝑦⋅ ( A(∇𝑤𝑗_𝑦+⃗𝑒_𝑗 )) =0 forally∈𝑌, 𝑤𝑗_isperiodicin_𝑌_.

Here,{⃗𝑒_𝑗}d_𝑗=1isthecanonicalbasisofdimension2.Toguarantee the uniquenessofthesolutiononerequiresthat

∫Ω_𝑙𝜔

𝑗₌₀_𝑑_y_,_for_all_Ω 𝑙.

Thehomogenizedpermeabilityintheauxiliaryproblemisobtained byconsideringthetermsoforder 𝜖0_,_in_which _̂𝑢

2 appears.Averaging

thesetermsandusingtheperiodicityofthefunctionsontherighthand sideof (17),oneobtains thattheauxiliaryunknown ̂𝑢0(x)solvesthe

homogenizedproblem ∇⋅ (𝐊𝑙⋅ ∇̂𝑢0)=0 inΩ,

̂𝑢0=0 on𝜕Ω. (19)

Herethematrixvaluedfunction𝐊𝑙_∶_Ω_{→ ℝ}2×2_has_the_elements

𝐊𝑙 𝑖,𝑗||_|Ω_𝑙=∫_Ω

𝑙

(

𝐾(⃗𝑒𝑗+∇𝑦𝜔𝑗))⋅ ⃗𝑒𝑖𝑑y.

Notethatthesestepsarecarriedoutonlytodeterminetheeﬀective permeabilitytensorKl_._The_solution_̂𝑢

0oftheeﬀectiveproblem(19)is

notofinteresthere.Inotherwords,weusetheauxiliaryproblemfor thesolepurposeofdeﬁninganeﬀectiveparameterthatcouldbeused ateachlevelofthedynamicmultilevelalgorithm.

Supplementarymaterial

Supplementarymaterialassociatedwiththisarticlecanbefound,in theonlineversion,atdoi:10.1016/j.advwatres.2020.103674.

References

Aarnes, J. , Hou, T.Y. , 2002. Multiscale domain decomposition methods for elliptic problems with high aspect ratios. Acta Math. Appl. Sinica 18 (1), 63–76 .

Abdulle, A. , Nonnenmacher, A. , 2009. A short and versatile ﬁnite element multiscale code for homogenization problems. Comput. Method. Appl. Mech. Eng. 198, 2839–2859 . Abdulle, A. , Weinan, E. , 2003. Finite diﬀerence heterogeneous multi-scale method for

homogenization problems. J. Comput. Phys. 191 (1), 18–39 .

Abdulle, A. , Weinan, E. , Engquist, B. , Vanden-Eijnden, E. , 2012. The heterogeneous multiscale method. Acta Numer. 21, 1–87 .

Allaire, G. , 1992. Homogenization and two-scale convergence. SIAM J. Math. Anal. 23 (6), 1482–1518 .

Amanbek, Y. , Singh, G. , Wheeler, M.F. , van Duijn, H. , 2019. Adaptive numerical homogenization for upscaling single phase flow and transport. J. Comput. Phys. 387, 117–133 . Amanbek, Y. , Singh, G. , Wheeler, M.F. , van Duijn, H. , 2019. Adaptive numerical homogenization for upscaling single phase flow and transport. J. Comput. Phys. 387, 117–133 . Amaziane, B. , Bourgeat, A. , Koebbe, J. , 1991. Numerical simulation and homogenization of two-phase flow in heterogeneous porous media. In: Mathematical Modeling for Flow and Transport Through Porous Media. Springer, pp. 519–547 .

Amaziane, B., Pankratov, L., Piatnitski, A., 2017. An improved homogenization result for immiscible compressible two phase ﬂow in porous media. Netw. Heterogen. Media 12 (1), 147–171. https://doi.org/10.3934/nhm.2017006 .

Arbogast, T. , Xiao, H. , 2013. A multiscale mortar mixed space based on homogenization for heterogeneous elliptic problems. SIAM J. Numer. Anal. 51 (1), 377–399 . Aziz, K. , Settari, A. , 2002. Petroleum Reservoir Simulation. Blitzprint Ltd., Cagary, Alberta . Bastidas, M. , Bringedal, C. , Pop, S. , Radu, F. , 2019. Adaptive numerical homogenization

of non-linear diﬀusion problems. arXiv preprint arXiv:1904.10665 .

Berger, M. , Oliger, J. , 1984. Adaptive mesh reﬁnement for hyperbolic partial diﬀerential equations. J. Comput. Phys. 53, 484–512 .

Bosma, S., Hajibeygi, H., Tene, M., Tchelepi, H.A., 2017. Multiscale ﬁnite volume method for discrete fracture modeling on unstructured grids (ms-dfm). J. Comput. Phys. 351, 145–164. https://doi.org/10.1016/j.jcp.2017.09.032 .

(13)

Bourgeat, A., Luckhaus, S., Mikelic, A., 1996. Convergence of the homogenization process for a double-porosity model of immiscible two-phase ﬂow. SIAM J. Math. Anal. 27 (6), 1520–1543. https://doi.org/10.1137/S0036141094276457 .

Brown, D.L. , Efendiev, Y. , Hoang, V.H. , 2013. An eﬃcient hierarchical multiscale ﬁnite element method for stokes equations in slowly varying media. Multiscale Model. Sim- ulat. 11 (1), 30–58 .

Carciopolo, L.D., Cusini, M., Formaggia, L., Hajibeygi, H., 2020. Adaptive multilevel space- time-stepping scheme for transport in heterogeneous porous media (adm-lts). J. Com- put. Phys. 6, 100052. https://doi.org/10.1016/j.jcpx.2020.100052 .

Celia, M.A. , Bouloutas, E.T. , Zarba, R.L. , 1990. A general mass-conservative numerical solution for the unsaturated ﬂow equation. Water Resour. Res. 26 (7), 1483–1496 . Christie, M.A. , Blunt, M. , 2001. Tenth SPE comparative solution project: a comparison

of upscaling techniques. SPE reservoir simulation symposium. Society of Petroleum Engineers .

Chung, E.T. , Efendiev, Y. , Lee, C.S. , 2015. Mixed generalized multiscale ﬁnite element methods and applications. Multiscale Model. Simulat. 13 (1), 338–366 .

Cortinovis, D. , Jenny, P. , 2014. Iterative galerkin-enriched multiscale ﬁnite-volume method. J. Comp. Phys. 277, 248–267 .

Cortinovis, D. , Jenny, P. , 2017. Zonal multiscale ﬁnite-volume framework. J. Comp. Phys. 337, 84–97 .

Cusini, M. , Fryer, B. , van Kruijsdijk, C. , Hajibeygi, H. , 2018. Algebraic dynamic multilevel method for compositional ﬂow in heterogeneous porous media. J. Comput. Phys. 354, 593–612 .

Cusini, M. , Gielisse, R. , Groot, H. , van Kruijsdijk, C. , Hajibeygi, H. , 2019. Incomplete mixing in porous media: todd-longstaﬀ upscaling approach versus a dynamic local grid reﬁnement method. Comput. Geosci. 23 (2), 373–397 .

Cusini, M. , van Kruijsdijk, C. , Hajibeygi, H. , 2016. Algebraic dynamic multilevel (adm) method for fully implicit simulations of multiphase ﬂow in porous media. J. Comput. Phys. 314, 60–79 .

van Duijn, C. , Eichel, H. , Helmig, R. , Pop, I. , 2007. Effective equations for two-phase flow in porous media: the effect of trapping at the micro-scale. Transp. Porous. Media 9, 411–428 .

Edwards, M.G. , 1996. A higher-order godunov scheme coupled with dynamic local grid reﬁnement for ﬂow in a porous medium. Comput. Method. Appl. Mech. Eng. 131 (3–4), 287–308 .

Efendiev, Y. , Hou, T.Y. , 2009. Multiscale Finite Element Methods: Theory and Applica- tions, 4. Springer Science & Business Media .

Faigle, B. , Helmig, R. , Aavatsmark, I. , Flemisch, B. , 2014. Eﬃcient multiphysics modelling with adaptive grid reﬁnement using a mpfa method. Comput. Geosci. 18 (5), 625–636 .

Gulbransen, A.F., Hauge, V.L., Lie, K.-A., 2010. A multiscale mixed ﬁnite element method for vuggy and naturally fractured reservoirs. SPE J. 15 (02), 395–403.

https://doi.org/10.2118/119104-PA .

Hajibeygi, H. , Bonﬁgli, G. , Hesse, M.A. , Jenny, P. , 2008. Iterative multiscale ﬁnite-volume method. J. Comput. Phys. 227 (19), 8604–8621 .

Hajibeygi, H. , Lee, S.H. , Lunati, I. , 2012. Accurate and efficient simulation of multiphase flow in a heterogeneous reservoir with error estimate and control in the multiscale finite-volume framework. SPE J. 17 (04), 1–071 .

Henning, P., Ohlberger, M., Schweizer, B., 2013. Homogenization of the degenerate two-phase ﬂow equations. Math. Model. Method. Appl. Sci. 23 (12), 2323–2352.

https://doi.org/10.1142/S0218202513500334 .

Henning, P. , Ohlberger, M. , Schweizer, B. , 2015. Adaptive heterogeneous multiscale methods for immiscible two-phase ﬂow in porous media. Comput. Geosci. 19 (1), 99–114 . Hornung, U. , 1997. Homogenization and porous media, 6. Springer Science & Business

Media .

HosseiniMehr, M. , Cusini, M. , Vuik, C. , Hajibeygi, H. , 2018. Algebraic dynamic multilevel method for embedded discrete fracture model (f-adm). J. Comput. Phys. 373, 324–345 .

Hou, T.Y. , Wu, X.-H. , 1997. A multiscale ﬁnite element method for elliptic problems in composite materials and porous media. J. Comput. Phys. 134 (1), 169–189 .

Jansen, J.-D., Brouwer, R., Douma, S.G., 2009. Closed loop reservoir manage- ment. SPE Reservoir Simulation Symposium. Society of Petroleum Engineers.

https://doi.org/10.2118/119098-MS .

Jenny, P. , Lee, S. , Tchelepi, H.A. , 2003. Multi-scale ﬁnite-volume method for elliptic problems in subsurface ﬂow simulation. J. Comput. Phys. 187 (1), 47–67 .

Klemetsdal, Ø.S., Lie, K.-A., 2020. Dynamic coarsening and local reordered nonlinear solvers for simulating transport in porous media. SPE J. Preprint (Preprint), 20.

https://doi.org/10.2118/201089-PA .

Klemetsdal, Ø.S. , Møyner, O. , Lie, K.-A. , 2020. Accelerating multiscale simulation of complex geomodels by use of dynamically adapted basis functions. Comput. Geosci. 24 (2), 459 —–476 .

Li, H., Leung, W.T., Wheeler, M.F., 2020. Sequential local mesh refinement solver with separate temporal and spatial adaptivity for non-linear two-phase flow problems. J. Comput. Phys. 403, 109074. https://doi.org/10.1016/j.jcp.2019.109074 . Mehrdoost, Z., 2019. Unstructured grid adaptation for multiscale finite volume method.

Comput. Geosci. 23 (6), 1293–1316. https://doi.org/10.1007/s10596-019-09878-9 .

Mehrdoost, Z. , Bahrainian, S.S. , 2016. A multilevel tabu search algorithm for balanced partitioning of unstructured grids. International Journal for Numerical Methods in Engineering 105 (9), 678–692 . 10.1002/nme.5003

Møyner, O. , Lie, K.-A. , 2016. A multiscale restriction-smoothed basis method for high contrast porous media represented on unstructured grids. J. Comput. Phys. 304, 46–71 . Parramore, E. , Edwards, M.G. , Pal, M. , Lamine, S. , 2016. Multiscale ﬁnite-volume

cvd-mpfa formulations on structured and unstructured grids. SIAM Multiscale Model. Simul. 14, 559–594 .

Pau, G.S. , Almgren, A.S. , Bell, J.B. , Lijewski, M.J. , 2009. A parallel second-order adaptive mesh algorithm for incompressible ﬂow in porous media. Philos. Trans. R. Soc. 367 (1907), 4633–4654 .

Pau, G.S.H. , Bell, J.B. , Almgren, A.S. , Fagnan, K.M. , Lijewski, M.J. , 2012. An adaptive mesh reﬁnement algorithm for compressible two-phase ﬂow in porous media. Com- put. Geosci. 16 (3), 577–592 .

Radu, F.A. , Kumar, K. , Nordbotten, J.M. , Pop, I.S. , 2017. A robust, mass conservative scheme for two-phase ﬂow in porous media including hölder continuous nonlineari- ties. IMA J. Numer. Anal. 38 (2), 884–920 .

Sammon, P.H. , 2003. Dynamic grid reﬁnement and amalgamation for compositional simulation. SPE reservoir simulation symposium. Society of Petroleum Engineers . Schmidt, G. , Jacobs, F. , 1988. Adaptive local grid reﬁnement and multi-grid in numerical

reservoir simulation. J. Comput. Phys. 77 (1), 140–165 .

Shah, S. , Moyner, O. , Tene, M. , Lie, K.-A. , Hajibeygi, H. , 2016. The multiscale restriction smoothed basis method for fractured porous media (f-msrsb). J. Comput. Phys. 318, 36–57 .

Singh, G. , Leung, W. , Wheeler, M.F. , 2019. Multiscale methods for model order reduction of non-linear multiphase ﬂow problems. Comput. Geosci. 23 (2), 305–323 . Singh, G. , Leung, W. , Wheeler, M.F. , 2019. Multiscale methods for model order reduction

of non-linear multiphase ﬂow problems. Comput. Geosci. 23, 305–323 .

Singh, G., Wheeler, M.F., 2018. A space time domain decomposition approach using enhanced velocity mixed ﬁnite element method. J. Comput. Phys. 374, 893–911.

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

Szymkiewicz, A. , Helmig, R. , Kuhnke, H. , 2011. Two-phase ﬂow in heterogeneous porous media with non-wetting phase trapping. Transp. Porous Media 86, 27–47 . Tene, M. , Al Kobaisi, M.S. , Hajibeygi, H. , 2016. Algebraic multiscale method for ﬂow in

heterogeneous porous media with embedded discrete fractures (f-ams). J. Comput. Phys. 321, 819–845 .

Tene, M. , Wang, Y. , Hajibeygi, H. , 2015. Adaptive algebraic multiscale solver for compressible ﬂow in heterogeneous porous media. J. Comput. Phys. 300, 679–694 .

Vasilyeva, M., Leung, W.T., Chung, E.T., Efendiev, Y., Wheeler, M., 2020. Learn- ing macroscopic parameters in nonlinear multiscale simulations using non- local multicontinua upscaling techniques. J. Comput. Phys. 412, 109323.

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

Wachspress, E.L. , 1966. Iterative Solution of Elliptic Systems: and Applications to the Neutron Diﬀusion Equations of Reactor Physics. Prentice-Hall .

Wang, Y. , Hajibeygi, H. , Tchelepi, H.A. , 2014. Algebraic multiscale solver for ﬂow in heterogeneous porous media. J. Comput. Phys. 259, 284–303 .

Weinan, E. , 2011. Principles of Multiscale Modeling. Cambridge University Press . Weinan, E. , Engquist, B. , Li, X. , Ren, W. , Vanden-Eijnden, E. , 2007. Heterogeneous multi-

scale methods: a review. Commun. Comput. Phys. 2 (3), 367–450 .

Weinan, E. , Yue, X.Y. , 2004. Heterogenous multiscale method for locally self-similar problems. Commun. Math. Sci. 2 (1), 137–144 .

Zhou, H. , Tchelepi, H.A. , 2012. Two-stage algebraic multiscale linear solver for highly heterogeneous reservoir models. SPE J. 17 (02), 523–539 .