Multiscale gradient computation for flow in heterogeneous porous media

Multiscale gradient computation for flow in heterogeneous porous media

Jesus de Moraes, Rafael; Rodrigues, José R P; Hajibeygi, Hadi; Jansen, Jan Dirk

DOI

10.1016/j.jcp.2017.02.024

Publication date

2017

Document Version

Final published version

Published in

Journal of Computational Physics

Citation (APA)

Jesus de Moraes, R., Rodrigues, J. R. P., Hajibeygi, H., & Jansen, J. D. (2017). Multiscale gradient

computation for flow in heterogeneous porous media. Journal of Computational Physics, 336, 644-663.

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

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.

‘You share, we take care!’ – Taverne project

https://www.openaccess.nl/en/you-share-we-take-care

Otherwise as indicated in the copyright section: the publisher

is the copyright holder of this work and the author uses the

Dutch legislation to make this work public.

Contents lists available atScienceDirect

Journal

of

Computational

Physics

www.elsevier.com/locate/jcp

Multiscale

gradient

computation

for

flow

in

heterogeneous

porous

media

Rafael J. de Moraes

a

,

b

,

∗

,

José R.P. Rodrigues

b

,

Hadi Hajibeygi

a

,

Jan Dirk Jansen

a

a_{DepartmentofGeoscienceandEngineering,FacultyofCivilEngineeringandGeosciences,DelftUniversityofTechnology,P.O. Box5048,2600,}

Netherlands

b_{PetrobrasResearchandDevelopmentCenter–CENPES,Av.HorácioMacedo950,CidadeUniversitária,RiodeJaneiro,RJ21941-915,Brazil}

a

r

t

i

c

l

e

i

n

f

o

a

b

s

t

r

a

c

t

Articlehistory: Received13July2016

Receivedinrevisedform6February2017 Accepted7February2017

Availableonline14February2017 Keywords: Gradient-basedoptimization Multiscalemethods Directmethod Adjointmethod Automaticdifferentiation

Anefficientmultiscale(MS)gradientcomputationmethodforsubsurfaceflowmanagement andoptimizationisintroduced.Thegeneral,algebraicframeworkallowsforthecalculation ofgradients usingboth theDirect and Adjoint derivative methods. Theframework also allows for the utilization of any MS formulation that can be algebraically expressed in termsof a restriction and a prolongation operator. This is achieved via an implicit differentiationformulation.Theapproachfavorsalgorithmsformultiplyingthesensitivity matrixanditstransposewitharbitraryvectors.Thisprovidesaflexiblewayofcomputing gradients in a form suitable for any given gradient-based optimization algorithm. No assumption w.r.t. the nature of the problem or specific optimization parameters is made. Therefore, the framework can be applied to any gradient-based study. In the implementation,extrapartialderivativeinformationrequiredbythegradientcomputation iscomputed viaautomaticdifferentiation.A detailedutilization oftheframework using theMSFiniteVolume(MSFV) simulationtechniqueispresented.Numerical experiments are performed to demonstrate the accuracy of the method compared to a fine-scale simulator.Inaddition,an asymptoticanalysisis presented toprovidean estimateofits computational complexity. The investigationsshow thatthe presented methodcasts an accurateand efficientMSgradientcomputationstrategythatcanbesuccessfullyutilized innext-generationreservoirmanagementstudies.

©2017ElsevierInc.Allrightsreserved.

1. Introduction

Model-based reservoirmanagement techniquestypically relyon theinformationprovided by derivatives.Forinstance, in sensitivityanalysisstudies, derivativescanbe directlyusedtoidentifythe mostinfluentialparameters inthereservoir response. Also, derivativeinformation canbe utilizedinhistorymatching [1]andcontroloptimization [2]studies,where gradient-basedoptimizationtechniquesareemployedintheminimization/maximizationofanobjectivefunction.

Thesetypesofmodel-basedreservoirmanagementstudiesare computationallydemanding.Theyrequiremultiple eval-uations of the reservoir modelin order to compute its response underthe influence ofdifferent inputs. Reduced-order modeling(ROM)techniqueshavebeenemployedtoreducethecomputationalcostofthereservoirresponseevaluation. In

*

Correspondingauthorat:DepartmentofGeoscienceandEngineering,Faculty ofCivilEngineeringandGeosciences,DelftUniversityofTechnology, P.O. Box5048,2600,Netherlands.

E-mailaddresses:r.moraes@tudelft.nl(R.J.deMoraes),jrprodrigues@petrobras.com.br(J.R.P. Rodrigues),h.hajibeygi@tudelft.nl(H. Hajibeygi),

j.d.jansen@tudelft.nl(J.D. Jansen).

http://dx.doi.org/10.1016/j.jcp.2017.02.024

sensitivityanalysisstudies,responsesurfacemodelsanddesignofexperimentsareoftenusedtoreducethecomputational costs (see,e.g.,[3]).Inhistorymatchingandoptimizationstudies,techniqueslikeupscaling[4],streamlinesimulation[5], and proper orthogonal decomposition [6] are employed to create reservoir models that are faster to evaluate. However, ROMandupscalingmethods usuallydonot provideaccurateenough systemresponsesduetoexcessivesimplifications of thefluid-rockphysicsandheterogeneousgeologicalproperties.Toresolvethischallenge,Multiscale(MS)methodshavebeen developed[7–9].

MS methodssolvea coarsersimulationmodel,thusincreasing thecomputational speed,whileresolving thefinescale heterogeneities. Note that the specific multiscale methods addressed here, map between fine and coarse grids that are bothatcontinuum(Darcy)scale,butwithdifferentcomputationalgridresolutions.Moreover,themapbetweenthenested fine andcoarse grids is developed by using multiscale basis functions[7].The basis functions are localsolutions of the fine-scale equation, which are adaptively updated [10,11] andallow the MS coarse systemto account forthe fine-scale heterogeneities(which typicallydonothaveseparationofscale). NotethatincontrastwithMultiGrid(MG) methods,MS methodsarenotdevelopedaslinearsolvers,butaremostefficientiftheyareused–similarasinthispaper–toprovide approximateconservativefine-scalesolutions(crucialformultiphasesystems).MSmethodsarefoundefficientandaccurate forlarge-scale reservoir models [12,13].Important inthis classof MS methods (comparedto upscaling methods)is that thecoarse-scale solutionsare mappedontothe originalfinescale, usingthesamebasis functions.As such,errorscan be calculated againstthe fine-scalereferencesystems (andnot upscaledaveraged ones).This allows forthedevelopment of convergentiterativeMSprocedures[14–16].RecentdevelopmentsincludeMSformulationsforfracturedmedia[17,18]with compositional effects[19,12,13]andcomplex wellconfigurations [20] andgravitationaleffects[21].Inaddition, algebraic formulation of the method has made it convenient to be integratedwithin existing simulators using structured [22,23] and unstructured [24–26] grids. The method hasbeen also extended to fully-implicit formulations where all unknowns crossmultipledynamically-definedscales[27].Althoughthesedevelopmentsarefoundefficient,theyaremainlylimitedto forwardsimulationmodeling.InthispaperthefocusisontheuseofMSmethodswithinreservoirmanagementworkflows. Reservoirmanagementtechniques includeoptimizationalgorithms,inwhich calculationofderivativesplays an impor-tant role. The classical approaches for calculation of derivatives are either computationally expensiveor inaccurate. For instance,numericaldifferentiation(see,e.g., [28,29])suffersfromdiscretizationapproximationsandtruncationerrors,and is impractical when the number of parameters is large. Alternatively, analytical methods – Direct Methods (or Gradient Simulators)[30]andAdjointMethods[31,32]–canprovideaccurateandefficientderivativesunderappropriate conditions (tobefurtherdiscussedintheSection2).However,theuseofanalyticalmethodshasnotbeenextensivelyadoptedmainly becausetheyarecode-intrusiveandrequireasubstantialimplementationeffort.Onthisissue,automaticdifferentiationcan partlyalleviate theburdenofcomputingderivativeinformation[33].Additionally,mostcommercialsimulatorsdonot pro-videanalyticalderivativecapabilitiesnordotheyprovideaccesstoextendtheirfunctionalityinthisdirection.Partiallydue tothesedrawbacks,ensemblemethods suchastheEnsembleKalmanFilter(EnKF)havebecomeverypopularinthedata assimilationcommunity[34].Similarly,stochasticapproximategradienttechniquessuchasensembleoptimization(EnOpt) andthe stochastic simplex approximategradient (StoSAG) methodare increasingly beingused forlife cycle optimization [35,36].Thesemethods,however,byconstruction,provideanapproximationofthegradient.

Multiscalegradientcomputation hasbeenstudiedinthe past.AMultiscalefinite-volumeAdjoint (MSADJ)methodhas been applied to sensitivity computation [37], where the global adjoint problem is solved via a set of coupled sub-grid problemsdescribed ata coarserscale. Thecoarse-scale sensitivitiesare then interpolatedtothe localfinegrid by recon-structingthelocalvariabilityofthemodelparameterswiththeaidofsolvingembeddedadjointsub-problems.Inafollow uppaper[38],theMSADJmethodwasefficientlyappliedtoinverseproblemsofsingle-phaseflowinheterogeneousporous media.Also,an efficientMultiscaleMixedFiniteElementmethodhasbeendevelopedformultiphaseadjointformulations, wherebothpressureandsaturationsaresolvedatthecoarsescale[39].IncontrasttoMSADJ,thismethoddidnotrequire fine-scalequantities.

The present development introduces a mathematical framework to compute sensitivities (gradients) in a multiscale strategy. The framework enables thesame computational efficiencyasexisting multiscale methods [37–39]. However, its formulationallows forfull flexibility withrespectto the typesofgradient informationthat are requiredby thedifferent model-basedreservoirmanagementstudies.Thisisachievedviaanimplicitdifferentiationstrategy,asopposedtothemore traditionalLagrangemultiplier formulation.Also,theformulationnaturallyprovides notonlytheAdjointMethod,butalso theDirectMethod.Itisimportanttonotethat,althoughinthisworkthemultiscalefinitevolume(MSFV)isbeingstudied, theproposedMS-gradientmethodisnotrestrictedtoaspecificMSmethod.Instead,itcanbeutilizedincombinationwith anyMS(andmulti-level)strategywhichisexpressedintermsofrestrictionandprolongationoperators.

Thispaperisstructuredasfollows.First,themultiscalegradientcomputationmethodisderivedbasedontheMS reser-voirmodelequationsandthe respectivemodelresponses. Thecomputation ofthe requiredprolongation(matrix ofbasis functions)operatorderivativesisdevelopedwithin theMultiscaleFiniteVolume(MSFV)formulation. Computational com-plexityofthemethodisalsodiscussedfromatheoreticalpointofviewviaasymptoticanalysisofthealgorithms.Thereafter, theNumerical Experimentssection describesa systematicinvestigationofthe validation,robustness, andaccuracyofthe MS-gradientmethodfortestcasesofincreasingcomplexity.Becausetheproposedmethodisquitefundamental,the exper-imentsareaimedatevaluatingthegradientcomputationitself,ratherthananyspecificapplication.

2. Derivationofthemultiscalegradientcomputationmathematicalframework

2.1. Forwardsimulationmodel

Thederivationoftheforwardsimulationmodel’sresponsewithrespecttotheparametersstartsbyfirstlydescribingits equationsinageneric,purely algebraicform.Theanalysisisrestrictedtoquasi-steadystateproblems.Inthatcase,theset ofequationsthatdescribestheforwardsimulationatthefinescalecanbealgebraicallyexpressed,withoutanyassumption regardingtheunderlyingphysicalmodel,as

gF

(

x

,

θ) =

0

,

(1)

where gF

: R

NF

× R

Nθ

→ R

NF representsthe setofalgebraic simulatorequations, x

∈ R

NF is thestate vector (which,for single-phase flow, contains the grid block pressures),

θ ∈ R

Nθ _{is the vector of parameters (which typically contains grid}

blockpermeabilitiesandporositiesbutpossiblyalsootherparameterssuchasfaulttransmissibilitymultipliersorstructural parameters),andthe subscript F refersto‘fine scale’.Thereare NF fine-scalecells and Nθ parameters. Eq.(1)implicitly assumesadependencyofthestatevectorx ontheparameters

θ

,i.e.

x

=

x

(

θ) .

(2)

Oncethemodelstateisdetermined,theobservable responsesoftheforwardmodelare computed.Theforwardmodel responses (typicallywellborepressuresand/or rates)maynotonlydependonthemodelstate, butalsoontheparameters themselves,andcanbeexpressedas

yF

=

hF

(

x

,

θ) ,

(3)

wherehF

: R

NF

× R

Nθ

→ R

Ny representstheoutputequations[40].ItisassumedthatgF canbedescribedas

gF

(

x

,

θ) =

A

(

θ)

x

−

q

(

θ) ,

(4)

whereA(

θ)

isan NF

×

NF matrixandq(

θ)

isan NF vector.

2.2. Algebraicmultiscaleformulationofflowinheterogeneousporousmedia

MSmethodsprovideaccurateandefficientsolutionstotheflowequationsbyincorporatingthefine-scaleheterogeneities in a coarse-scale operator [8]. Thisis achieved by basis functions, which are local solutions of the governing equations withoutright-hand-side(RHS)terms,subjecttoapproximatelocalboundaryconditions.Theselocalbasisfunctionsconstruct the prolongation(interpolation)operator, P

=

P(

θ)

,whichisan NF

×

NC matrix thatmaps(interpolates)the coarse-scale solutiontothefine-scaleresolution.Forthepurposeofthedevelopmentpresentedhere,therequiredbasicknowledgeabout themultiscalestrategy isthatacoarse-scalesystemcanbealgebraicallydescribed oncetherestrictionoperator,R

=

R(

θ)

, isdefinedasan NC

×

NF matrixwhichmapsthefinescaletothecoarsescale(moreinformationcanbefoundin[41,22]). Let x

˘

∈ R

NC _{bethe coarsescale solution(N}

C



NF),R

=

R(

θ)

be an NC

×

NF matrixmapping fromthe finescaleto the coarsescaleandP

=

P(

θ)

bean NF

×

NC matrixmappingfromthecoarsescaletothefinescale.x is

˘

obtainedbysolving

˘

g

= (

RAP

)

x

˘

− (

Rq

)

= ˘

Ax

˘

− ˘

q

= ˘

0

.

(5) Finally,theapproximatedfine-scalesolutionx isobtainedbyinterpolatingthecoarsescalesolutionx,

˘

i.e.,

x

=

Px

˘

.

(6)

2.3. Derivativecalculationofsimulatorresponses

Inordertoderiveanexpressiontocomputethemultiscalederivativeswithrespecttotheparameters,animplicit differ-entiationschemeisfollowed[42,43].TheimplicitdifferentiationschemefacilitatesboththeDirectandtheAdjointmethods, incontrasttotheformulationsbasedonthestandardLagrange multiplier[2,28].Inaddition,throughitsspecificalgebraic form,itallowsforprovidinganygradientinformationrequiredbytheselectedoptimizationalgorithm.

BasedonEq.(6),thefunctiongisdefinedas

g

=

x

−

Px

˘

=

0

,

(7)

whichrepresentsthemultiscaleproceduretofindapproximateprimaryvariablesatthefinescale.Thestatevectorisnow describedasacombinationofbothsetsofprimaryvariablesatthefineandcoarsescales,i.e.,

x

=



˘

x x



,

(8)

g

(

x

,

θ) =



˘

g g



=

0

.

(9)

Thedefinitionofthestatevector asinEq.(8),asakeyaspect ofthisdevelopment,allows fordescribing thestatenot only atthe fine scale, butalso atthe coarse scale. The simulator responses y obtained from the multiscalemethod are representedas

y

=

h

(

x

,

θ) .

(10)

ThesensitivitymatrixG isthencomputedbyobtainingthederivativeofEq.(10)withrespectto

θ

as

G

=

dh d

θ

=

∂

h

∂

x dx d

θ

+

∂

h

∂

θ

.

(11)

Inorderto findarelationship thatdefinesthederivative ofthestate vector withrespect totheparameters, Eq.(9)is differentiatedwithrespectto

θ

∂

g

∂

x dx d

θ

+

∂

g

∂

θ

=

0

,

(12) sothat dx d

θ

= −



∂

g

∂

x



₋1

∂

g

∂

θ

.

(13)

SubstitutingEq.(13)inEq.(11)gives

G

= −

∂

h

∂

x



∂

g

∂

x



−1

_∂

_g

∂

θ

+

∂

h

∂

θ

.

(14)

FromEqs.(5),(7),(8)and(9),itfollowsthat

∂

g

∂

x

=



RAP 0

−

P I



.

(15) Notethat



∂

g

∂

x



₋1

=



(

RAP

)

−1 0 P

(

RAP

)

−1 I



(16) holds.Thus,Eq.(14)canberestatedas

G

=



∂

h

∂

x

˘

+

∂

h

∂

xP



(

RAP

)

−1

∂

g

˘

∂

θ

+

∂

h

∂

x

∂

g

∂

θ

+

∂

h

∂

θ

.

(17)

ByderivingEq.(5)withrespectto

θ

onefinds

∂

g

˘

∂

θ

=



∂

R

∂

θ

(

AP

)

+

R

∂

A

∂

θ

P

+ (

RA

)

∂

P

∂

θ



˘

x

−

∂

R

∂

θ

q

−

R

∂

q

∂

θ

(18)

andalso,fromEq.(7),

∂

g

∂

θ

= −

∂

P

∂

θ

x

˘

.

(19)

Notethatthe partialderivativesofthematricesA, R andP withrespectto thevectorofparameters

θ

resultin(third order)tensors.TheappropriateinterpretationoftensoroperationscanbefoundinAppendixA.ThesubstitutionofEqs.(18) and(19)inEq.(17)leadstoanexpressionthatwillserveasthebasisforthemultiscalegradientcomputation,i.e.,

G

=



∂

h

∂

x

˘

+

∂

h

∂

xP



(

RAP

)

−1



∂

R

∂

θ

(

AP

)

+

R

∂

A

∂

θ

P

+ (

RA

)

∂

P

∂

θ



˘

x

−

∂

R

∂

θ

q

−

R

∂

q

∂

θ



−

∂

h

∂

x

∂

P

∂

θ

x

˘

+

∂

h

∂

θ

.

(20)

Eq.(20)isquitegeneral,i.e.,itcanbeemployedtocompute gradientsforanyMS(andmulti-level)method,giventhat thepartialderivativesofR andP areprovided.

Forthesakeofsimplicityandtobecoherentwiththenumericalexperimentspresentedbelow,fromnowonthe depen-dencyoftherestrictionoperatorR andtheright-hand-sidevectorq ontheparametersisneglected.Thisisconsistentwith

theMSFV method(where therestrictionoperatorisbasedonafinitevolumeintegrationoperatoratcoarsescale(seee.g. [22]).Afterthesesimplifications,Eq.(20)canberewrittenas

G

=



∂

h

∂

x

+

∂

h

∂

xP



(

RAP

)

−1R



∂

A

∂

θ

P

+

A

∂

P

∂

θ



 x

−

∂

h

∂

x

∂

P

∂

θ

 x

+

∂

h

∂

θ

.

(21)

The orderoftheoperationsinvolvingthe

(RAP)

−1 _{termiscrucial todeterminethe computationalperformance ofthe}

sensitivity matrixcomputation. Thisorderdefines twodifferentalgorithms known intheliterature asthe Direct Method andtheAdjoint Method.The followingsectionsdiscusshowthetwomethods arederivedfortheMS scenariodefinedby Eq.(21),aswellasthe(dis)advantagesofutilizingeachofthem.

3. Computationofgradientinformation:generalizationoftheframework

Accordingto thetypeofstudyone wantstoperform,differentderivative informationhastobeprovided.Forinstance, foroptimizationmethods,Quasi-Newtonmethodsrequirethegradientoftheobjectivefunction(andpossiblyconstraints). Inhistorymatchingalgorithms,aswellasinsensitivityanalysisstudies,thesensitivitymatrixG,definedasthematrixof derivativesofallresponses withrespecttoall parameters,playsan importantrole. InGauss–Newtonmethods,thematrix product GTG isdirectlyused,whileconjugategradientmethodsrequireproductsofG andGT witharbitraryvectors.More detailedinformationonthedifferentoptimizationalgorithmscanfoundin[44].

To preserve the general applicability of our method, the general problem of multiplying G by arbitrary matrices W (of order m

×

NY) and V (of order Nθ

×

n) from the left andthe right, respectively, is considered. When n

=

1, GV is simply theproductofG witha vector,whereas,whenm

=

1,

(WG)

T

₌

_GT_WT _{givestheproduct ofG}T _{withthe(column)} vector WT_{.Those examples show how, by defining algorithms to calculate GV andWG for arbitraryV and W, different} typesofderivativeinformationcanbeaccommodatedinasingleframework.

3.1. Directmethod

ThederivationstartsbyconsideringthecalculationofGV foragivenNθ

×

n matrixV.FromEq.(21),bydefining

Z

=



(

RAP

)

−1R



∂

A

∂

θ

P

+

A

∂

P

∂

θ



 x



V

,

(22)

theproductGV canberewrittenas

GV

=



∂

h

∂

x

+

∂

h

∂

xP



Z

−



∂

h

∂

x

∂

P

∂

θ

 x



V

+

∂

h

∂

θ

V

.

(23)

Here,Z isobtainedasthesolutiontothefollowinglinearsystemofequations

(

RAP

)

Z

=



R



∂

A

∂

θ

P  x

+

A

∂

P

∂

θ

 x



V

.

(24)

ThisisknownastheForwardMethod[42],GradientSimulator [30],orDirectMethod[28].NotethatZ hasdimensionsof

NC

×

n and,therefore, itrequires n linearsystemsto be solved. Hence,the cost ofcomputing GV isproportional to the numberofcolumnsinV.Inparticular,toobtainthefullsensitivitymatrixG withtheDirectMethod,inwhichcaseonehas tosetV equaltotheidentitymatrixoforderNθ,thecostwillbeproportionaltothenumberofparameters.

The algorithm to compute the product of the sensitivitymatrix with a matrix via the Direct Method is depicted in Algorithm 1.

Algorithm 1: RightmultiplyingthesensitivitymatrixbyanarbitrarymatrixviatheDirectMethod.

Input : R,A,P,∂∂Aθ,x˘,∂∂hx˘, ∂h ∂x, ∂h ∂θ,V Output: GV 1 Computeα =∂h ∂x˘+ ∂h ∂xP andβ =∂A ∂θPx˘ 2 foreach j=1,2,...,n do 3 Computeγ =∂P ∂θx˘ V.,j;// Algorithm 3, where m=V.,j 4 Computeδ =RβV.,j+Aγ 5 Solvez= (RAP)−1δ 6 Compute(GV).,j= αz−∂∂xhγ + ∂h ∂θV.,j

InAlgorithm 1,thenotation ‘

.,

j’representsthe j-thcolumnofa matrix.Algorithm 1requiresthecomputation ofthe productof

∂P/∂

θ˘

x byacolumnvector.ThisisdiscussedinthenextsectionwhenAlgorithm 3ispresented.

3.2. Adjointmethod

NextthecalculationofWG foragivenm

×

NY matrixW isconsidered.FromEq.(21),defining

ZT

=

W



∂

h

∂

x

˘

+

∂

h

∂

xP



(

RAP

)

−1

,

(25) oneobtains WG

=

ZTR



∂

A

∂

θ

P

+

A

∂

P

∂

θ



˘

x

−

W

∂

h

∂

x

∂

P

∂

θ

x

˘

+

W

∂

h

∂

θ

,

(26)

whichcanberearrangedas

WG

=



ZTR

∂

A

∂

θ

Px

˘

+



ZTRA

−

W

∂

h

∂

x



∂

P

∂

θ

x

˘



+

W

∂

h

∂

θ

.

(27)

MultiplyingEq.(25)byRAP fromtherightandtransposing,Z isobtainedasthesolutiontothefollowinglinearsystem ofequations

(

RAP

)

TZ

=



∂

h

∂

x

+

∂

h

∂

xP



T WT

.

(28)

Thisisknown intheliterature astheAdjoint (orBackward) Method[32].Notethat nowZ has dimensionsof NC

×

m, henceitrequiresm linearsystemstobesolved.Assuch,thecostofcomputingWG isproportionaltothenumberofrows inW.Inparticular,toobtainthefullsensitivitymatrixG withtheAdjointMethod,inwhichcaseonehastosetW equalto theidentitymatrixoforderNY,thecostwillbeproportionaltothenumberofresponses.

Thealgorithmtocomputetheproductofthesensitivitymatrixwithamatrixfromtheleftviathebackwardmethodis depictedinAlgorithm 2.

Algorithm 2: LeftmultiplyingthesensitivitymatrixbyanarbitrarymatrixviatheAdjointMethod.

Input : R,A,P,∂A ∂θ,x˘,∂∂hx˘, ∂h ∂x, ∂h ∂θ,W Output: WG 1 Computeα =  ∂h ∂x˘+ ∂h ∂xP andβ =R∂∂Aθx˘ 2 foreach i=1,2,...,m do 3 Solvez= (RAP)−TαWT i,. 4 ComputemT₌_zT_RA−W i,.∂∂xh 5 Computeγ =mT∂P ∂θx˘ ;// Algorithm 4 6 Compute(WG)i,.=zTβ + γ +Wi,.∂∂hθ

Algorithm 2 requiresa left multiplication of

∂P/∂

θ˘

x by a row vector mT_{.This will be further discussed in the next} section,whereAlgorithm 4ispresented.

3.3. Remarksabouttheframework

Thispaperpresentsanoveltechniqueforcomputationofthegradientsusingmultiscalemethods.Itentailssome impor-tantfeatureswhicharesummarizedinthissection.

Firstofall,thegeneralformulationofthemethodaccommodatesboththeDirectandAdjointapproaches.Notethatthe previously-developedmultiscalegradientcalculationswereallapplicabletoAdjointformulationonly[37–39].

It is important to also note that the algebraic multiscale-gradient formulation, presented here, does not make any assumptionwithrespecttoneitherthenatureoftheproblemnorthespecificoptimizationparameters.Thisflexible frame-work provides anygradient informationwhich is requiredby the chosen optimization algorithm.Note that the previous methodshavebeendevelopedforspecifictypesofparameters (e.g.,permeabilityin[37]). Assuch,they cannotbereadily used toprovide gradient informationifthe requiredparameters by the chosen optimizationalgorithm arenot those the methodswere developedfor.Forinstance,theworkpresentedin[37]and[38]addressedthecomputationofsensitivities where the parameters were specifically the grid-block permeabilityvalues. Thisspecific parameterleads to the sensitiv-ityofthecoarse-scale quantities(transmissibility) beingrelatedto thefine-scalepermeabilityvia additionallocaladjoint problemsforthebasisfunctions.Assuch,themethodpresentedin[37,38]isnotapplicableif,e.g.,wellcontrolsarebeing usedasoptimizationparameters.Ontheotherhand,theworkpresentedin[39]developsanalgorithmthatallowsthewell controlsbeingused.However,itdoesnotaccountforthesensitivitycalculationsofthecoarse-scalequantitieswithrespect

Fig. 1. IllustrationofMSFVcoarsegridsfora2Ddomain.Givenafine-scalegrid(showninlightsolidblacklines),thecoarsegrid(showninsolidbold black)isimposedasanon-overlappingpartitionofthecomputationaldomain.Thecoarsenodes(vertices)arethenselected(redcells).Connectingcoarse nodesconstructsthedual-coarsegrid(bluecells)wherebasisfunctionsaresolved.(Forinterpretationofthereferencestocolorinthisfigurelegend,the readerisreferredtothewebversionofthisarticle.)

to the fine-scale parameters. Instead, only a coarse-scale adjoint problemwas solved to compute the required gradient information.

Theimplicitdifferentiationstrategyofthepresentmethodmakesitmoreflexible,comparedtothosemethodsdeveloped basedontheLagrangemultiplierformulation(asin[37–39]).NotethattheLagrangemultiplierformulationsaredeveloped basedonanobjectivefunction,whichrequirestobepre-definedforeachspecificapplication[2,28].Instead,thepresented implicitdifferentiationstrategy[42,43]isbuiltonamoregenericsystem,i.e.,thesensitivitymatrix.Moreimportantly,the pre- and post-multiplicationofthis sensitivitymatrix by arbitrarymatricesprovides a flexible wayto conform it tothe specifictypeofgradientinformationneeded(withoutrequiringthepre-computationofthefullsensitivitymatrix).

Finally, the multiscale-gradientformulation is developedthrough an algebraic formulation, which is convenientto be implemented inthenext-generationoptimizationframeworks,andalsobenefitsallmultilevelapproachesthatcan be de-scribedthroughrestrictionandprolongationoperationssuchasmultigrid[45]anddomain-decomposition[46,47]methods. 4. PartialderivativeofMSFVprolongationoperatorwithrespecttotheparameters

Aparticularaspectofthemethodologyisthecomputationofpartialderivativesoftheprolongationoperatorwithrespect to the parameters, i.e.,

∂P/∂

θ

in Eq. (21), and operations with this tensor. This is particularly challenging because the computationofP mightinvolvecomplexoperations,e.g.,thesolutionoflinearsystemsinthecaseofMSFVmethods.

In thiswork,theMSFV methodis considered; extensionstoother MS methodologiescanbe obtainedalongthe same lines.Inthiscase,P istheassemblyofallbasisfunctionsobtainedviathesolutionoflocalflowproblems,i.e.,

−∇ · (λ · ∇

ϕ

)

=

0

,

(29)

where

λ

isthemobilityand

ϕ

isthebasis functionvalue.Thedual-gridsub domainswherethebasisfunctionsare com-putedaredefinedasfollows.Aprimalcoarsegrid(onwhichtheconservativecoarse-scalesystemisconstructed)andadual coarse grid,whichisobtainedby connectingcoarsenodes,aredefinedona givenfine-scale grid.Acoarsenode isafine cellinside(typicallyatthecenter of)eachcoarsecell.Thebasisfunctionsaresolvedlocallyonthesedualcoarsecells.Such overlappingcoarseanddualcoarsegridsarecrucialforconservativesolutionsinMSFV.AnillustrationoftheMSFVgridsis providedinFig. 1.Theprolongationoperatorcanbeexpressedintermsofthebasisfunctionscorrespondingtoeachcoarse cell j

=

1

,

. . . ,

NC as

P

=



ϕ

1

ϕ

2

· · · ϕ

NC−1

ϕ

NC

,

(30)

wherethebasisfunction belongingtocell j,

ϕj

,isatcolumn j,beingavector ofdimension NF [22].Theconstructionof thebasisfunctions

ϕj

isbasedontheFiniteVolume(FV)discretizationofEq.(29)onthedualcoarsegridcells.Inorderto compute thebasisfunctions,reduced-dimensionalproblemsarefirstsolvedattheedgecellstoprovideDirichletboundary conditionsforinternal cells.Forthesakeofclarity, letusassumeaCartesian two-dimensionalsolutiondomain(although extension to3Disstraightforward).Morespecifically,let



j bethesetofalledgesemanatingfromcoarse node j (in 2D, #



j

=

4).Foreachedgee

∈ j

,let

A_eE_,j

ϕ

_eE_,j

−

b_eE_,j

=

0 (31)

bethereduced-dimensiondiscreteformofEq.(29)alongtheedgee,where

ϕ

E

e,jisthebasisfunctionatthatedge,andbeE,j istheRHSwhichresultsfromspecifyingcornervaluesof

ϕ

E

e,j

=

1 atthevertex j and

ϕ

E

e,j

=

0 attheoppositevertexofe. If



j isthesetofallfaceswhichhavecoarsenode j asacommonvertex(#



j

=

4 in2D),foreachface f

∈ 

j,according toEq.(29),onecansolve

AF_f,j

ϕ

F_f,j

−



e∈j

for

ϕ

F_f,j,which is the basis function at theinternal (face)cells. Also, E_eF_,f is a transformation matrix that appropriately assemblesa vectorrepresentedintheedgetopology ofe intothefacetopology of f .Notethat EF

e,f iszerowhen e does notbelongtotheboundaryof f andthatEq.(32)implicitlydefinestheboundaryconditionaszeroforanyboundarycells of f which donotbelongto anyoftheedgesin



j.Finally,

ϕ

j istheresultofassemblingthecontributionofeach

ϕ

Ff,j,

f

∈ 

j,intotheoverallfinemeshtopology:

ϕ

j

=



f∈j

Ef

ϕ

Ff,j

,

(33)

whereEf isa transformationmatrixthatassembles avector representedintheface topologyof f intoasize NF vector followingthefinegridtopology.

Thederivativeof

ϕj

w.r.t.theparametersisobtainedfromEq.(33),i.e.,

∂

ϕ

j

∂

θ

=



f∈j Ef

∂

ϕ

F_f,j

∂

θ

.

(34)

DifferentiatingEq.(32),onecanwrite

∂

AF_f,j

∂

θ

ϕ

F f,j

+

AFf,j

∂

ϕ

F_f,j

∂

θ

−



e∈j E_eF_,f

∂

ϕ

E e,j

∂

θ

=

0

,

(35)

fromwhichitfollowsthat

∂

ϕ

F_f,j

∂

θ

=



AF_f,j

−1

⎛

⎝−

∂

AFf,j

∂

θ

ϕ

F f,j

+



e∈j E_eF_,f

∂

ϕ

E e,j

∂

θ

⎞

⎠ .

(36)

Similarly,fromEq.(31),

∂

ϕ

_eE_,j

∂

θ

= −



A_eE_,j

−1

∂

A E e,j

∂

θ

ϕ

E e,j (37)

holds.CombiningEqs.(36)and(37)into(34)resultsin

∂

ϕ

_j

∂

θ

=



f∈j Ef



AF_f,j

₋1

⎛

⎝−

∂

AFf,j

∂

θ

ϕ

F f,j

−



e∈j E_eF_,f



A_eE_,j

₋1

∂

AE e,j

∂

θ

ϕ

E e,j

⎞

⎠ .

(38)

ThethirdordertensordefiningthederivativeofP w.r.t.theparametersisobtainedbygroupingall

∂

ϕj

/∂

θ

atitsslices,

∂

P

∂

θ

=



∂ϕ1 ∂θ ∂∂ϕθ2

· · ·

∂ϕ_NC−1 ∂θ ∂ϕ_NC ∂θ



NF×NC×Nθ

.

(39)

SeeAppendixAfora discussiononthisnotation. AsdiscussedinAlgorithm 1,theproduct

(∂

P/∂

θ˘

x)m,wherex

˘

∈ R

NC

andm

∈ R

Nθ_{,isrequired.Since}

∂

P

∂

θ

x

˘

=

NC



j=1

∂

ϕ

_j

∂

θ

x

˘

j (40)

wherex

˘

jdenotesthe j-thentryofvector



x,itfollows,fromEq.(38),that



∂

P

∂

θ˘

x



m

=

NC



j=1

˘

xj



f∈j Ef



AF_f,j

₋1

⎛

⎝−



∂

AF_f,j

∂

θ

ϕ

F f,j



m

−



e∈j EF_e,f



A_eE_,j

₋1



∂

AE e,j

∂

θ

ϕ

E e,j



m

⎞

⎠ .

(41)

Algorithm 3depictsthesolutionproceduretocalculate

(∂P/∂

θ˘

x)m basedonEq.(41).

Ina similar manner,Algorithm 2 requirestheproduct mT

(∂P/∂

θ˘

x),where x

˘

∈ R

NC _andm

∈ R

NF_{. FromEqs. (38)and}

(40),oneobtains mT



∂

P

∂

θ

x

˘



=

NC



j=1

˘

xj



f∈j mTEf



AF_f,j

−1

⎛

⎝−

∂

AFf,j

∂

θ

ϕ

F f,j

−



e∈j E_eF_,f



A_eE_,j

−1



∂

A_eE_,j

∂

θ

ϕ

E e,j

⎞

⎠ .

(42)

Algorithm 3: Computationoftheproduct

(∂P



∂

θ˘

x)m,wherex

˘

∈ R

NC _andm

∈ R

Nθ_. Input :x˘,m,ϕE e,j,AeE,j, ∂AE e,j ∂θ ,j=1,...,NC,e∈ j,ϕ F f,j,AFf,j, ∂AF f,j ∂θ ,j=1,...,NC,f∈ j Output:∂∂Pθx˘ m 1 Set:∂P ∂θx˘ m=0

2 foreach j=1,2,...,NC(primalcoarsenodes) do

3 foreach f∈ jdo 4 Computeα = −  ∂AF f,j ∂θ ϕFf,j  m 5 foreach e∈ jdo 6 Computeβ =  ∂AE e,j ∂θ ϕeE,j  m 7 Solvelocaledgesystemγ =AE

e,j −1

β 8 Accumulateresultinα:α = α −EF

e,fγ

9 Solvelocalfacesystemδ =AF f,j −1 α 10 Accumulateresultin∂P ∂θx˘ m :∂P ∂θx˘ m=∂P ∂θx˘ m+ ˘xjEfδ Bydefining

γ

T f,j

=

m T_E f



AF_f,j

₋1 (43) and

δ

T e,f,j

= γ

T_EF e,f



A_eE_,j

₋1

,

(44)

Eq.(42)canberewrittenas

mT



∂

P

∂

θ

x

˘



=

NC



j=1

˘

xj



f∈j

⎛

⎝−γ

T f,j



∂

AF_f,j

∂

θ

ϕ

F f,j



−



e∈j

δ

T e,f,j



∂

AE_e,j

∂

θ

ϕ

E e,j

⎞

⎠ .

(45)

RightmultiplyingEq.(43)byAF_f,jandtransposing,onecanseethat

γ

f,j isobtainedas



AF_f,j

T

γ

f,j

=

ETfm

.

(46)

NotethatET

f isthetransformationmatrixthatassemblesavectorfollowingthefinegridtopologyintoitsproper restric-tiontothefacetopologyof f .Analogously,

δe

,f,jisobtainedbysolving



A_eE_,j

T

δ

e,f,j

=



E_eF_,f

T

γ.

(47)

SimilarasforET_f,

(E

_eF_,f

)

T _{isthetransformationmatrixthatassembles avectorfollowingthefacetopologyof f intoits} restrictiontotheedgetopologyofe.Equations(45),(46)and(47)arethebasisforthealgorithmtocalculatemT

(∂P/∂

θ˘

x), asdepictedinAlgorithm 4.Notethatthefacesystemsaresolvedbeforetheedgeones,incontrasttotheorderoftheMS calculation.

4.1. Prolongationanditsderivativeinthepresenceofwells

Wellsandotherfine-scalesourcetermsareconsideredinthemultiscaleformulationbysupplementarywellbasis func-tions [20,17].WellfunctionsarelocalsolutionstoEq.(29)subjecttoDirichletconditionsof1atthewellcelland0atthe coarsenodes.Hence,consideringwellfunctions,theprolongationoperator(forporousrock)isenrichedandreads

P

=



ϕ

₁

· · · ϕ

_NC

 ψ

1

· · · ψ

NW

,

(48)

where each well function

ψw

adds a column vector to the porous rockprolongation operator. Note that there are NW wellsinthedomain.Forrate-constrainedwells,onehastoadd additionalrowstotheprolongation(andconsequently,the coarse-scale system)inordertosolveforthewellpressures asadditionalunknowns[17].The derivativeof(48)w.r.t.the parametersisgivenby

Algorithm 4: ComputationoftheproductmT

_(∂P



_∂

_θ˘

_x

₎

_,wherex

_˘

_{∈ R}

NC _andm

∈ R

NF_. Input :x˘,m,ϕE e,j,AEe,j, ∂AE e,j ∂θ ,j=1,...,NC,e∈ j,ϕ F f,j,AFf,j, ∂AF f,j ∂θ ,j=1,...,NC,f∈ j Output: mT∂P ∂θx˘ 1 Set:mT∂P ∂θx˘ =0 2 foreach j=1,2,...,NC do 3 foreach f∈ jdo

4 Solvetransposelocalfacesystemγ =AF f,j −T ET fm 5 Computeα= −γT  ∂AF f,i ∂θ ϕFf,i  6 foreach e∈ jdo

7 Solvetransposelocaledgesystemδ=AE e,j −T EF e,f T γ 8 Accumulateresultinα:α= α− δT _∂AE e,j ∂θ ϕeE,j  9 AccumulateresultinmT∂P ∂θ˘x :mT∂P ∂θx˘ =mT∂P ∂θx˘ + ˘xjα

∂

P

∂

θ

=



∂ϕ1 ∂θ

· · ·

∂ϕ_NC ∂θ



∂ψ1 ∂θ

· · ·

∂ψ_NW ∂θ



.

(49)

ComputationofthepartialderivativesofthewellbasisfunctionsfollowssimilarstrategyasdiscussedforEq.(38). 5. ComputationalaspectsoftheMS-gradientmethod

5.1. Partialderivativecomputationandautomaticdifferentiation

Inthecomputationofanalyticalgradients,asignificantpartoftheoverallimplementationeffortisassociatedwiththe computation ofpartial derivativematrices. Thiscomputation stepoftenrequiresaccessto thesourcecodebecausenotall partialderivativesarereadilyavailable,asopposedtothestatevariablederivativesviatheJacobianmatrix(see,e.g.,inputs required by Algorithm 1 andAlgorithm 2). Tocalculatethe partial derivative matrices an AutomaticDifferentiation (AD) approach isused,because ofits flexibility andaccuracy.A review ofdifferentAD approachesis providedby [33].Inour work,anOperatorOverloadingADtechnique(BendtsenandStauning1996)basedonExpressionTemplates[48] isapplied. Thisfacilitatesthecomputationandassemblageofallthepartialderivativematricesbecausetheyareobtainedatthesame timeaswhenthesystemmatrixandthesimulatorresponsearecomputed.

5.2. Computationalefficiency

AnasymptoticanalysisisperformedtoassesstheefficiencyoftheMS-gradientinformation,comparedwiththefine-scale gradientcomputation.Notethatthecomputationalcostoflinearsystemassemblyandmatrix-vectorproductsarenegligible whencomparedtothecostinvolvedinsolvingthelinearsystemofequations.Thecomplexityofsolvingalinearsystemof sizeN isassumedtobe

O(

aNb

₎

_{,wherea andb areconstantsassociatedwiththespecificlinearsolveremployed.}

FollowingAlgorithm 1,n (numberofcolumnsintheV matrix)linearsystemsofsize NC mustbesolvedtofullydefine thegradientinformation,henceits computationalcomplexityreads

O(

n

(

aM SNCbM S

))

.Additionally,Algorithm 3consistsof solving NL

×

n

×

NC linearsystems ofsize NR

=

NF

/

NC toprovidethepartialderivative ofP w.r.t.theparameters.Here,

NR isthemultiscalecoarseningratio,andNL isthenumberoflocalproblemsthatmustbesolvedpercoarsegridvertex(4 in2Dand8in3Dproblems).Thus,thecomputationalcomplexityofAlgorithm 3is

O(

NLnNC

(

aM SNRbM S

))

,whichresultsin acomplexityof

O

_{M S}D

(

n

(

aM SNCbM S

+

NLNC

(

aM SNRbM S

)))

fortheMSDirectMethod.

AsimilaranalysiscanbeperformedforAlgorithm 2andAlgorithm 4,resultinginanestimateofthecomplexity ofthe MSAdjointMethod,i.e.,

O

A_{M S}

(

m

(

aM SNCbM S

+

NLNC

(

aM SNRbM S

)))

,wherem isthenumberofcolumnsofW.

Thecomplexity ofthefine-scalegradient computationfortheDirect andAdjoint Methodsis givenby

O(

n

(

aF SNCbF S

))

and

O(

m

(

aF SNCbF S

))

,respectively.Therefore,thecostratiobetweentheMSandfine-scalecomputationalefficiencycanbe definedas

O

M S

O

F S

=

aM S aF S



NCbM S NFbF S

+

NL NFbM S NFbF S NC(1−bM S)



,

(50)

whichholdsforboththeDirectandtheAdjointMethod.

For the sake of simplicity, it is assumedthat the solver employed to the MS system is equally efficient to the one employed to the fine-scale system, i.e., aM S

=

aF S and bM S

=

bF S

=

b. As such, the expression defining the cost ratio simplifiesto

Fig. 2. CostratiobetweenMSandfine-scalegradientcomputationmethods,asafunctionofthemultiscalecoarseningratio,NR,fora2D(blue)and3D

(red)domainwithNF=107fine-scalecells.(Forinterpretationofthereferencestocolorinthisfigurelegend,thereaderisreferredtothewebversionof

thisarticle.)

O

M S

O

F S

=

1 Nb_R

+

NL Nb_R−1 Nb_F−1

.

(51)

For a fixed numberoffine grid cells NF, it ismainly the coarseningratio NR that determines the complexity of the MS-gradientalgorithm.Toillustratethispoint,adomainwithNF

=

107fine-scalegridcellsisconsidered.TheMS-gradient speed-up,

O

M S

/

O

F S,fordifferentcoarseningratios NR forboth2Dand3DcasesispresentedinFig. 2,whereb

=

1

.

3. 6. Numericalexperiments

Inthissection,performanceoftheMS-gradientmethodisstudiedforsingle-phaseincompressibleflowinheterogeneous porous media. The following numerical experiments are presented to first validate andthen assessthe accuracy of the gradientinformationcomputedbythemethod.Forthispurpose,amisfitobjectivefunctionwithnoregularizationterm

O

(

θ) =

1 2

(

h

(

x

,

θ) −

dobs

)

T_C−1 D

(

h

(

x

,

θ) −

dobs

) ,

(52) withagradient

∇

θO

=

GTC−_D1

(

h

(

x

,

θ) −

dobs

) ,

(53)

isconsidered[28].Inallexperiments,thefittingparametersarecell-centeredpermeabilities.Theobservedquantity,dobs,is thefinescalepressureatthelocationof(non-flowing)observationwells,therefore

∂

h

∂

x

=

I

,

(54)

and

∂

h

∂

x

˘

=

0

.

(55)

Thetensor

∂A/∂

θ

representspartialderivativesofthesystem(transmissibility)matrixwithrespecttopermeability.Note thattheresultsareexpressedintermsofanon-dimensionalpressure,i.e.,

pD

=

p

−

pprod pinj

−

pprod

,

(56)

where pinj andpprodaretheinjectionandproductionpressures,respectively.Inalltheexperiments,pinj

=

1

.

0 and pprod

=

0

.

0,thegrid-block dimensionsare

x

=

y

=

z

=

1 m andthefluid viscosityis1

.

0

×

10−3 _{Pas.Inaddition,inall the}

followingtestcases,wellbasisfunctionsareincluded.

6.1. Validationexperiments

TheMS-gradientmethodisvalidatedagainstthenumericaldifferentiationmethodwithahigher-order,two-sidedTaylor approximation

∇

θOi

=

1 2

δθ

i

O

θ

i

, . . . , θ

i−1

, θ

i

+ δθ

i

, θ

i+i

, . . . , θ

Nθ



−

O

θ

i

, . . . , θ

i−1

, θ

i

− δθ

i

, θ

i+i

, . . . , θ

Nθ



,

(57)

Fig. 3. Fineand coarsegridsandwellssetupfor1Dnumericalexperiments.Thesolidthinlinesrepresentthefinegrid-blocks.Thebold dashedlines representtheprimal-coarsegridblocks.Vertexcellsidentified withredcircles.Thecrossedcirclerepresentstheinjectionwell,thedottedcirclethe productionwell,andthesolidcircletheobservationwell.(Forinterpretationofthereferencestocolorinthisfigurelegend,thereaderisreferredtothe webversionofthisarticle.)

Fig. 4. (a) Fine, primal and dual coarse grids for 2D validation test cases and (b) reference permeability field.

where

δ

isamultiplicativeparameterperturbation.Therelativeerrorcanbedefinedas

ε

=

∇

θOF D

− ∇

θOAN



2

∇

θOAN



2

,

(58)

where

∇

_θOF D isobtainedby performingthe appropriatenumberofmultiscale reservoirsimulations requiredtoevaluate Eq.(57)and

∇

θOANisobtainedbyeitheremployingtheDirectortheAdjointMethodtoevaluateEq.(53).Notethat

m

=

C−_D1

(

h

(

x

,

θ) −

dobs

) ,

(59)

where m is an auxiliary vector, sothe gradientof O canbe written as

∇

_θO

= (

mT_G)T _{andAlgorithm 1, withW}

₌

_mT_, calculates

∇

_θO with a cost proportional to one extraMS simulation,making the Adjoint Method a much moreefficient waytocalculatethegradientthan theDirect Method.Moreover, forallcases,unlessstatedotherwise,forsimplicity,itis assumedCD

=

I.

Inordertovalidatetheproposedderivativecalculationmethods,aswellastheirimplementation,thelineardecreaseof theerror

ε

bydecreasingtheperturbationvalue

δ

(seechapter8of[29])from10−1_to10−4 _{(therangewithinwhichonly}

discretizationerrorsareobserved)isinvestigated.

Theinvestigationiscarriedoutinthreeexamplesofincreasingcomplexity.Thefirstcaseisaone-dimensional, homoge-neousmediumwith45gridblocks.Aprimalcoarsegridofjust3gridblocksisemployed(coarseningratioof15).Injection andproductionwellsarelocatedat,respectively,cells1and45.Oneobservationwellislocatedatthecenter ofthedomain, i.e.atgridblock23.Thepressuremeasuredintheobservationwellistakenfromareferencecasewitharandomlysampled permeabilityfield.Thenon-dimensionalinjectionandproductionpressuresareoneandzero,respectively.Fig. 3illustrates thesetupforthisexperiment.

Intheothertwoexperiments,theaccuracyofthemethodisassessedfor2Dtestcases,one homogeneousandanother oneheterogeneous.Inbothcasesthefinegridsizeis21

×

21,whilethecoarsegridsizeis3

×

3 (coarseningratioof7

×

7). Theprimal- anddual-coarsegridsareillustratedinFig. 4(a).

Thepermeabilityfield isextractedfrom1000geological realizations,asshowninFig. 4(b),andservesto computethe observedpressure.Fortheheterogeneouscase,anothergeologicalrealizationischosenfromtheensemble.Theensembleis generatedviathedecompositionofareferencepermeability“image”usingPrincipalComponentAnalysisparameterization. Fig. 5illustrates4differentpermeabilityrealizationsfromtheensemble.See[49]formoredetails.

Aquarterfive-spot wellconfigurationis considered,withtwoobservationwellscloseto theoperatingwells.The well positionsandoperatingpressuresaredescribedinTable 1.

Fig. 6showsthattheMS-gradientandfine-scalemethodshavethesameorderofaccuracywithrespecttothe perturba-tion

δ

,forallcases.Theexpectedbehavior oflinearlydecreasingerrorvaluesastheperturbationsizedecreasesisobserved inallexperiments.OnecanalsonotethattheDirect andAdjointMethods areequallyaccurate,i.e.,theybothprovidethe analyticalgradientatthesameaccuracy.Also,itisimportanttonoticethatthelocalizationassumptionsinvolvedin2Dare consistentlyrepresentedby theanalyticalmethods.Theresultofthethirdexperiment(Fig. 6(c))indicatesthecorrectness ofthemethodwhenappliedtocomputegradientsforheterogeneousmedia.

Fig. 5. Four different permeability realizations from the ensemble of 1000 members used in the 2D numerical experiments.

Table 1

Wellconfigurationforthehomogeneous,two-dimensionalcase. Well Fine scale position (I, J) Well type

INJE (1, 1) Injection

PROD (21, 21) Production

OBSWELL1 (3, 3) Observation OBSWELL2 (19, 19) Observation

Fig. 6. Validation ofMS gradient computation method via comparisonwith numerical differentiation.(a) One-dimensional, homogeneous, (b) two-dimensionalhomogeneous and(c)two-dimensionalheterogeneoussteady-stateflowtestcases.

6.2. Gradientaccuracy

InordertoassessthequalityoftheMSgradient,theanglebetweenfine-scaleandMSnormalizedgradients,i.e.,

α

=

cos−1



∇

T θO

ˆ

F S

∇

θO

ˆ

M S

,

(60) ismeasured.Here,

∇

θO

ˆ

F S

=

∇

θ OF S

∇

θOF S



2 (61) and

∇

θO

ˆ

M S

=

∇

θ OM S

∇

θOM S



2

.

(62)

Also,

∇

_θOF S and

∇

θOM S denotethe fine-scale andMS analytical gradients, respectively. As a minimumrequirement, acceptableMSgradientsareobtainedif

α

ismuchsmallerthan90◦[50].

Firstly, forthe 1D test case, both methods result in

α

=

0◦,indicating that finescale and MS gradients are perfectly alignedinthiscase.Thisisduetothefactthat,in1D,noapproximations(duetolocalization)aremadeintheMSsolution, andthus,intheMSgradientcomputation.

Forthe2D homogeneous case,thefinescaleandMS gradientsresultin

α

=

10

.

97◦,usingthesamesetupdepictedin Fig. 4(a).Althoughthegradientsarepracticallypointinginthesamedirection,adeviationbetweenthetwoisobserved.To better investigatethisdifference,Fig. 7(a)andFig. 7(b)presentthefine-scaleandMS pressuresolutions,respectively,and Fig. 7(c)showsthedifferencebetweenthesetwopressuresolutions.Fig. 7(d)andFig. 7(e)presenttheabsolute,normalized

Fig. 7. Finescale(p)pressuresolution(a)andMSFV(p)pressuresolution(b).Differencebetweenfine-scaleandMSFVpressuresolutions(c).Absolute, normalizedfinescalegradient(d)andMSFV(e).Differencebetweenfine-scaleandMSabsolute,normalizedgradients(f).Inallfiguresthedual-coarsegrid isshown.Alsoshown,intheredboxes,arethecoarsenodes.(Forinterpretationofthereferencestocolorinthisfigurelegend,thereaderisreferredto thewebversionofthisarticle.)

Table 2

Wellconfigurationforthehomogeneous,two-dimensionalcase. Well Fine scale position (I, J) Well type

INJE (11, 11) Injection PROD1 (1, 1) Production PROD2 (21, 1) Production PROD3 (1, 21) Production PROD4 (21, 21) Production OBSWELL1 (3, 3) Observation OBSWELL2 (19, 3) Observation OBSWELL3 (3, 19) Observation OBSWELL4 (19, 19) Observation

gradientoftheOFcomputedviafine-scaleandMS gradientmethods,respectively.Fig. 7(f)showsthedifference between theabsolute,normalizedgradientcomputedbythetwomethods.

The differencebetweenthe MSandfine-scale methodsis dueto thelocalizationassumption inthecalculation ofthe basisfunctions[14].Notethat,thankstothewellfunctions,themultiscalesolutionsareaccuratelycapturingthewells.

6.3. Effectofheterogeneitydistributionandcoarseningratio

Inthe firstcase, OFgradients arecomputedforthe wholeensemble ofheterogeneouspermeabilityfields(see Fig. 4). Moreover,thesamerealizationdepictedinFig. 4(b)isconsideredasthereferencefromwhichtheobservedpressuresare computed.

Aninvertedfive-spotwell patternisemployed,whiled fourobservationwellsareplacedclosetotheproductionwells. ThewellconfigurationisdepictedinTable 2.

Twodifferentcoarseningratiosareapplied,oneresultinginacoarsegridof3

×

3 (asillustratedinFig. 4(a))andanother oneresultinginacoarsegridof7

×

7.TheanglesbetweenthefinescalegradientandtheMSgradientforeachrealization arecomputed,usingEq.(60).AhistogramillustratingtheangledistributionforeachcoarseningratioispresentedinFig. 8. Notethat,forthiscase,thequalityofthegradientsissignificantlyimprovedoncethecoarsegridsizeof3

×

3 isincreased to7

×

7.Importanttonote isthat,asshowninFig. 9,theMSgradientsareleastaccuratewhenthenormofthegradient vectorissmall.Therefore,theyarenotexpectedtohaveamajoreffectontheoptimizationprocedure.

Forthe3

×

3 coarsegridcase,9.71%oftheanglesaregreaterthan90◦,indicatingthatsome MSgradientspointinthe oppositedirectionofthedecreasingOFdirection.Ontheotherhand,ifa7

×

7 coarsegridisusednoangleisgreaterthan 90◦.

Fig. 8. Histogramsofanglesbetweenfinescaleandmultiscalegradientsfor3×3 (a)and7×7 (b)coarsegrids.Thefine-scalecomputationaldomain contains21×21 gridcells.Theverticalreddashedlineindicatesthe90◦limit.(Forinterpretationofthereferencestocolorinthisfigurelegend,the readerisreferredtothewebversionofthisarticle.)

Fig. 9. Cross-plotsoffinescalegradientnormvs.αforthe1000realizationsintheheterogeneousensemblefor(a)3×3 and(b) 7×7 coarsegrid.Note thatMSgradientsareaccurateforthecaseswithlargevaluesofgradientnorms.Theirrelativelyinaccurateestimates(largeα)happenmostlywhenthe normsofgradientsaresmall.Theverticalreddashedlineindicatesthe90◦limit.(Forinterpretationofthereferencestocolorinthisfigurelegend,the readerisreferredtothewebversionofthisarticle.)

Fig. 10. Permeabilitydistributionoffourdifferentrealizationstakenfromthesetsof20geostatisticallyequiprobablepermeabilityfieldswith0◦(a),15◦ (b),and45◦(c)correlationangles.Also,apatchyfield(d)withasmallcorrelationlengthisconsidered.

Inordertofurtherexplorethispoint,foursetsof20equiprobable realizationsoflog-normallydistributedpermeability fields with a spherical variogram anddimensionless correlation lengthsof

1

=

0

.

5 and

2

=

0

.

02 are generated using

sequentialGaussiansimulations[51].Foreachset,thevarianceandthemeanofln

(

k

)

are2.0and3.0,respectively,wherek

isthegridblockpermeability.AsdepictedinFig. 10,fortherealizationswithalongcorrelationlength,theanglesbetween the permeability layers and thehorizontal axis are 0◦, 15◦, and 45◦.A patchy (small correlation length) patternis also considered(Fig. 10(d)).Comparedwiththepreviousset,thepermeabilitycontrastismuchhigherinthiscase.

Thefine-scaleandcoarsegridscontain100

×

100 and20

×

20 cells,respectively.Thewellconfigurationutilizedinthis numericalexperimentisdepictedinTable 3.

From Fig. 11,one can observe that all caseswith

α

>

50o _{are associated withsmallgradient norms andthat, forall}

geological sets,theMS gradientprovidesgradientdirectionsthatareveryaccurate,comparedwiththefinescalesolution, whenthegradientnormsarelarge.

Notealsothat, inthecaseofasmallgradientnorm,theOFhasaweakdependency ontheparametersinthevicinity ofthepoint wherethegradientisbeingcalculated.As such,theoverallperformance ofthe optimizationalgorithmisnot