Feature-based estimation for applications in geosciences

(1)

(2)

(3)

for appli ations in geos ien es

Proefs hrift

ter verkrijging van de graad van do tor

aan de Te hnis he Universiteit Delft,

op gezag van de Re tor Magni us prof. ir. K.C.A.M. Luyben,

voorzitter van het College voor Promoties,

in het openbaar te verdedigen

op maandag 10 de ember 2012 om 10.00 uur

door

Wiktoria AWNICZAK

Master of S ien e in Applied Mathemati s,

Delft University of Te hnology, The Netherlands

(4)

Prof.dr.ir. A.W.Heemink

Prof.dr.ir. J.D. Jansen

SamenstellingPromotie ommissie:

Re tor Magni us voorzitter

Prof.dr.ir. A.W.Heemink Te hnis he Universiteit Delft, promotor

Prof.dr.ir. J.D. Jansen Te hnis he Universiteit Delft, promotor

Prof. dr. W. R.Rossen Te hnis he Universiteit Delft

Prof.dr.ir. J.Biemond Te hnis he Universiteit Delft

Prof. dr. ir. F.C. van Geer Utre ht University

Dr. P. J. van den Hoek Shell InternationalE&P

Dr. R.G.Hanea Te hnis he Universiteit Delft

Copyright

2012 by Wiktoria awni zak

reprodu ed or utilized in any form or by any means, ele troni or me hani al, in luding

photo opying, re ording, or by any information storage and retrieval system, without

written permissionfromthe author.

Cover design: Proefs hriftmaken.nl ||UitgeverijBOXPress

Printed & Lay Out by: Proefs hriftmaken.nl ||Uitgeverij BOXPress

(5)

(6)

(7)

1 Introdu tion 3

2 Numeri al modeling of reservoirs 9

2.1 Modelingof reservoirs . . . 9

2.2 simsim . . . 12

2.3 State-spa e representation of the model . . . 16

3 Data assimilation and parameter estimation 17 3.1 Problem setup . . . 17

3.2 The lassi al Kalmanlter . . . 18

3.3 The Ensemble KalmanFilter - EnKF . . . 21

3.3.1 The Ensemble Square-RootFilter -EnSRF . . . 22

3.4 Implementation issues of the Ensemble Kalman Filter . . . 24

3.5 Parameter estimation problem . . . 28

4 Ensemble Multis ale Filter - EnMSF 31 4.1 Multis ale ensemble ltering formeasurementupdate . . . 31

4.1.1 Introdu tion . . . 31

4.1.2 EnKF and EnMSF - theoreti al ba kground . . . 32

4.1.3 Appli ation . . . 40

4.1.4 Con lusions . . . 43

4.2 Towards the use of the ensemble multis alelter for history mat hing . . . 44

4.2.2 Pixel numbering s heme . . . 44

4.2.3 Interpolationproblem . . . 47

4.2.4 Coarsening parameters in the EnMSF. . . 47

4.2.5 The designof the tree and itsimpa t . . . 49

4.2.6 Covarian ematrix approximation . . . 53

4.2.7 History mat hing usingEnMSF . . . 54

4.3 Con lusions . . . 58

5 Feature-based methods 61 5.1 Grid distortionfor a 2D reservoir model . . . 61

5.1.2 Basi notions . . . 63

5.1.3 Field Alignment . . . 65

(8)

5.1.6 Twin Experiment . . . 76

5.1.7 Dis ussion and Con lusions . . . 91

5.2 Grid distortionfor a3D groundwater ow model. . . 91

5.2.1 The 3D grid distortionmethod . . . 91

5.2.2 3D groundwater owmodel . . . 92

5.2.3 Experimentsetup . . . 93

5.2.4 Dis ussion and results . . . 98

5.3 Grid distortionfor a3D reservoir model . . . 101

5.3.1 Dis ussion and on lusions . . . 108

6 Con lusions and re ommendations 109

Summary 119

Samenvatting 121

(9)

Thisresear hwas arriedoutwithinthe ontextoftheISAPPKnowledgeCentre. ISAPP

(Integrated SystemsApproa htoPetroleum Produ tion) isajointproje t ofthe

Nether-lands Organizationfor Applied S ienti Resear h TNO, ShellGlobal Solutions

Interna-tional,and Delft University of Te hnology.

Bana h[41℄isdeadbutthereareothergreatmathemati iansImetonmywaythrough

s hools. I would liketo thankthem fortea hing me rigorousthinkingand I hope thatby

writingthisthesis Iproved thattheir timewasnot ompletelywasted. Espe ially,I want

toname: Ali ja Gande ka, JanSzajkowski, Dorota Krassowska, Mie zyslaw Trad, Anna

Kubi ka and Jolanta Misiewi z. All of them really s ared me at the time but now I am

gratefulfor the solid mathbasis I obtained.

I onsider myself lu ky to have four (4!) supervisors during the onstru tion of this

do toratethesis. I would like tothank Dr. RemusHanea for giving me the opportunity

to start this PhD program, and his friendliness; Prof. Arnold Heemink for good sense

of humor, great ideas and ontagious pea efulness; Prof. Jan Dirk Jansen for being fair,

supportiveand onsistent; andProf. Dennis M Laughlinforallthe timehe took towork

withme, I greatly appre iatethe trips toBoston for multiplereasons. I amvery grateful

to Prof. Heemink for translating the propositions into Dut h and for all the help and

understanding duringthe last part of gettingthis degree.

The toughest job was to make all the administration, papers, omputers work. It

would have never been done if it was not for: Dorothee Engering, Evelyn Sharabi, Carl

S hneider, JimLong and Anke Dahlmann. Thank you forbeing sokind.

Forme the se rets ofgoodtea hing [44℄ are stilltobe dis overed but I amglad I had

good ourse oordinators: Roelof Koekoek, Harry Kneppers and Henri Corstens.

I realize that there are many great books to learn from but I think there are even

more great people to learn from. The literature meetings organized by Martin Verlaan

and dis ussions within our littlegroup are an irrepla eablesour e of knowledge. Thank

you all!

Always the o e mates make aday atthe university fun, atleast mine did. All sorts

ofproblems ouldbesolvedor reated withSvetlana,MariyaandCristi;Igreatlyenjoyed

my timewith Elena, Luis, Gift and Raf. Thanks guys! I alsowanttothank Angéline for

avery pleasant ooperation,and olleagues and friends at DIAM.

I have made great friends during my studies that substituted for my family, thank

you: mylovelyroomieswho havebeen makingevery-daylifeinteresting andentertaining,

Karolinafor allthe onversations,trips, her greatsense ofhumor andemergen y

transla-tion of the summary, and Angela for being full of positive inspiring energy; my Mexi an

family who always makes me feel wel ome: Sandra, Camila y OswalDOS; my Ameri an

(10)

(with Olek and Zosia); the always-there, like brother to me: Mar in (with Monika and

Natalia); the hilled and warm Paisas: Alejandra, Manuela y Sebastian. At the same

time, away from Poland, I ould also keep my dear quarter-of-a- entury friends. I want

to thank: Gosia (with Tomek and Klaudia) for her onden e and view on life, and Ola

(with Daniel and u ja) for always being in a good mood. And Adrian for being so

fun and easygoing. An additionalthank you goesto: Eka, Israel, Christiaan, Manuelito,

Maryam, Diego and Nora, Gosia, Eduardo, Isabel, Tatiana and Dirk, Ivan, Asia, An a,

Joaquin and Marisa, Negar.

I would alsoliketothank Ri hardand the small Commer ialteam, forthe great new

beginnings.

Most of all I would like to thank my family inPoland. My grandparents who always

miss me a lot;my parentsthat are alwaysworried but keep it all ool; and my aunt who

likes alling me.

Mar o&F., I am glad we met just at the right time. Thank you for the good times

together, youare irrepla eable.

Thank you. Dziekuje. Gra ias. Bedankt.

(11)

Introdu tion

The energy requirements around the world are expe ted to grow. Even with blooming

and expanding green energy R&D, the developing so iety will still need petroleum for

another tens of years. A on ern might be that the oil and gas reserves that are easily

a essible have, to a large extent, been produ ed already. The sour es that are left are

eitheringeologi ally ompli atedareasor ontain heavy oilthatisdi ulttoextra tand

pro ess. Therefore, new te hnologies and any possible improvements to all stages of the

oilre overy an ontribute tothe world energy supply.

Mathemati al modeling of reservoirs

A omputer model an help fore ast some better known phenomena: in meteorology

it improves weather predi tion, in y lone tra king it an save lives when people are

eva uated on time from the endangered areas, in reservoir engineering it helps optimize

e onomi s. Nevertheless, omputersimulations anneverpre iselypredi tthepro essesin

thenature whether itisreservoir engineering,weatherpredi tion, spe ies spread,market

u tuations;therealworldisfartoo omplex. Therefore,dataassimilationmethodsneed

tobe employed.

Reservoir modeling plays an important role in enhan ed oil re overy planning. It

is a basis for eld development and involves a signi ant number of resear hers around

the world in the industry as well as in the a ademia. The basi physi s of uid moving

through a porous medium has been known sin e Henry Dar y in 1856 published his

resear h. Modern omputers allow to e iently simulate the behavior of the uids in

dierent reservoir onditions and with several hemi al omponentspresent.

Even though the models are not perfe t due to simpli ationsand approximations,a

bigger issue in reservoir engineering is reservoir ondition spe i ation. Unlike in some

otherdis iplines,likemeteorologyforinstan e, itisnotpossibletoobservethe domainof

interest that is deep underground or, even more hallenging, under the sea bottom. For

that reason, the ru ial reservoir parameters like permeability or porosity are unknown.

Therefore, initially the reservoir stru tural information omes from the geologists and

geophysi ists.

Thewaythesedimentsformandwhatpro esses reatethelands apethroughthe ages

determine whi h stru tures are to beexpe ted inthe underground formations. Basedon

(12)

a ording totheir respe tive experien e and interpretations. These property models are

sometimesverydetailedandneedtobeups aledtomakethemsuitableforowsimulation.

Observations of the behavior of the reservoirs

Whenanoileldisdis overedthere anbeexplorationdata olle ted. Themost ommon

larges aleobservationsinthis asewouldbeseismi observations. Togatherseismi data,

a ontrolledexplosionisoftenperformedand thesignalofthe ree tedwavesisre orded.

The seismi waves ree tfrom subsurfa e stru tures in adierent way depending onthe

stru ture's density, and that reatesa seismi image. This image an be interpreted and

inverted intosubsurfa e stru tures.

However, today there are not many new elds being dis overed. The industry highly

relies on oilelds that are at various stages of development. Then, one wants to olle t

dataforprodu tionimprovementbutatthelowestpossible ost, andmostdatatypesare

expensive to obtain.

Figure 1.1: Simulated time-lapse gravity

[mi ro Gal℄measurements ofa water-driven gas

reservoirprodu tionpro essat(fromthetop)

4

,

7

and

10

years. Courtesyof M. Glegola. Additionally to seismi data, there are

other large s ale observation te hniques

available,for example,gravity(Figure1.1)

or magneti observations. While the wells

are drilled inthe eld, the ro k that is

be-ing removed an be a sour e of valuable

information about the underground

stru -ture. It is a type of hard data (data that

does not need to be inverted) that gives

dire t values of permeability and porosity

throughlabtestsonaverysmalls ale.

Ad-ditionally, there an be wireline log data

olle ted from instruments lowered down

an open borehole. These an provide very

a urate measurements in a small vi inity

of the well.

There are dierent types of wells, the

more ompli atedthe moreexpensive they

get, rea hing up to tens of millions of

dol-lars per single well. The simplest wells

are the verti al monitoring wellsthat only

gather data and do not intera t with the

surroundings ( ommon in hydrogeology). There are inje tion/produ tion verti al wells

that an also olle tpressureand/oruidowdata. Themostte hni allyadvan edwells

aretheso- alledsmarthorizontalwellsthat anrea hareasthataredi ulttoa essand

an be kilometersaway. Thesewells an sometimes be opened/ losed toowin separate

(13)

In the ase that a model is well alibrated and its parameters are known, it would have

good predi tive apa ities, and it would be good enough to use the model in de ision

making on erning eld development. When the data are ri h and pre ise, they might

be su ient for su essful operation of the eld without using a model. Unfortunately,

in most ases neither the model nor the data provide good-enough information due to

simpli ationsand un ertainties, and the ombinationof the two an help to redu e the

s ar ity of information.

Data assimilation is a on ept that ombines the model with observed data in order

to minimize un ertainties and improve predi tions. The model as a theoreti al

repre-sentation of a natural pro ess is not perfe t and its un ertain parameters an be better

estimated given a quired data. An improved model is expe ted tohave better predi tive

apabilities, and in reservoir engineering an help planning the eld development and

designing produ tionstrategies.

Figure 1.2: Map of Muskauer Park at the

German-Polishborder. Ameandering river uts throughthepark.

The park's south-east border ar is the biggest terminal

morainein the world. Formally, data assimilation

is a mathemati al optimization

problem. Typi ally, one is

inter-estedtominimizeamismat h

be-tween the data observed and the

data modeled,and this mismat h

is a basis for reating an

obje -tive(or ost)fun tion. The

obje -tivefun tionis ahighlynonlinear

fun tion of its variables in

large-s ale appli ations. Then, the

problemsaredi ulttosolveand,

therefore, a variety of data

as-similationte hniques exists, none

of whi h is universally

appli a-ble. Additionally, the problem

is non-unique, i.e. there exist

its many solutions that minimize

the given obje tive. In other

words,real-lifeoptimization

prob-lems are severely ill-posed.

A variety of methodsfor data

assimilation has been developed

for the oil industry, [62℄. Espe ially, within the last years emphasis has been put to

produ e geologi ally onsistent permeability and porosity elds. Geologi al onsisten y

requires the hara teristi s of the geometri stru ture of the parameters to be preserved

duringthe dataassimilationpro ess, [10℄,[11℄. These hara teristi sarepropertiesof the

eldsthatgobeyondthemeanandthe ovarian e. Naturalformationsshowawiderange

of patterns like river meanders, salt domes, layer- ake formations, faults, wind or water

eddies, or moraines (Figure 1.2). Dealing with patterns and handling images is learly

(14)

The image pro essing dis iplineoers a variety of methodsfor dealingwith images,

pat-terns,featuresandotherstru tures. Ourgoalisto ombinedataassimilationwithpattern

mat hingideas. Thetwoareasarenotfarapart,andamoredetailedinsightrevealsmany

similarities. Imagepro essingproblemsdealwithphenomenawithoutanunderlying

phys-i al model, like a sequen e of movie snapshots. The onse utive states (snapshots) are

known andthe goalistodeneanautomati warpingpro essfromoneimagetoanother.

The warping should be done in a natural way, i.e. in a way a human observer would do

it, whi h rea hes as far as ognitive s ien e. An example ould be fa ere ognition when

seeing one side of somebody's fa e would make us re ognize a person in real life. Our

brains interpolate that partial informationin aninstant but this every-day experien eis

avery omplextaskfor omputerfa ere ognitionsoftware. Patternre ognitionproblems

in lude also ommon nger-print mat hing, as present nowadays at airports, in laptops,

or a poli e database. An e ologi ally oriented appli ation is, for example, salamander

mat hing,[95℄, whereadigitalpi tureofthe animalneedstobe omparedtoa olle tion

of previously gathered images.

Appli ations like y lonepredi tion,[2℄, movingre fronts, [54℄,[6℄, pre ipitationand

thunderstorms,[42℄, [100℄,[26℄,orepidemi spreading, [53℄,deal withstru tured,

feature-driven elds. These approa hes tou h the area of image pro essingbut atthe same time

ontain dynami models of the underlying phenomena. There, data assimilationis

per-formed with a large amount of observations, namely, the parameter eld of interest is

observed but it is dierent from the model output. The two images present the same

feature that might be displa ed or aligned dierently. The di ulty lies in making

on-sistent orre tions to the model predi tions that have to be approa hed globally where

underlyingfeatures are taken intoa ount.

In this thesis we will adopt ideas from the image pro essing area within data

assim-ilation methods for reservoir engineering appli ations. A typi al problem in reservoir

engineering is that the observations are usually very s ar e and nonlinearly related to

the variables that need to be estimated. We oftendeal with ill-posedness due to a large

amount of unknown variables ompared to the number of observations. The spatial

un- ertainty wouldusuallyberepresented withasmallensemble ofnodes,mu hsmallerthan

the number of variables. This representation is too simple to pi ture the omplex

spa-tial un ertainty. This way spurious orrelations arise between physi ally or ideologi ally

distant states.

We will see that in reservoir engineering the unknown variables to be estimated are

ommonlypermeability eldsthat are nothingelse but images. Often, theseimages have

a predened stru ture that disappears during the data assimilation and, therefore, the

whole pro ess loses the geologi alba kground and realism.

Obje tive of the thesis

Wepropose tworesear h dire tions to resolvethe aforementionedproblems:

•

diminishingthe spurious orrelations between thesparse dataand distantlylo ated variables through ups aled treatment of the ovarian es,

(15)

•

taking into a ount features in the domain and redu ing the state size through an ee tivereparametrization.

First, a type of ensemble Kalman lter is applied, namely, an ensemble multis ale

lter. In reservoir engineering,ensemble sequentiallteringstartedwith [59℄ in2002 and

sin e then various modi ations were implemented, [20℄. The ensemble multis ale lter

has originated from multis aletrees in image pro essing ([22℄, [92℄, pyramidal mat hing

[2℄), and it was rst implemented for an in ompressibleow model, [99℄. It ups ales the

ovarian e matrix on a tree stru ture, extra ting stronger dependen ies and introdu ing

lo alizedimprovements. Wedevelop itfor the estimation of thepermeabilityinreservoir

models and test it thoroughlyfor several ases in reservoir engineering.

Se ond, a feature modeling te hnique alled grid distortion is proposed. In reservoir

engineering a number of feature-based methods have been implemented already: level

sets ([58℄, [14℄), Karhunen-Loève expansion ([72℄), dis rete osine transformation ([34℄,

[35℄, [36℄), hannel parametrization ([97℄, [84℄), training-image based sampling ([49℄ and

referen es therein) and elasti gridding ([76℄, [75℄). The grid distortionmethod is based

onsmooth gridgenerationte hniques, [80℄,but it anequivalentlybeviewedasanimage

warping te hnique, [67℄. It provides smooth distortions of images adjusting the features

inthedomain,and anbeexpressed usingveryfewparameters,whi hmakestheproblem

better-posed. We pair grid distortionwith the ensemble Kalman lter algorithmas well

aswith deterministi sear h methods.

Overview of the thesis

Before presenting the resear h results, an introdu tion to the modeling of reservoirs is

provided. InChapter2wedes ribethebasi reservoirequations,showasmallexampleof

aforwardrunofanin-housereservoirsimulator,andoutlinea on eptofaforwardmodel.

Data assimilationwith (ensemble) Kalmanlters and optimization methods used in the

thesis omprise Chapter 3. There, we derive basi (ensemble) Kalman lter equations,

point out some improvements, and present other parameter estimation methods. The

ensemblemultis alelterand itsappli ationare presented inChapter 4,rstforasimple

ase and later as a full data assimilation s heme in the ontext of reservoir engineering

appli ation. The ontent of this hapter has been published in [46℄ and [45℄. F

eature-basedmethodsaredis ussed inChapter 5;here,thefo usismainlyonthegrid distortion

methoddeveloped inthisthesis. Theleadingequationsforgriddistortionarederived and

multipleimplementationsin2Dand 3D ases areshown. Thereisanintentiontopublish

(16)

(17)

Numeri al modeling of reservoirs

Inthis hapter,abrief des riptionof reservoir modelingispresented, followed by

spe i- ationsof the main model used in this thesis. A detailed overview of the subje t an be

found in[37℄, [4℄, [5℄.

2.1 Modeling of reservoirs

LetuspresentequationsforathreedimensionalprobleminCartesianspa e, where

x

and

y

are horizontal oordinates,

z

is the verti al oordinate, and

t

is time.

Several phenomena play a role in reating a ow model of a reservoir with oil and

water,namely,themassbalan e on ept,Dar y'slaw,equationsofstate,andthe apillary

pressure on ept. Weexplain them further here.

flux out

flux in

V

Figure 2.1: Representation of avolume with a uxgoingthrough it.

The mass onservation lawgivesrise toone mass balan eequationper hemi al

om-ponent

c

. Let

C

c

be the mass on entration of omponent

c

in a unit volume

V

. Then, thereexists aux(arateof hangeofmass),

F

,transferringthemassthroughthevolume (Figure2.1). The amount of omponent

c

inthe volume an be expressed as:

Z ZZ

V

C

c

dV.

The mass on entration

C

c

an be expressed as

C

c

= φ

P

l

ρ

l

S

l

y

c,l

with

φ

- porosity,

ρ

l

-thedensity of phase

l

,

S

l

-thesaturation ofphase

l

,

y

c,l

-the massfra tion of omponent

c

inphase

l

. Saturationofphase

l

isafra tion ofthevolumeoftheporesthatiso upied by phase

l

; naturally,it isa value in the interval

[0, 1]

and the sum over allphases is

X

l

(18)

The hange in time of mass of omponent

c

is equal to

∂

∂t





Z ZZ

V

C

c

dV



 =

uxin

−

uxout

.

The dieren e inuxes is a surfa e integral of the ux over a boundary of

V

,

∂V

: uxin

−

ux out

= −

ZZ

∂V

F · n dA,

where

n

denotes a normalve tor (in the dire tion out of the volume), and the symbol

·

indi ates a dot produ t. This expression, a ording to the Divergen e Theorem, is equal

to:

−

ZZ

∂V

F · n dA = −

Z ZZ

V

∇ · F dV,

where

∇ · F

is the divergen e of

F

. This leadsto

∂

∂t





ZZZ

V

C

c

dV



 = −

Z ZZ

V

∇ · F dV,

and hen e

∂

∂t

C

c

= −∇ · F .

Finally,allowing additionalinput/outputsour es

q

leads to:

∂

∂t

C

c

= −∇ · F + q.

(2.1)

Assuming dead omponents, i.e., omponentsthat donot travel between phases, and

assumingthat thevarious hemi al omponentsinthe oil an belumped intoone 'pseudo

omponent', two basi equations arise: one for oiland one for water (subs ript

w

stands for water,

o

for oil):

C

w

= φρ

w

S

w

, C

o

= φρ

o

S

o

.

(2.2) We know that for density

ρ

and velo ity

v

, ux is equal to

F = ρv

, and we want to work out the velo ity. Dar y's law is based on a proportionality between the pressure

gradient,

∇p

, and the negative velo ity:

∇p ∼ −v

or,more pre isely,

−

_µ

k

∇p = v,

where

k

and

µ

denote the permeability tensor and vis osity, respe tively. This is an equation for the ase of one phase present. The two-phase Dar y's law dierentiates

between the water and oil equations:

−

_µ

k

w

(19)

−

_µ

k

o

k

ro

∇p

o

= v

o

.

Relativepermeabilities,

k

rw

∈ [0, 1]

forwaterand

k

ro

∈ [0, 1]

foroil,arefun tionsofwater saturation. These relations are found experimentally. Additionally, in 3D domains the

impa tof gravityhas tobe taken into a ountand the two-phase Dar y's lawbe omes:

−

_µ

k

w

k

rw

(∇p

w

− ρ

w

g∇d) = v

w

,

(2.3)

−

k

µ

o

k

ro

(∇p

o

− ρ

o

g∇d) = v

o

,

(2.4) where

g

isthegravitationala eleration,

ρ

w

and

ρ

o

arewaterandoildensities,respe tively, and

d

isdepth pointing verti ally downwards.

The ompressibilitiesfor water

c

w

,ro k

c

r

and oil

c

o

are dened as:

c

w

=

1 ρ

w

∂ρ

w

∂p

, c

r

=

1 φ

∂φ

∂p

, c

o

=

1 ρ

o

∂ρ

o

∂p

.

The apillarypressure (

p

c

) onstraintis relatedto the interfa ialtension between the phasesandthewettingpropertiesofthero kwhi hneedtobedeterminedexperimentally:

p

c

(S

w

) = p

w

− p

o

,

where

p

w

and

p

o

are water and oil pressures, respe tively.

Forea hphaseweinsertinDar y'slaw(2.3)and(2.4),and on entrationdenitions (2.2),

intoEquation (2.1), given that

F

l

= ρ

l

v

l

for

l ∈ {w, o}

:

∂

∂t

(φρ

l

S

l

) = −∇ ·

ρ

l

−

k

µ

l

k

rl

(∇p

l

− ρ

l

g∇d)

+ q

l

.

Usingthe hain rule gives:

φ

1 φ

∂φ

∂p

| {z }

c

r

∂p

∂t

ρ

l

S

l

+ φρ

l

1 ρ

l

∂ρ

l

∂p

| {z }

c

l

∂p

∂t

S

l

+ φρ

l

∂S

l

∂t

= ρ

l

∇ ·

k

µ

l

k

rl

(∇p

l

− ρ

l

g∇d)

+ q

l

.

It leads totwo oupled equationsthat need to besolved for

S

w

and

p

, [37℄, [5℄:

φ

S

w

(c

w

+ c

r

)

∂p

∂t

+

∂S

w

∂t

=

∂

∂x

k

µ

w

k

rw

∂p

∂x

− ρ

w

g

∂d

∂x

+

∂

∂y

k

µ

w

k

rw

∂p

∂y

− ρ

w

g

∂d

∂y

+

∂

∂z

k

µ

w

k

rw

∂p

∂z

− ρ

w

g

∂d

∂z

+ q

w

,

(2.5)

φ

(1 − S

w

)(c

o

+ c

r

)

∂p

∂t

−

∂S

w

∂t

=

∂

∂x

k

µ

o

k

ro

∂p

∂x

− ρ

o

g

∂d

∂x

+

∂

∂y

k

µ

o

k

ro

∂p

∂y

− ρ

o

g

∂d

∂y

+

∂

∂z

k

µ

o

k

ro

∂p

∂z

− ρ

o

g

∂d

∂z

+ q

o

,

(2.6)

(20)

where

p = p

w

= p

o

when apillaryfor esarenegle ted. Theinitial onditionsfor

p(x, y, z)

and

S

w

(x, y, z)

need tobespe ied and typi allythere isa no-ow onditionassumed at the boundaries.

These equations an be spatially dis retized and solved numeri ally with a nite

dif-feren e s heme. Typi ally, aregular omputationalgrid is used but advan ed simulators

might allow non-standard unstru tured grids. The equations shown here are 3D,

two-phase equations; gas an be added as a third phase. The presen e of aquifers, faults,

tra ers an alsobe a ounted for.

2.2 simsim

In this thesis a horizontaltwo-dimensional, two-phase (oil-water)reservoirmodel isused

where the gravity ee t an be negle ted. Therefore, Equations (2.5) and (2.6) be ome:

∂

∂x

k

µ

w

k

rw

∂p

∂x

+

∂

∂y

k

µ

w

k

rw

∂p

∂y

+ q

w

= φ

S

w

(c

w

+ c

r

)

∂p

∂t

+

∂S

w

∂t

,

∂

∂x

k

µ

o

k

ro

∂p

∂x

+

∂

∂y

k

µ

o

k

ro

∂p

∂y

+ q

o

= φ

(1 − S

w

)(c

o

+ c

r

)

∂p

∂t

−

∂S

w

∂t

.

Reservoirsimulatorsimsim (simplesimulator,[37℄)has been implementedand

devel-oped at Delft University of Te hnology.

After the equations are dis retized in spa e for the two-phase ow, four methods for

the time integration are available,[37℄, [5℄, [4℄: expli it Euler, impli itEuler with Pi ard

iteration, impli it Euler with Newton iteration, and IMPES (IMpli it Pressure Expli it

Saturation).

The reservoirdomainis2Dwith no-owthrough the boundaries. The wells(inje tors

or produ ers) an be implemented with pres ribed initialpressures and rate onstraints,

or with pres ribed initialrates and pressure onstraints.

As an example of how simsim operates, a

49 × 49

reservoir is simulated for a year. The hannelized permeability eld and homogeneous porosity eld are shown in the top

row of Figure 2.2. There, the inje tion (I) and produ tion (P) wells are indi ated as

white dots. The initial ondition for reservoir pressure

p

0

(Pa) and water saturation

s

0

is

p

T

₀

s

T

₀

T

=

_{3 · 10}

7 _{· 1}

T

49

2 ×

₁

0.2 · 1

T

49

2 ×

₁

T

, where

1

is a olumn ve tor of ones. The streamlines,[37℄,betweenthethreeinje torsattheleftboundaryandthethreeprodu ers

at the rightboundaryare shown in Figure2.3. The wells operate under pres ribed rates

(

±0.001 m

3 _/s

) and nopressure onstraints.

After one year, the pressure and saturation elds look like in the bottom row of

Figure2.2. The gradual saturation hange through time is shown in Figure2.4.

Thesimulateddatafromea hwell anbe olle tedandplottedagainsttime. Figure2.5

shows the bottom hole andthe well gridblo kpressures together with the owrates and

water saturation inwells. Cumulativeprodu tion data over the whole eld are shown in

(21)

Permeability field (log10 k , m

2 )

−13.5

−13

−12.5

Porosity field (

φ

, −)

0

0.5

1 Final pressure field (Pa)

2

4

6

8

10

12

14 x 10

7 Final oil saturation field (−)

0.2

0.4

0.6

0.8 I

P

I

Figure2.2: Permeability and porosityelds - top row, pressure and saturation afterone year

-bottom row.

Streamlines

• I

• P

(22)

0

100

200

300

400

0

0.5

1

1.5

2 x 10

8 Time t , days

Bottom hole pressure

p

wf

,

Pa

0

100

200

300

0

0.5

1

1.5 x 10

−3

Time t , days

Well flow rates

|q

o

|, |q

w

| ,

m

3 /s

0

100

200

300

400

0

5

10

15 x 10

7 Time t , days

Well grid block pressure

p ,

Pa

0

100

200

300

400

0

0.5

1 Time t , days

Well grid block water saturation

S

w

,

−

(23)

0

100

200

300

0

1

2

3

4 x 10

−3

Time t , days

Field production rates

−q

o,f

, −q

w,f

,

m

3 /s

oil

water

0

100

200

300

400

0

2

4

6

8 x 10

4 Time t , days

Cum. field pr. rates

−q

o,c,f

, −q

w,c,f

,

m

3 oil

water

0

100

200

300

0

1

2

3

4 x 10

−3

Time t , days

Field injection rate

q

w,f

,

m

3 /s

0

100

200

300

400

0

5

10 x 10

4 Time t , days

Cum. field inj. rate

q

w,c,f

,

m

3

(24)

Let

x

k

be a ve tor ontaining all the dynami model variables at a dis rete time

k

. In reservoir engineering it would typi ally be dis retized pressure and saturation per

grid blo k:

x

k

=

p

T

_k

s

T

_k

T

. Stati variables, dis retized permeability and porosity per grid blo k, are stored in ve tor

m

=

k

T

φ

T

. The deterministi nonlinear model representation an be formulated as:

x

_k+1

= f

k→k+1

(m, x

k

).

(2.7)

Here,

f

k→k+1

(m, x

k

)

isamodeloperatorthatdependsonspe iedstati parameters

m

and propagatesdynami variables

x

k

fromtime

k

to

k + 1

. Aninitial ondition

x

0 =

p

T

₀

s

T

₀

T

(25)

Data assimilation and parameter

estimation

This hapter presents sequential ensemble data assimilation and parameter estimation

methods that have been implemented for appli ationsin this thesis. The des riptions of

Kalman, ensemble Kalman and ensemble square root algorithms are presented.

Imple-mentation improvementsand parameter estimation methods on lude the hapter.

3.1 Problem setup

Thetaskofdataassimilationistoimprovetheresultsofanumeri almodelofthe

phenom-enaone isinterested inwith theavailableobservations. Even thoughthe physi s ofmany

pro esses may be well-understood, the numeri al models are often un ertain due to

dis- retization errors and modeling approximations. Therefore, the model from Se tion2.3,

Equation (2.7), an in general ontain a sto hasti term expressing the model

un er-tainty. Moreover, ameasurementequation isadded for afull des riptionof a state-spa e

representation as given by:

x

_k+1

= f

k→k+1

(m, x

k

) + ε

k+1

,

y

k+1

= h(x

k+1

) + ν

k+1

.

(3.1)

Theinitialstate

x

0

andthemodelparameters

m

needtobespe ied;

h

isameasurement operatorexpressing the relationbetween the state

x

k+1

and observations

y

k+1

ata given time; and

ε

and

ν

are model noise and observation noise, respe tively. The noise in the model isassumed tobenormally distributed

ε

∼ N(0, Q)

with a model error ovarian e matrix

Q

; the noise in the data is assumed to be normally distributed

ν

∼ N(0, R)

withanobservation error ovarian ematrix

R

;and the modeland observation errorsare independent. The error ovarian e matri esmightdepend on time

k

.

The measurementoperator

h

an depend ontimewhen dierent datatypesare avail-able atdierent times. Large s ale measurements like seismi observationsare expensive

and not sensitive to small s ale hanges, therefore, they are olle ted relatively less

fre-quently than well observations, for example. In real-life appli ations data ome from

instruments and an be ontaminated with noise due tohuman or devi e error.

Upon observation arrival we an assess how well the model predi ts just a quired

(26)

dierent ways. A ommonstartingpointisthe spe i ationof anobje tivefun tion. Let

usdene a ve tor

x

[k]

that ontains states

x

k

for alltime steps

k

, and thenthe obje tive fun tion tobeminimizedis

J x

[k]

=

1

2 X

k

(y

k

− h(x

k

))

T

R

−

₁

(y

k

− h(x

k

)) +

1

2 X

k

(x

b

_{− x}

k

)

T

B

−

₁

(x

b

_{− x}

k

) +

1

2 X

k

(x

k+1

− f(m, x

k

))

T

Q

−

₁

(x

k+1

− f(m, x

k

)).

(3.2)

We are looking for

x

∗

[k]

= argmin

x

[k]

J x

[k]

, where

J

is a nonnegative s alar fun tion

J : R

n

x[k]

_{−→ R}

+

_{∪ {0}}

with

n

x

[k]

thesizeofve tor

x

[k]

. Thefun tionin ludesaquadrati data mismat h weighted by the noise ovarian e matrix

R

, a quadrati ba kground mis-mat h measuring how dierent the state

x

k

is fromsome ba kground state

x

b

, weighted

by a given ovarian e matrix

B

, and a quadrati model mismat h weighted by a given ovarian e matrix

Q

.

3.2 The lassi al Kalman lter

Kalman ltering is a sequential optimization pro edure that omprises of two steps: a

fore ast step and an update step, alternatingly applieduntil the last update time. First,

the fore ast step integrates the model to the rst update time to olle t the predi ted

observations, then the update is performed with the aid of the data, and the model

ontinues with the new updatedparameters untilthe next update time.

Let us onsider a sto hasti system like in Equation (3.1) but with linear model and

observation operators:

x

_k+1

= F

k

x

k

+ ε

k+1

,

y

k+1

= H

k

x

k+1

+ ν

k+1

,

(3.3)

where, as before,

ε

∼ N(0, Q)

and

ν

∼ N(0, R)

, and model parameters have been omitted. For simpli ity of notation, we use

F

k

and

H

k

without the subs ripts from here on. Then, the obje tive fun tion inEquation (3.2) be omes:

J x

[k]

=

1

2 X

k

(y

k

− Hx

k

)

T

R

−

₁

(y

k

− Hx

k

) +

1

2 X

k

(x

b

_{− x}

k

)

T

B

−

₁

(x

b

_{− x}

k

) +

1

2 X

k

(x

k+1

− Fx

k

)

T

Q

−

₁

(x

k+1

− Fx

k

).

(3.4)

(27)

At ea h updatetime wewantto ompute adistributionof state

x

k

giventhe data

y

k

. This distribution an be shown tobeGaussian,[38℄. Then, the estimateof the mean

x

¯

k

, andthe ovarian e

P

k

ompletelydes ribeitsshape. Wewant toderivethe equationsfor

¯

x

k

and

P

k

that are time and measurement updates, respe tively, where

x

¯

k

is anoptimal state estimateand

P

k

is its orrespondingminimizederror attime

k

.

We will dierentiate between fore asted and updated variables using supers ripts

f

and

a

(for analyzed), respe tively. In this way,

x

f

k

denotes fore asted dynami variables attime

k

and

x

a

k

denotes analyzed dynami variablesat the same time after the update step;

¯

x

f

k

and

x

¯

a

k

denotefore asted and analyzedmean estimates at time

k

,respe tively. Sin e the model error mean is

0

, the mean estimate

x

¯

k

propagates in time as given by:

¯

x

k

time update

¯

x

k+1

= F¯

x

k

.

(3.5)

Sin ethedenitionofa ovarian ematrix anbewrittenas

P

k

= E[(x

k

− ¯x

k

)(x

k

− ¯x

k

)

T

_]

,

letusdene the fore asted and updated ovarian e matri es:

P

f

_k

= E[(x

k

− ¯x

f

k

)(x

k

− ¯x

f

k

)

T

],

P

a

_k

= E[(x

k

− ¯x

a

k

)(x

k

− ¯x

a

k

)

T

],

forthe given mean estimates

x

¯

f

k

and

x

¯

a

k

,where

E[·]

denotes the expe ted value. To estimatethe forward ovarian epropagation we ompute:

x

_k

_{− ¯x}

_k

(3.3),(3.5)

=

Fx

_k−1

+ ε

k

− F¯x

k−1

= F(x

k−1

− ¯x

k−1

) + ε

k

,

and this gives

P

k

time update

P

k

= FP

k−1

F

T

+ Q.

Let

k

beaxed observation time. A Kalmanlteraims atprovidinganoptimalstate estimate given the obje tive fun tion at the observation time

k

. The obje tive fun tion inEquation (3.4) for the

x

¯

a

k

estimatebe omes:

J(¯

x

a

_k

) =

1

2 (y

k

− H¯x

a

k

)

T

R

−

₁

(y

k

− H¯x

a

k

) +

1

2 (x

f

k

− ¯x

a

k

)

T

(P

f

k

)

−

₁

(x

f

_k

_{− ¯x}

a

_k

),

wherethe model mismat hterm disappears due toEquation (3.5). Here the ba kground

mismat hterm from Equation (3.4) is a distan e to the fore ast

x

f

k

, and

P

f

k

denotes the errorin this fore ast mismat h.

A ne essary ondition for nding the obje tive fun tion's minimum is that the

fun -tion'sgradientvanishes:

∇

¯

x

a

k

J = 0.

We ompute the estimate

¯

x

a

k

. Following[38℄, we get:

−H

T

_R

−

₁

(y

k

− H¯x

a

k

) + (P

f

k

)

−

₁

(¯

x

a

_k

_{− ¯x}

f

_k

) = 0,

(28)

¯

x

a

_k

=

H

T

R

−

1 H

+ (P

f

_k

)

−

1 −

₁

H

T

R

−

1 y

_k

+ (P

_k

f

)

−

1 x

¯

f

_k

=

RH

−

T

H

T

R

−

1 H

+ RH

−

_T

(P

f

_k

)

−

₁

−

₁

y

k

+

P

f

_k

H

T

R

−

1 H

+ P

f

_k

(P

f

_k

)

−

₁

−

₁

¯

x

f

_k

=

H

+ RH

−

_T

(P

f

_k

)

−

₁

−

₁

y

k

+

P

f

_k

H

T

R

−

1 H

+ I

−

₁

¯

x

f

_k

.

(3.6)

Wewant to workout the terms infrontof

y

k

and

x

¯

f

k

in Equation(3.6). Wehave

H

+ RH

−

_T

(P

f

_k

)

−

₁

−

₁

=

h

H(H

−

_T

(P

f

_k

)

−

₁

)

−

₁

+ R

H

−

T

(P

f

_k

)

−

₁

i

−

₁

= P

f

_k

H

T

(HP

f

_k

H

T

+R)

−

₁

,

(3.7)

and we want toshow that

(P

f

_k

H

T

R

−

1 H

+ I)

−

₁

= I − P

f

k

H

T

(HP

f

k

H

T

+ R)

−

₁

H

(3.8) by simple multipli ation:

P

f

_k

H

T

R

−

1 H

+ I

I

_{− P}

f

_k

H

T

(HP

f

_k

H

T

+ R)

−

₁

H

=

I

+ P

f

_k

H

T

R

−

1 H

_{− P}

f

_k

H

T

(HP

f

_k

H

T

+ R)

−

₁

H

_{− P}

f

_k

H

T

R

−

1 HP

_k

f

H

T

(HP

f

_k

H

T

+ R)

−

₁

H

=

I

+ P

f

_k

H

T

R

−

1 H

_{− P}

f

_k

H

T

R

−

1 R(HP

f

_k

H

T

+ R)

−

₁

H

_{− P}

f

_k

H

T

R

−

1 HP

_k

f

H

T

(HP

f

_k

H

T

+ R)

−

₁

H

=

I

+ P

f

_k

H

T

R

−

1 H

_{− P}

f

_k

H

T

R

−

1 (R + HP

f

_k

H

T

)(HP

f

_k

H

T

+ R)

−

₁

H

= I.

Substituting Equations(3.7) and (3.8) in Equation(3.6), we obtain anexpression for

the update equation inthe Kalmanlter, [38℄:

¯

x

k

measurement update

¯

x

a

_k

=

H

+ RH

−

T

(P

f

_k

)

−

1 −

₁

y

_k

+

P

f

_k

H

T

R

−

1 H

+ I

−

₁

¯

x

f

_k

=

P

f

_k

H

T

HP

f

_k

H

T

+ R

−

₁

y

k

+

I

_{− P}

f

_k

H

T

(HP

f

_k

H

T

+ R)

−

₁

H

¯

x

f

_k

=

¯

x

f

_k

+ P

f

_k

H

T

(HP

f

_k

H

T

+ R)

−

₁

(y

k

− H¯x

f

_k

).

(3.9) The term

K

= P

f

_k

H

T

(HP

f

_k

H

T

+ R)

−

1

(3.10) is alledthe Kalmangain.

The updated ovarian e matrix

P

a

k

= E[(x

k

− ¯x

a

k

)(x

k

− ¯x

a

k

)

T

]

for the given mean estimate

x

¯

a

k

an nowbe omputed from

x

k

− ¯x

a

k

(3.9)

=

x

_k

_{− ¯x}

f

_k

_{− K(y}

_k

_{− H¯x}

f

_k

)

(3.3)

=

x

k

− ¯x

f

_k

− K(Hx

k

+ ν

k

− H¯x

f

_k

) =

(I − KH)(x

k

− ¯x

f

k

) − Kν

k

,

whi hgives

(29)

P

k

measurement update

P

a

_k

=

(I − KH)P

f

k

(I − KH)

T

− KRK

T

=

(I − KH)P

f

k

− (I − KH)P

f

k

H

T

K

T

− KRK

T

=

(I − KH)P

f

k

− K(K

−

₁

P

f

_k

_{− HP}

f

_k

+ RH

−

_T

)H

T

K

T

(3.10)

=

(I − KH)P

f

k

− K

(HP

f

_k

H

T

+ R)H

−

_T

(P

f

_k

)

−

₁

P

f

_k

_{− HP}

f

_k

+ RH

−

_T

H

T

K

T

=

(I − KH)P

f

k

.

Kalmanltering,[40℄, was originallydeveloped for linear,Gaussian problems. It

pro-videsthe optimalestimateof thestateof thesystem andthe ovarian eofthe estimation

error. It is also able to propagate these statisti s in time. Models are, however, rarely

linear and variables are Gaussian only in some spe i ases. Lately, for omputational

reasonsandtoallowfor nonlinearmodels, the ensemble Kalmanlter(EnKF) was

intro-du ed, [16℄, [20℄,and be amea new standard for sequential data assimilation.

3.3 The Ensemble Kalman Filter - EnKF

Wepresentasequentialensemble-basedalgorithmfornonlinearGaussianpro esses,re all

Equations(3.1):

x

k+1

= f

k→k+1

(m, x

k

) + ε

k+1

,

y

k+1

= h(x

k+1

) + ν

k+1

.

TheensembleKalmanlter,[16℄,representsthedistributionofthestateve tor

x

k

∈ R

n

s

×

1

by a sample (a sample from the distribution of interest), i.e., a olle tion of possible

realizations,alsoknown asan ensemble, with

n

e

members:

X

= [x

1 x

2 ... x

n

e

] ∈ R

n

s

×

n

e

_.

The time index has been omitted for larity sin e all the variables onsidered here are

taken at the same time step, namely,a dis rete update step.

Let

Y

be a matrix holding

n

e

opies of the observation ve tor:

Y

_{= [y y ... y] ∈ R}

n

o

×

n

e

_,

where

n

o

isthe numberof observation points. The ve tor ofobservations

y

isaninputto the data assimilationalgorithm.

Let

Y

ˆ

beamatrixholdinganensembleofpredi tedmeasurementsfromea hrepli ate:

ˆ

Y

= [h(x

1 _{) h(x}

2 _{) ... h(x}

n

e

_{)] ∈ R}

n

o

×

n

e

_.

This olle tion of ve tors is a result of integrating the given model to the urrent

up-datetimeforea hensemble member,

n

e

-times,and omputingthe fore asted observation values. For ompli ated models this fore ast step an be very time onsuming for large

(30)

X

a

= X

f

_{+ K(Y − ˆ}

_Y),

and the Kalman gain

K

∈ R

n

s

×

n

o

isequal to:

K

= Cov(X

f

, ˆ

Y)[Cov( ˆ

Y, ˆ

Y) + R]

−

₁

,

where

Cov(·, ·)

denotes a sample( ross-) ovarian e matrix.

This algorithm was initially des ribed in [16℄ in 1994. The original formulation was

wrongandits orre tionwaspublishedin[8℄whereitisexplainedthattopreserve orre t

statisti s the observations need to be perturbed. Sin e then a number of improvements

havebeen introdu edof whi hmany were applied toreservoir engineering,[1℄. Let

Y

be anensemble of perturbed measurements:

Y

= [y + ν

1 y

+ ν

2 ... y + ν

n

e

] ∈ R

n

o

×

n

e

_,

where

ν

isameasurementerrorsampledfromanormalzeromeandistributionwithagiven ovarian e matrix

R

that expresses our belief inthe un ertainty of the measurements.

The setup of the ensemble Kalman lter, where the distribution of a variable is

rep-resented by a sample, is easy to implement. Even though it has been used su essfully

inmany appli ations,it isknown to ause ompli ations andsome of them are dis ussed

here further.

The square rootimplementationof the ensemble Kalmanlter is presented sin e itis

the version used throughout this work. The basi idea of the ensemble Kalman lter is

anintrodu tiontoalmost everymajorse tion forthe sake of ompleteness. It mightalso

be presented fromdierent (but equivalent)angles.

3.3.1 The Ensemble Square-Root Filter - EnSRF

In the lassi al ensemble Kalman lter as formulated in the previous se tion, the

on-stru tion of the ensemble of perturbed measurements

Y

ontains noise that an be an additionalsour e of sampling error for small ensemble sizes, [19℄. Therefore, square root

algorithms were developed to avoid the use of perturbed measurements. Note that the

versionpresented here in ludesamatrix of perturbationsas asamplerepresentation of a

ovarian e matrix, [18℄, whi h lowers the omputationalburden.

Square root algorithms use an ensemble representation like the EnKF, but update

the mean and deviations from the mean separately as in the traditional Kalman lter.

In fa t, both equations, the update of the mean and the update of the perturbations,

ould be dire tly derived fromthe lassi alKalman lter equationsbut with the

modi- ationsregarding theensemble representations. Asquare rootversionof thelter isused

throughout this thesis. The basis of the algorithm was taken from [18℄, Se tion7.4.2 1

therein, and improved as des ribed in [70℄ and [51℄. Other versions of square rootlters

exist, see for example[83℄.

Let usagain dene:

• X = [x

1 _x

2 _{... x}

n

e

_{] ∈ R}

n

s

×

n

e

- anensemble of state ve tors,

1

(31)

• Y = [y + ν

1 _y

_{+ ν}

2 _{... y + ν}

n

e

_{] ∈ R}

n

o

×

n

e

-anensemble of perturbed measurements,

• ˆ

Y

= [h(x

1 _{) h(x}

2 _{) ... h(x}

n

e

)] ∈ R

n

o

×

n

e

- anensemble of predi ted measurements.

Additionally,we need:

• 1

n

∈ R

n×n

- a matrix whereea h element isequal to

1 n

,

• I

- identity matrix of a propersize,

• ¯

X

= X1

n

e

∈ R

n

s

×

n

e

- a matrix storingthe ensemble mean

n

e

-times,

• X

′

= X − X1

n

e

= X − ¯

X

- anensemble-perturbation matrix,

• E = [ν

1 _ν

2 _{... ν}

n

e

_]

- an ensemble of measurement perturbations,

• S = ˆ

Y

_{− ˆ}

Y1

n

e

-a matrix holdingperturbationsof the predi ted measurements. We want the above variablestobeunderstood as(derived from)fore asted variables,

for the ease of notation omitting the supers ript

f

. The EnSRF denes the updated ensemble asa sum ofupdated ensemble meanand updated ensemble perturbations, that

is:

X

a

= ( ¯

X)

a

+ (X

′

)

a

and

( ¯

X)

a

= X1

n

e

+ X

′

S

T

X

1 (I + Σ

2

1 )

−

₁

X

T

₁

(Y1

n

e

− ˆ

Y1

n

e

),

(X

′

)

a

= X

′

V

2 q

I

_{− Σ}

T

₂

Σ

2 V

T

2 .

The last multipli ation by

V

T

2

provides the unbiasedness of the lter, [70℄, [51℄, [19℄, sin e it was shown that the original algorithm without the last multipli ation does not

onserve the mean. To ompute these equationswe need toknow matri es:

Σ

1

,

X

1

,

Σ

2

, and

V

2

. Let the operator

tSVD

=

denote the thin singular value de omposition, [27℄, then the following sequen e of equations ompletes the EnSRFupdate step:

1.

S

tSVD

= U

0 Σ

0 V

0 T

,

2.

X

0 = Σ

−

₁

0 U

T

0 E,

3.

X

0

tSVD

= U

1 Σ

1 V

1 T

,

4.

X

1 = U

0 (Σ

−

₁

0 )

T

U

1 ,

5.

X

2 = (I + Σ

2

1 )

−

1

2 X

T

1 S,

6.

X

2

tSVD

= U

2 Σ

2 V

2 T

.

Note that the equation for the perturbationupdate

(X

′

)

a

omes fromthe ovarian e

measurementupdateinthe lassi alKalmanlter

P

a

_{= (I−KH)P}

f

sin e

P

a

_{= (X}

′

)

a

_[(X

′

)

a

_]

T

,

P

f

= X

′

(X

′

)

T

(32)

P

a

=

_{(I − KH)P}

f

P

a

=

P

f

_{− P}

f

H

T

(HP

f

_H

T

_{+ R)}

−

₁

HP

f

(X

′

)

a

_[(X

′

)

a

_]

T

₌

_X

′

(X

′

)

T

_{− X}

′

(X

′

)

T

_H

T

_(HX

′

(X

′

)

T

_H

T

_{+ R)}

−

₁

HX

′

(X

′

)

T

(X

′

)

a

_[(X

′

)

a

_]

T

₌

_X

′

I

_{− (X}

′

)

T

_H

T

_(HX

′

(X

′

)

T

_H

T

_{+ R)}

−

₁

HX

′

(X

′

)

T

where we used Equation (3.10)for the Kalman gain. Then, if we write

I

_{− (X}

′

)

T

H

T

(HX

′

(X

′

)

T

H

T

+ R)

−

1 HX

′

= TT

T

we obtain

(X

′

)

a

[(X

′

)

a

]

T

= X

′

TT

T

(X

′

)

T

= X

′

T(X

′

T)

T

.

Finally,the square rootperturbation update in ageneral formis:

(X

′

)

a

= X

′

T.

Matrix

T

is alleda square roottransformation matrix and itis not unique, see [51℄ and [70℄ for a detaileddis ussion.

3.4 Implementation issues of the Ensemble Kalman

Fil-ter

The EnKF or EnSRF even in the most e ient form will almost never give satisfa tory

results when implemented for the rst time ina new appli ation. To be able to apply it

su essfully,dierentmodi ations,[1℄,[3℄,[19℄, anbeusedthatdependontheproblems

en ountered inagiven setup. Anumberof improvements usedin our implementationsis

presented inthis se tion.

In today's resear h we want not only the dynami model variables

x

k

to be esti-matedbut alsoadditionalun ertain parameters

m

. When

x

k

follows the modelequation

x

k+1

= f

k→k+1

(m, x

k

)

,

m

doesnot hangeintime,i.e.,

m

k+1

= Im

k

where

I

isanidentity matrix.

Updating all variables

Typi ally, all variables, dynami and stati , are updated and they omprise the state