Adaptive multilevel space-time-stepping scheme for transport in heterogeneous porous media (ADM-LTS)

Delft University of Technology

Adaptive multilevel space-time-stepping scheme for transport in heterogeneous porous

media (ADM-LTS)

Carciopolo, Ludovica Delpopolo; Cusini, Matteo; Formaggia, Luca; Hajibeygi, Hadi

DOI

10.1016/j.jcpx.2020.100052

Publication date

2020

Document Version

Final published version

Published in

Journal of Computational Physics: X

Citation (APA)

Carciopolo, L. D., Cusini, M., Formaggia, L., & Hajibeygi, H. (2020). Adaptive multilevel space-time-stepping

scheme for transport in heterogeneous porous media (ADM-LTS). Journal of Computational Physics: X, 6,

[100052]. https://doi.org/10.1016/j.jcpx.2020.100052

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

multilevel

space-time-stepping

scheme

for

transport

in

heterogeneous

porous

media

(ADM-LTS)

Ludovica Delpopolo Carciopolo

a

,

∗

,

Matteo Cusini

b

,

Luca Formaggia

a

,

Hadi Hajibeygi

c

a_{MOX, Dipartimento di Matematica, Politecnico di Milano, Via Bonardi 9, 20133 Milano, Italy}

b_{Atmospheric, Earth, & Energy Division, Lawrence Livermore National Laboratory, 7000 East Ave., Livermore, CA 94550-9234, United States of} America

c_{Department of Geoscience and Engineering, Faculty of Civil Engineering and Geosciences, Delft University of Technology, Stevinweg 1, 2628} CN, Delft, the Netherlands

a

r

t

i

c

l

e

i

n

f

o

a

b

s

t

r

a

c

t

Article history:

Received 11 June 2019

Received in revised form 10 February 2020 Accepted 19 February 2020

Available online 3 March 2020 Keywords:

Local time-stepping strategies Conservative multirate methods Algebraic multilevel methods Multiphase flow

Porous media

We present ADM-LTS, an adaptive multilevel space-time-stepping schemefor transport in heterogeneous porous media. At each time step, firstly, the flow (pressure) solution is obtained. Then, the transport equationis solved using the ADM-LTS method, which consists oftwostages.Inthe firststage,an initialsolutionis obtainedbyimposingthe coarsest space-time grid. Thisinitialsolution is thenimproved, inthe secondstage,by imposing aspace-time adaptive grid on the cells where the solution does not satisfy the desired quality. The quality control is based onerror estimators with user-defined threshold values. The time-integration procedure, in which the coarsest-scale solution provideslocal flux boundary conditions forsub-domains with local timerefinement, is strictlymassconservative.Inaddition,themethodemploysspace-timefinegridcellsonly atthemoving saturationfronts.Inordertoensurelocal massconservationatalllevels, finite-volumerestrictionoperatorsandunityprolongationoperatorsaredeveloped.Several numerical experiments have been performed to analyzethe efficiency and accuracy of the proposed ADM-LTS method for bothhomogeneous and heterogeneous permeability fieldsontwoandthreedimensionaldomains.Theresultsshowthatthemethodprovides accuratesolutions,atthesametimeitmaintainsthecomputationalefficiency.The ADM-LTSimplementationispubliclyavailableathttps://gitlab.com/darsim2simulator.

©2020TheAuthors.PublishedbyElsevierInc.ThisisanopenaccessarticleundertheCC BY-NC-NDlicense(http://creativecommons.org/licenses/by-nc-nd/4.0/).

1. Introduction

Simulationoftransportinnaturalporousmediaischallenging duetothevarietyoftimeandlengthscalesinvolvedin theprocess. Infact, geologicalformations extendseveralhundreds ofmeterswhereas physicalandchemical phenomena, whicharerelevantinmanyapplications(e.g.,energyandgreenhousegasstorageandgeo-energyproduction),occuratmuch smallerscales(cmandbelow).Additionally,atthecontinuum(orDarcy)scale,porousmediapresenthighlyheterogeneous discreteproperties(e.g.,permeability)withnoseparationofscales.Thisleads tohavingfastandslowtransportprocesses

*

Corresponding author.

E-mail addresses:ludovica.delpopolo@polimi.it(L. Delpopolo Carciopolo), cusini1@llnl.gov(M. Cusini), luca.formaggia@polimi.it(L. Formaggia), H.Hajibeygi@tudelft.nl(H. Hajibeygi).

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

2590-0552/©2020 The Authors. Published by Elsevier Inc. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).

coexistinonecomputationaldomain.Accuratenumericalmodels,thus,requireveryhighresolutiongridsbothinspaceand time tocaptureall relevantphysics.However,duetothebigsizeofthedomains andthelargenumberofsimulationsfor uncertaintyreduction[1,2],field-scalesimulationsonsuchhighresolutiongridsareimpractical.

The globalsystemcomplexitycannotbe reducedbyupscaling approaches,astheheterogeneouspropertiesoftentimes donotentailclearseparationofscales[3].Alternatively,onecanmapscalablepartofthesolution(e.g.pressure)incoarser spatial resolutions,asinthemultiscalemethods[4–8].Forthetransport (e.g.saturation)equation,duetoitslocalnature, Adaptive MeshRefinement (AMR) schemes providea promising framework forthe reduction ofthe computational com-plexity[9,10].Bothapproachesaddressresolutionchallengesinspaceonly,whichcanbeintegratedintoanoperator-based strategybytheAlgebraicDynamicMultilevel(ADM)method[11,12].

ADMmapsthespatialfine-scalesystemintoadynamicmultilevelsystembyusingsequencesoflocaloperators (restric-tion andprolongation).Itcan beappliedtocoupledflow-transport, oronlytransport equation.Inadditiontoexplicitand implicittimeintegrationschemes,bothsequentialandfully-implicitcouplingstrategiescanbetreatedinastraightforward manner. Thefocusofthepresentstudy isonthetransportequation,whichissequentiallycoupledwiththepressureand velocitysolverateachtimestep.

SimilartoAMR,ADMdevelopsanadaptivespatial-gridsteppingscheme.IncontrasttoAMR,however,ADMisapplicable toheterogeneoussystemswithoutdependencyonupscalingparameters.

Althoughimportant,spatialgridisonlypartofthefulltransportsimulationadvancement.Tofurtheradvancetheoverall simulation efficiency,one mayneedtotakelarge timestepsizes.However, theexcessivenumericaldispersionintroduced by theuseoflargetime stepscan significantlyimpacttheaccuracy ofthe solution,e.g., bysmearing theadvancing satu-rationfront.Itisthereforecruciallyimportantto,notonlyallowforlargetime stepsizes,butalsotocontroltheaccuracy by resolving the sharp fronts using smaller time steps. Such a development is made possible–for static spatial grids–by Local TimeStepping (LTS)ormultiratemethods[13].RecentLTS methodsallowforself-adjustingtime-steppingstrategies based ona posteriorierrorestimator[14] ina strictly conservativeimplicittime integrationframework [15]. TheLTS for conservativetransportsimulationonfine-scalespatialgridhasbeenalsosuccessfullyintegratedwiththepressuresolverof themultiscalefinitevolumemethod[16].ThisideawasalsofollowedtodevelopanexplicitLTSscheme[17].Notethatall ofthesementionedmethodsdependonagivenstaticsingle-levelspatialgridfortransport.

Ofparticularinterestistodevelopan adaptivespace-timesteppingschemefortransport simulation.Thisdevelopment wouldallow notonlyto adjustthe spatial,butalsotemporalgridinorderto achievetheoptimum efficientandaccurate transportsolutions.

Inthispaper,theadaptivespace-timesteppingscheme(ADM-LTS)isdevised.Thisnewdevelopmentmakes itpossible tomarchadaptivelybothinspaceandtime(withanarbitrarynumberoflevels)toobtainaccurateandefficientsolutions forcomplextransportinheterogeneousporousmedia.Additionally,theproposedADM-LTSmethodpreservesmass conser-vation,eventhoughittakesdifferenttimestepsizesindifferentsubsetsofthedomain.Thisismadepossiblebyobtaining an initialsolutionestimatoronthecoarsestspace-timegrid,andimposeitasNeumann(fixedflux)boundaryconditionto localregionswithsmallerspace-timesteps.

NotethattheADM-LTSisdevelopedforthetransportequation,withinthecoupledflow-transportmultiphasesimulation. Inparticular,itisemphasizedthatthepressureequationsolutionisassumedtobeavailablebyanyadvancedsolver,which staysoutofthescopeofthispaper.Oncethepressureandvelocityfieldsareobtained,ADM-LTSisimposedonthetransport equation,withadaptivespace-timegridbasedonspace-timeerrorestimators.Asaconsequence,themethodisabletouse a fine-gridresolution(in spaceandtime)onlyatthelocationofthemoving saturationfronts.Moreprecisely,themethod firstappliesa globaltimesteponthecoarsestpossiblegridresolution.Then,itdefinesamultilevelgridresolutiononthe basisontheerrorestimators.Onthisnewgrid,thesolutionisrecomputedwithsmallertimestepsonlyinafractionofthe domain.

ADM-LTS is different, yet complementary, to other advanced methods for transport equation, specially the Adaptive Implicit Method (AIM) [18–20]. AIM develops an adaptive implicit-explicit time-integration scheme, but takes the same space-time gridseverywhere. ADM-LTS,ontheother hand,is developedtoallowdifferentspaceandtimestep sizesina domain,irrespectiveofthechoiceofthetime-integrationscheme.

TheADM-LTSmethodisappliedtohomogeneousandheterogeneous2Dand3Dtestcasesincludingcomplexnonlinear transportphysics(i.e.,gravitationalforcesandcapillarityheterogeneity).Systematicstudiesoftheperformance(i.e.,accuracy andsystemcomplexity)arepresented.NumericaltestcasesshowthattheADM-LTSapproachprovidesanaccuratesolution reducingthenumberofactivecellsbothinspaceandintime.

Thepaperisorganizedasfollows.Theequationsdescribingmultiphaseflowinporousmediaarepresentedinsection2 along withthefine-scale discretesystemsfortransportequation.TheADM-LTSmethodisexplainedindetailinsection 3 whereasnumericalexperimentsarepresentedinsection4.Finally,conclusionsareprovidedinsection5.

2. Sequential-implicitsimulation

MassconservationequationsfortheflowofNp incompressiblephasesinad-dimensionalporousdomain



⊂ R

d read

∂

Here, Sα,

ρ

α andqα arethe

α

-phasesaturation,densityandsourceterm(i.e.,wells).The Darcyvelocityfor

α

-phasecan bewrittenas

uα

= −K

λ

α

(

∇

pα

+

ρ

αg

∇

z

)

∀

α

∈ {

1

, . . . ,

Np

}.

(2)

Here,the phasemobilityisdenoted by

λ

α ,and

λ

α

=

krα

/

μ

α holds,wherekrα and

μ

α arethephase relative permeabil-ityandviscosity,respectively. Furthermore,K isthe rockabsolutepermeability(definedatdiscrete fine-scaleresolution). Additionally, pα arethephasepressures,whicharerelatedtothecapillarypressure Pcas

pα

−

pβ

= (

1

− δ

α,β

)

Pcα,β

∀

α

, β

∈ {

1

, . . . ,

Np

}.

(3)

Here,

δ

α,β is the Kronecker delta, which is equal to 1 if

α

= β

and0 otherwise, while Pcα,β is a nonlinear function of wetting phase saturation. Notethat the constraint



_αN₌p₁Sα

=

1 allows foreliminating one of the unknowns. Assuming incompressiblefluids androck, ina two-phase systemone can use thenon wetting (nw) pressureand the wetting (w) saturationasprimaryvariables[21,22].Wewillthensetp

=

pnw

,

S

=

Sw anddefineatotalvelocityas

ut

=

unw

+

uw

.

(4)

Sequential-implicitsimulation requiresadecoupledpressure-saturationsets ofequations.Theequation forpressureis ob-tainedbycombiningthemass-balanceequations,tohave

∇ ·

ut

=

qw

ρ

w

+

qnw

ρ

nw in

,

(5)

wherethetotalvelocityisexpressedas ut

=

K



−λt

∇

p

+ λw∇

Pc

− (λw

ρ

w

+ λnw

ρ

nw

)

g

∇

z



.

(6)

Ontheotherhand,the(wetting-phase)saturationequation isre-writtenbyintroducingthefractional flowfunction fw

=

λ

w

/λ

t,thusobtaining

∂(φ

S

)

∂

t

+ ∇ ·



K fw

λ

nw

(

∇

Pc

+ (

ρ

nw

−

ρ

w

)

g

∇

z

)

+

fwut



=

qw

ρ

w in

.

(7)

Equations (5), (6) and (7) are coupled by the total velocity, the (nonlinear) phase relative permeability krα and the capillarypressure Pc.

Sequentialimplicitsimulation(SIM)consistsofdecouplingthepressureandtransportequationsatnumericallevel,and solving each ofthem implicitly intime. Giventhe state ata current time tn_{,the solution attime t}n+1 _{isfound by first}

solvingEq.(5),freezingallsaturationdependencies,i.e.,

−∇ ·



K

λ

nt

∇

pn+1



=

qt

− ∇ ·



K

(λ

nw∇Pn_c

− (λ

n_w

ρ

w

+ λ

nnw

ρ

nw

)

g

∇

z

)



,

(8)

whereqt

=

qw

/

ρ

w

+

qnw

/

ρ

nw.Then,thetotalvelocityiscomputedas

un_t+1

= −

K



λ

n_t

∇

pn+1

− λ

nw∇P_cn

+ (λ

n_w

ρ

w

+ λ

nnw

ρ

nw

)

g

∇

z



.

(9)

Finallythetransportequationissolved,i.e.,

φ

S n+1

₋

_Sn

t

+

∇ ·

K f_wn+1

λ

nnw+1



∇

P_cn+1

+ (

ρ

nw

−

ρ

w

)

g

∇

z



+

f_wn+1un_t+1

−

qw

ρ

w n+1

=

0

.

(10)

Thesaturationequationisanonlinearfunction.Thus,weemployaNewton-Raphson’smethodforEq.(10),whichleads toa sequenceofsystemsof theformJν

δ

xν+1

= −

rν.Here,

δ

x isthevector ofincrementforsaturation,J istheJacobian matrix,ristheresidual,and

ν

theiterationindex.

ForeachcellK ,onecanwritetheresidualas rn_K+1

= φ(

S_Kn+1

−

Sn_K

)

−

1

|

K

|



eK L∈EK eK LFn_{K L}+1

−

qw

ρ

w n+1

=

0

.

(11)

Here,

|

K

|

isthevolume ofelement K ,eK L thearea ofthe interfacebetweencells K and L and

E

K theset offaces e of

theelementK.Finally,Fn_{K L}+1 isanumericalflux,whichcanbedecomposedintothreepartstohaveaseparatetreatmentof eachmobilityterm.Namely,

Fn_{K L}+1

=

V_{K L}n+1

+

Gn_{K L}+1

+

Cn_{K L}+1

,

(12) whereVn_{K L}+1,Gn_{K L}+1 andC_{K L}n+1aretheviscous,buoyancyandcapillarynumericalfluxes,respectively.Intheviscousnumerical flux,themobilityiscalculatedusinganupwindmethodbasedonthesignofthetotalvelocity,i.e.,

Vn_{K L}+1

=



t fw

(

Sn_U+1

)

uTK L if uTK L

>

0

t fw

(

Sn_D+1

)

uTK L otherwise

.

(13) Here, SU and SD denote theupstream anddownstreamsaturation values,respectively. Inthegravity numericalflux, the

mobilityiscomputedbasedonthedensitydifferences,i.e.,

Gn_{K L}+1

=



t TK Lfw

(

SUn+1

)λ

nw

(

SnU+1

)

ρ

g

(

zL

−

zK

)

if

ρ

g

(

zL

−

zK

) >

0

t TK Lfw

(

SnD+1

)λ

nw

(

SDn+1

)

ρ

g

(

zL

−

zK

)

otherwise

.

(14) Here,

ρ

= (

ρ

nw

−

ρ

w

)

,andTK Listheinterfacenormalizedtransmissibility(normalizedbythegridcellcross-sectionarea),

whichiscomputedastheharmonicaverageoftheneighboringcellparameters,i.e.,

TK L

=



dK e KK

+

deL KL



−1

,

(15)

wheredK eanddeL arethedistancesbetweenthefaceeK LandthecentersofcellsK andL,respectively.Finally,thecapillary

numericalfluxiscalculatedas

Cn_{K L}+1

=



t TK Lfw

(

Sn_U+1

)λ

nw

(

Sn_U+1

)(

PncL+1

−

P n+1 cK

)

if uTK L

>

0

t TK Lfw

(

Sn_D+1

)λ

nw

(

Sn_D+1

)(

PncL+1

−

P n+1 cK

)

otherwise

.

(16) ThequalityofthesolutionofEq.(10) ishighlyinfluencedbytheresolutionofthespatialandtimediscretizationscheme. The objectiveof thiswork isto developa space-time multileveladaptivemethod (ADM-LTS)to solveEq.(10) accurately andefficiently.

3. ADM-LTSmethod

In thissection, first,the originalADMmethod [11] isreviewed,then, thenewly proposed ADM-LTS algorithmis pre-sentedindetail.

3.1. TheADMmethod

TheADMmethodisemployedtoreducethecomputationalcostassociatedwiththesolutionofthelinearsystemarising fromthelinearizedproblemofEq.(10).

Let usconsider a domain discretized with a high resolutiongrid which is assumed to be fine enough to capture all relevantphysicsandtohonortheheterogeneousdistributionofthegeologicalproperties.Giventhisfine-scalediscretization, a hierarchyofnl nestedcoarsegrids isconstructed.Eachgridis formedby Nl

=

Nlx

×

Nly

×

Nlz grid cells,wherel isthe

resolutionindexandl

=

0 representsthefinegridresolution.

The setofall gridcellsbelongingtoresolutionlevell iscalled



l.AteachtimestepADMdefinesamultilevelgrid by combininggridcellsbelongingtothehierarchyofgridspreviouslydefined.GivenamultilevelADMgrid,letusdefine



l as thesetofgridcellsbelongingtoalllevelsfrom0tol whichare presentintheADMgrid.Additionally,itisconvenientto definetheset

l as

l

= 

l

∩ 

l.

GivenanADMgridformedbythesetofgridcells



l,ADMassumesthatthefinescalesolutioncanbeapproximatedby employingasequenceofprolongationoperators,i.e.

δ

xf

≈ δ

x

= ˆ

P10

. . . ˆ

Pll−1

δ

xA D M

.

(17)

Here, operator P

ˆ

i_i−1 interpolates the solution atlevel i to the finer resolution level

(

i

−

1

)

and

δ

xA D M is the vector of incrementfortheADMsolutionontheadaptivemultilevelgrid.Thefine-scaleJacobiansystemismappedtotheADMgrid by

ˆ

Rl_l−1

. . . ˆ

R0₁JP

ˆ

1₀

. . . ˆ

Pl_l−1

δ

xA D M

= − ˆ

R_ll−1

. . . ˆ

R0₁rf

,

(18)

whereR

ˆ

i_i−1 istherestrictionoperatorwhichmapsthesolutionfromresolutionatlevel i tocoarserlevel

(

i

−

1

)

.Inorder to ensure massconservation atall level,a finitevolume restriction operatoris considered [5]. Thus,the entry

(

i

,

j

)

of a restrictionoperatorreads

ˆ

Rl_l−1

(

i

,

j

)

=



1 if cell i

∈

l_{and cell j}

_∈

l−1

_,

δ

i j otherwise.

(19) Additionally,constantinterpolationisconsideredforsaturation,

ˆ

Pl_l−1

=

ˆ

Rl_l−1

T

.

(20)

Consequently,nosolutionstepisrequiredtocomputetheprolongationoperator. 3.2. ADMmethodwithlocaltime-stepping(ADM-LTS)

Ateachtimestepn,havingsolvedthepressureequationandcomputedthetotalvelocityfield,weaddressthetransport equationemployingtheADM-LTSalgorithm.

First,Eq.(18) issolvedwithtimestep

t overthewholedomainonthecoarsestgridresolution(lmax)formedbycells

belongingto



lmax _{andlimitingrefinementonlyaroundthewells.Then,basedonthecoarsesolutionobtained,theproper} ADMgridresolutionischosenaccordingtoafront-trackingcriterion.Twoalternativefront-trackingstrategiesareconsidered inthiswork:

•

acriterionbasedonthesaturationdifferencebetweenneighboringcells.Acelli belongingtolevell isrefinedwhenever thesaturationdifference,asdefinedin[11],betweeni andoneofitsneighborsexceedsauser-definedtolerance

x

.

•

atime-dependent criterion combinedwiththeprevious oneto determinewhethercells belongingto



0 shouldstay

fine. Letus define

ψ

S

=

Sn+1

−

Sn. A fine cell i is kept atthe fine resolution only if

ψ

Si

>



t, where

t

is a user-definedtolerance.Asimilartime-basedcoarseningcriterionhassuccessfullybeenusedintheliteratureforchannelized heterogeneousproblemswherestationarygradientsarepresent[23].

Once theADM gridresolution has beendefined, the solution isrecomputed forall cells belonging to



lmax−1 _witha time step

tlmax−1

=

t

/

η

byimposinglocalboundary conditionsasdescribed indetailinthefollowingsubsection.Here

η

isthetimerefinementratio.Then,thesameoperationisrepeatedforallresolutionlevelsl untill

=

0 hasbeenreached. Thus,eachresolutionlevell (formedbythesetofgridblocks



l)issolvedwithatimestep

tl

=

tl+1

/

η

.Forthefinest

level(l

=

0)ADM-LTS onlyrecomputesthesolution,withtime step

t0

=

t

/

η

lmax,forasubset,defined



0A,ofthecells

belongingto



0.Infact,onlyfinecellsforwhich

ψ

S

=

Sn+t1

−

Sn

>



t arepartoftheset



0A.

Themethodadvancesintimefortheactivecellsin



0_A untiltheyreacht

=

tn

+

t1.Oncethey aresynchronized,cells

in



1 advancein time.At thispoint, a newsetof cells



0_A isselected andthesecells are advanced by

t1 performing

η

t0 time steps. Once all cells in



1 have reachedtime t

=

tn

+

t2 another time step

t2 can be performed forall

cells belonging to



2. This is a recursive procedure which is performedfor all levels until all cells have reached time tn+1

=

tn

+

t.

Fig.1illustratesaschematicoverviewoftheADM-LTSmethodwhere

η

andlmax arebothtakenequalto2.Fig.2shows

anexampleoftheADMgridateachstepandtherefiningarea.Attheglobaltimestep

t,thesolutioniscomputedonthe coarsestresolutionlmax.Attheintermediate timesteptheADMgridresolutionisdefinedandthesolutionisrecomputed

withtheintermediatetime stepeverywhereexceptatthecoarsestregion (middlefigure).At theend,the methodchecks theerrors anddefinesthe setofactivecells



0_A (pinkregionon theright),thesolutionis recomputedwiththe smallest timestep.Notethatthesolutionateachlocaltimestepisperformedatadifferentspace-timeresolution.Furthermore,all cellsbelongingtothecoarsestlevelmarchwiththebiggesttimestepsize.Ontheotherhand,thesaturationatcellsin



0_A (i.e.,wherehigherresolutionisneeded)naturallyadvancewithsmallertimestepsizes.

WeremarkthatthequalityofthesolutionandtheefficiencyoftheADM-LTSalgorithmwillbesignificantlyaffectedby the choiceof theuser-defined tolerances

x

and

t

. Eventhough,a generalstrategy cannot be defined, simpleanalytical solutions(e.g.,Buckley-Leverett)canbe usedtoestimate thevalue ofthe saturationfrontanddeduceeffectivevaluesfor thetolerancesthatwillhavetobedeterminedtoobtainthedesiredtrade-offbetweenaccuracyandefficiency.

3.2.1. Localsystemsandlocalboundaryconditions

Foreachresolutionlevell,thesetofgridcells



l issolvedwiththecorrespondingtimestep

tl

=

t

η

(lmax−l).Thenumber ofactivecellscontainedby



l isdenotedby NlA.

When solving forthecells belonging to



l, thenumerical fluxat theinterface betweentwo cells K and L such that K

∈ 

l

∧

L

∈

l+1 isapproximatedby Fn_{K L}+β(i)

=

F n+γ(j) K L

η

(21) where

Fig. 1. Schematic overview of a time step for the ADM-LTS strategy withη=2 and lmax=2.

Fig. 2. Example of ADM grid and active regions for the refinement time steps withη=2 and lmax=2.

β(

i

)

=

i

η

(lmax−l)

,

i

= {

1

· · ·

η

(lmax−l)

}

₍₂₂₎

γ

(

j

)

=

j

η

(lmax−l−1)

,

j

= {

1

· · ·

η

(lmax−l−1)

}.

₍₂₃₎

Fig. 3. Test case 1[99×99]- CFL values for different global time steps at time t=500 days.

Thus,Eq.(11) canbemodifiedtoaccountforthepresenceofdifferenttimelevelsas rn_K+β(i)

= φ(

Sn_K+β(i)

−

Sn_K+β(i)−1

)

−

1

|

K

|



eK L∈EK A eK LFnK L+β(i)

−

qw

ρ

w n+β(i)

−

1

|

K

|



eK L∈EK L eK L Fn_{K L}+γ(j)

η

.

(24)

Here,

E

K A isthesetofinterfacefluxesexchangedbetweentwocells K and L bothbelongingto



l.Additionally,

E

K Lis

thesetoffluxesattheinterface betweentwocells K and L where K

∈ 

l and L

∈

l+1.Notethat,forl

=

0 theresidual fortheactivecellsisthesamedescribedinEq.(11),but

E

K L wouldbethesetoffluxesattheinterfacebetween



0A and



1

\ 

0_A.

Weremarkthat,foreachlevell,thelinearsystemtobesolvedhassize Nl_A

×

Nl_A,whichissignificantlysmallerthanthe fullfine-scalesystem.Ateachtimestep,onlyasub-setofthecellsofeachlevell isrestrictedandprolonged.

The following section presentsa study of theperformance of the ADM-LTSmethod forvarious 2D and3D test cases involvingdifferentfluidphysics.

4. Numericalresults

TheperformanceofthenewlydevelopedADM-LTSstrategyisthoroughlyinvestigatedforseveralchallenging testcases. Forallthecasespresented,quadraticrelativepermeabilitycurvesareconsidered.Gravitational andcapillaryforcesareonly considered intest cases 2and7, respectively.Furthermore, allerrors are computedwithrespect to areference solution, obtainedbyemployingahighresolutiondiscretizationbothinspaceandintime.

4.1. Testcase1:2Dhomogeneousreservoir

Thefirsttestcaseisa100

×

100

[

m2

_]

_{homogeneousreservoir,withisotropicpermeabilityof5}

_×

₁₀−15

_[

_m2

_]

_.A

pressure-constrainedwetting-phaseinjector well ispositioned inthe bottom-left cornerofthedomain withapressure pinj

=

108

[

Pa

]

, whereas a production well is presentin the top-right cornerwith a relative pressure of pprod

=

0 [Pa]. The phase

viscosityvalues are

μw

=

10−3

_[

_Pa

_·

_s

_]

_and

_μnw

₌

₁₀−2

_[

_Pa

_·

_s

_]

_{forthe wettingandnon wetting phase,respectively. The}

finalsimulatedtimeis600 [days]afterinjectionhasstarted.

A fine-scalegrid with 99

×

99 cells is imposed on thedomain. ADM-LTS employs, a time refining ratio

η

=

2 and a spacecoarseningratioequal to3 inalldirections. Theuser-defined tolerances forthecoarseningandrefinementcriteria are

x

=

0

.

07 and

t

=

5

×

10−2_.

Simulationsarerunemployingthreedifferentglobaltimestepsizes:5,10and20days.Fig.3reportstheCFLvaluesat timet

=

500 daysforthethreedifferenttimestepsforfine-scaleinspacesimulations.

Fig.4showsacomparisonoftheADM-LTSsolutionwiththereferencesolution(ontheleft)attimet

=

500 [days]using threedifferentsizesoftheglobaltimesteps.

Fig.5reportstheerrorforthesaturationattimet

=

500 daysbetweenareferencesolutionandtheADMmethodwith finetimesteps(firstcolumn)withtheLTSapproach(secondcolumn)andwiththecoarsetimesteps(thirdcolumn)forthe threedifferenttimestepssizes

t

=

5 (firstrow),

t

=

10 (secondrow)and

t

=

20 (thirdrow).InallcasestheAMD-LTS approachimprovestheerrorswithrespecttothecoarsetimestepapproach.

ThecomplexityofthealgorithmisshowninFig.6.Inparticular,eachcolumnrepresentsthetotalamountofactivecells multipliedbythenumberofNewtoniterationsinvolvedtocomputethesolution,forthethreeapproachesandforthethree differentglobaltime stepsizes. Notethat,toobtainthesolutionattime t

=

600 [days], 120

,

60 and30 globaltime steps havebeenperformedusingthethreeanalyzedtime steps.We remarkthattheerrors obtainedby employing theoriginal ADMmethodwithafinetimesteparecomparabletothoseobtainedwithADM-LTSintermsofaccuracy.

Fig. 4. Test

case 1

[99 ×99]- Reference solution (first column) and ADM-LTS solution using a global time step size: t=5, 10 and 20 [days] at time t=500 days for the second, third and fourth column, respectively.

Fig. 5. Test

case 1

[99 ×99]- Saturation errors for the ADM method with fine time steps (first column), ADM-LTS method (second column) and ADM coarse time steps method (third column) for the three different global time step sizes.

Fig.7showsthecomplexity ofasingleglobaltimestep.FortheADMmethodwithfinetimesteps,thelocalstepsare justthesmallstepsappliedatthewholedomain.AttheendofthelocalstepsboththeADM-LTSmethodandtheADMfine stepmethodreachthesametime.FortheADM-LTSmethod,localstep1indicatestheglobalsteponthecoarsestgrid,step 2and5aretheintermediatetimestepsperformedonlevel0and1oftheADMgrid,andtheotherlocalstepsarethesmall time steps forthe activecells detected by theerror estimatorintime. Inparticular, wecan noticethat the intermediate time stepshavealmostthesamecomplexity ofthesmalltimestepsoftheADMfinemethod,evenifthesize ofthetime stepistwotimesbiggerwithalmostthesamenumberofactivecells.Thisisduetotheimprovementoftheinitialguessfor theNewtonloop.Intheintermediatetimestepsweuseasinitialguessalinearcombinationofthesolutionoftheprevious timetnandthesolutionobtainedonthecoarsestgridatthenewglobaltimetn+1.Inthesmalltimestepsisnotnecessary toperformthistechniquesinceasmallstepisusedtoadvanceintime.

The same test caseis analyzedafter performing a 2

×

2 refinementof thespace fine-scalegrid. Inorder to obtaina reasonablesolution,usingaglobal

t equalto20[days],weneedtocomputemorelocaltimestepsinsidetheglobalone, soarefiningratioequalto4hasbeentakenintoaccount.

Fig. 8reports thecomplexity forthe entiresimulation usingADM-LTSmethod andtheADM withfinetime steps. To obtain thesolutionatfinaltimet

=

600 [days]withaglobaltimestepequalto20days,thesamenumberofglobaltime

Fig. 6. Test

case 1

[99 ×99]- Total amount of active cells multiplied by number of Newton iterations for the three different time step sizes. On the top of each bar the mean in time of the averaged absolute difference respect to the reference solution for the saturation is displayed Es=mean|S(tf) −Sre f(tf)|

where tf is the final time 600 days.

Fig. 7. Test

case 1

[99 ×99]- Computational complexity history at each local times step within a global step. The computation complexity is the number of active cells multiplied by the number of Newton iterations.

Fig. 8. Test

case 1

[198 ×198]- Total amount of active cells multiplied by number of Newton iterations for the ADM with fine time steps and the ADM LTS method.

stepare involved(30 timesteps intotal). Ofcoursethenumberoflocal timestep forboththeLTS methodandthefine timestepsapproachhasincreased;buttheratiobetweenactivecellsandtotalcellsdecreases.

Fig.9showstheaveragednumberofactivecellstimesthenumberofNewtoniterationsforeachlocaltimestepwithin aglobaltimestep.

InFig.10wecanseethattheADM-LTSapproachreducestheerrorsobtainedusingacoarsegridintime. 4.2. Testcase2:3Dhomogeneousreservoir

A 3D 108

×

108

×

108

[

m3

_]

_{homogeneous reservoir is considered in thistest case. The domain is discretized, atthe}

Fig. 9. Test

case 1

[198 ×198]- Computational complexity history at each local times step within a global step. The computation complexity is the number of active cells multiplied by the number of Newton iterations.

Fig. 10. Test

case 1

[198 ×198]- Saturation errors at time t=540 [days] for the ADM method with fine grid in time (left), ADM-LTS method (center) and the ADM method with coarse grid in time (right).

Fig. 11. Test

case 2 - Saturation profile (top row) at time t

=1500 days (left) and at time t=8750 days (right) for the case without gravity. Active cells for

the level lre f=0at time t=1500 (bottom-left) and saturation profile inside the domain at time t=8750 (bottom-right).

testcase.Thesizeoftheglobalstepsisequalto125days.Thesimulationendsafter70globalsteps.Thetolerancesforthe coarseningcriteriainspaceandtimearesettobe

x

=

0

.

2 and

t

=

5

×

10−2.

Fig.11reportsthesaturationmapsattwodifferentsimulationtimes(on thetop)neglecting thegravityeffects.Italso displays thesetofactivecells



0_A attimet

=

1500 days (left,bottom)andasectionofthesolutionatfinal timet

=

8750 days (right,bottom). Notethat ADM-LTSautomaticallyemploys finecells only aroundtheadvancing saturationfront and thattheactivecellsintimeareonlyafractionofthem.

Fig. 12showsthe saturationmapsandtheactive cells



0_A atdifferenttimesconsideringthe samescenariobut intro-ducingadensityratioof

ρw

/

ρ

nw

=

5

/

4.Asexpected,theheavierfluidoccupiesthebottomofthereservoir.

Fig. 12. Test

case 2 - Saturation profile (top row) at time t

=1500 days (left), t=6375 (center) and t=8750 days (right) for the case with gravity. Active cells for the level lre f=0at time t=1500 (bottom-left), active cells for the level lre f=0at time t=6375 (bottom-center) and saturation profile inside the domain at time t=8750 (bottom-right).

Fig. 13. Test

case 2 - Computational complexity – total amount of active cells multiplied by the count of Newton iterations for different scenarios (top).

Computational complexity history for ADM-LTS and ADM fine-scale time are also presented (bottom).

Fig.13showsthetotalandthemeancomputationalcomplexityperlocaltimestepforbothADM-LTSandADMmethods. Thisfigureillustratestheresultswithfinetimesteps,withandwithoutgravityeffects.Thetransportequation,inthe pres-enceofgravitationalforces,becomeshighlynonlinear.As aconsequence,moreNewtoniterationsarerequiredtoconverge comparedwiththecasewithoutgravitationaleffects.Notethat,asshowninFig.14,thenumberofADM-LTSactive cells forthe twodifferent scenariosis thesame during theentiresimulation.The numberofactive cells increasesduring the progressofthesimulationtime,resultinginreductionoftherelativeerrors(Fig.14).TheerrorsofADM-LTSforthecases withandwithoutgravityarecomparable.

4.3. Testcase3:2Dhomogeneousreservoirwithbarrier

A2Dhomogeneous reservoirwithlowpermeabilitybarriersisconsidered,asshowninFig.15.Thesamepermeability field was presentedin [23]. The domaindimensionsandthe physicalparameters are thesame ofthe firsttest case, the same99

×

99 finescalegridisimposed.Theglobaltimestepisequalto50 [days]andthesimulationendsafter100 global timesteps(t

=

5000 days).

Simulationsare carriedoutboth withthe originalADMmethodemployinga globalfine time-stepandwithADM-LTS. Thecoarseningandthetime-refinementcriteriatolerancesaresetto

x

=

0

.

05 and

t

=

0

.

005.

Fig. 14. Test

case 2 - Number of ADM-LTS active cells for the test cases with and without gravity (left). Also shown (right) are the saturation l

1_{-norm errors.}

Fig. 15. Test case 3 - Absolute permeability field.

Fig. 16. Test

case 3 - Saturation profile and ADM grid at different time steps (columns) for ADM with coarse time steps and classical ADM grid resolution

(first row) and for ADM-LTS method with the new ADM grid resolution (second row).

Fig.16showsacomparisonofthesaturationprofileandthegridresolutionforthetwodifferentstrategies.Theoriginal ADMmethod withasaturation difference-basedcoarseningcriterion (toprow)employs a largenumberoffine gridcells whereversaturationgradientsarepresenteveniftheyarestationary.Ontheotherhandthenewlyproposedgridresolution criterion (bottomrow) for the ADM-LTSapproach uses finecells only in those regions where the saturation gradient is moving,reducingthenumberofactivecells.

Fig. 17. Test

case 3 - number of active cells employed in ADM with fine grid in time and ADM-LTS simulations expressed as percentage of fine grid cells

(left) and the saturation relative errors in l1_{-norm for the ADM fine and ADM-LTS method (right).}

Fig. 18. Test

case 3 - Total amount of active cells multiplied by number of Newton iterations (left) and computational complexity history at each local times

step within a global step (right) for the ADM with fine time steps approach and for the ADM-LTS method.

Fig. 19. Test case 4 - Natural logarithm of the permeability.

Fig.17reportstheevolutionoftheactivegridcellspercentageforthetwodifferentapproaches(left)andtheevolution oftherelativesaturationerrorinl1_{-norm(right).Intheearly steps,wecanseethat theADMfinewithjust thegradient}

inspace criterion approachemploys almost the samenumberof active gridcells used by theADM-LTS method.Forthe ADMwithfinetimestepsateverysmalllocaltimestepwesolveboththeflowandthetransportequations,insteadforthe ADM-LTSapproach onlythetransport equation issolvedforthe localsteps.Thisis thereasonwhyinthefirst fivesteps thesaturationerrorsfortheADM-LTSapproacharelargerwithrespecttotheADMfinestepsapproach.Instead,inthelast stepstheerrorsincreasebecausealowernumberoffinegridcellshasbeenused.

Fig.18showsthetotalcomplexity(numberofactivecellsmultipliedbythenumberofNewtoniterations)fortheADM withfinestepsandtheADM-LTSapproach.NotethatthelocaltimestepsofADM-LTSmethodreducethecomplexityofthe systemcomparedtotheclassicalADMapproach.

4.4. Testcase4:heterogeneousreservoir(SPE10toplayer)

Inthistestcasea heterogeneousreservoirisconsidered.The permeabilitymapisthe toplayeroftheSPE10test case [24] anditispresented,inlogarithmic scale,inFig.19.The sizeofthereservoiris2200

×

600

[

m2

_]

_anda216

_×

_{54 grid}

isemployedatthefinestlevel.Theinjectorisatthetopleftcornerandhasaconstrainedpressure107

_[

_Pa

_]

_{.Aproduceris,}

instead,locatedatthebottomrightcornerofthedomainwithapressureequalto0

[

Pa

]

.Theporosity ofthereservoir

φ

isequalto0.2.Theviscosityforthewettingphaseis10−5

[

Pa

·

s

]

,whereas,forthenon-wettingphase,is10−4

[

Pa

·

s

]

.The coarseningratioforthespacegridisequalto2aswellasthetimerefiningratio.Theerrortoleranceforthetimeestimator isequalto5

×

10−2_.

Fig. 20. Test

case 4 - Saturation map and ADM grid for the ADM with fine time step approach with classical grid criterion for different values of the

threshold x=0.05, 0.1, 0.2 (row 1, 2, and 3) and for the ADM-LTS method with the new grid criterion x=0.05 and t=0.05 (row 4) at time t=1200

days (first column), t=15000 days (second column) and t=20000 days (third column).

Fig. 21. Test case 4 - Active cells for the refinement in time, at time t=1200 days (left), t=15000 days (center) and t=20000 days (right).

Fig. 22. Test

case 4 - Number of active cells expressed as percentage of fine grid cells (left) and saturation relative errors in l

1_{-norm (right) for the ADM} with fine grid in time with different values of the threshold and for the ADM-LTS simulation.

Fig.20reportsthesaturationmap andtheADMgridfordifferentthresholdvaluesoftheADMgridresolutioncriterion usingthe classicalADMapproachwithfine timesteps,andtheADM-LTSapproachwiththenewgridresolutionstrategy. Theclassicalapproachuses,forsmallthresholdvalues,alargenumberoffinegridcells.Ifwerelaxthethresholdparameter the methodis notabletocapturethefronts.Thanks tothenewADM-LTSapproach,themethodisabletoapply thefine gridcellsonlywherethefrontismovingfast(highpermeabilityregions).

Fig.21showstheactivecellsintimeatthefinestlevellre f

=

2 for differentglobaltime steps.Themethodrecomputes

thesolutionwithsmalltimestepsonlyforafewpercentageofcellswherethefrontcrosseshighpermeabilityregions.In fact,inthelastsnapshot,thesaturationprofileisalmostdevelopedeverywhereandso,thesetofactivecellsisverysmall.

In Fig. 22 we compare the number of active cells and the saturation errors for the different simulations. Using the classical ADMapproachwithsmallvaluesofthetolerancealotofactivegridcellsare employedgivingverysmallerrors. The classicalADMapproachwithlargertolerancevalue andtheADM-LTSmethodarecomparable intermsofactivecells duringallthesimulationbuttheADM-LTSapproachgivesbetterresultsintermoferrors.

Fig. 23reportsthe complexity ofthefour simulations.The ADMapproachwith finegrid intime andsmallthreshold valuesisreallyexpensive. TheADM-LTSapproachiscomparabletotheADMwithfinetimestepapproachandlargevalue of

x

but,asshownpreviously,thesolutionoftheclassicADM,inthiscase,isnotasaccurate.

Fig. 23. Test

case 4 - Total amount of active cells multiplied by number of Newton iterations (top) and computational complexity history at each local times

step within a global step (bottom) for the ADM approach.

Fig. 24. Test case 5 - Natural logarithm of the permeability.

4.5. Testcase5:heterogeneousreservoir(SPE10bottomlayer)

The permeabilityofSPE10bottomlayer isused forthistest case,asshown in Fig.24.This layer,withrespect tothe previouscasehashighercontrastsandmoreachannelizeddistribution.

Theglobaltime stepisequalto 10daysandthe simulationendsafter50globaltime steps.Theinput parametersfor thewellsandthephysicalpropertiesareidenticaltoTestCase4.

Thetop2rowsofFig.25showthesaturationdistributionatsimulationtime of150,250 and350 daysobtainedwith



x

=

0

.

15 and

x

=

0

.

2, respectively. Thebottom rows, instead, show thesaturation map, atthesame simulations times,

obtainedbyemployingtheADM-LTSmethodwith

x

=

0

.

05 and

t

=

5

×

10−2and

t

=

5

×

10−3.TheclassicADMapproach employsalargenumberofactivecellduringthesimulation,insteadtheADM-LTSmethodisabletoselectafinescalegrid onlywherethefrontsaremoving.

Fig.26showstheactivecellsintimeforlre f

=

2attime 150,250 and 350 days.Asexpected,forsmallervalueofthe

thresholdmorecellsareinvolvedintherefinementstep.

Thehistoryofthepercentageofactivecellsemployedbythedifferentsimulationstrategiesforthevarioustolerancesis showninFig.27(left),alongwiththel1 _{normofthesaturationerror(right). Forboth theADM-LTStolerancevaluesless}

activecellsareinvolvedrespecttotheclassicalADMapproach.Sinceasmallernumberofcellsisemployed,thesaturation errorsarehigherbutstillofthesameorderofmagnitude.

Fig.28reportsthecomplexityofthefoursimulationsfordifferenttolerancevalues. 4.6. Testcase6:heterogeneousreservoirswithdifferentlayeringorientations

A500

×

500 m2 2Dreservoirisconsideredonwhicha99

×

99 gridisimposed.Thefluidproperties,thelocationofthe wellsandtheirconstraintsarethesameasintheprevious testcases.Fivesetsofpermeabilityfields,withdifferent layer-ing orientationandcreatedusingsequentialGaussiansimulations withsphericalvariogramanddimensionlesscorrelation lengths0

.

5 and0

.

02 asproposedin[25],areconsidered.Eachsetconsistsof20statisticallyidenticalrealizations.

Fig. 25. Test case 5 - Saturation map and ADM grid at 150, 250 and 350 days for the ADM approach with fine time steps and the ADM-LTS approach.

Fig. 26. Test case 5 - Active cells for the refinement level lre f =2, at 150 (left), 250 (center) and 350 (right) days for the two threshold values.

Fig. 27. Test

case 5 - Number of active cells expressed as percentage of fine grid cells (left) and saturation relative errors in l

1_{-norm (right) for the ADM} with fine grid in time and for the ADM-LTS simulations.

Fig.29showsonerealizationforeachset.Injectionofthewettingphase,for560days,issimulatedforeachrealization. Simulations are runwiththe ADM-LTSmethod.Forallruns, the spatialcoarseningcriterion tolerance is

x

=

0

.

008.Two differentvaluesareinsteadconsideredforthetime-basedcriteriontolerance,

t

:5

×

10−2 _and5

_×

₁₀−3_.

Fig.30showsacomparison,foronepermeabilityrealizationofeachset,ofthesaturationmapattheendofthe simu-lationobtainedwithfine-scale(timeandspace)simulation(toprow),ADM-LTSemployingafixedrefinedtime-step.

Fig. 31 displays the active cells in time for the last refinement level ofthe last globaltime step.As expected, using a bigger value of the tolerance forthe time error estimator, just few cells need to be computedwith smalltime steps. Moreover,thespacegridchangesandallowstousecoarsergridcells.

Fig. 32 represents the mean andthe standard deviations of the complexity for the ADM-LTS method using the two differenttime-basedcriterion tolerancesandforthesolutioncomputedwiththefinegridresolutionbothinspaceandin time.Notethatthey-axisscaleforthetwopicturesaredifferent.

Fig. 33showsthe meanandthestandard deviations ofthesaturation errorsrespect to thereferencesolutionfor the ADM-LTS method usingthe two differenttime-based criterion tolerances. From thesestudies,one can concludethat the ADM-LTSperformsrobustlywhenseveralequiprobablerealizationsareconsidered.Inotherwords,theerrorand computa-tionalcomplexitiesforall20realizationsarenotmuchdifferentcomparedwiththeaveragevalues.

Fig. 28. Test

case 5 - Total amount of active cells multiplied by number of Newton iterations (top) and computational complexity history at each local times

step within a global step (bottom) for the ADM approach.

Fig. 29. Test

case 6 - One of the 20 realization of each of the 5 sets of permeability fields with different angles (0 deg, 15 deg, 30 deg, 45 deg and patchy

from left to right).

Table 1

Test case 7 - Wells coordinates and constraints.

Well x y Pressure [bar] Prod (W1) 1 1 120 Prod (W2) 99 1 120 Prod (W3) 99 99 100 Prod (W4) 1 99 100 Inj (W5) 50 50 150

4.7. Testcase7:capillaryforces

A500

×

500

[

m2

_]

_{heterogeneousreservoirisconsidered.The fine-scalegridcontains99}

_×

_{99 cells.Thephase viscosity}

valuesare

μw

=

1e

−

4

[

Pa

·

s

]

and

μnw

=

1e

−

3

[

Pa

·

s

]

.Fivepressure-constrainedwellsarepresent,asshowninFig. 34 alongwiththeheterogeneouspermeabilityfield.ThewelllocationsandthepressurevaluesarepresentedinTable1.

Thesimulationisrununtil2500

[

days

]

afterinjectionisreached.Atimestepsizeequalto10

[

days

]

isused.The ADM-LTSemploysatimerefiningratioof

η

=

2 andaspacecoarseningratioof3.Theuser-definedtolerancesaresetas

x

=

0

.

25 and

t

=

1

×

10−2_{.Thesamecapillarypressurefunctionasin[12,26] isconsidered,i.e.,}

Pc

=

σ

cos

(θ )

φ

K J

(

S

),

(25) where J

(

S

)

=

0

.

05



1 S



₋0.5

.

(26)

Fig. 30. Test

case 6 - Comparison of the saturation profile, for one realization of each set of permeability fields at time t

=560 days. Two different threshold

values for the time error estimator are employed for the ADM-LTS simulation (center row and bottom row), the fine scale solution are also shown (top row).

Fig. 31. Test

case 6 - Active calls at the last refinement level for the last global time step using two different threshold values for the error estimator in

time.

Here, S isthesaturationofthewettingphase,

σ

=

4

.

361

×

10−2

_[

_Pa

_·

_m

_]

_{thesurfacetensionand}

_θ

₌

_{0 thecontactangle.}

Fig. 35compares the referencesolution(fine-scale) andthe ADM-LTSapproach withandwithout includingthe effectof capillarypressure.TheADM-LTSapproach,asshowninFig.36,inpresenceofcapillarypressureemploysmoreactivecells, since thesaturationmap ismorecomplexcomparedwiththecasewithoutcapillaryeffects.Inthistestcase,thecomplex physics inducedbycapillaryheterogeneityispresentinalmosttheentirecomputational domain.Therefore,theadvancing saturationfrontispresentinalargesub-regionofthedomain.Astheresult,inthiscase,70% oftheactivecellsareusedto obtainaccurate solutions.Fig.36showstheerrorsinl1 _{normfortheADM-LTSmethodwithandwithoutcapillaryeffects.}

After the globaltimestep 150,the saturationerrorincreasesdueto coarseningsome finecells. Thenafter200days,the errorisstabilized.

Fig. 32. Test

case 6 - Mean and standard deviation of complexity over 20 realization for the ADM-LTS method (left) and for the reference solution computed

with fine grid resolution both in space and time (right).

Fig. 33. Test

case 6 - Mean and standard deviation errors of the saturation errors over 20 realization for the ADM-LTS method with different time threshold

values respect to the reference solution ES=meantN=t1



mean|Sf(t)−S(t)|



.

Fig. 34. Test case 7 - Base 10 logarithm of the permeability field.

5. Conclusions

Adynamiclocalspace-timesteppingschemefortransportinporousmedia(ADM-LTS)wasdevisedandintegratedwithin thesequentially-coupledmultiphaseflow-transportsystem.Thestudiedcoupledmultiphaseflowprocessedobservepressure andvelocitychangeateverytime step.TheADM-LTSmethodoperatesintwostages:first thesolutionisobtainedatthe coarsestspace-timegrid.Thenthisinitialsolutionisimprovedbyimposinganadaptivemulti-resolutiongridinspaceand time.The resolutionofthespace-time gridisdefinedbasedonan errorcriteria,withuser-defined thresholdvalues.This strategyalsoguaranteeslocalmassconservation.

Fig. 35. Test case 7 - Reference solutions (left) and ADM-LTS saturation maps (right) without (top) and with (bottom) capillary heterogeneity effects.

Fig. 36. Test

case 7 - Number of active cells expressed as percentage of fine grid cells (left) and saturation relative errors in l

1_{-norm (right) for the ADM-LTS} simulations.

ComparedwiththeclassicalADMapproach,theADM-LTSmethodemployscoarsegridswherehighsaturationgradients exist, butthephase velocityis nearlyzero.This happenswhenthesaturation frontfaces impermeable zonesorbarriers. The ADM-LTSgridresolutionisdecidedimplicitlybythementionedtwo-stagestrategy,whereinitiallythecoarsest space-time gridisusedeverywhere.Then,basedontheintroducederrorcriteria,thetransportsolutionislocallyimprovedwith a multirate multilevelstrategy onan adaptivespace-timegrid. Thismethodisfound promisingto reducethe sizeofthe systeminthenonlinearloopwithoutlossofaccuracy.

Several numericaltest caseswithchallenging transport physics were considered. Resultsshowed thatthe ADM-LTS is promisingforreal-fieldapplications.Theresearchsimulator(withallnumericaltestcases) ismadeavailabletothepublic (https://gitlab.com/darsim2simulator).

Future work includes development of enhanced prolongation operators to further improve the saturation quality at coarser levels (specially for slow-moving fronts). Anotherongoing work is the implementation ofADM-LTS in the DAR-Sim1C++researchsimulatorforlarge-scale3Dtestcases.

Declarationofcompetinginterest

The authors declare that they have noknown competingfinancial interests orpersonal relationships that could have appearedtoinfluencetheworkreportedinthispaper.

Acknowledgements

HadiHajibeygiwasfundedbythe2019DutchNationalScienceFoundation(NWO)grant,underVidischeme,withaward number of 17509(Project ADMIRE). Matteo Cusini was sponsoredby the auspicesof the U.S. Departmentof Energy by Lawrence LivermoreNationalLaboratory(LLNL)underContractDE-AC52-07NA27344.Authors thankDelftAdvanced Reser-voir Simulation (DARSim) research group members for the fruitful discussions during the development of the ADM-LTS method.

References

[1] D.R. Brouwer, J.D. Jansen, Dynamic optimization of water flooding with smart wells using optimal control theory, SPE J. 4 (2004) 391–402, https:// doi.org/10.2118/78278-MS.

[2] R.J. de Moraes, J.R. Rodrigues, H. Hajibeygi, J.D. Jansen, Multiscale gradient computation for flow in heterogeneous porous media, J. Comput. Phys. 336 (2017) 644–663, https://doi.org/10.1016/j.jcp.2017.02.024.

[3] M. Cusini, R. Gielisse, H. Groot, C. van Kruijsdijk, H. Hajibeygi, Incomplete mixing in porous media: Todd-longstaff upscaling approach versus a dynamic local grid refinement method, Comput. Geosci. 23 (2) (2019) 373–397, https://doi.org/10.1007/s10596-018-9802-0.

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

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

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

[7]H.Zhou,H.A.Tchelepi,etal.,Two-stagealgebraicmultiscalelinearsolverforhighlyheterogeneousreservoirmodels,SPEJ.17 (02)(2012)523–539. [8] O. Møyner, K.-A. Lie, A multiscale restriction-smoothed basis method for high contrast porous media represented on unstructured grids, J. Comput.

Phys. 304 (2016) 46–71, https://doi.org/10.1016/j.jcp.2015.10.010.

[9] M.J. Berger, J. Oliger, Adaptive mesh refinement for hyperbolic partial differential equations, J. Comput. Phys. 53 (3) (1984) 484–512, https://doi.org/10. 1016/0021-9991(84)90073-1.

[10] B. Faigle, R. Helmig, I. Aavatsmark, B. Flemisch, Efficient multiphysics modelling with adaptive grid refinement using a mpfa method, Comput. Geosci. 18 (5) (2014) 625–636, https://doi.org/10.1007/s10596-014-9407-1.

[11] M. Cusini, C. van Kruijsdijk, H. Hajibeygi, Algebraic dynamic multilevel (adm) method for fully implicit simulations of multiphase flow in porous media, J. Comput. Phys. 314 (2016) 60–79, https://doi.org/10.1016/j.jcp.2016.03.007.

[12] M. Cusini, B. Fryer, C. van Kruijsdijk, H. Hajibeygi, Algebraic dynamic multilevel method for compositional flow in heterogeneous porous media, J. Comput. Phys. 354 (2018) 593–612, https://doi.org/10.1016/j.jcp.2017.10.052.

[13]J.Andrus,Numericalsolutionofsystemsofordinarydifferentialequationsseparatedintosubsystems,SIAMJ.Numer.Anal.16(1979)605–611. [14]V.Savcenco,W.Hundsdorfer,J.Verwer,Amultiratetimesteppingstrategyforstiffordinarydifferentialequations,BITNumer.Math.47(2007)137–155. [15] L. Delpopolo Carciopolo, L. Bonaventura, A. Scotti, L. Formaggia, A conservative implicit multirate method for hyperbolic problems, Comput. Geosci.

23 (4) (2019) 647–664, https://doi.org/10.1007/s10596-018-9764-2.

[16] L. Delpopolo Carciopolo, L. Formaggia, S. Anna, H. Hajibeygi, Conservative multirate multiscale simulation of multiphase flow in heterogeneous porous media, J. Comput. Phys. 404 (2020) 109134, https://doi.org/10.1016/j.jcp.2019.109134.

[17] P. Jenny, Time adaptive conservative finite volume method, J. Comput. Phys. 403 (2020), https://doi.org/10.1016/j.jcp.2019.109067.

[18] G.W. Thomas, D.H. Thurnau, Reservoir simulation using an adaptive implicit method, Soc. Pet. Eng. J. 23 (05) (1983) 759–768, https://doi.org/10.2118/ 10120-PA.

[19] D.A. Collins, L.X. Nghiem, Y.-K. Li, J.E. Grabonstotter, An efficient approach to adaptive-implicit compositional simulation with an equation of state, SPE Reserv. Eng. 7 (02) (1992) 259–264, https://doi.org/10.2118/15133-PA.

[20] J. Maes, A. Moncorgé, H. Tchelepi, Thermal adaptive implicit method: time step selection, J. Pet. Sci. Eng. 106 (2013) 34–45, https://doi.org/10.1016/j. petrol.2013.03.019.

[21]R.Helmig,MultiphaseFlowandTransportProcessesintheSubsurface:AContributiontotheModelingofHydrosystems,EnvironmentalEngineering, Springer,1997.

[22]Z.Chen,G.Huan,Y.Ma,ComputationalMethodsforMultiphaseFlowsinPorousMedia,vol.2,Siam,2006.

[23]M.Cusini,H.Hajibeygi,Algebraicdynamicmultilevel(adm)methodforsimulationsofmultiphaseflowwithanadaptivesaturationinterpolator,in: ECMORXVI-16thEuropeanConferenceontheMathematicsofOilRecovery,2018.

[24] M. Christie, M. Blunt, Tenth SPE comparative solution project: a comparison of upscaling techniques, in: SPE Reservoir Simulation Symposium, 11-14 February 2001.

[25]N.Remy,A.Boucher,J.Wu,AppliedGeostatisticswithSGeMS:AUser’sGuide,CambridgeUniversityPress,2009.

[26] B. Li, S.M. Benson, Influence of small-scale heterogeneity on upward CO2plume migration in storage aquifers, Adv. Water Resour. 83 (2015) 389–404,