A locally conservative multiphase level set method for capillary-controlled displacements in porous media

Contents lists available atScienceDirect

Journal of Computational Physics

www.elsevier.com/locate/jcp

A locally conservative multiphase level set method for capillary-controlled displacements in porous media

Espen Jettestuen

^a,∗

, Helmer André Friis

^b

, Johan Olav Helland

^b

aNORCENorwegianResearchCentre,Essendropsgate3,N-0368Oslo,Norway

bNORCENorwegianResearchCentre,Prof.OlavHanssensvei15,N-4021Stavanger,Norway

a rt i c l e i nf o a b s t ra c t

Articlehistory:

Received28January2020

Receivedinrevisedform27October2020 Accepted28October2020

Availableonline31October2020

Keywords:

Levelsetmethod Volumeconservation Three-phasedisplacements Gangliasplittingandcoalescence Porescale

Capillarypressure

We presentamultiphaselevel setmethodwithlocalvolumeconservationforcapillary- controlled displacement in porous structures. Standard numerical formulations of the level set method for capillary-controlled (or, curvature-driven) motions assume phase pressures and interface properties are spatially uniform and disregard the fact that separatephasegangliatypically havedistinct pressures.Thisisamajorproblemforthe suitabilityofsuchmethodstosimulatecapillarytrappinginporousrocksas itwilllead to severemass loss. The methodpresented herepreserves volumes ofindividual phase ganglia, whileitpredictscapillarypressuresbetweengangliaand surroundingphases. A conservativevolumeredistributionalgorithmhandlesgangliabreakupandcoalescence.The method distinguishes between three-phase systems, where separate level set functions describe the different phases, and two-phase systems, where one level set function representsinterfaces.We presentsequentialandparallelalgorithms forthenewmethod andemphasizeimportantaspectsspecifictothepatch-basedparallelimplementation.

Wevalidatethemethodnumericallybyapplyinglocalvolumeconservationtosimulations oftwoandthreephasesystemsinbothtwoandthreespatialdimensions.Themodel is testedforbothsaturationandpressurecontrolledsystemsandhandlesbothmergingand splittingofphaseganglia.

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

1. Introduction

Understanding themechanisms leading tofluid phase trappingandmobilization inthree-phase flowin porousmedia areessentialtostrategiesforimprovedoilrecoveryandcarbondioxide(CO2)storageinmatureoilreservoirs.Thisrequires knowledge of howfluid gangliabehave, withcoalescence andsplitting, inthe presenceof other phasesin thereservoir, aswell astheir impact oncapillary pressureandrelative permeability relationsused in reservoirsimulation. Inmost of the reservoir, away from wells and fractures, capillary forces typically dominate over viscous and gravity forces during multiphasedisplacementsattheporescale[1].Therefore,gangliaofthenon-wettingphase(e.g.,gasoroil),surroundedby acontinuouswettingphase(e.g.,water),likelygettrappedbycapillaryforces[2].ThisisadesiredmechanismforsafeCO2 storage[3].However,inathree-phasesystem,theactionofcapillarypressurecanmobilizegangliabydoubledisplacement, inwhichacontinuousphase(e.g.,gas)displacesanisolatedclusterofthesecondphase(e.g.,oil)thatinturndisplacesthe continuous thirdphase (e.g.,water)[4,5]. Acomplexitywithsuch three-phasedisplacements istheir differentbehaviours

*

Correspondingauthor.

E-mailaddress:[email protected](E. Jettestuen).

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

0021-9991/©^{2020TheAuthors.PublishedbyElsevierInc.ThisisanopenaccessarticleundertheCCBYlicense} (http://creativecommons.org/licenses/by/4.0/).

in the presence of triple lines, where all three phases meet, and inthe presence ofspreading oil layers separating the other twophases.Pore- andcore-scaleexperimentshaveshownthatspreadingoillayers,whichmaketheoilphase more continuous, anddouble displacements, can mobilize oil and lead to low residual oil saturation after gas invasion [6–8].

An extension of the double-displacementmechanism refers to the situationwhere the invading phase displaces a chain ofseveralneighbouringfluidclustersthatdisplacesanothercontinuousphase.Thiscanalsohappenintwo-phasesystems underweakly-wetstates whenthefluid phase present onopposite sidesofa ganglionare disconnectedfromeach other andassociatedwithdifferentpressures.Multipledisplacementsplayanimportantroleasmechanismsforoilrecoveryover multiplewater-alternate-gasinjectioncycles[9–11,8,12].

X-raymicrotomographyimagingoffluiddistributionsduring multiphaseflowexperimentsinporousmediahasallowed the analysisof shapeand topology offluid ganglia[7,12], capillarypressure statesof gangliacompared withcontinuous phases[13,14], wetting state anddisplacement mechanisms[15,16,12]. The useofdigital rockphysics, specifically direct numericalsimulationonsegmentedrockimages,ispotentiallyatoolforinvestigatingthree-phaseprocesses.Fundamental tosimulatingmultiphaseflowwithisolatedgangliaonrealisticporestructuresisthenumericalmethod’sabilitytoconserve mass(orvolume,forincompressiblefluidflow)inacorrectmanner.

There are many methods forconservation ofvolume inmultiphase flow. Lattice Boltzmann [17] andphase field [18, 19] methodsadd interaction betweenphasesanduse theemerging phaseseparation asthe basis formultiphasesolvers, wherethephaseinteractionpotentialsetsurface tensionsandcontactangles.Intheseschemesfluidmassesareexplicitly conservedbutundesirednumericaleffectssuchasmutualpollutionofnonmisciblephases, largediffuseinterfaceregions, andspuriouscurrents,canoccur.

For capillary-dominated multiphase flow regimes, a reasonable assumption is that capillary forces alone control the interface displacement. Examples ofsuch quasi-static, capillary-controlled methods include pore-network modelling that typicallyusesanidealizationoftheporegeometry[e.g.,20] andpore-morphologymethods[21].Intrinsically,thesemethods requireadditionalrulesortechniquesforvolumeconservation.

Implicit interface-tracking techniques are direct numerical simulation methods suitable to capillary-controlled (or, curvature-driven)displacements.Withinthiscategory,thetwomostwidelyusedapproachesarethelevelset(LS) method [22–25] andthevolumeoffluid(VOF) method [26].VOFmethodshavetheadvantagethat theyconservemassbutcalcu- lationofinterfacialpropertieslikenormalvectorsandcurvaturesarecumbersome.Forthelevelsetmethoditistheother wayaround.Calculatinginterfacialpropertiesisstraightforward,butthemethodhasknownissueswithmassconservation.

However,therearemethodsforimprovingmassconservationinLSmethods.Forinstance,SussmanandPuckett [27] devel- opedahybridmethodusingbothVOFandLSapproaches, Enrightetal. [28] usedahybridparticle-LSmethodtopreserve volume,andLuoetal. devisedamethodbasedonspectrallyrefined interfaceapproaches [29].SayeandSethian[30] pre- servedvolumesofeachphasewithintheirLS-basedVoronoiimplicitinterfacemethod,bymodifyingthephasepressuresto controlinterfacemotioninnormaldirection.

Inthiswork, we developa multiphaselevel setmethodwithlocalvolume conservation(MLS-LVC)tostudy capillary- controlledthree-phasedisplacementsattheporescale.Wehavepreviouslyimplementedglobalvolumeconservationinour LS[31] and MLS[32] methodsusingthe approachofSayeandSethian [30],which allowssimulations withpredictionof capillarypressureforspecifiedphasesaturations[33,32,34].Here,weextendthisapproachfurthertoconservetheindividual volumes of all isolated ganglia of one or more phases, while treating ganglia coalescence andsplitting through a new conservative volume redistributionalgorithm. The resultingMLS-LVC methodpredicts surface area, volume, andcapillary pressureforindividualfluidgangliaatequilibriumstates,thusprovidingcomplementaryinsightstopore-scaleexperiments addressingganglionbehaviour[e.g.,14,7].

Thepaperisorganizedasfollows:First,Section2presentsthemainfeaturesoftheMLSmethod.Then,wedescribethe LVC methodinthe MLSsetting(Section3) andexplain themodifications neededtoimplementit fortwo-phasesystems usinga standardLSapproach(Section 4). Section 5presentsimportantaspectsofour parallelimplementationusingLVC.

Section 6presentsthree-phaseMLS-LVC simulationsofdoubleandmultipledisplacements througha 2-Dsinusoidalpore channel whereseveralgangliacoalesce andsplit. Weinvestigateboth theeffects ofapplyingLVC toone andtwo phases aswellasthebehaviour ofvolumeconservationandcapillarypressureconvergenceofdisconnectedphasegangliaasgrid resolutionincreases.Thereafter,we applytheLVC algorithmtotheLSandMLSmethods tocarryout respectivetwo- and three-phasesimulationsofprimarydrainageandsubsequentgasinvasionsinasynthetic2-Dporousmediumandonasmall subvolumeof3-Dsandstone.Section7summarizesourworkwithmainconclusions.Futureworkwilltargetsimulationof three-phase residual saturations afterwater-alternate-gasinvasion cyclesforcomparison withexperimentsandempirical correlations.

2. Multiphaselevelset(MLS)method

Thelevelsetmethod [22–24] isa numericalschemeforpropagationofcurvesandsurfacesinaccordancewithagiven velocityfunction, F.Thisisaccomplishedbydefiningacurve orsurface implicitlyasthezerocontourofascalarfunction φ (^x),whereφ (^x)isthepositionvector,thatevolvesaccordingtotheequation:

φ

t

+

^F

· ∇φ =

⁰

.

(1)

Fig. 1.MLSrepresentation offluid/fluid/solidcontactlinesinaporethroatcontainingfluidphaseαandβbefore(left)andafter(right)enforcingthe projectionstepofLosassoetal. [35].(Forinterpretationofthecoloursinthefigures,thereaderisreferredtothewebversionofthisarticle.)

It iscustomarytosplit thevelocitytermintonormalandadvectivecomponents, F_N and F_A,respectively, whichleads to thefollowingequation:

φ

t

+

FN

|∇φ| +

FA

· ∇φ =

0

.

(2)

Astrengthofthelevelsetmethodisthatnormalvectornofaninterfaceanditscurvature

κ

iscalculatedfromderivatives ofthelevelsetfunctionsothatn= ∇φ/|∇φ|^and

κ

_{= ∇ ·} n.

For systems containing more than two fluid phases, we will use a multiphase level set (MLS) approach [34,32] that assignsonelevelsetfunctionandoneevolutionequationperfluidphase

α

.Thisyields

(φ

_α

)

t

+

^F_N^α

|∇φ| +

^Fα_A

· ∇φ

α

=

⁰

.

(3)

Weidentifydomainswithφ_α<0 aspartoffluidphase

α

.Thestaticsolidphaseisalsodefinedbyascalarfieldψ,where ψ <0 isidentified aspartofthesolid andψ >0 is porespace. Onechallenge withusingmultiplelevel setfunctionsis thepresenceofoverlapregions,wheretwoormorephasesarepresent,andvacuumregions,wherenophasesarepresent.

Toavoidtheseissuesweuseaprojection stepmethoddescribed byLosassoetal. [35],whichamountstoidentifying the twosmallestlevelsetfunctionsateachpointofthedomainandthensubtracttheaverageofthesetwo fromalllevelset functions.

We will investigate quasi-static, capillary-controlled multiphase systems [36,34,32], with well defined contactangles betweenfluid/fluidinterfacesandporewalls [31,34].ThebulkpartofthesystemisgovernedbytheYoung-Laplaceequation forthefluid/fluidinterfaces,

P_αβ

=

^pα

−

^pβ

= σ

_αβC_αβ

,

(4)

where pα andp_β arephasepressures, Pαβ isthecapillarypressure,

σ

αβ istheinterfacialtension,andCαβ istheinterface curvature,forinterfacesbetweenphases

α

andβ.

Thecontactangleθαβ,measuredthroughphaseβ,istheanglebetweenthe

α

β-interfaceandthesolidporewall.Itis givenbythefluid/fluidinterfacialtension,

σ

αβ,andthesolid/fluidinterfacialtensions

σ

sαand

σ

sβthroughYoung’sequation:

σ

_αβcos

(θ

_αβ

) = σ

sα

− σ

sβ

.

(5)

The MLSmethodassigns a surfacetension contribution

γ

_α andasolid/fluidintersectionangle β_α toeach fluidphase [34,32].The

γ

α’sarerelatedtothefluid-fluidinterfacialtensionsthroughtheidentity[36]

σ

αβ

= γ

α

+ γ

β

.

(6)

Wedefinethesolid-fluidinterfacial tensionsas

σ

sα=

γ

s−

γ

αcosβα,where

γ

sischosensothatall

σ

sα>0.Further,βα isaspecifiedanglebetweentheporewallandtheboundaryofphase

α

,definedbycosβα= ⁿα· ⁿs,whereⁿα andⁿsare thenormalvectorsderivedfromφ_α andψ,respectively,seeFig.1.ConsistentwithEq. (5),thecontactangleθ_αβ isgiven by

cos

θ

_αβ

=

ⁿβ

·

ⁿs

= γ

_βcos

β

_β

− γ

_αcos

β

_α

σ

_αβ

^.

⁽⁷⁾

Note that thisexpression is independentof

γ

s. Forgiven sets of

σ

_αβ andθ_αβ we use

γ

- and β-parameters that satisfy Eqs. (6) and(7) asinputtotheMLSsimulations.Intwo-phasesimulations,usingtwolevelsetfunctions,apossiblechoice fortheintersectionanglesisββ=θαβ andβα=^180◦−θαβ.

Thealgorithmusedtofindthelocalequilibriumofamultiphasesystem, whereeachphaseisgivenauniformpressure, isbasedonassigninganenergycost

γ

α foreachphaseboundarysothatthetotalenergyofthesystembecomesthesum ofbulk andsurfacecontributions foreachphase separately.Wehaveexplored thisapproachpreviously [32,34].Here,we

alsowanttoconsiderthatthepressureineachisolatedregionofeachphasecanbedifferent,sothatweneedtokeeptrack ofindividual regionsthat canchangeshapeandtopology intime.However, weshallassume thatall regionsofthesame phasestillhavethesamesurfacetensionandcontactangle.Thetotalenergyofthesystem [32] cannowbewrittenas

E

=

α

⎛

⎝

p_αH

(−φ

α

)

H

(ψ )

dV

+ γ

_α

δ(φ

_α

)|∇φ

α

|

H

(ψ )

dV

+ σ

sα

δ(ψ )|∇ψ|

H

(−φ

α

)

dV

⎞

⎠ ,

(8)

where is the computationaldomain, H(· · ·) ^{is theHeaviside stepfunction and} δ(· · ·)is the Diracdelta function.The velocityfunctionsthatareusedtominimizethisenergyare [32]

F^α_N

=

^H

(ψ )(

p_α

− γ

_α

κ

_α

) −

^H

(−ψ )

x S

(ψ ) γ

_α

|∇ψ |

^cos

β

_α

,

(9)

F^α_A

=

^H

(−ψ )

x S

(ψ ) γ

_α

_∇ ψ,

(10)

where

κ

_{α= ∇ ·}(∇φ_α/|∇φ_α|) is the curvature field calculated from phase

α

’s level set function, and S(· · ·) is the sign function.Wenotethatthepressurepα,andconsequently

κ

α,isnotonlydependentonphasebutalso,possibly,ondomain aswillbeexplainedinthenextsection.

WeuseanexplicitnumericalschemebasedonfinitedifferencestosolveEqs. (3),(9) and(10),asdescribedelsewhere [32,24].WiththeprojectionmethodofLosassoetal. [35] enforcedattheendofeach iteration-timestep,theequilibrium solutionssatisfy[32]

(

p_α

− γ

α

κ

α

) |∇φ

α

| =

p_β

− γ

β

κ

β

|∇φ

β

|

⁽¹¹⁾

inporespace,and

γ

_α

( ∇ φ

_α

· ∇ ψ −

^cos

β

_α

|∇ φ

_α

||∇ ψ | ) = γ

_β

_∇ φ

_β

· ∇ ψ −

^cos

β

_β

|∇ φ

_β

||∇ ψ |

⁽¹²⁾

atthe porewalls. Around

α

β-interfaces, we havethat φ_α= −φ_β,whichimplies ∇φ_α= −∇φ_β, |∇φ_α|= |∇φ_β|^and

κ

α=

−

κ

β.Thus,Eqs. (11) and (12) areconsistentwithEqs. (4) and (7).

3. Localvolumeconservation(LVC)method

Our goalis to develop a simple andlightweight algorithm that, to a goodapproximation, conserves the volume ofa phase andeach isolated region ofthat phase. When evolving thissystem, regions cancoalesce andsplit, sowe need to deviseanalgorithmthatcanrobustlyredistributevolumesbetweenregionsthatchangebothshapeandtopology.

Weintroducethefollowingnotation:Foragivenfluidphase

α

,letIα(t)⊂Nbetheenumerationofallisolatedregions ofthatphaseatiterationtimet,similartoan indicatorfunction.Then,a(t)⊂,wherea∈^Iα(t),isanisolatedregionin space. Toimplementvolumeconservationoftheisolated regionsweusethemethodintroducedbySayeandSethian [30]

wherethepressure p_a(t)ofregiona,a∈^Iα(t),atiteration-timet,isgivenas

p_a

(

t

) =

^V

(0)

a

(

t

) −

Va

(

t

)

t A_a

(

t

) .

(13)

Here, Va⁽⁰⁾(t)isthetarget volumeofregiona, Va(t) istheregionvolumeatiteration-timet, Aa(t)istheregionsurface area,andt istheiteration-timestep.Thepressure p_aonlyequalsthepressureinregiona whentheevolutionoftheLS fieldshasreachedits equilibriumstate.Inthisstationarystate, Va willdeviatefrom V_a⁽⁰⁾ byasmallamountgovernedby t.Hence,thismethodworksbyassigninganartificialnumericalcompressibilitytoeveryisolatedregion.Wenotethatt isnotaphysicaltimesothatitsdimensionsarechosenaccordingtoEq. (13),alsoinagreementwiththeiterativetimeunit inEq. (3).

WeidentifyisolatedregionswiththegrassfirealgorithmdescribedbyProdanovi ´cetal. [37].Toeachregionweassigna volume V_a⁽⁰⁾ sothatthetotaltargetvolumeofphase

α

isgivenas

V_α⁽⁰⁾

=

a∈^Iα

V_a⁽⁰⁾

.

(14)

Ourmethod requiresa strategyto assign volumestoregions atthe next time basedon thevolume ofregions atthe previous time step. For this purpose,we define extendedregions ^ext_a (t)⊃a(t),a∈^Iα(t), that completely partition the computationaldomain,i.e.,= ∪a∈^Iα(t)_a^ext(t).Weconstructtheseextendedregions^ext_a (t)byincludingallpointsnearest to a given region a(t) using a distance transform. In thiswork we haveused the method described by Rosenfeld and Pfaltz [38].Usingthisalgorithm,wecancalculatetheintegerdistanced(^x,a)fromavoxel(or,gridcell)atposition^xto theboundaryofeachregiona,a∈^Iα.Foreachvoxelx∈/a,a∈^Iα,wemarkitaspartoftheextendedregion_a^ext,given

Fig. 2.Illustrationofthevolumeredistributionalgorithmwithanexampleofdomaincoalescence(t2>t1)andsplitting(t1>t2).Calculationoftarget volumesfollowsEq. (16).Theverticallinesindicatethebordersbetweenextendedregions,whilethegrey areasshowtheregionsofconservedphase.

Thestippledlinesshowthepositionofinterfacesintheprevious/nexttimestep.Forthecasewithdomaincoalescence(t2>t1),weobtainthefollowing targetvolumesatt2:V₁⁽⁰⁾(t2)=^V1⁽⁰⁾(t1)and V₂⁽⁰⁾(t2)=^V2⁽⁰⁾(t1)+^V(3⁰⁾(t1).Forthecasewithdomainsplitting(t1>t2),weobtainthefollowingtarget volumes att1 (byswapping aand b, andt1 and t2,inEq. (16)): V₁⁽⁰⁾(t1)=^V1⁽⁰⁾(t2), V₂⁽⁰⁾(t1)=^f22V(2⁰⁾(t2)and V₃⁽⁰⁾(t1)=^f23V(2⁰⁾(t2).Here, f22=

|22(t2,t1)|/|2(t2)|^{and f23}= |23(t2,t1)|/|2(t2)|^.

thatd(x,a)≤^d(x,b)forallb∈^Iα (wherewehavethatd(x,)≥^{0).Incaseswherethesmallestdistancetoseveralother} domainsisequal,wechoose^{xtobepartoftheextendeddomainwithsmallestdomainnumber.}

Nowassumewehaveidentifiedthesets{a(t1):^a∈^Iα(t1)}^and{^exta (t1):^a∈^Iα(t1)}^attimet1,andthecorresponding sets{b(t2):^b∈^Iα(t2)}^and{^ext_b (t2):^b∈^Iα(t2)}^attimet2.Wedefinetheintersectionbetweenregionaattimet1 and theextendedregion^ext_b attimet₂>t₁as

ab

(

t1

,

t2

) =

a

(

t1

) ∩

^ext_b

(

t2

).

(15)

Then,weusetherelativevolumeofoverlapbetweenthet1 andt2regionsasaweightfortheredistributionofvolume.The targetvolumeassignedtoregionb(t2)isthusgivenas

V_b⁽⁰⁾

(

t2

) =

a∈^Iα(t1)

|

_ab

(

t₁

,

t₂

) |

|

a

(

t1

) |

^V

(0)

a

(

t1

),

(16)

where|· · · |^{denotesthevolumeofaregion,measuredinthenumberofvoxelsitcontains.Inourschemeweknowthatthe ext}regionscoverthewholecomputationaldomainsothatnovolumeislostinthevolumereassignment.Notethatwhen coupled withMLS iterationscheme, restrictions ontime steplead to relatively smallchanges inthe fluid configurations withinaniterationstep.ThisensuresthattargetdomainvolumesareredistributedcorrectlywithLVC.Fig.2illustratesthe volumeredistributionalgorithmforsimpleexampleswithdomaincoalescenceandsplitting.

ThepressuretermgiveninEq. (13) forregionb(t₂),b∈^Iα(t₂),isobtainedbyusingV_b⁽⁰⁾(t₂)asthetargetvolumeand byusingthefollowingvolumeintegralstoevaluatethevolume,V_b(t2),andarea, A_b(t2):

V_b

(

t2

) =

_b^ext(t2)

H

(ψ )

H

( − φ

_α

)

dV

,

(17)

Ab

(

t2

) =

_b^ext(t2)

H

(ψ )δ(φ

_α

)|∇φ

α

|

dV

.

(18)

Althoughtheaccuracyofthesevolumeandareacalculationscandecreasewhendomainsaresufficientlyclosetoeachother, typicallyonthelengthscaleoftheneighbourhoodusedbythenumericalstencil,thesumoftargetvolumes, Vα⁽⁰⁾,willstill be conserved.The region phasepressures p_b areassigned tothecorresponding extended regions ^ext_b , b∈^Iα, forusein the MLSiteration scheme.Thisensures that the numericalstencils usethe correctpressures inthe vicinityof fluid/fluid interfaces.

Toallowformassenteringand/or leavingthesystem,wetagall regionsofthesamephaseincontactwithinlet/outlet boundariesofthecomputationaldomain.Weassumeallsuchregionsareconnectedtothesamefluidreservoirandbelong to the continuous phase withthe same pressure.This phase pressure iseither assignedor calculatedby controlling the phasesaturation[32] (seealsoSection6).Regionstaggedasboundariesareexcludedfromthelocalconservationschemeas thepressureisgivenbythecontinuous phase.Oncenewisolatedregions detachfromthecontinuousphase,we calculate

their target volumes withEq. (17). Then, Eq. (13) implies that thephase pressures forsuch regions are zero inthe first iteration afterregiondetachment.However, insubsequentiterationsthesephasepressures eventually convergeto correct valuesasthesystemrelaxestowardequilibrium.

Forsimulationsonrealisticporegeometriesitisinevitable thattinyphasedomains containingonlyafewvoxelsarise.

This can lead to highandunrealistic region phase pressures calculated fromEq. (13), which in turncan make the MLS iterationschemeunstable.Inthiswork,weresolvethisissuebysettingthephasepressureofregionswithtargetvolumes smaller than a few grid-cell volumesequal to the pressure of the corresponding continuous phase. Specifically, we use 3(x)² (in2-D)and20(x)³ (in3-D)asvolumelimits.Althoughthisimpliesthat tinyregionscanvanishduring theMLS iterationsteps, oursimulations(seeSection 6)indicatethat theeffectonthevolume redistributionisnegligible.Wenote thatanalternativeapproach,thatwehavenotinvestigatedhere,wouldbetoexcludetinyregionscompletelyfromtheLVC algorithm.

3.1. Formulationofthealgorithm

WenowgivethemainstepsinvolvedintheimplementationofourLVCmethodintheMLSsetting.

AlgorithmSeq

1. Initialize the system: (i)Construct theset{a:^a∈^Iα(t0)} ^{(usinggrassfirealgorithm) withunique tagsa}∈^Iα(t0)as- signed to each isolated fluid region of the conserved phase

α

, and (ii) generate the accompanying set ofextended regions{^ext_a :^a∈^Iα(t0)}^.Settn=^t0.

Loop over iteration time (fromt_n tot_n+1):

2. Set pressures of conserved phase

α

equal to the continuous-phase pressure in regions connected to inlet or outlet boundaries. For the other isolated phase

α

regions, set pressures to pa, a∈^Iα(tn),using Eq. (13), enforcingvolume conservation.

3. Iterate the multiphase system fromtn totn+¹ withEqs. (3), (9) and (10),usingthephase pressuresfromstep2,and thenenforcetheprojectionstep[35].

4. Repartition the system attimestept_n+1,thatis,constructtheset{b:^b∈^Iα(t_n+1)}^{.Tagregionsasboundariesifthey} belongtothephaseconnectedtoinletoroutletboundaries.

5. Extend the partitioning{b(tn+¹)}^{toobtaintheset}{^ext_b (tn+¹)}^.

6. Redistribute volumes according to Eq. (16), usingt_n andt_n+1 astime step t₁ andt₂,respectively. Regions taggedas boundarybelongtothecontinuousphaseandareexcludedfromthevolumeredistribution.

7. Identify new regions that detached fromthe continuous phase and calculate their region target volumes based on Eq. (17).

8. Go to step 2 withn←ⁿ+^1.

Afterrepeatingthisschemeaspecifiednumberoftimes,wereinitializetheφα-functionstosigneddistancefunctionsby solving[39]

(φ

α

)

t

+

^S

(φ

α

)(|∇φ

α

| −

¹

) =

⁰

.

(19)

Theiteration-loopendswhenthesystemisinastationarystategivenby[32]

maxα

H

(ψ + )

H

(φ

ⁿ_α

+ )

H

( − φ

ⁿ_α

+ ) | φ

ⁿ_α

− φ

_α^m

|

H

(ψ + )

H

(φ

_αⁿ

+ )

H

(−φ

ⁿ_α

+ )

<

c

x

,

(20)

whereweusec=⁰.001 and =¹.5x,andnandmarethetwolastiterationstepswherereinitializationofφ_α occurred.

Thisendstateistreatedasthequasi-staticequilibriumofthesystem.Byincrementalchangesinboundaryconditionsand/or constraintsthesystemapproximatesthequasi-staticevolutionofamultiphasesystem.

4. Two-phasesystemasaspecialcase

A two-phase systemrequires onlya single level set functionφ that evolvesaccordingto Eq. (2).In thiscasewe can rescalethevelocitiesfromEqs. (9) and(10) bymeansofinterfacialtensiontoobtain[31]

FN

=

H

(ψ )(

C

− κ ) −

^H

( − ψ )

x S

(ψ )|∇ψ|

cos

β,

F_A

=

^H

(−ψ )

x S

(ψ )∇ψ,

(21)

whereβ=^180◦−θ.InEq. (21),C isa prescribedvalue forcapillarypressure dividedbyinterfacial tension,representing interfacecurvaturebyEq. (4).Withinterfacialtensioneliminatedfromtheevolutionequation,capillaryequilibriuminthe porespaceoccurswhenthecomputedcurvaturefield

κ

hasbecomeequaltotheprescribedvalueC aroundtheinterfaces

(i.e.,C=

κ

).NotethatintheevolutionequationsoftheMLSapproachwehavechosentousephasepressures ratherthan interface curvaturesC becausewe cannoteliminateinterfacial tensionsfromsimulationsofsystemscontainingmorethan twofluidphases.

Compared tothe MLS approach,applying LVC algorithm in a single level setframework requiresonly a few changes relatedtothesignofφ.Assumingφ <0 inphase

α

andφ >0 inphaseβ,wenowupdateinterfacecurvatureCaofisolated domainsasfollows:

Ca

(

t

) =

^sV

(0)

a

(

t

) −

^Va

(

t

)

t Aa

(

t

) ,

(22)

whereitis crucialthat s=1 fora∈Iα and s= −1 fora∈Iβ.Domain volumeandareacalculationsfollowEqs. (17) and (18),exceptfora∈^Iβ whereφ >0 sothatthevolumecalculationis

H(ψ)H(φ)dV. 5. ParallelaspectsoftheLVCalgorithm

We haveimplemented the levelset methodsdescribed in the previous sectionswithin the StructuredAdaptive Mesh Refinement Application Infrastructure (SAMRAI) software system [40]. SAMRAI is an object-oriented MPI andC++-based programmingframeworkthatprovidesanextensibleandscalableinfrastructureforparalleladaptivemeshrefinement(AMR) applications [41–43].

TheimplementationstrategyforstructuredAMRemployedinSAMRAIistheso-calledpatch-basedapproach,asopposed to the quad/octreeapproach which is the other typical alternative. In SAMRAI a givenprocessor can be associated with one orseveralpatches.WehaveimplementedAMR inourprevious levelsetmethodsfortwo- andthree-phasecapillary- controlled displacementsinporousmediaasdemonstratedinourrecentpapers [44,32].

While thepatchbased parallelapproach inSAMRAI can beeasily exploitedforthe levelset methodsin ourprevious works, itismore challengingwhen includingthe aboveLVC algorithm.Thepurposeofthissection isthusto explainthe implementation of certain important aspects of the parallel LVC algorithm. In the present paper we will only consider uniformgrids.TheinclusionofAMRintroducesadditionalchallengesandiscurrentlyworkinprogress.

5.1. Parallelalgorithm

We nowgive themain stepsinvolved intheimplementationofourLVC methodin theparallelMLS setting.This will closelyfollow thesteps fromAlgorithmSeq (Section 3.1), butwiththe inclusion ofadditionalsubsteps andcommentsto explain importantextensionsneededforparallelism.Wewilldiscusssome ofthesestepsinmoredetailinSection 5.2as wellasinAppendixA.

AlgorithmPar

1. Initialize the system:

(a) Oneachpatch:Usethegrassfirealgorithmtoidentifyfluidregionsoftheconservedphase

α

withlocaltagsunique toeachprocessor.

(b) Oneachprocessor:Gathertheinformationin(a)fromthecorrespondingpatchesandcommunicateinformationto themasterprocessor.

(c) On the masterprocessor: Constructunique globaltags a∈^Iα(t0)assignedto each isolated fluid regionusing the algorithmdescribed inAppendixA.Construct theset{a:^a∈^Iα(t0)}^{withthehelp ofthe}informationfromstep (a)above.Communicaterelevantinformationbacktoallprocessorsandpatches.

(d) Oneachpatch:Generatethesetofextendedregions{^exta :^a∈^Iα(t₀)}^.

(e) Oneachpatch:Modifythesetofextendedregions{^ext_a :^a∈^Iα(t0)}^{byusingthedistancefunctiond}(x,)together withpatchneighbourhoodinformation.

(f) Sett_n=^t0

Loop over iteration time (fromtn totn+¹):

2. Set pressures of conserved phase

α

equal to the continuous-phase pressure in regions connected to inlet or outlet boundaries. For theother isolated phase

α

regions, setpressures to pa, a∈^Iα(tn), usingEq. (13) (enforcingvolume conservation)asfollows.

(a) Oneachpatch:PerformcomputationsaccordingtoEq. (17) andEq. (18).

(b) Oneachprocessor:Similarapproachasinstep1(b).

(c) Onthemasterprocessor:PerformEq. (13) andcommunicaterelevantinformationbacktoallprocessorsandpatches.

3. Iterate the multiphase system fromt_n tot_n+1 withEqs. (3), (9) and (10),usingthephase pressuresfromstep2,and thenenforcetheprojectionstep[35].Thisisdoneinparallelusingthepatch-basedparallelisminSAMRAI.

4. Repartition the system attimesteptn+¹,thatis,constructtheset{b:^b∈^Iα(tn+¹)}^{.Tagregionsasboundariesifthey} belong tothephase connectedtoinlet oroutletboundaries. Inparallel thisis doneina similarmanner asinstep 1 (a)-(c)above.

5. Extend the partitioning{b(t_n+1)}^{toobtaintheset}{_b^ext(t_n+1)}^{.Inparallelthisisdoneinasimilarmannerasinstep} 1(d)-(e)above.

Fig. 3.Anexampleofthedecompositionofthecomputationaldomainwheretheleft,middleandrightcolumnsshowtheresultsfor1,8and16processors, respectively.Eachprocessorcontainsonlyonepatchintheseexamples.Thetoprowshowseachprocessor’slocaltagnumbersforthedomains.Themiddle rowshowsthedistancemapcalculatedforindividualpatches.Thedistancemapconsistsoftheglobaltagsa∈Iαinsidetheregionsaandintheexterior itcontainsthenegativedistancetothevariousregions.Thebottomrowshowstheextendedglobalmap,wheregreenistheextensionofdomain1,redis theextensionofdomain2andblueshowspatcheswithoutdomains.

6. Redistribute volumes

(a) Oneachpatch:PerformcomputationsaccordingtoEq. (15) andEq. (16) (denominator).

(b) Oneachprocessor:Similarapproachasinstep1(b).

(c) Onthemasterprocessor:PerformEq. (16),usingtnandtn+¹ astimestept1andt2,respectively.Regionstaggedas boundarybelong tothe continuous phaseandare excluded fromthevolume redistribution.Finally,communicate relevantinformationbacktoallprocessorsandpatches.

7. Restofalgorithmresemblessteps7and8inAlgorithmSeq.

5.2. Simpleillustrationofsomeaspectsoftheparallelalgorithm

Here,weelaborateonessentialpartsoftheparallel algorithmgivenaboveusingthefollowingsimpletwo-dimensional example containing two isolated regions (see Fig. 3). Twodomains need tobe identifiedin thissimple example,and in thetop rowofthefigurewe showlocallabels(tags)aswellasthevariouspatchesforone,eightandsixteenprocessors, respectively. Obviously,local andgloballabelscoincide inthe caseofa single processor.Note that inparallel computing localtagsuniquetoeachprocessorareidentifiedbyperformingthegrassfirealgorithmoneachpatch.Aftercommunicating thevariousresultstothemasterprocessorglobaltagshavetobeassociatedwitheachisolatedregion,sinceagivenregion cancoverseveralpatches(andprocessors).ThisisdonebyusingthealgorithmdescribedinAppendixA.

InthemiddlerowofFig.3theregionsaredenotedby1and2,withcorrespondinggloballabels1 and2,respectively.

Outsideofthetwodomainsinthesethreesubfiguresweshowthedistancemap,whichisusedforcomputingtheextended regions. Thisisjustthenegative integerdistancecomputedwithrespecttothenearest domain(i.e.−^d(x,)).Inpatches

Table 1

Weakscalingresultsforvariouscasesofconservedphasesandnumberofisolateddomains.

Absolute timing (s) for weak scalability tests

# cores Without LVC LVC1phase (32domains)

LVC2phases (64domains)

LVC1phase (108domains)

LVC2phases (216domains)

1 11.97 13.77 15.61 14.36 15.85

8 12.78 14.48 16.15 16.33 19.8

64 16.41 19.72 23.07 22.48 28.07

512 17.33 21.14 24.83 27.02 36.19

without domains thisnegative integer distanceis simply putto the negative sum ofthe numberof grid points in each coordinatedirection.

The corresponding extended regions are shown in Fig. 3, bottom row. In the parallel computing process we iterate through the patches connected to a processor and perform the “extension algorithm” (described in Section 3) for each patch. However, when thisis completed we must yet againiterate through the patches andpossibly correct the results from the extension based on patch neighbourhood information. In essence we compare the “negative distance map” of neighbouringpatchesandgivepreferencetotheshortestnegativedistancewhenupdatingthe“extensionmap”.

AsseeninSection3,volumeandareacomputationsofthedifferentregionsarecentraltotheLVCalgorithm.Weperform suchcomputationsoneachpatchinparallel.Differentcontributionsarethengatheredoneachprocessorandcommunicated to themasterprocessor wheretheycan be puttogetherwiththe helpofthe globaldomain identificationtagsdiscussed above.

Thevolume redistributionalgorithmdescribed intheprevioussection isalsoatfirstdoneoneach patchinparallel in ordertofindthevarious domainsizesandcorrespondingintersections (“hits”),bothrepresentedasintegers,betweenold andnewdomains.Communicationsareagainneededtogatherthevariousdatastructuresonthemasterprocessor.Aftera somewhatinvolvedprocessinvolvingmatrix-likestructuresthenewdomainvolumescanbecomputedasinthesequential algorithm (seeSection 3).Finally,thepressures inthedifferentregions canbe computed,andthen communicatedto the variousprocessorsforuseinthelevelsetalgorithm.

5.3. ParallelscalabilityoftheLVCalgorithm

TheparallelversionoftheLVCalgorithmintroducesadditionalMPIcommunicationintheLS/MLScodes.Inthissection we thusinvestigatetheparallelperformance oftheLVC algorithmusingtheMLSmethod.Wepresenttestsforbothweak andstrongscalabilityoftheparallelcodebasedonsimulationsofcurvature-drivenevolutionofisolateddomainstowarda steadystatewherethedomainsattainsphericalshapes.Inmostrealisticapplicationsthesizeofthedatastructuresinvolved intheMPIcommunicationwillbe relativelymodest,butit willgrowwiththenumberofisolateddomainsthat needsto beidentified.Sincetheparallelscalabilitythusclearlywilldependonthisnumber,weshow resultswherethenumberof isolateddomainsrangesfrom0 to216,representingtypicalsimulationsinrealisticporegeometries.

In the scaling teststhe distributions of thephases subjectedto LVC (thatis, the oiland gasphases) constitute fixed numbers ofisolated cubic domains, whereas theremaining phase (water)is continuous.We consider two such cases,as showninFig.4(a).Thefirstcasecontains64 (=^{43)domains(32ofeachphase)andthesecond casecontains 216 (}=⁶³⁾ domains (108ofeach phase). Weexplore scalingbehaviour ineach ofthesecases forthesettingswithoutLVC andwith LVC applied to one (i.e.,oil) andtwo (i.e., oil andgas) phases. To maintain control of the number ofisolated domains presentinthesystemthecomputationaldomainconsistsentirelyofporespace.Intheweakscalingtestswefirstsimulate acomputationaldomainwith40³ gridcellsusingasingleprocessor.Wethenrefinethisconfigurationthreetimes(keeping thenumberofisolateddomainsconstant)sothattheproblemsizeandnumberofprocessorsusedinthesimulationsboth increase witha factorof8foreach refinement.Wethusendupwithaproblemcontaining320³ grid cellsinthecaseof 512processors.Inthestrongscalingtestsweperformsimulationsonthelargestproblemsize(320³ gridcells)using128, 256,512 and1024 processors.For512 processors,thesimulationsareidenticaltothecorrespondingsimulationsfromthe weakscalabilitytest.

Intheweakandstrongscalingtestseachsimulationexecutes20iterationstepswitheachofthethreefluidLSevolution equations, during which three reinitialization solves (each consisting of 10 iterations) occur after the 10th and20th LS iterationsteps,respectively.Inaddition,attheendofeachLSiteration-timestep,thesimulationsalsoincludetheprojection step[35] and,ifapplicable,theLVCalgorithm.Overall,thethree-phasesimulationscarryout120iterations,ofwhich60are main LSiterationsand60arereinitializationiterations. Thesefractionsarerepresentative fora typicalusage ofourcode.

All computationswereperformedonthe supercomputerFramprovided byUNINETTSigma2–the NationalInfrastructure forHighPerformanceComputingandDataStorageinNorway.

Results forthe weak scalingare presented inTable 1.With aperfect (ideal)parallel speedupthe absolutetiming for each ofthesesimulations shouldequal thetiming obtainedforthecorresponding simulation(i.e., withanequal number ofconservedphasesandisolateddomains)onasingleprocessor. InthecasewithoutLVC weobserveavery goodparallel performance.However,inouropiniontheresultsarealsosatisfactorywhenweincludeLVCononeortwophases,evenfor

Table 2

Strongscalingresultsforvariouscasesofconservedphasesandnumberofisolateddomains.

Absolute timing (s) for strong scalability tests

# cores Without LVC LVC1phase (32domains)

LVC2phases (64domains)

LVC1phase (108domains)

LVC2phases (216domains)

128 61.26 71.28 80.99 74.94 87.81

256 31.88 38.4 44.43 43.21 54.14

512 17.33 21.14 24.83 27.02 36.19

1024 11.04 14.0 17.2 19.77 27.85

Fig. 4.(a)Configurationsofthedisconnectedphases(oilinredandgasingreen)usedinthescalingtests.Eachphasecontains32 domains(top)and 108 domains(bottom).Theintermediatecontinuousphaseistransparent.SimulationswithLVCappliedtoonephaseconservethe reddomains,while simulationswithLVCappliedtotwophasesconservebothredandgreendomains.(b)Speedupnormalizedto128computingcores(processors)asa functionofnumberofcoresfromthestrongscalingtests.

ahighnumberofisolateddomains.Eventhoughthetimingsobtainedwith216domainsmayseemabitfarfromidealfor 64 and512 processors,theystillrepresentefficientandusefulhighperformancecomputations.

TheresultsfromthestrongscalingtestsarepresentedinTable2andFig.4(b).Weobservethattheparallelspeedupis very goodforthecasewithoutLVC.ForthecaseswithLVC theparallelspeedupdecreases withanincreasing numberof isolated domainspresentinthesystem. Thisisduetoparallelcommunicationintroducedby theLVCalgorithm.However, the parallel speedupobtained isstill absolutely acceptable anddefinitely demonstrates that our codesrun efficientlyon modernsupercomputers.

6. Simulations

In the numericalmodel given by Eqs. (3), (9) and (10), we introduce thedimensionless variables

γ

ˆα =

γ

α/

γ

^∗, ^pˆα= pα/p^∗, φˆα=φα/L^∗, and ^xˆ=^x/L^∗, where L^∗ and

γ

^∗ are the characteristic length and surface tension, respectively, and p^∗=

γ

^∗/L^∗ isthecharacteristicpressurefollowingYoung-Laplace’sEq. (4).It thenfollowsthatthedimensionlessformof Eq. (3) is unchangedwiththeiterativetime unitchosen as(L^∗)²/

γ

^∗.Clearly, thiswillnot affectthecapillaryequilibrium solutions.Thus,hereafterwetreat

γ

α,pα,φα andxasdimensionlessvariables.

WedemonstratetheLVC algorithmondifferentporegeometriesusinga fluidsystemofgas(g), oil(o)andwater(w) withphysicallyrealisticinterfacialtensions

σ

g w=⁰.03,

σ

ow=⁰.02 and

σ

go=⁰.011 N/m.Thecharacteristicsurfacetension

γ

^∗ is0.01 N/m,whichtogether withEq. (6) impliesthat

γ

w=¹.95,

γ

o=⁰.05 and

γ

g=¹.05.Wewillalsoensurethatall nonzero

γ

α aresufficientlylarger thanx,toobtaincorrectbehaviour inthelimit

γ

α→^{0.Foradiscussion ofvanishing} viscositysolutionswithlevelsetmethodswereferto[36].

Forathree-phasefluidsysteminthermodynamicequilibriumwithasolidthethreecontactanglesθ_αβ cannotbespeci- fiedindependently.Instead,Eq. (5) (or,equivalentlyEq. (7))impliesthattheinterfacialtensionsandcontactanglessatisfy

σ

g wcos

θ

g w

= σ

gocos

θ

go

+ σ

owcos

θ

ow

.

(23) In this work, we specify θow and calculate θgo and θg w from linear cosθgo(cosθow)- and cosθg w(cosθow)-relationships [45] that satisfy Eq. (23). Ourmodel uses a corresponding set of intersection angles β_α that are consistent with these

relationships[32]:

β

w

= θ

ow

, β

o

=

^180◦

− θ

ow

,

and

β

g

=

^180◦

.

(24)

6.1. 2-Dsinusoidalpore

Weperformsimulationsonanidealized2-Dporethroatgeometrytoinvestigateeffectsofgridresolutiononconservation ofvolume (areain2-D)andpressureconvergence ofisolateddomains. Thecharacteristicquantities are L^∗=¹×^{10−4 m} andp^∗=

γ

^∗/L^∗=^{100 N/m2.In}dimensionlessvariables,weconsiderthesinusoidalporethroat

_P

= { (

x

,

y

) :

⁰

≤

^x

≤

²

,

y₂

≤

^y

≤

^y1

} ,

(25)

where

y1

(

x

) =

¹ 100

52

+

^{5 sin}

25

9

π

x

+

^{8 sin}

50

9

π

x

,

y₂

(

x

) =

¹ 100

17

−

^{5 sin}

25

9

π

x

−

^{8 sin}

50

9

π

x

.

(26)

An analytic signed distance function for ψ describes the idealized pore geometry in the simulations. The initial φα- configurationisgivenbyfivecircularoildropssurroundedbywaterwithdifferentradiiandcentresatdifferentpositionson theaxisofsymmetryalongthex-directionoftheporethroat.Someofthesecirclesoverlapwithsolid,allowingforpossible contactangleformationattheporewalls.Initially,gasispresentinalayerwiththickness3xattheleftboundary(x=^0).

Asaturation-controlledapproachsimulatesquasi-staticgasinvasionfromx=^{0 whilethepressureofthecontinuousoil} andwaterphases, p_o and p_w,arekept constantduringthesimulation.Thisamountstosolving Eq. (3) foraseriesofgas saturations Sg,k=^Sg,1+(k−¹)Sg,k=¹,. . . ,m, thatcorresponds to target gasvolumes V_g^target_,k =^Sg,kVP, where VP is porevolume.Foreachsaturation S_g,k wesolveEq. (3) whileupdating p_gineachlevelsetiterationstepbasedon

p_g

(

t

) =

^V

target g,k

−

^Vg

(

t

)

t A_g

(

t

) .

(27)

We useinitialgassaturation Sg,1=⁰.05 andsaturationstepSg=⁰.02.The convergedfluidconfigurationfromthelast saturation step is the initial configurationfor the next. Notethat the saturation-controlled approach mimics quasi-static invasionatconstantlow flowrate,whileitpredicts capillarypressure ateach stationarystate [32,33].Pressure-controlled displacement,whichisthe otheralternative,predicts saturations foreach stationarystategivenaseriesoffixed capillary pressures.Fluidscanenterandleavethesystemattheinletandoutletboundaries(thatis,x=^{0 orx}=^2).

We performsimulationswithx=⁰.02, 0.01,0.005,and0.0025.Thecorresponding reinitializationfrequencies ofφ_α are 5,10,25and50time-iterationsteps.Forx=⁰.01,weuse

γ

α asstatedaboveandcontinuousphasepressurespw=^1, po=^{10,and p}g obtainedfromEq. (27).Fortheothergridresolutions,we multiplyEq. (3) by100xsothatthesmallest

γ

-value(

γ

o) remainsequalto5x.Notethat westillsolvethedimensionlessequation,thistimewithiterativetime unit 100x(L^∗)²/

γ

^∗.Whilethiswillimpactthespeedofgradientdescent (i.e.,numberofiterations),thecapillaryequilibrium solutionsremainunaffected.

AtstationarystatesweestimatetherelativeerrorE[^Va]^{ofsimulatedvolumesV}awithrespecttotargetvolumesV_a⁽⁰⁾for allindividualdomainsa∈^Iα,aswellastherelativeerrorE[^Vα]^{ofsimulatedtotalvolume V}α=

a∈^IαV_a withrespectto totaltargetvolumeVα⁽⁰⁾givenbyEq. (14).Fortherelativeerrorsofthetotalvolumewealsoestimatetheaverageabsolute value overthe numberofequilibriumstatesm, E_av[^Vα]=_m

i=1|^E[^Vα]i|/^{m.Weevaluate domainpressures, p}a,a∈^Iα,by means of their pressure difference to the other two continuous phases (i.e., local capillary pressure), which represents interfacecurvaturewhenscaledbyinterfacialtension,seeEq. (4).

6.1.1. Doubledisplacementwithdomainfragmentationandcoalescence

First we investigateLVC of oilonly, while waterpressure isconstant. Hence, we assume sub-resolution wetting films along thepore wallsconnect thewaterto inlet/outletboundaries. Forthispurpose weuse θow=^20◦.Calculationof the other contact angles basedon Eqs. (7) and (24) yieldsθgo=²⁴.2^◦ andθg w=¹⁶.1^◦. We wish to investigatedouble dis- placements where a continuous invading phase (e.g., gas) displaces a disconnected phase cluster (e.g.,oil) that in turn displaces acontinuous thirdphase(e.g.,water).Pore-scale experimentshaveshownthat thismechanismcanmobilizeoil incapillary-dominatedthree-phaseflow[4,5,15,16].

Inoursimulation example,gasdisplacesan oildropthat coalesceswithstaticoildropsaheadofthe frontaloil/water interfaceduringthemotion,whileoil-dropfragmentationformstrappedoillayersinporecornersbehindthetrailinggas/oil interface.Fig.5showsonlysmallgridresolutioneffectsonthestationaryfluidconfigurationsfromthesegas-invasionsimu- lations,exceptforoildropcoalescenceeventsthatinsomecasesoccuratslightlydifferentgassaturations.Asexpected,the capillarypressures P_aw,a∈^Io,betweenstaticoildropsandsurroundingwaterareconstant anddisplay excellentconver- gencebehaviourwithdecreasingx,seeFig.6(a).Fortheoildropdisplacedbygas,thelocalcapillarypressures Paw and