Metoder for lokal data-assimilering relatert til varsling av turbulens ved flyplasser

1 Innledning 2

2 Problemstilling 2

3 Lokal data-assimilering 3

4 Invers problemstilling 6

4.1 Invers løsning basertpå observasjonsdata . . . 6

4.2 Invers løsning basertpå forenklet 3d-Var . . . 8

5 Varsling basert på forenklet 3d-Var 10 6 Variasjonsmetode for estimering av startbetingelse og tidsvarierende randbetingelser 12 6.1 Analyse iforkantav prediksjon . . . 12

6.1.1 Eksempelmedsinusformetbakgrunnsløsningogenikke-systematisk amplitudefeil . . . 13

6.1.2 Eksempel med sinusformet bakgrunnsløsning og en systematisk faseforskyvning . . . 14

6.1.3 Hvor mange målingertrengs? Bakgrunnsløsning med både fase- og amplitudefeil . . . 15

6.1.4 Eksempel med sinusformet bakgrunnsløsning med variabel ampli- tude og systematisk fasefeil . . . 16

6.1.5 Lineær, kvadratiskog kubisk inuens-prol. . . 16

6.2 Ekstrapolering av analysekurven inn i prediksjonsintervallet . . . 17

Referanser 26 A Forenklet 3d-Var: Lokal data-assimilering 27 A.1 Generell kostfunksjon . . . 27

A.2 Eksempel med ett observasjonspunkt . . . 27

A.2.1 Fordeling av høyde-korreksjon . . . 28

A.3 Tidsfaktoren ved korreksjon av randverdier . . . 29

Referanser 29 B Sensitivitetsanalyse - forenklet 2D-analyse 30 B.1 Navier-Stokes ligningerfor plan strømning . . . 30

B.2 Førsteordens sensitivitetsligningerfor

c

^{. . . . . . . . . . . . . . . . . . . . 32}

Problemetmedbrukavdata-assimileringimodellsimuleringkanseessomtodelt:Etinvers

problem hvor observasjonsdata benyttes alene eller sammen med en bakgrunnsløsning

for å rekonstruere initial- og randbetingelser, og et varslings/prediksjonsproblem hvor

observasjonsdata bare er kjent i starten, men bakgrunnsløsningen kan være gitt også i

varslingsintervallet.

Det spesielle med lokal data-assimilering er at området er relativt lite, slik at trans-

porttiden gjennom området blir kort ved sterk vind. Dette gjør at randbetingelser blir

relativtviktigere enn initialbetingelser,ogkonvensjonelle assimilasjons-metoderermind-

re egnet. Et tilleggsmoment ved lokal data-assimilering er relativt få data-punkter i de

estetilfeller.Med bakgrunnidennetypen vurderingerharvivalgtåtesteenformulering

basertpå forenklet 3d-Var,og gjennomført noen elementære numeriske testsimuleringer.

En konklusjon fra denne analysen erat korreksjon av randbetingelser i prediksjonsin-

tervallet er et hovedproblem for en eektiv bruk av lokal data-assimilering. I tråd med

ideen fra 4d-Var kan man benytte et tidsvindu for å analysere trenden i avvik mellom

data og tidligereprediksjon. Dette kanså benyttes i den videre korreksjonen av randbe-

tingelsene i prediksjonsperioden. Rapporten inneholder en analyse av slike vurderinger,

inkludertnoenstilisertetesteksempler.Dette arbeidettenkesgeneralisertfortestingimer

realistiske situasjoner.

2 Problemstilling

Vibenytter en enkel tre-dimensjonalgeometri som representerer en symmetriskfjelltopp

(Hunt & Snyder 1980). Det er valgt en relativt grovmasket gridoppløsning (40 x 40 x

35), men god nok til å representere hovedstrukturen i strømningen, og beregningene er

gjortmed denserielle versjonen av SIMRA.Alleberegningeneergjennomførtfor nøytral

strømning (uten stratikasjon). For den aktuelle testingen antas følgende størrelser for

geometriogvindstyrke:

Fjellhøyde H= 500 m

Høydevind

U ∞ = 20

^m/s

Vindretningen antas å variere sinus-formet med en tilleggsforstyrrelse som gir inhomo-

gen fristrømning med vindskjær. Hovedretningen ligger i området mellom sør-vestlig og

nord-vestlig, og perioden på dreiningen er 3 timer, noe som også tilsvarer hele predik-

sjonsperioden. Med gitt hastighetsamplitude

u a

^{gir dette følgende vindfelt i} uforstyrret del av feltet (dvs. høydevind og randverdierlangt borte fra fjellet):

u = u a cos θ; v = u a sin θ

⁽¹⁾

hvor retningen er

θ = (π/4) sin ωt + 0.01(ǫ x + ǫ y ),

⁽²⁾

frekvensen

ω = 2π/T

^{, og}perturbasjonen

(ǫ x , ǫ y )

^er proporsjonal med hhv.koordinatene i de to horisontale retningene. Hastighetsamplitudenlangs innløpsrendene antas å variere

logaritmiskmed høyden:

u a (z) = min

U ∞ , u ∗

0.4 ln z z o

(3)

Herrepresenterer

u ^∗

friksjonshastighet og

z o

overateruhet. Perioden på vindepisoden er

T = 3

^timer

= 10800

^s

Figur 1: Øyeblikksbilde av roterende vind- og turbulensfelt. Rotasjonsretningen for

fristrømningenerfra nordvestlig tilsørvestligvind. Horisontalsnitt i z= 0.4 H overbak-

ken.

3 Lokal data-assimilering

I det følgende gjennomføres simuleringstester for å undersøke en enkel form for lokal

data-assimilering.Globale størrelser som inngår er:

Eksakt løsning:

x t

Observasjonsverdier:

y

x

I de numeriske eksperimentene antas den 'eksakte' løsningen

x t

^{å være lik løsning av}

RANS-modellen(SIMRA) medbruk av initial-ograndbetingelserspesisertfraligninge-

ne(1-3)foran.Figur1illustrererdenneløsningenforettidspunkthvorhøydevindenligger

i sørvestlig retning, oger iferd med å dreie mot sørlig retning.Hele simulerings-syklusen

erogså dokumentert som animasjon.

Bakgrunnsløsningen

x b

^{beregnes her som en} perturbasjon av

x t

^{: Vi har valgten tilfeldig}

løsning med fase-forsinkelse på 15 min og et noe mindre utsving sammenlignet med

x _t

^.

Forskjellen mellomde to løsningene indikeres i gur 2, som viser beregnete høydeverdier

for

x t

^og

x b

ⁱ høydepunktet (3.2, 0.0, 2.2)h.

Fra fysiske vurderinger antas at det primært er høydedata som er interessant for bruk i

den lokaledata-assimileringen.Den enkleste tilnærmingenerdaåanta atvi har bareett

observasjonspunkt(f.eks.fraypåveimotlanding),eventuelt noenfå observasjoner som

slåssammen.Idisseeksperimentenevelgesetobservasjonspunktsomliggerietuforstyrret

områdeover fjelltoppen, med koordinatverdier

(x, y, z) = (3.2, 0.0, 2.2)h.

Fradetteobservasjonspunktet tasdetutdatafradeneksakte løsningen

x _t

^{,dvs.en tidsse-}

rieavdataforhastighetsfeltetidettepunktet,ogdettegirdermedobservasjonsvektoren

y

^.

For å estimere en korrigert løsning benytter vi en forenklet versjon av variasjonsformu-

leringen 3d-Var, med bruk av bare ett observasjonspunkt. Dette gir følgende korrigerte

løsningeller analyse-løsning (se Appendiks):

x a = x b + y k − x b,k

σ ² _b + σ _o ² [B 1k , B 2k , ..., B N k ] ^T

⁽⁴⁾

hvor

σ b

representerer en bakgrunns-kovarians og

y _k

^er observasjonsdata i gridpunkt 'k'.

Vi antar i tillegg athøydedata sprer informasjontilnærmet uniformt isamme høyde i et

lokaltområde.Dermed kanden korrigerteløsningen forenkles til formen

x _a = x _b + (y k − x _b,k )f (z/H ) σ _b ²

σ ² _b + σ _o ² [1, 1, ..., 1 h ] ^T

⁽⁵⁾

hvor

f (z/H )

^{eren antatt vertikal fordeling(jfr.} Appendiks).

Ethovedpunktvedlokaldata-assimileringeratrandbetingelseneerrelativtmyeviktigere

enn initialbetingelsene. Typisk tar det bare 10 - 20 min. før en partikkel transporteres

gjennom hele det lokale området, og 'hukommelsen' om det opprinnelige initialfeltet er

like kort. Derimot styrer randbetingelsene løsningen gjennom hele prediksjonsperioden.

Problemet blir da hvordan vi kan korrigere randbetingelsene i selve varslingsintervallet.

Formuleringen ovenfor (og alle klassiske assimilerings-formuleringer) gjelder i prinsippet

bareforetstarttidspunkt,ogvierdermedforeløpighenvisttilgroveantakelser,jfr.seksjon

5for et numerisk eksperiment.

0.4 0.6 0.8 1 1.2 1.4 1.6 1.8

0 20 40 60 80 100 120 140 160 180

Solutions, u/U_inf

Prediction time (min)

exact solution, u-comp.

background solution, u-comp.

-1 -0.5 0 0.5 1

0 20 40 60 80 100 120 140 160 180

Solutions, v/U_inf

Prediction time (min)

exact solution, v-comp.

background solution, v-comp.

Figur 2: Bakgrunnsløsning

x _b

^{og eksakt løsning}

x _t

^. Hastighetskomponenter i høydepkt.

(3.2,0.0,2.2)som funksjon av varslingstid.

Somet første punkt ianalysen ser vi på følgendeproblem: Antaat observasjonsdata iet

høydepunkterkjentforhele tidsperiodensom skalsimuleres,sammenmedbakgrunnsløs-

ningen.Vi glemmerinitial-ograndbetingelsenefor deneksakte løsningen, ogkonstruerer

tilnærmetebetingelser vhabakgrunnsløsningen ogobservasjonsdataene.

Denne problemstillingen skillerseg fra situasjonen med varsling/prediksjon ved atvi her

antar at observasjonsdata er gitt for hele tidsperioden, mens dette naturligvis ikke er

tilfelle i en varslingssituasjon, hvor slike data bare er gitt fram mot tidspunktet hvor

varslingenbegynner.

4.1 Invers løsning basert på observasjonsdata

Idette tilfelletlegges allvektpåobservasjonsdatanår vi skal konstruere initial- ogrand-

betingelser, og bakgrunnsløsningen neglisjeres. Vi antar altså at observasjonsdata er til-

gjengeligfor hele simuleringsperioden.Istedenfor ligning(5) fåsdermedfølgende relasjon

forinitial- ograndbetingelser:

x _a = y _k f(z/H )[1, 1, ..., 1 h ] ^T

Initial- og randverdier i samme høydenivå som observasjonspunktet 'k' får uniformt

sammeverdi, mens den vertikale fordelingen

f (z/H )

^antas logaritmisk.

Figur 3 illustrerer denne løsningen i et kontrollpunktet (2.2 h, 0.4 h, 0.4 h), og viser

rimelig overensstemmelse med den eksakte løsningen som også er inkludert i guren.

Inhomogeniteten i høydevinden for den 'eksakte' løsningen er den viktigste årsaken til

forskjellen mellom disse løsningene. (Tilsvarende beregninger med en uniform, homogen

høydevindgir løsningersom ligger betydelignærmere hverandre.)

0 0.5 1 1.5 2

0 20 40 60 80 100 120 140 160 180

Solutions, u/U_inf

Prediction time (min)

obs.based solution, u-comp exact solution, u-comp.

-1 -0.5 0 0.5 1 1.5

0 20 40 60 80 100 120 140 160 180

Solutions, v/U_inf

Prediction time (min)

obs.based solution, v-comp exact solution, v-comp.

Figur3: Sammenligningmellomløsning basert på observasjonsdata

x a

^{ogeksakt løsning}

x _t

^. Hastighetskomponenter i kontrollpkt (2.2,0.4,0.4)som funksjon av tid.

I dette tilfellet måvi antanoe om forholdet

σ b /σ o

^{, dvs. forholdet mellom}standardavvik (eller varians) for hhv. bakgrunnsløsning og observasjonsdata. Altså gjøre en antakelse

om hvor stor vekt vi vil legge på observasjonsdata sammenlignet med bakgrunnsløsnin-

gen, og observasjonsdata antas gitt for hele simuleringsperioden. Dersom all vekt legges

på observasjonsdata, blir løsningen som gitt i forrige avsnitt. Dersom all vekt legges på

bakgrunnsløsningen,er løsningen gittav

x _b

^{(jfr. g. 2).}

Her velger vi å legge samme vekt på bakgrunnsløsning og observasjonsdata, samtidig

somden vertikalefordelingenav korreksjonenvelges logaritmisksom tidligere.Initial-og

randbetingelser erda gittfra (5) på formen

x a = x b + 1

2 (y k − x b,k )f(z/H )[1, 1, ..., 1 h ] ^T

Det bør merkes atvi her også har valgt å beholde samme vekt på beregning av randbe-

tingelseneoverhele prediksjonsperioden. Dette eropplagtet åpentspørsmål,men idette

spesielle tilfellethar vi ingen tilleggsopplysninger som kan gi andre preferanser. Med an-

tatt samme vertikale fordeling for

x _a

^,

x _b

^og

y _k

^{, gir dette følgende for} høydenivå-verdier (markert med

x ˆ

^):

ˆ

x _a = ˆ x _b + 1

2 (y k − x ˆ _b,k )[1, 1, ..., 1 h ] ^T

Dersom

x ˆ b

^{er uniformtfordelt, forenkles denne relasjonen til}

ˆ x _a = 1

2 (y k + ˆ x _b )

Mankannaturligvisogsåjustere vektforholdet

σ b /σ o

^{slik atdennerelasjonenoppnåsselv}

med ikke-uniform fordeling av

(ˆ x _a , x ˆ _b )

^.

Figur 4 viser den resulterede analyse-løsningen med bruk av sistnevnte uttrykk.

Analyse-løsningen ligger omtrent midt mellom eksakt og bakgrunnsløsning, og er derfor

en klarforbedringav bakgrunnsløsningen.

0 0.5 1 1.5 2

0 20 40 60 80 100 120 140 160 180

Solutions, u/U_inf

Prediction time (min)

analysis solution, u-comp.

background solution, u-comp.

exact solution, u-comp.

-1 -0.5 0 0.5 1 1.5

0 20 40 60 80 100 120 140 160 180

Solutions, v/U_inf

Prediction time (min)

analysis solution, v-comp.

background solution, v-comp.

exact solution, v-comp.

Figur4:Sammenligningmellomanalyse-løsning

x _a

^,bakgrunnsløsning

x _b

^{ogeksaktløsning}

x _t

^. Hastighetskomponenter i kontrollpkt. (2.2,0.4,0.4)som funksjon av tid.

Idettetilfelleterobservasjonsdatabarekjentinntilstartenav simuleringen,ogproblemet

ersomtidligerenevnthvordanvikanbenyttedennekunnskapen tilåforbedrevarslingen.

Den klassiske situasjonen er klar nok, vi kan korrigere initialfeltet vha en metodikk som

den vi alleredehar vist i det foregående. Spørsmålet erså hvordan randbetingelsene skal

korrigeres framover itid.

Vitar utgangspunkt irelasjonen (5):

x a = x b + (y k − x b,k )f (z/H ) σ _b ²

σ ² _b + σ _o ² [1, 1, ..., 1 h ] ^T

⁽⁶⁾

Korreksjonsleddet

∆ k = y _k − x _b,k

ergenerelt tidsavhengig,menproblemetnåeratobservasjonsdata

y _k

^{ikkekanværekjent}

førvarslingsperioden!

Forå illustrere problemet har vi gjennomført to forskjelligetyper korreksjoner, begge så

enkle som mulig: Vi korrigerer bare det som har å gjøre med vindretning, og beholder

alt annet fra bakgrunnsløsningen. Korreksjonen holdes kostant gjennom hele varslings-

perioden. I begge tilfeller antas at den reelle varslingsperioden starter fra tidspunkt

t = 20

^{min., hvor nyeste data er kjent. (De første 20 min. kan her betraktes som en}

initialiseringsfase.)

(1) Her benyttes den korreksjonensom kanmåles i starttidspunktet

t = 20

^min.

I vårt tilfelletilsvarer det en dreining

∆θ ≈ 20

^{grader mellom} observasjon og bakgrunnsløsning i høydepunkt 'k'.

(2) Her antas at vi overtidsintervallet

0 ≤ t ≤ 20

^{har observert en}

faseforskyvning itid på20 min., oginnfører en tilsvarendekorreksjon.

Resultatene av disse korreksjonene er vist i gur 5 sammen med bakgrunnsløsnin-

gen og den eksakte løsningen. Det er klart at korreksjonen (1) har vært mislykket i

dettetilfelle,menskorreksjon(2)girengevinst(somventet utfravårforhåndskunnskap).

Fradettekanvi konkludereatkorreksjon av randbetingelser i varslingsintervallet

er hovedproblemet for en eektiv bruk av lokal data-assimilering. Eksempel (1) ovenfor

illustrerer at en konstant korreksjon basert på data ved et enkelt tidspunkt kan være

risikabelt.Eksempel(2) viserderimotatmer relevantinformasjonkanfåsved bruk av et

tidsvinduforutforvarslingsperioden.Detteindikereratenmetodikkbasertpåanalyseav

forutgående trend i observasjonsdata (over et tidsintervall) er interessant i fortsettelsen

av denne problematikken.

Transport-tiden gjennom det lokale området er så kort at vi ikke vil ha nok ydata

i løpet av denne tiden til å kunne benytte 4d-var. Men en mulig tilnærming kan være å

benytte en variasjons-formulering basert på data og bakgrunnsløsningen over siste vars-

lingsperiode, og bruke dette som korreksjon. En slik metode analyseres nærmere i neste

seksjon.

-1 -0.5 0 0.5 1 1.5

0 20 40 60 80 100 120 140 160 180

Solutions, v/U_inf

Prediction time (min)

analysis solution, v-comp.

background solution, v-comp.

exact solution, v-comp.

-1 -0.5 0 0.5 1 1.5

0 20 40 60 80 100 120 140 160 180

Solutions, v/U_inf

Prediction time (min)

analysis solution, v-comp.

background solution, v-comp.

exact solution, v-comp.

Figur 5: Sammenligning mellom eksakt løsning

x _t

^, bakgrunnsløsning

x _b

^{og to analyse-}

løsningeri kontrollpkt. (2.2,0.4,0.4) for en hastighetskomponent. Øverst: Korreksjon ba-

sert på forskjell i vindretninger på tidspunkt

t = 20

^{min. Nederst: Korreksjon basert på}

faseforskyvningi tid.

og tidsvarierende randbetingelser

Ethovedproblemhvaangår lokalværvarslingogdata-assimileringer åbestemme etgodt

estimatav startbetingelse ograndbetingelser for det tidsintervallet vi betrakter, som ty-

pisk er 2-3 timer. For lokal varsling er randbetingelser spesielt viktig fordi endring av

disseiprediksjonsintervallet ofteforplanterseg tilhele områdetisammetidsintervall.Ut-

gangspunkteteratvi kunhar målingerframtilstarttidspunktet fordettetidsintervallet.

Derimot har vi en bakgrunnsløsning (fra en 'grovere' og større modell) for det aktuelle

tidsintervallet. Vår antagelse er at det sannsynligvis vil være fordelaktig å benytte må-

linger (f.eks. fra y) i forkant av dette starttidspunktet for derved å kunne ekstrahere

eventuellekarakteristikkasom bakgrunnsløsningenhariforkantav starttidspunktet.Det-

te kan f.eks. være systematiske feil som faseforskyvning og amplitudedempning. Det vil

davære rimeligå anta atbakgrunnsløsningen fortsetter åha en slik karakter for predik-

sjonsintervalletogså. Visetterossderfor sommålåtilpasseen 'analyse'iettidsintervalli

forkantav prediksjonsintervallet basert påtilgjengeligemålingerogbakgrunnsløsning for

derved å ekstrapolere forskjellsforløpetmellom bakgrunnsløsning oganalyseløsning inn i

prediksjonsintervallet.

6.1 Analyse i forkant av prediksjon

En vanlig måte å bestemme en kurvetilpasning for et gitt datasett er gjennom bruk

av en eller annen form for 'minste kvadraters tilpasning' og antagelse av en polynomisk

kurveformmedgittgrad.Ivårttilfellebestårdatasettetavbådemålingerogverdierfraen

bakgrunnsløsning.Imidlertiderdetvanskeligåutsinoeomkurveformenpåforhånd.Forå

bestemmeetveietmiddelmellombakgrunnsløsningogmåledatabenytterviderforistedet

en variasjonsmetode. Viser til[2℄og referansene her. Den ergittved en 'kostfunksjon'

J (x) = 1

2 (x − x _b ) ^T B ^{− 1} (x − x _b ) + 1

2 (Hx − y) ^T R ^{− 1} (Hx − y)

⁽⁷⁾

medtobegrensningersomrepresentererh.h.v. bakgrunnsløsning ogmåledata. Måletvårt

er å minimere

J (x)

^{. Den}

x

^{som gjør det kaller vi} 'analyseløsningen' og benevner vi

x _a

^.

I ligning (7) er

x _b

bakgrunnsløsningen,

B

^erbakgrunnsløsningens kovariansmatrise,

y

^er

måledatavektoren,

H

^eren 'avbildningsmatrise'som bringer

x

^overimåledataposisjonog

R

^{ermålingenes}kovariansmatrise.Sistnevntematrisekanantasåværeendiagonalmatrise

medkonstant verdi

ρ

^pådiagonalen.Denne verdien sier noeomusikkerheten imålingene eller alternativt hvor mye vekt vi legger på målingene. Den andre kovariansmatrisen,

B

^{, bestemmer hvordan målingene spres til} nærtliggende (her i tid) punkter. Vi antar

at

B

^{er en symmetrisk} båndmatrise med en gitt båndbredde som representeres ved en 'inuensfaktor'in. Verdieneidette båndet larvi først variere lineærtmellom1 og0 fra

diagonalelementettilyttersteelementibåndet,påbeggesideravdiagonalen.Senerevilvi

benytteandreinuensprologså(jfr.Figur14).Viantarat

x _b

^(og

x _a

^)harNkomponenter,

y

^{har M} komponenter og at

M ≪ N

^{. Det betyr at}

B

^og

R

^er

N × N

^{-matriser og}

H

^er

en

M × N

^-matrise.

Gradienten (vektor) til

J (x)

^ergittved

∇ J(x) = B ^{− 1} (x − x _b ) + H ^T R ^{− 1} (Hx − y).

⁽⁸⁾

Forå bestemme analyseløsningen

x = x a

^{setter vi gradienten lik} nullvektoren

B ^{− 1} (x a − x _b ) + H ^T R ^{− 1} (Hx a − y) = 0.

⁽⁹⁾

Multipliserervi ligningen med

B

^{, setter}

R = ρI

^{ogløser m.h.p.}

x _a

^{får vi}

( I + ρ ^{− 1} BH ^T H ) x _a = x _b + ρ ^{− 1} BH ^T y

⁽¹⁰⁾

Ligning 10 er kun en nødvendig betingelse for minimum, men vi ser også at dette er en

tilstrekkeligbetingelsesiden Hess-matrisen, gittved

∇ ² J (x) = B ^{− 1} + ρ ^{− 1} H ^T H,

er positiv denitt (symmetri ok, positive egenverdier ikke innlysende men er sjekket for

mindreeksempler).

Observerat

x _a → x _b

^når

ρ → ∞

^(nullvektpåmåledataene)ogatligningenkonverger mot

Hx a = y

når

ρ → 0

^(allvektpåmålingene).Begge'ekstremaltilfellene'ersomde børvære.Ligning (10) kanskrives som

Ax a = b,

der

A = I + ρ ^{− 1} BH ^T H, b = x _b + ρ ^{− 1} BH ^T y .

Observer at

A

^{har samme båndbredde som}

B

^{hvis alle}målepunktene sammenfallermed komponenter ibakgrunnsvektoren. Hvisminst etmålepunktliggermellomto komponen-

ter i

x _b

^{og vi anvender en lineært} interpolerende avbildning, vil båndbredden øke med 1. For store systemer vil det selvsagt være lurt å benytte denne matrisestrukturen når

ligningssystemet (10) skal løses. For eksemplene i det etterfølgende har vi benyttet en

klassisk Gauss-eliminasjon uten hensyntagen til koesientmatrisens båndstruktur fordi

de resulterendeligningssystemene errelativt små.

6.1.1 Eksempel med sinusformet bakgrunnsløsning og en ikke-systematisk

amplitudefeil

Vi ser først på et eksempel med en sinusformet bakgrunnsløsning med 41 tidsverdier

(

N = 41

^{).Deretterantarviatvihar4målinger(}

M = 4

^)idettetidsintervallet.Målingene erkonstruert slik at de representerer en usystematisk korreksjon til bakgrunnsløsningen,

-6 -5 -4 -3 -2 -1 0 1 2 3 4 5

0 5 10 15 20 25 30 35 40

t

x

xb y

Figur 6:Eksempel med hovedsaklig amplitudefeili bakgrunnsløsning.

sterkerevektpåmålingene(

ρ = 0.1

^{)måvibenytteenstørre}inuensfaktorforåfremskae en 'glatt' analyseløsningenn ellers.

Vi har med hensiktbenyttet en usystematisk amplitudefeilfor å teste metoden påen

realistisksituasjon.Detbetyrselvsagtatdet ikkeermuligellerønskeligåbeholdeformen

påbakgrunnsløsningenianalyseløsningen.Dette gjenspeilesogsåiresultatenesomervist

iFigur 7.Mangel på målingeri den første fjerdedelen av tidsintervallet erhovedårsak til

atvi ikke eri stand tilå nærmeoss formen tilbakgrunnsløsningen her.

6.1.2 Eksempel med sinusformet bakgrunnsløsning og en systematisk fase-

forskyvning

Vi ser nå på et eksempel som har den samme sinusformede bakgrunnsløsningen som i

forrige eksempel, men nå antar vi at de 4 målingene (

M = 4

^{) er} faseforskjøvet

π/8

radianer. Dette er illustrerti Figur8.

Idetteeksempletharvivariert

ρ

^{oginpåsammemåtesomieksempel1.}Målsettingen eråminimereavstandentilmålingenesamtidigsomanalyseløsningenharsammeformsom

bakgrunnsløsningen,spesieltforsistedelavintervallet(sidensisteverdipåanalysekurven

ertenktbenyttetietterfølgendeprediksjon).FraFigur9serviatvalget

ρ = 0.1, inf l = 9

girdet beste resultatet. I dette kurveplottet vises også resultatet for

ρ = 0.1, inf l = 18

som grønn kurve. Vi har også lagt inn den 'eksakte' kurven, vist som orange kurve. Vi

seratvitilfelletmed ekstrastor inuensfaktor(in=18)giren kurvesom liggernærmere

deneksakte kurven enntilfelletin=9idenførstefjerdelenavtidsintervallet.Grunnentil

dettemåvære atmålingi punkt 20begynner åfå virkningogså her når inuensfaktoren

erså stor (18).

I Figure10er vistresultatet med en ekstra målingi startenav tidsintervalletog med

'tunede' verdier for

ρ

^{ogin med det som formålå nneut hvor nærtvi kankommeden}

ρ =1.0, infl=3

-6 -5 -4 -3 -2 -1 0 1 2 3 4 5

0 5 10 15 20 25 30 35 40

t

x

xb xa y

ρ =1.0, infl=6

-6 -5 -4 -3 -2 -1 0 1 2 3 4 5

0 5 10 15 20 25 30 35 40

t

x

xb xa y

ρ =1.0, infl=9

-6 -5 -4 -3 -2 -1 0 1 2 3 4 5

0 5 10 15 20 25 30 35 40

t

x

xb xa y

ρ =0.1, infl=3

-6 -5 -4 -3 -2 -1 0 1 2 3 4 5

0 5 10 15 20 25 30 35 40

t

x

xb xa y

ρ =0.1, infl=6

-6 -5 -4 -3 -2 -1 0 1 2 3 4 5

0 5 10 15 20 25 30 35 40

t

x

xb xa y

ρ =0.1, infl=9

-6 -5 -4 -3 -2 -1 0 1 2 3 4 5

0 5 10 15 20 25 30 35 40

t

x

xb xa y

Figur7:Usystematiskamplitudefeil:Analyse forforskjelligeverdierav

ρ

^oginuensfaktor in.

eksakte løsningen i dette faseforskyvningseksemplet. Som guren viser ligger analyseløs-

ningen'svært nær' den eksakte løsningen (grønn kurve).

6.1.3 Hvor mange målinger trengs? Bakgrunnsløsning med både fase- og

amplitudefeil

Viser igjenpåeteksempelmedsinusformetbakgrunnsøsningienenkeltperiode(

0 − 2π

^).

Denneantarvinåharenkombinertsystematiskfase-ogamplitudefeil.Etutvalgavtestene

-5 -4 -3 -2 -1 0 1 2 3 4 5

0 5 10 15 20 25 30 35 40

t

x

xb y

Figur8: Eksempelmed systematisk faseforskyvning.

vendepunkt idetspesielleeksempletvihartestet.I praksisvildetvære muligåanalysere

en gitt bakgrunnsløsning for etgitt tidsintervall m.h.p. de samme karakteristikka. Dette

kan da være til hjelp med tanke på valg av målinger eller vurdering av tilgjengelige

målingerstilstrekkelighet.

6.1.4 Eksempel med sinusformet bakgrunnsløsning med variabel amplitude

og systematisk fasefeil

Vi ser på et eksempel med sinusformet bakgrunnsøsning i tre perioder (

0 − 6π

^{) med}

variabelamplitude. I Figur12 er vistresultatene med 11 og 13målinger.Siste alternativ

samsvarer med antall endepunkt, ekstremalpunkt og vendepunkt som viser (også her)

at det antallet målinger erakkurat nok til å generere en glatt analyseløsning som ligger

megetnær den eksakte løsningen.

6.1.5 Lineær, kvadratisk og kubisk inuens-prol

Vi ser fra en rekke av plottene knyttet til tidligere tester at vi mangler

C ¹

-kontinuitet i analyseløsningen i måletidspunktene. Det er naturlig å anta at dette er p.g.a. 'hatte-

funksjonsprolet' på inuensfaktoren i kovariansmatrisen

B

^{, m.a.o. har den en knekk i}

midten.Hvis vi istedenf.eks. benytter ethøyereorden polynomisk prol, erdet forventet

atviskalkunne få en glattereanalyseløsning.Med tilpasningav båndbredden fårvi også

detforventederesultatmedbådeetkvadratiskogetkubiskprolsom Figur13viser. Det

er verdt å merke seg at båndbredden for det 'kvadratiske tilfellet' måvære mindre eller

likavstandenmellommålepunktene.Foretkubiskprolmåvibruke en båndbredde som

er 2 ganger avstanden mellom målepunktene. Dette krever en nærmere analyse. De tre

ρ =1.0, infl=3

-5 -4 -3 -2 -1 0 1 2 3 4 5

0 5 10 15 20 25 30 35 40

t

x

xb xa y

ρ =1.0, infl=6

-5 -4 -3 -2 -1 0 1 2 3 4 5

0 5 10 15 20 25 30 35 40

t

x

xb xa y

ρ =1.0, infl=9

-5 -4 -3 -2 -1 0 1 2 3 4 5

0 5 10 15 20 25 30 35 40

t

x

xb xa y

ρ =0.1, infl=3

-5 -4 -3 -2 -1 0 1 2 3 4 5

0 5 10 15 20 25 30 35 40

t

x

xb xa y

ρ =0.1, infl=6

-5 -4 -3 -2 -1 0 1 2 3 4 5

0 5 10 15 20 25 30 35 40

t

x

xb xa y

ρ =0.1, infl=9

-5 -4 -3 -2 -1 0 1 2 3 4 5

0 5 10 15 20 25 30 35 40

t

x

xb xa y

Figur9:Systematiskfaseforskyvning:Analyseforforskjelligeverdierav

ρ

^oginuensfaktor in.

forskjelligeprolene er visti Figur 14oger gittved

β 1 = 1 − | i − j |

inf l

⁽¹¹⁾

β 2 = 1 − | i − j | ²

inf l ²

⁽¹²⁾

β 3 = 1 − | i − j | ³

inf l ³ × (3 inf l − 2abs(i − j))

⁽¹³⁾

ρ=0.02, infl=15

-5 -4 -3 -2 -1 0 1 2 3 4 5

0 5 10 15 20 25 30 35 40

t

x xb

xa y eksakt

Figur 10: Eksempel med tilleggsmålingi startenav tidsintervallet.

iguren er kurven til dieransen

x b (t) − x a (t)

^{en glatt kurve som viseren} karakteristisk 'bølgeaktig'variasjon med en gitt fase og amplitude-endring.Dette vil det være muligå

utnytte i prediksjonsintervallet. En konservativt estimatville være å benytte en middel-

verdi for amplitude ogperiode samtfaseforskyvning, lest utfra dieransen

x b (t) − x a (t)

^.

En 'polynomisk parametertilpasning' er vurdert, men ikke funnet å inneha nødvendige

egneskaperforvårtformål,nemligåfange spesikkekarakteristikkasomvi kanekstrapo-

lere inn iprediksjonsintervallet. I det etterfølgendevil vi isteden gåutifra atdieransen

haren spesikk variasjonder både amplitudevariasjon (

a 0 + a 1 t

^),faseforskyvning(

ω 0

^)og

periode eller bølgelengde (

ω 1

^{) inngåreksplisitt}

dx ab (t) ≡ x b (t) − x a (t) := (a 0 + a 1 t) sin(ω 0 + ω 1 t),

⁽¹⁴⁾

Anvender vi en slik tilpasning, kan vi anvende en 'minste-kvadraters metode' med

E i = x b (t i ) − x a (t i )

a 0 ,a min 1 ,ω 0 ,ω 1

F (a 0 , a 1 , ω 0 , ω 1 ) ≡ 1 2

N

X

i=1

[(a 0 + a 1 t i ) sin(ω 0 + ω 1 t i ) − E i ] ²

! .

Løsningen måtilfredsstillekravene

∂F

∂a 0

= ∂F

∂a 1

= ∂F

∂ω 0

= ∂F

∂ω 1

= 0

ρ =0.01, infl=15, M=4

-1.5 -1 -0.5 0 0.5 1 1.5

0 5 10 15 20 25 30 35 40

t

x

xb xa xt y ρ =0.01, infl=15, M=3

-1.5 -1 -0.5 0 0.5 1 1.5

0 5 10 15 20 25 30 35 40

t

x

xb xa xt y

ρ =0.01, infl=15, M=5

-1.5 -1 -0.5 0 0.5 1 1.5

0 5 10 15 20 25 30 35 40

t

x

xb xa xt y

Figur 11: Eksempelmed forskjellig antallmålinger.

som gir

g 1 (a 0 , a 1 , ω 0 , ω 1 ) ≡

N

X

i=1

[(a 0 + a 1 t i ) sin(ω 0 + ω 1 t i ) − E i ] sin(ω 0 + ω 1 t i ) = 0

⁽¹⁵⁾

g 2 (a 0 , a 1 , ω 0 , ω 1 ) ≡

N

X

i=1

[(a 0 + a 1 t i ) sin(ω 0 + ω 1 t i ) − E i ] t i sin(ω 0 + ω 1 t i ) = 0

⁽¹⁶⁾

g 3 ≡

N

X

i=1

[(a 0 + a 1 t i ) sin(ω 0 + ω 1 t i ) − E i ] (a 0 + a 1 t i ) cos(ω 0 + ω 1 t i ) = 0

⁽¹⁷⁾

g 4 ≡

N

X [(a 0 + a 1 t i ) sin(ω 0 + ω 1 t i ) − E i ] (a 0 + a 1 t i )t i cos(ω 0 + ω 1 t i ) = 0

⁽¹⁸⁾

La

r = 0, 1, 2, . . .

^være iterasjonsindeksen.For dette systemet fårda metoden formen

a ^(r+1) ₀ = a ^(r) ₀ + da ^(r+1) ₀ a ^(r+1) ₁ = a ^(r) ₁ + da ^(r+1) ₁ ω ₀ ^(r+1) = ω ₀ ^(r) + dω ₀ ^(r+1) ω ₁ ^(r+1) = ω ₁ ^(r) + dω ₁ ^(r+1)

der

a ^(r+1) ₀ , a ^(r+1) ₁ , ω ₀ ^(r+1) , ω ₁ ^(r+1)

^{ergittved det lineære}ligningssystemet









 ( ^∂g _∂a ¹

0 ) ^(r) ( ^∂g _∂a ¹

1 ) ^(r) ( _∂ω ^{∂g 1}

0 ) ^(r) ( _∂ω ^{∂g 1}

1 ) ^(r) ( ^∂g _∂a ²

0 ) ^(r) ( ^∂g _∂a ²

1 ) ^(r) ( _∂ω ^{∂g 2}

0 ) ^(r) ( _∂ω ^{∂g 2}

1 ) ^(r) ( ^∂g _∂a ³

0 ) ^(r) ( ^∂g _∂a ³

1 ) ^(r) ( _∂ω ^{∂g 3}

0 ) ^(r) ( _∂ω ^{∂g 3}

1 ) ^(r) ( ^∂g _∂a ⁴

0 ) ^(r) ( ^∂g _∂a ⁴

1 ) ^(r) ( _∂ω ^{∂g 4}

0 ) ^(r) ( _∂ω ^{∂g 4}

1 ) ^(r)





















da ^(r+1) ₀ da ^(r+1) ₁ dω ^(r+1) ₀ dω ^(r+1) ₁











= −









 g ₁ ^(r) g ₂ ^(r) g ₃ ^(r) g ₄ ^(r)











(19)

Koesientene i matrisen ovenfor får vi ved å dierensiere

g 1

^,

g 2

^,

g 3

^og

g 4

^m.h.p.

a 0

^,

a 1

^,

ω 0

^og

ω 1

^. Ligningssystemet(19) blirsymmetrisk og ergittved ligningene (15-18) samt

∂g 1

∂a ₀ =

N

X

i=1

sin ² (ω 0 + ω 1 t i )

⁽²⁰⁾

∂g 2

∂a 1

=

N

X

i=1

t ² _i sin ² (ω 0 + ω 1 t i )

⁽²¹⁾

∂g 3

∂ω 0

=

N

X

i=1

(a 0 + a 1 t i ) [(a 0 + a 1 t i ) cos[2(ω o + ω 1 t i )] + E i sin(ω 0 + ω 1 t i )]

⁽²²⁾

∂g 4

∂ω 1

=

N

X

i=1

t ² _i (a 0 + a 1 t i ) [(a 0 + a 1 t i ) cos[2(ω o + ω 1 t i )] + E i sin(ω 0 + ω 1 t i )]

⁽²³⁾

∂g 1

∂a 1

= ∂g 2

∂a 0

=

N

X

i=1

t i sin ² (ω 0 + ω 1 t i )

⁽²⁴⁾

∂g 1

∂ω 0

= ∂g 3

∂a 0

=

N

X

i=1

(a 0 + a 1 t i ) sin[2(ω o + ω 1 t i )] + E i cos(ω 0 + ω 1 t i )

⁽²⁵⁾

∂g ₂

∂ω 1

= ∂g ₄

∂a 1

=

N

X

i=1

t ² _i (a 0 + a 1 t i ) sin[2(ω o + ω 1 t i )] + E i cos(ω 0 + ω 1 t i )

⁽²⁶⁾

∂g ₃

∂ω 1

= ∂g ₄

∂ω 0

=

N

X

i=1

t i (a ₀ + a ₁ t i ) [(a ₀ + a ₁ t i ) cos[2(ω o + ω ₁ t i )] + E i sin(ω ₀ + ω ₁ t i )]

⁽²⁷⁾

∂g 1

∂ω 1

= ∂g 4

∂a 0

= ∂g 2

∂ω 0

= ∂g 3

∂a 1

=

N

X

i=1

t i (a 0 + a 1 t i ) sin[2(ω o + ω 1 t i )] + E i cos(ω 0 + ω 1 t i )

⁽²⁸⁾

Nå gjenstår kun korreksjon av bakgrunnsløsningenen i prediksjonsintervallet ved eks-

trapolasjon av

dx ab (t)

x a (t) = x b (t) − dx ab (t).

betingelse for hele prediksjonsintervallet i disse punktene. Som tidligere nevnt (og argu-

mentert for) bør disse målepunktene være lokalisert i et viss høyde over terrenget slik

at terrengeekter i størst mulig grad er elliminert fra måleverdiene; m.a.o. ønsker vi å

fange oppdynamiskeogtrendsettende vindforhold(eventuelt også trykk- ogtemperatur-

forhold) i en periode før prediksjonsintervallet som i størst mulig grad er representative

forheleprediksjonsområdet.Enmetodeforåekstrapolerede estimertepunktvariasjonene

iprediksjonsintervallet tilrendene erbeskrevet tidligere,men i hovedsak går metoden ut

påførst åekstrapolere i horisontalplanetder målingenerforetatt og deretter benytte en

'prol-antagelse'for variasjonenned mot bakken.

Konklusjon

Figur 15 viser en bemerkelsesverdig god analyseløsning i prediksjonsintervallet (til

høyreforden vertikalerødelinjeniFigur15).Grunnen tildet ernokatvårantagelseom

'avviksprol' på formen gitt ved ligning (14) er meget god, bedre enn det vi antagelig i

alminneligheteristand tilå anta idet generelle tilfellet.

Disse testene er meget lovende, men det gjenstår å gjøre mer omfattende tester på

andreeksempler.Det vilikkeminstvære viktigåtestemetodenpårealistisketilfellerder

bakgrunnsløsningenhar andresystematiskeavviki forhold tilmålinger.

Siden

x b (t) − x a (t)

^kanhavidtforskjelligformi intervallet der viharmålinger(og der vi gjør parametertilpasning)er det naturligå tenke seg anvendt en endelig Fourier-rekke

istedenforuttrykket i ligning(14). Dette vil bliundersøkt idet videre arbeid.

-5 -4 -3 -2 -1 0 1 2 3 4

0 20 40 60 80 100 120

t

x

xb xa xt y -5

-4 -3 -2 -1 0 1 2 3 4

0 20 40 60 80 100 120

t

x

xb xa xt y

Figur 12: Eksempel med varierende amplitude.

ρ =0.01, infl=10, M=5, lineær

-1.5 -1 -0.5 0 0.5 1 1.5

0 5 10 15 20 25 30 35 40

t

x

xb xa xt y

ρ =0.01, infl=10, M=5, kvadratisk

-1.5 -1 -0.5 0 0.5 1 1.5

0 5 10 15 20 25 30 35 40

t

x

xb xa xt y

ρ =0.01, infl=20, M=5, kubisk

-0.5 0 0.5 1 1.5

0 5 10 15 20 25 30 35 40

x

xb

xa

xt

y

-10 10 0.1

0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

-10 10

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

-10 0 10

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

Figur14: Lineært, kvadratisk og kubiskinuensprol.

-5 -4 -3 -2 -1 0 1 2 3 4

0 12 24 36 48 60 72 84 96 108 120

t

x

xb xa xt y

-2 -1.5

-1 -0.5

0 0.5 1 1.5 2 2.5

0 12 24 36 48

xb xa xt dx1 (a) xa beregnet kun med 5 målinger i intervallet [0,48]

-4 -3 -2 -1 0 1 2 3 4

0 12 24 36 48 60 72 84 96 108 120

t

x

xb xa xt dx1 dx2 (b) dx1=xb-xa i intervallet [0,48]

(c) xa=xb-dx1(ekstrapolert) i intervallet [48,120]

[1℄ Hunt, J.C.R., Snyder, W.H. (1980): Experiments on stably and neutrally stratied

ow overa model three-dimensional hill.J.Fluid Meh. (1980),96, 671704.

[2℄ Lund, B., Eidsvik, K., Sørli, K., Utnes, T., Berntsen, H., ALMOD - Summary of

Leture Series onData AssimilationTehniques, SINTEF Report A4205, 2007.

A Forenklet 3d-Var: Lokal data-assimilering

A.1 Generell kostfunksjon

Kostfunksjonen for analysen kan ifølge standard 3d-Var variasjonsformulering uttrykkes

påformen (Bouttier &Courtier [1℄):

J (x) = (x − x b ) ^T B ^{− 1} (x − x b ) + (y − Hx) ^T R ^{− 1} (y − Hx)

^(A-1)

Her representerer første ledd på høyre side eekten fra bakgrunnsløsningen

x b

^{, og siste}

ledd gir eekten fra observasjonsdata

y

^. Tilstandsvektoren

x

^{inneholder i vårt tilfelle}

variablesom hastighet, temperaturm.m.,

H

^eren transformasjonsmatrise, mens

B

^og

R

erkovariansmatriser for hhv. bakgrunns- ogobservasjonsdata.

Minimumsverdien,minJ(x), nnesved åsette gradienten

∇ J(x) = 0

^for

x = x a

^:

∇ J(x a ) = 0 = 2B ^{− 1} (x a − x b ) − 2H ^T R ^{− 1} (y − Hx a )

^(A-2)

som girfølgende uttrykk:

x a − x b = BH ^T R ^{− 1} (y − Hx a )

^(A-3)

Vektoren

x a

^betegnes 'analysevektoren', dvs. den søkte verdien av tilstandsvektoren som minimaliserer kostfunksjonen, og som gir korrigerte initial- og randverdier. Relasjonen

(A-3)kan også uttrykkes eksplisitt påformen ([1℄)

x a − x b = BH ^T (HBH ^T + R) ^{− 1} (y − Hx b )

^(A-4)

A.2 Eksempel med ett observasjonspunkt

Med tanke pårelativt liten datatilgangeni det lokaleområdet, erdet interessant å sepå

følgendeforenklinger:

•

^{Anta at} tilstandsvektoren

x

^{er gitt i N} gridpunkter,

og at vi har bare én observasjon, spesisert i gridpunkt k.

Dette girfølgende forenkling:

Observasjons-kovariansmatrisen reduseres tilen skalarvarians denert ved

R ≡ σ _o ²

^.

Videre forenkles matriseproduktene

BH ^T = [B 1k , B 2k , ..., B N k ] ^T

dvs.at bare den k-tekolonnen ibakgrunns-kovariansmatrisengjenstår,og

H(BH ^T ) = B kk ≡ σ _b ²

som representerer den skalare bakgrunns-variansen.

Løsningen (A-4)reduseres dermed tilfølgende uttrykk:

x a − x b = y k − x b,k

σ _b ² + σ ² _o [B 1k , B 2k , ..., B N k ] ^T

^(A-5)

hvor

(y k − x b,k )

^{angirden skalaredieransen mellom}observasjonsverdi ogbakgrunnsverdi i gridpunkt k. Vi ser dermed at bakgrunns-kovariansmatrisen er avgjørende for hvordan

informasjonenspres utfra punkt 'k' (hvorobservasjonspunktet benner seg).

A.2.1 Fordeling av høyde-korreksjon

Vigjør nåfølgende tilleggsantakelse ut frafysiske betraktninger:

•

^{Anta at høydedata sprer} informasjon tilnærmet uniformt i høydenivå langt over lokale fjellformasjoner.

Dvs. atbakgrunns-kovariansmatrisenfor høydedata reduseres til

B h,k ≈ σ ² _b [1, 1, ..., 1 h ] ^T

hvor indeks 'h,k' refererertil allegridpunkter isamme høydesom 'k'.For dettehøydeni-

vået blirdermed analysen forenklet til

x a,h = x b,h + (y k − x b,k ) σ _b ²

σ _b ² + σ ² _o [1, 1, ..., 1 h ] ^T

^(A-6)

Med gitthorisontal fordelingav høyde-korreksjonen gjenstår dermedå nne en fornuftig

vertikalfordelingnedmotbakkenivået.Påbakkenerhastighetskorreksjonenliknull,ogvi

måantaenvertikalfordelingfordelingutfrafysiske vurderinger,hvorkorreksjonseekten

avtar mot null i bakkenivå. Formelt antas dermed at fordelingen kan faktoriseres i en

horisontalog vertikaldel, og den totale distribusjonen blirdaav typen

x a = x b + (y k − x b,k )f (z/H ) σ _b ²

σ ² _b + σ _o ² [1, 1, ..., 1 h ] ^T

^(A-7)

hvor

f (z/H )

^{er en antatt vertikal fordeling,}

H

^{er høyden til øverste} modellnivået hvor høydedatamåles,

z

^{erdenaktuellehøyden, ogallehøydermålesfrabakkenivå.Likevekts-}

prolet ved nøytral strømning er tilnærmet logaritmisk over att terreng, og dette kan

derfor være etmuligvalgav vertikalefordeling i sliketilfeller.Mer generelt kandet være

realistisk å velge samme vertikale fordeling som allerede gitt fra bakgrunnsløsningen, og

vi fårdermed

f (z/H ) ∝

( ln(z/z o ),

^{for konstant horisontalposisjon}

x b (z/H ),

^d.s.

(A-8)

Ligning(A-7)representerer korrigerte verdier (

x a

^{) for de variableomkring det tidspunkt}

t o

^hvor observasjonen

y k

^{er gjort. Dette girkorrigerte} initialbetingelserfor modellen ved samme tidspunkt (

t o

^{). Men problemet er nå hvordan} randbetingelsene skal korrigeres framover itid iløpetav varslingsperiodensom starterved

t o

^.

Problematikken erfølgende: I løpet av varslingsperioden som starter ved tidspunktet

t o

^,

oppdateres bakgrunnsverdien

x b

fortløpende i tid med verdier fra en storskala modell (UM1).Itilleggharvikorreksjonenfraobservasjonen

y k

^ved

t = t o

^,slikatanalyseverdien

x a

^{gir oss det beste estimatet både for initial- og} randverdier ved starttidspunktet

t o

^.

Men hvordan skal data-korreksjonen

(y k − x b,k )

^{behandles for}

t > t o

^?

For en eektiv bruk av lokal data-assimilering er dette kanskje både det vanskeligste

og viktigste problemet å nne ut av, spesielt fordi randbetingelsene styrer den lokale

løsningen i mye større grad enn initialbetingelsene. Typisk tar det bare 10 - 20 min. før

en partikkel transporteres gjennom hele det lokale området, og 'hukommelsen' om det

opprinneligeinitialfelteter likekort. Derimot styrer randbetingelsene løsningen gjennom

hele varslingsperioden, som typisk kan være noen få timer.

Referanser

[1℄ BouttierF.,CourtierP.:Dataassimilationoneptsandmethods.MeteorologialTrai-

ning Course LetureSeries, ECMWF, 2002.

B Sensitivitetsanalyse - forenklet 2D-analyse

I dette appendiks rapporterer vi resultater vedrørende sensitiviteten til strømningsfeltet

m.h.t. variasjoner i vindstyrke, på de segmentene av beregningsområdetder luften kom-

mer inn i området. For enkelhets skyld larvi beregningsområdetvære 2-dimensjonaltog

enhetskvadratetmedensirkulærhindringmedgittdiameterogplassering.Detteerillust-

rert i Figur 16. Vindretningen antas enten å være konstant sør-vestlig"(Figure17) eller

kontinuerligvarierendemellomsørogvest,somresultereriatden sørligeogvestligedelen

av randaer"innstrømsegmenter". Spørsmåletvistiller oss er hvilkedeler av strømnings-

områdetsom ermest sensitiv tilendringer i vindstyrken pådisse innstrømsegmentene.

a b

c d

e f

g

h

i

0 1

1

Figur16: Strømningsområde med et sirkulært hinder. Også punktplassering for spesielle

sensitivitetskurver(

a, . . . , i

^{)er vist.}

B.1 Navier-Stokes ligninger for plan strømning

Ikke-stasjonær, inkompressibel strømning av en viskøs væske i et to-dimensjonalt områ-

de

Ω

^er fullstendig beskrevet av væskens hastighet

u

^{med dens horisonale og vertikale}

komponenter

u(x, y)

^,

v(x, y)

^ogtrykk

p(x, y)

^.

Funksjonene

u

^,

v

^and

p

^blirkalttilstandsfunksjoner, og dette settet av tilstandsfunk- sjoner blir kalt en tilstandsløsning. Dette betyr at hastighet og trykk i hvert punkt i

Ω

fullstendigbeskriver tilstanden til det fysiske systemet som betraktes.

Gitt våre antagelser om problemet så antar vi at strømningsfunksjonene

u

^,

v

^og

p

tilfredsstillerNavier-Stokesligningerforikke-stasjonær, inkompressibel,viskøs strømning

ihvert punkt

(x, y)

ⁱ strømningsområdet

Ω

^{. Disse ligningene kanskrivessom:}

∂u

∂t + u ∂u

∂x + v ∂u

∂y − ν ∇ ² u + ∂p

∂x = 0

^(A-9)

∂v

∂t + u ∂v

∂x + v ∂v

∂y − ν ∇ ² v + ∂p

∂y = 0

^(A-10)

∂u

∂x + ∂v

∂y = 0

^(A-11)

med følgenderandbetingelser

u = v = c

√ 2

for

x = 0

^&

y = 0 (

^hastighet

| u | ≡ √

u ² + v ² = c)

^(A-12)

u = v = 0

^for

(x, y) ∈ Γ ₀ (

^heft

)

^(A-13)

∂p

∂x + ∂p

∂y = 0

^for

x = 1

^&

y = 1 (

uforstyrret

)

^(A-14)

∂p

∂n = 0

^for

(x, y) ∈ Γ 0 (

^heft

)

^(A-15)

hvor

c

^,

r

^og

(x 0 , y 0 )

^{er h.h.v.} vindhastighet, radius og senter til strømningshinderet, og randatil strømningshinderetkan dermedbeskrivesved

Γ ₀ = { (x, y) | x = x ₀ + r cos(s), y = y ₀ + r sin(s), 0 ≤ s < 2π } .

ogen kontrollkurverundthovedparten av virvelgaten.

Merk at ligningene (A-9) og(A-10) er symmetriske på den måten at begge ligninger

kantransformeres tilden andreved åbenyttesubstitusjonene

x ↔ y

^og

u ↔ v

^.Detbetyr

at påstander om begge ligninger kan forenkles til påstander om bare den horisontale

ligningen.

B.2 Førsteordens sensitivitetsligninger for

c

Foråutledeførsteordenssensitivitetsligningerantarviatdeteksistererenspesiellløsning

(u ⁰ , v ⁰ , p ⁰ )

^av Navier-Stokesligninger(A-9)-(A-11),med randbetingelsene (A-12)-(A-15).

Vi er interessert i de førstederiverte av tilstandsfunksjonene m.h.p. parameteren

c

^{. For}

åmnne ligningersom sier noeomdissefunksjonene oppførerseg dierensierer viganske

enkeltNavier-Stokesligninger(medrandbetingelser)m.h.p.

c

^ogbytterpådierensierings- rekkefølgen slik at dierensiering m.h.p. denne parameteren blirutført først. Vi innfører

nåkortnotasjonen:

u c ≡ ∂u

∂c v c ≡ ∂v

∂c p c ≡ ∂p

∂c

^(A-16)

De dierensierterandbetingelsenehar sammeformsom de opprinneligerandbetingel-

sene,menmedhomogenehøresider(unntaketerførsteordenssensitivietforkomponentene

u

^og

v

^{).Slikvildet være, uansetthvormangegangervigjentar} dierensieringsprosessen.

Brukervi kortnotasjonen ovenfor blirde resulterendedierensierte Navier-Stokes lig-

∂u c

∂t + u ∂u c

∂x + v ∂u c

∂y − ν ∇ ² u c + ∂p c

∂x = − u x u c − u y v c

^(A-17)

∂v c

∂t + u ∂v c

∂x + v ∂v c

∂y − ν ∇ ² v c + ∂p c

∂y = − v x u c − v y v c

^(A-18)

∂u c

∂x + ∂v c

∂y = 0

^(A-19)

med følgenderandbetingelser

u c = v c = 1

√ 2

for

x = 0

^&

y = 0

^(A-20)

u c = v c = 0

^for

(x, y) ∈ Γ 0

^(A-21)

∂p c

∂x + ∂p c

∂y = 0

^for

x = 1

^&

y = 1

^(A-22)

∂p c

∂n = 0

^for

(x, y) ∈ Γ 0

^(A-23)

B.2.1 Iterasjon

En viktig ting å merke seg er at venstresiden av ligningene (A-17) og (A-18) er identisk

i form med strømningsligningene. Derfor er det naturlig å formulere følgende iterative

prosedyre

∂u ^(k) c

∂t + u ∂u ^(k) c

∂x + v ∂u ^(k) c

∂y − ν ∇ ² u ^(k) _c + ∂p ^(k) c

∂x = − u x u ^(k _c ^{− 1)} − u y v ^(k _c ^{− 1)}

^(A-24)

∂v c ^(k)

∂t + u ∂v c ^(k)

∂x + v ∂v c ^(k)

∂y − ν ∇ ² v _c ^(k) + ∂p ^(k) c

∂y = − v x u ^(k _c ^{− 1)} − v y v _c ^{(k − 1)}

^(A-25)

∂u ^(k) c

∂x + ∂v ^(k) c

∂y = 0

^(A-26)

sidenvidakangjenbrukedensammekoesientmatrisentildediskretiserte systemetsom

forstrømningsligningene.Eksperimenterviseratkunetfåtalliterasjonerernødvendigfor

åoppnå "akseptabel konvergens".

0.5786 0.4962 0.4137 0.3312 0.2487 0.1663 0.0838 0.0013 -0.0811 -0.1636 -0.2461 -0.3285 -0.4110 -0.4935 -0.5760 -0.6585 -0.7409 -0.8234 -0.9059 -0.9884

Figur 19: Hastighet og trykk ved

t = 2.0

^. Beregningene er utført med FreeFem+

(http://www.freefem.org)

0.2 0.4 0.6 0.8 1 1.2 1.4

time

s e n s it iv it y

a (0.5,0.5) b (0.2,0.8) c (0.8,0.2) d (0.5,0.8) e (0.8,0.5) f (0.2,0.5) g (0.5,0.2) h (0.8,0.8) i (0.2,0.2)

a b

c d

e f

g h

i 0 1

1

0.2

0.1 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

0

Figur20:Sensitivitet

u c (x P , y P )

^{som funksjonavtidi}forskjelligepunkter

P = a, b, . . . , i

^.