Het model van Vallis voor El Niño

faculteit Wiskunde en Natuurwetenschappen

Het model van Vallis voor El Niño

Bacheloronderzoek Wiskunde

Juni 2014

Student: C. van der Lei

Eerste Begeleider: Dr. A.E. Sterk Tweede Begeleider: Dr. Ir. F.W. Wubs

Samenvatting

El Ni˜no is een natuurverschijnsel dat rond de evenaar eens in de 3 tot 7 jaar optreedt.

De ogenschijnlijke stabiele situatie, genaamd zuidelijke oscillatie, is dan omgeslagen.

Dit is waarneembaar door de sterke opwarming van het normaal koele zeewater voor de kust van Zuid-Amerika. Wetenschappers zijn al lange tijd op zoek naar een goed model om El Ni˜no’s te kunnen voorspellen. Inmiddels is duidelijk dat het een sterk oceaan-atmosfeer gekoppeld systeem is. Vallis’ model is een simpel, conceptueel model waarmee enkel interessante eigenschappen van El Ni˜no onder de loep kunnen worden genomen, waaronder chaos. We beginnen met een basis model en introduceren vervol- gens een aandrijfwind en een seizoenscyclus om het model realistischer te maken.

Inhoudsopgave

1 Inleiding 3

1.1 ENSO . . . 3

1.2 Gekoppeld systeem . . . 4

1.2.1 Indices . . . 4

1.2.2 Samenhang . . . 6

2 Een El Ni˜no 7 2.1 Aspecten die een rol spelen . . . 7

2.1.1 Walker circulatie . . . 7

2.1.2 Opwelling . . . 7

2.1.3 Thermocline . . . 8

2.1.4 Oost-west asymmetrie . . . 8

2.2 Begin van een El Ni˜no . . . 9

2.3 Herstel van El Ni˜no . . . 9

3 Modellen 10 3.1 Gekoppeld systeem . . . 10

3.2 Modellen . . . 11

4 Het model van Vallis 12 4.1 Analyse u^∗ = 0 . . . 13

4.2 Analyse u^∗ = u₀ . . . 16

4.2.1 Lyapunov exponenten . . . 18

4.3 Analyse u^∗ = u₀(1 + u_ssin ωt) . . . 20

5 Conclusie en discussie 22 APPENDICES 25 A Routh-Hurwitz criterium 26 B Mathematica notebook codes 28 B.1 ut-plot, uy-plot en 3D plot . . . 28

B.2 Bifurcatie plot . . . 29

C Matlab codes 30 C.1 Vallis system . . . 30

C.2 Lyapunov exponenten . . . 31

C.3 Run script . . . 33

D Data NINO3.4 en SOI 34

D.1 NINO 3.4 . . . 34 D.2 SOI . . . 36

Hoofdstuk 1 Inleiding

1.1 ENSO

El Ni˜no ofwel het ‘kerstjongetje’ is de Peruaanse naam voor een natuurverschijnsel dat al eeuwen lang zo nu en dan optreed aan de westkust van Zuid-Amerika. De Peruaanse vissers constateren dit fenomeen om de 3 tot 7 jaar rond kerst; een aantal maanden waarin de visvangst beduidend minder is. Deze kleinere vangst is een gevolg van El Ni˜no. Langs de evenaar in de Grote Oceaan is dan het normaal koele oceaanwater sterk opgewarmd. Hierdoor verandert de luchtdruk en zal de westpassaat zwakker worden. Dit heeft tot gevolg dat het zeewater minder snel wegstroomt en dus nog verder opwarmt. De passaatwind wordt daardoor nog zwakker en het regengebied wat normaal boven Indonesi¨e ligt, is nu verschoven naar het midden van de Grote Oceaan.

Uitlopers hiervan bereiken ook de kust van Zuid-Amerika, waar het veel schade aan kan richten.

In plaats van El Ni˜no kunnen we beter spreken van ENSO, El Ni˜no Southern Oscillation. Dit is een samenstelling van El Ni˜no, La Ni˜na en de zuidelijke oscillatie.

Deze drie verschijnselen hebben samen een sterk onderling verband. In de Grote Oceaan is altijd een van de drie situaties aanwezig. Om te weten welk verschijnsel zich voordoet, wordt er gekeken naar bijvoorbeeld de zogenaamde NINO3.4-index. Deze index wordt bepaald aan de hand van de afwijking van de zeewatertemperatuur langs de evenaar.

Bij zuidelijke oscillatie komt koud zeewater omhoog bij de kust van Zuid-Amerika.

De oostpassaat brengt dit naar het westen. Wanneer het Azi¨e bereikt heeft, is het zo’n vijf graden opgewarmd. Boven dit warme zeewater, stijgt de lucht en regent het vaker dan bij het koude zeewater van Zuid-Amerika. Door stijgende lucht, wordt de luchtdruk lager, wat op zijn beurt weer de oostelijke passaat versterkt. Het lijkt erop dat dit patroon zichzelf in stand houdt. Maar nu is de grote vraag: Waarom gaat een ogenschijnlijke stabiele situatie toch steeds over in een andere situatie?

Het ’zusje’ van het El Ni˜no fenomeen is de La Ni˜na, ofwel het ’meisje’. Dit ver- schijnsel komt vaak voor na een El Ni˜no. Het warme water voor de Zuid-Amerikaanse kust wordt weggedreven en maakt plaats voor een koude stroom. Hierdoor wordt de oostpassaat versterkt en valt er veel regen in Australi¨e en Azi¨e [8, 18].

Figuur 1.1: Een El Ni˜no: Relatief warm water voor de kust van Peru (rood) [3].

1.2 Gekoppeld systeem

1.2.1 Indices

Zoals ook al in de introductie genoemd, is er tijdens een El Ni˜no sprake van warmer oppervlakte water rond de evenaar aan de oostkant van de Grote Oceaan, dus bij de kunst van Peru. Er zijn een aantal verschillende indices die gebasseerd zijn op deze observaties. Ze worden vaak gedefineerd door de gemiddelde zee temperatuur aan de oppervlakte (sea surface temperature, SST) te meten over een bepaald deel van oostelijke Grote Oceaan over een bepaalde periode.

Voor nu gaan we gebruik maken van de index NINO3.4. Dit is de gemiddelde SST in het gebied begrensd door 5^◦N tot 5^◦Z en van 170^◦W tot 120^◦W. Naast de NINO3.4 zullen we gaan kijken naar de zogenaamde SOI (Southern Oscillation Index). Deze index is gedefinieerd als de maandelijkse afwijking van de luchtdruk bij Tahiti gedeeld door zijn standaard deviatie minus de maandelijkse afwijking van de luchtdruk bij Darwin gedeeld door zijn standaard deviatie. Zoals in figuur 1.2 zijn er verschillende indices: NINO1+2, NINO3, NINO4 en dus NINO3.4. Al deze indices zijn begrensd door verschillende lengte en breedte graden.

Wat zeggen deze indices nu eigenlijk? We beginnen weer bij de index NINO3.4, zoals hierboven staat, zegt deze index wat over de temperatuur van het oppervlakte water. Je zou zeggen als deze waarde een bepaald afwijking heeft van de gemiddelde SST, dan is er sprake van een El Ni˜no of La Ni˜na.

De tweede index heeft te maken met drukverschillen. Zoals bekend stroomt lucht van hoge naar lage druk en zullen er dus pasaatwinden ontstaan bij drukverschillen. De drukverschillen in de normale situatie (zuidelijke oscillatie) zorgen voor passaatwinden die waaien van oost naar west [11].

Figuur 1.2: Darwin en Tahiti aangegeven op de kaart en de verschillende NINO indices [28].

Figuur 1.3: NINO3.4 en SOI, periode 1900-2000. Hierbij zijn voor zowel SOI als NINO3.4 de standaardafwijkingen genomen en vergeleken [18].

1.2.2 Samenhang

Nu we bekend zijn met deze indices kunnen we gaan kijken naar de correlatie tussen de twee verschijnselen: Zuidelijke oscillatie en El Ni˜no. De zuidelijke oscillatie wordt namelijk gemeten met de SOI index en een El Ni˜no met de NINO3.4 index. Kijkend naar de data vermoed je al snel een verband. Over het algemeen kan men zeggen dat een hoge NINO3.4 index samenhangt met een lage SOI index en vice versa, zie figuur 1.3. Ook cijfers wijzen dit uit. Met de data van 1951 tot en met 2013 zijn de volgende correlatieco¨effici¨enten bepaald [12].

Maand Jan Feb Mrt Apr Mei Jun

Correlatieco¨effici¨ent -0,80 -0,75 -0,71 -0,65 -0,56 -0,60

Jul Aug Sep Okt Nov Dec

-0,66 -0,82 -0,80 -0,76 -0,76 -0,74 Tabel 1.1: Correlatieco¨effci¨enten per maand.

Dit duidt op een samenhang tussen de twee verschijnselen en het is dus een oceaan- atmosfeer gekoppeld systeem.

Hoofdstuk 2 Een El Ni˜ no

2.1 Aspecten die een rol spelen

2.1.1 Walker circulatie

De situaties die hierboven beschreven staan volgen allemaal de zogenaamde Walker circulatie. Deze circulatie ziet er als volgt uit: de passaatwind brengt de vochtige lucht vanuit het oosten naar het westen. Boven Port Darwin stijgt deze op, waardoor er condensatie optreedt en daarmee gepaard wolkenvorming. Vervolgens stroomt de lucht op grote hoogte weer terug naar het oosten en daalt weer boven Tahiti, waardoor daar een hogedrukgebied ontstaat.

De walker circulatie tijdens de normale situatie/zuidelijke oscillatie ziet er als volgt uit:

1. Passaatwinden waaien vanuit het oosten naar het westen over de Grote Oceaan.

2. De passaten stuwen het oppervlakte water op, waardoor het water voor Azi¨e ongeveer 60 centimeter hoger staat dan voor de kust van Zuid-Amerika.

3. De passaten nemen waterdamp op als ze over de oceaan waaien en raakt dit vervolgens weer kwijt als moessons boven Indonesi¨e.

4. Het oppervlakte water beweegt mee naar westelijke richting. Bij de evenaar wordt deze afgebogen naar de Zuidpool als gevolg van de aardse rotatie.

5. Deze afbuiging zorgt voor een opwelling van kouder water in het oostelijke ge- deelte van de Grote Oceaan.

6. En dan nu het punt waar de Peruaanse vissers bij gebaat zijn: Als het koelere water dichterbij het zonlicht komt, voedt het plankton zich met voedingsstoffen.

Dit plankton zorgt voor een grote hoeveelheid zeeleven voor de kusten van Peru en Chili [1, 14].

2.1.2 Opwelling

Oceaanstromen worden aangedreven door de wind en de kracht van de rotatie van de aarde, genaamd de Corioliskracht. Op het noordelijk halfrond werkt deze kracht met de klok mee en op het zuidelijk halfrond tegen de klok in. Deze opwelling wordt gedre- ven door sterke wind en deze Corioliskracht. Opwelling bij de kust wordt veroorzaakt door het feit dat de passaatwinden het water van de kust wegduwen. Om deze ‘lege’

plaats op te vullen welt er water op naar de oppervlakte. Dit koude water is vaak

Figuur 2.1: Kust opwelling bij Zuid-Amerika. Relatief koud water komt omhoog wellen [37].

rijk aan voedingstoffen, waardoor er de oceaanop die plaatsen rijk is aan leven. Dit fenomeen, de zogenaamde kust opwelling, zie je bij de kust rond de evenaar, bijvoor- beeld de kust van Peru. In figuur 2.1 zie je deze kust opwelling. Een andere vorm van opwelling noemen we evenaar opwelling. Dit gaat als volgt: Bij de evenaar waaien de passaatwinden in westelijke richting door de rotatie van de aarden van oost naar west. Verder zorgt de Corioliskracht ervoor dat deze winden op het noordelijk half- rond worden afgebogen naar rechts en op het zuidelijk halfrond naar links. Waardoor het water dus van de evenaar vandaan wordt gestuwd. Ook hier zal koud, vruchtbaar water opwellen [37].

2.1.3 Thermocline

Een ander begrip wat we hier tegenkomen is de zogenaamde thermocline. De oceaan is ruim 4, 2 kilometer diep. Echter is voor ons slechts 5% van de bovenste laag interessant.

Hier bevindt zich namelijk de 20^◦C thermocline. Een thermocline is een overgangslaag in het water, waarbinnen de temperatuur in verticale richting over een korte afstand sterk verandert. Het temperatuurverschil tussen boven en onder de thermocline kan ongeveer geschat worden op 10^◦C.

Deze thermocline ligt in het westen dieper dan in het oosten, omdat door de pas- saatwind een lichte stroming in de oceaan is ontstaan. In het westen, voor de kust van Australi¨e en Indonesi¨e, daalt het warme water, terwijl in het oosten het koude water uit het diepe opwelt. Naar schatting ligt deze themocline in het westen ruim 200 meter onder de waterspiegel en in het oosten maar rond de 50 meter [15, 37].

2.1.4 Oost-west asymmetrie

Door de opwelling bij de kust, zoals hierboven omschreven ontstaat er een oost-west asymmetrie in de oppervlakte water temperatuur. Door het relatief warme water in het westen, stijgt daar de lucht en ontstaat er een lagedrukgebied. Aan de andere

kant, boven het koude water ontstaat een hogedrukgebied. Doordat de wind waait van hoge naar lage druk, ontstaat er een versterkte wind van oost naar west [14, 37].

2.2 Begin van een El Ni˜ no

In de normale situatie is de kracht van de stuwing in evenwicht met de gravitatie kracht van het water, zoals te zien in figuur 2.2. Maar eens in de zoveel jaar zorgt een turbulentie ervoor dat deze situatie wordt verstoord. De normale westelijk winden moeten plaats maken voor hier en daar oostelijke winden die ervoor zorgen dat er Kelvin golven ontstaan rond de evenaar. Dit zijn grote golven die zich vanuit het westen richting Peru bewegen. Hierdoor wordt het evenwicht verstoord. De oost-west asymmetrie zal verdwijnen en de thermoline zal meer horizontaal komen te liggen.

Doordat de opwelling van relatief koud voedzaam water bij de kust van Peru zal stoppen en de visvangst van de Peruaanse vissers daalt [1, 15, 37].

Figuur 2.2: Oost-west asymmetrie in de oceaan tijdens de normale situatie en tijdens een El Ni˜no [37].

2.3 Herstel van El Ni˜ no

Wanneer de thermocline bijna horizontaal ligt, als in het geval van El Ni˜no, is dit een stabiele, maar ongebruikelijke situatie. Hier begint dan ook het herstel naar de normale toestand. Water aan de oost kant ten hoogte van de evenaar splitst langzaam weer naar noord en zuid. De hoogte van het wateroppervlakte daalt daardoor weer en er ontstaat weer opwelling en de thermocline zal in het oosten weer iets omhoog komen. Op deze manier komt de oost-west asymmetrie weer een stukje terug. Door de het koeler wordende water in het oosten, zullen de passaatwinden weer sterker worden.

Op deze manier versterken de processen elkaar en zal de situatie weer terugkeren in normale, onstabiele toestand. Wel is het zo dat na een sterke El Ni˜no een La Ni˜na volgt. Het terugkeren naar normale toestand ‘slaat door’. Het water bij de kust van Peru wordt kouder dan normaal, de therocline ligt schever dan normaal en de passaatwinden zijn daarmee sterker dan normaal [37].

Hoofdstuk 3 Modellen

3.1 Gekoppeld systeem

In figuur 3.1 zien we een overzicht van het hele gekoppelde systeem. In dit figuur is de standaard situatie geschetst met een koude tong in het oosten en een warm bad in het westen. Verder zie je passaatwinden die deel uit maken van grotere circulaties: de Walker circulatie en de Hadley circulatie. Daarnaast zie je Kelvin golven en Rossby golven, hier een korte toelichten van deze onderdelen.

Kelvin golven zijn golven met een hele grote golflengte die zich bewegen van west naar oost. Deze zijn waarneemnbaar tijdens een El Ni˜no. Rossby golven zijn juist een stuk kleiner en bewegen van oost naar west.

De Walker circulatie wordt uitgebreid genoemd in hoofdstuk 2. De Hadley cir- culatie is de circulatie die ontstaat door de draaiing van de aarde, ofwel de Coriolis schijnkracht. Dit verklaard dat stromen op het noordelijk halfrond afbuigen naar rechts en op het zuidelijk halfrond naar links.

Figuur 3.1: Overzicht gekoppelde systeem [16].

3.2 Modellen

Voor dit systeem zijn vele verschillende modellen gemaakt. In het artikel van McCreary [31] wordt onderscheid gemaakt tussen drie types modellen.

1. Conceptuele en simpele modellen 2. Modellen van gematigde complexiteit

3. Gekoppelde algemene circulatiemodellen (ofwel coupled general circulation mo- dels, CGCM)

Conceptuele modellen zijn sterk vereenvoudigde modellen. Hierbij is er een be- perkt aantal variabelen en zijn de gekoppelde vergelijkingen alleen tijdsafhankelijk.

De simpele modellen bestaan uit vereenvoudigde differentiaalvergelijkingen. Beide bo- vengenoemde types zijn nuttig om belangrijke processen te bekijken, zoals koppeling tussen oceaan en atmosfeer. Echter is het moeilijk deze modellen te vergelijken met de werkelijkheid. Deze modellen zijn dus moeilijk toetsbaar. Artikelen over deze mo- dellen zijn onder andere die van McCreary [29], Bjerkness [10], Graham en White [21], McCreary en Anderson [30], Suarez and Schopf [36], Vallis [42], Hirst [23, 22], Lau [26], Yamagata [43] en Philander et al.[32].

De tussenliggende modellen komen al een stuk dichter bij de werkelijkheid. De modellen bestaan uit een set van gekoppelde systemen die vele realistische verschijnse- len bevatten als thermodynamica. Modellen van dit type zijn goed toepasbaar en nog relatief eenvoudig te analyseren. Wetenschappers die hier over schreven waren onder meer Anderson en McCreary [7, 6], Hirst [22], Schopf en Suarez [34] en het bekende duo Zebiak en Cane [44].

Ten slotte komen we bij de CGCM’s. Dit zijn momenteel de beste, maar daarmee ook de lastigste voorspellingsmodellen voor El Ni˜no. Er is een realistische koppeling tussen de relatief langzame circulaties in de oceaan en de relatief snelle circulaties in de atmosfeer. Voorbeelden van CGCM’s zijn te vinden in de werken van Latif et al.

[25] en AchutaRao en Sperber [4, 5]

Wij zullen echter aan de slag gaan met een simpel model, waarbij we een aantal componenten even links laten liggen, maar toch een aantal belangrijke processen onder de loep kunnen nemen [39, 31, 38].

Hoofdstuk 4

Het model van Vallis

Geoffrey K. Vallis publiceerde in 1986 een simpel, maar realistisch model. Er be- staan vele andere modellen voor El Ni˜no, varie¨erend van zeer complexe, nauwkeurig uitgewerkte oceaan-atmosfeer gekoppelde modellen tot relatief simpele stochastische dynamische of probabilistische modellen. Om een El Ni˜no te voorspellen is een gekop- peld circulatie model misschien nodig, maar ook een simpelere theorie is van belang.

Het model wat we hier gaan behandelen verschilt van de meeste andere omdat er geen stochastische en tijdsgebonden krachten de aperiodiciteit veroorzaken. Hier is geen expliciete golfdynamica nodig om de tijdsschalen te verklaren en geen invloeden van gematigde breedte zijn nodig om de varieteit in intensiteit duidelijk te maken.

We beschouwen de Grote Oceaan als een box met water. De watertemperatuur aan de westkant noemen we T_wen aan de oostkant T_oen de stroming u. Deze stroom wordt veroorzaakt door de passaatwinden, die worden ontstaan door een het temperatuur verschil (T_o − T_w)/∆x. Dit maakt deel uit van de Walker circulatie zoals hierboven beschreven. Een koudere temperatuur in het oosten zorgt dus voor een westelijke passaatwind. We kunnen dit als volgt formuleren:

du

dt = B(T_o− T_w)/2∆x − C(u − u^∗). (4.1) Hierin zijn B en C constanten. De term B(T_o − T_w)/2∆x + Cu^∗ representeert de door de wind geproduceerde spanning. De term −Cu staat voor de demping. Als we aannemen dat het diepe oceaanwater een constante temperatuur heeft van ¯T . Dan is de eenvoudigste weergave van de aanpassing van temperatuur als volgt:

dT_w dt = u

2∆x( ¯T − T_o) − A(T_w − T^∗), (4.2) dTo

dt = u

2∆x(T_w− ¯T ) − A(T_o− T^∗). (4.3) Bij deze twee vergelijkingen staat de eerste term voor het horizontale energietrans- port en opwelling en de tweede term voor de thermische demping, waarbij A en T^∗ constanten zijn. T^∗ is de temperatuur waarbij er evenwicht zou zijn en er dus geen stromingen meer zouden zijn in het water. Wel zorgt deze parameter dan nog voor warmte-uitwisselingen met de atmosfeer.

De vergelijkingen 4.1-4.3 vormen nu samen een basis model. Bij dit model zijn vele dingen verwaarloosd. Zo is er bijvoorbeeld geen aandacht besteed aan onder andere

Kelvin golven. Toch gaan we kijken hoe dit model zicht gedraagt. Zonder verlies van algemeenheid kunnen we ¯T = 0 stellen. We meten nu alle temperaturen ten opzichte van deze ¯T .

T_w Temperatuur aan de west kant van de Grote Oceaan To Temperatuur aan de oost kant van de Grote Oceaan T^∗ Evenwichtstemperatuur van het oceaanwater

T¯ Temperatuur van het diepe oceaanwater (constant verondersteld) u Stroming in het oceaanwater door passaatwinden

∆x Afstand tussen de oost en westkust A, B, C Constanten (nader vast te stellen)

Tabel 4.1: Overzicht van alle variabelen

4.1 Analyse u

= 0

We stellen nu ook u^∗ = 0, op deze manier kunnen we inzicht verkrijgen in het gedrag van dit systeem. Voordat we iets zinnigs kunnen zeggen over dit systeem zal er eerst wat omschrijf- en substitutiewerk aan vooraf gaan. We beginnen met de volgende substituties:

ˆt = At, u =ˆ U

2A∆x, Tˆ_w = T_w

T^∗, Tˆ_o = T_o T^∗.

Verder introduceren we de constanten ˆB en ˆC zo dat ˆB = T^∗B/(∆xA²) en ˆC = C/A.

Hiermee krijgen we het volgende systeem:

dˆu dˆt =

Bˆ

2( ˆTo− ˆTw) − ˆC ˆu, d ˆT_w

dˆt = −ˆu ˆT_o− ( ˆT_w− 1), d ˆT_o

dˆt = −ˆu ˆT_w− ( ˆT_o− 1).

Wanneer we ook nog de volgende substitutie uitvoeren:

y = (To− Tw)/2, z = (T_o+ T_w)/2, krijgen we de volgende vergelijkingen:

dˆu

dˆt = ˆBy − ˆC ˆu, (4.4)

dy

dˆt = ˆuz − y, (4.5)

dz

dˆt = −ˆuy − (z − 1). (4.6)

Dit systeem lijkt enigzins op de Lorenz vergelijkingen:

du

dt = −σu + σy, du

dt = ru − uz − y, dz

dt = uy − bz.

Hierbij is u gesubstitueerd voor de normale x. Wanneer we Vallis systeem vergelijken met die van Lorenz, zien we dat r = 0, b = 1 en z⁰ = z − 1 de overeenkomsten duidelijker maken.

Eerst zullen we opzoek gaan naar de evenwichtspunten van het systeem 4.4-4.6 door steeds het rechterlid nul te stellen. Op deze manier vinden we 3 evenwichtspunten voor B 6= 0 en ˆˆ C 6= 0, namelijk de triviale (0, 0, 1) (namelijk z⁰ = z − 1) en

ˆ u = ±

v u u t

Bˆ

Cˆ 1 − Cˆ Bˆ

! ,

y = ± v u u t

Cˆ Bˆ 1 −

Cˆ Bˆ

! ,

z = Cˆ Bˆ.

Voor ˆB = 0 ´of ˆC = 0 krijgen we alleen de triviale oplossing (0, 0, 1). Wanneer zowel B als ˆˆ C gelijk is aan nul, krijgen we oneindig veel evenwichtspunten. In dit geval is er geen koppeling tussen ˆu, y en z en is dus fysisch niet relevant.

Het is makkelijk te zien dat voor ˆC ≥ ˆB alleen de triviale oplossing mogelijk is.

Verder staat ˆB voor de mate van feedback tussen de temperatuur van het oceaanwater en de oppervlakte winden. Wanneer deze dus 0 is of bijna 0 is er een rust situatie mogelijk.

Zoals we hierboven al kort hebben genoemd is alleen de triviale oplossing mogelijk als ˆC ≥ ˆB. Voor ˆC < ˆB vinden we alle 3 oplossingen. Met andere woorden er vindt een zogenaamde pitchfork bifurcatie plaats bij ˆC = ˆB. Door het systeem te lineariseren kunnen we ook de aard van deze oplossingen gaan bekijken. De linearisatie rond (0, 0, 1) gaat als volgt:



 ˆ u⁰ y⁰ z⁰



=





− ˆC Bˆ 0 z −1 uˆ

−y −ˆu −1







 ˆ u y z



.

Voor de triviale oplossing p₀ ≡ (0, 0, 1) vinden we de eigenwaarden λ₀ = −1, λ_1,2 = 1

2(−1 − ˆC ± q

( ˆC + 1)²− 4( ˆC − ˆB).

Hier zien we meteen dat λ₀ altijd ree¨el en negatief is. Over λ1,2 weten we dat het altijd re¨eel is, omdat ( ˆC + 1)²+ 4( ˆC − ˆB) ≥ 0. Wel zien we een bifurcatie bij ˆB = ˆC,

dan klapt namelijk het teken van λ2 = ¹₂(−1 − ˆC + q

( ˆC + 1)²− 4( ˆC − ˆB) om. Voor B > ˆˆ C wordt deze positief en daarmee onstabiel. Kortom voor ˆB < ˆC hebben we 3 oplossingen, waarvan de triviale onstabiel is en de andere 2 stabiel en voor ˆB > ˆC hebben we alleen de stabiele triviale oplossing. Deze overgang noemen we een pitchfork bifurcatie.

Naast een pitchfork bifurcatie is er ook een Hopf bifurcatie waar te nemen bij



± qBˆ

Cˆ(1 −^Cˆ_ˆ

B), ± qCˆ

Bˆ(1 −^Cˆ_ˆ

B),^Cˆ_ˆ

B



. Hiervoor bekijken we het karakteristieke poly- noom van deze evenwichtspunten:

λ³+ ( ˆC + 2)λ²+ ( ˆC + ˆB/ ˆC)λ + 2( ˆB − ˆC) = 0.

Het kritieke punt is wanneer ( ˆC − 2)( ˆC + ˆB/ ˆC) = ˆB − ˆC (Ofwel de co¨effici¨ent van λ² maal de co¨effici¨ent van λ gelijk aan de constante co¨effi¨ent) Hierbij is gebruik gemaakt van het Routh-Hurwitz criterium, voor nadere toelichting zie bijlage A. Dit kritieke punt noemen we ˆB_h en is gelijk aan

Bˆ_h = Cˆ²(4 + ˆC) C − 2ˆ .

Wanneer ˆB < ˆB_h dan zijn alle eigenwaarden λ complex, maar het re¨ele deel is negatief.

Hiermee is het systeem dus stabiel. Als ˆB > ˆB_h dan is het re¨ele deel van de eigen- waarden positief en daarmee is het systeem instabiel. Kortom bij ˆB = ˆBh ondergaat



± qBˆ

Cˆ(1 −^Cˆ_ˆ

B), ± qCˆ

Bˆ(1 −^Cˆ_ˆ

B),^Cˆ_ˆ

B



een Hopf bifurcatie. Om dit inzichtelijk te maken moeten we eerst schattingen maken van onze parameters. We gebruiken de volgende waarden:

∆x = 7500 km, T^∗ = 12 K, ˆB = 127, ˆC = 3, ˆB_h = 63.

Met de bovenstaande gegevens is figuur 4.1 opgesteld. Merk op dat ˆB > ˆB_h en dus in het onstabiele gebied ligt [27, 35, 41, 42].

Figuur 4.1: Bifurcatie diagram, u^∗ = 0 [35].

4.2 Analyse u

= u

In bovenstaande model hebben we u^∗ = 0 aangenomen om de analyse makkelijker te maken. Echter hebben we daarmee gezegd dat er complete symmetrie is tussen Tw en T_o. Dit is in werkelijk anders: er wordt steeds over gegaan van de situatie die hoort bij p₊ naar p− (El Ni˜no). Om deze complete symmetrie eruit te halen, moeten we een u^∗ introduceren. Echter zou dan de asymmetrie die ontstaat niet weer kunnen herstellen. Vandaar dat we u − u^∗ in het model stoppen. Verder correspondeert een westelijke passaatwind met een u^∗ < 0. Door deze verandering in het model wordt een analytische analyse erg complex, maar met een numerieke analyse kunnen we toch interessante dingen vinden over de bifurcaties. Om nu weer de evenwichtspunten te vinden moeten we de volgende vergelijkingen oplossen:

By − ˆˆ C(ˆu − u^∗) = 0, ˆ

uz − y = 0,

−ˆuy − z − 1 = 0.

Hieruit volgt de vergelijking voor het evenwichtspunt (ˆu − u^∗)(ˆu²+ 1) −^Bˆ_ˆ

Cu = 0. Merkˆ op dat dit overeenkomt met vergelijkingen uit de analyse hierboven als we u^∗ = 0 stellen.

In figuur 4.3 en 4.4 zien we voor ˆB = 127 en u^∗ = 0 en u^∗ = 0, 95 de (u, t) en (u, y)-plots. Verder zijn in figuur 4.2 de bifurcatiediagrammen te zien.

Als we nu goed kijken naar figuur 4.2, zien we dat de bifurcatiediagrammen van u^∗ = 0 en u^∗ = 0.95 erg van elkaar verschillen. Door deze u^∗ te introduceren verdwijnt de symmetrie in het diagram en daarmee de symmetrie tussen oost en west. Ook treed er een andere bifurcatie op. Bij u^∗ = 0 zien we een co-dimensie ´e´en pitchfork bifurcatie terwijl we bij u^∗ = 0.95 een co-dimensie twee bifurcatie zien. Verder zien we dat voor u^∗ = 0.95 geldt dat rond ˆB = 10, 6 (af te lezen) de twee nieuwe oplossingen p₀ en p₊ erbij komen. Daarnaast is p− stabiel tot ongeveer ˆB = 101, dan krijgen twee van de eigenwaarden een postief ree¨el deel. p+ is stabiel tussen ˆB = 10, 6 en ˆB = 36, daarna wordt het een zadelpunt. Ten slotte is p₀ steeds een zadelpunt. We merken meteen dat deze bifurcatie een stuk complexer is dan die van u^∗ = 0. Echter is er ´e´en belangrijke overeenkomst: Na ˆB = 101 zijn alle drie punten zadelpunten en hebben dus alsnog te maken met chaos, zoals ook te zien in figuur 4.4 [35].

Figuur 4.2: Bifurcatie diagrammen a)u^∗ = 0 en b) u^∗ = 0.95.

Figuur 4.3: (u, t)-plot, ˆB = 127 a) u^∗ = 0 en b) u^∗ = 0.95.

Figuur 4.4: (u, y)-plot, ˆB = 127 a) u^∗ = 0 en b) u^∗ = 0.95.

4.2.1 Lyapunov exponenten

Naast deze manier om chaos op te sporen, kunnen we dat ook doen met de zogenaamde Lyapunov karakteristieke exponenten. Deze methode is vernoemd naar de Russische wiskundige A.M. Lyapunov (1857-1918). Voor een systeem met n vergelijkingen zijn er ook n Lyapunov exponenten. We gaan eerst kijken naar 1 variabele en dus 1 exponent, vervolgens zullen we dit uitbouwen naar n variabelen. We beschouwen een systeem met twee beginvoorwaarden, die een heel klein beetje van elkaar verschillen, zeg ; dus we hebben x₀ en x₀+. We willen nu kijken wat er na n iteraties gebeurt is, met andere woorden we willen x_n van beide beginwaarden met elkaar vergelijken. De Lyapunov exponent is de co¨effici¨ent van gemiddelde exponenti¨ele groei per tijdseenheid tussen deze twee toestanden. Het verschil d_n na n iteraties tussen de twee waarden van x_n is te benaderen met:

d_n= e^nλ.

Aan deze vergelijking kunnen we meteen zien dat ze convergeren wanneer λ negatief is en ze divergeren als λ positief is. In dat laaste geval is er dus sprake van chaos.

Zoals eerder gezegd beginnen we met een ´e´en dimensionale afbeelding, namelijk x_n+1 = f (x_n). Het verschil in waarden x_n is in het begin d₀ =  en na 1 iteratie:

d₁ = f (x₀+ ) − f (x₀) ≈ df dx|_x0.

De laatste benadering aan de rechterkant kan je doen aangezien  erg klein is. Verder is dn, dus het verschil na n iteraties gegeven door

d_n= fⁿ(x + ) − fⁿ(x₀) = e^nλ.

Wanneer we deze vergelijking omschrijven, door te delen door  en de ln te nemen aan beide kanten krijgen we

ln fⁿ(x + ) − fⁿ(x₀))





= ln(e^nλ) = nλ.

en daarmee

λ = 1

n ln fⁿ(x +  − fⁿ(x₀))





= 1 n ln

dfⁿ(x) dx |_x0

    .

Met de notatie fⁿ(x₀) bedoelen we de n^eiteratie van f (x₀), dus fⁿ(x₀) = f (f (...(f (x₀))...)).

Vervolgens passen we de kettingregel toe dfⁿ(x)

dx |_x0 = df

dx|_xn−1df

dx|_xn−2...df dx|_x0. Als laatste stap nemen we de limiet van t → ∞

λ = lim

n→∞

1 n

n−1

X

i=0

ln    

df (xi) dx

    .

Voor ons systeem met drie vergelijkingen, zullen we dus ook 3 Lyapunov exponenten krijgen. Hiervoor gebruiken we in plaats van de afgeleide, de jacobiaan in de voor- afgaande uitwerking. We zullen kort uitleggen hoe deze waarden verkregen kunnen worden. Vervolgens zullen we Matlab het werk laten doen en met behulp van de 4^e orde Runge-Kutta (RK4) methode de gewenste waardes berekenen.

Voor de ontwikkeling van het verschil , hebben weer de linearisatie nodig.

d

dt =





− ˆC Bˆ 0

z −1 u

−y −u −1



.

We beginnen met drie orthogonale vectoren t1 = (1, 0, 0), t2 = (0, 1, 0) en t3 = (0, 0, 1) en gaan kijken hoe de componenten van elke vector zich ontwikkelen volgens 4.2.1. Na een aantal iteraties bepalen we de vergroting van t₁ en normaliseren deze. Vervolgens projecteren we t2 normaal ten opzichte van t1 en bepalen de vergroting van de vector die we hiermee verkijgen en normaliseren deze. Ten slotte projecteren we t₃ zo dat deze zowel loodrecht staat op t₁ als t₂. Ook deze vector wordt genormaliseerd. Wanneer we het product nemen van al deze vergrotingsfactoren van heel veel iteraties krijgen we e^λit. Verder is in het bijzonder op ons systeem van toepassing dat P λ_i = trace Jacobiaan = − ˆC − 2 [13].

Met behulp van de Matlab code die te vinden is in de appendix zijn de volgende resultaten verkregen:

λ₁ λ₂ λ₃ P λ_i

B = 60ˆ −0, 004 −0, 010 −4, 986 −5, 000 B = 65ˆ 0, 422 0, 001 −5, 423 −5, 000 B = 127ˆ 0, 573 0, 001 −5, 572 −4, 998

Tabel 4.2: De waarden van λ voor verschillende B’s met t = 100000 en u^∗ = 0

λ₁ λ₂ λ₃ P λ_i

B = 95ˆ −0, 018 −0, 019 −4, 963 −5, 000 B = 105ˆ 0, 388 0, 002 −5, 388 −4, 999 B = 127ˆ 0, 465 0, 002 −5, 465 −4, 998

Tabel 4.3: De waarden van λ voor verschillende B’s met t = 100000 en u^∗ = 0, 95 Hierbij hebben we t = 100000 genomen, deze waarde is groot genoeg om een goede benadering te krijgen van de Lyapunov exponenten λ_i. Zoals hierboven ook al kort genoemd leidt een positieve λ tot exponti¨ele divergentie en daarmee ontstaat chaos.

Dus wanneer ´e´en van de λ groter dan nul wordt, is er sprake van chaos.

In tabel 4.2 zien we dat inderdaad de som van de exponenten steeds uitkomt op

− ˆC − 2 = 5. Ook zien we dat inderdaad het omslagpunt van stabiel naar chaotisch ligt tussen ˆB = 60 en ˆB = 65. Dit komt overeen met de waarde ˆB_h = 63 die we in paragraaf 5.1 hebben bepaald. Ook de overgangswaarde van stabiel naar chaotische in tabel 4.3 (u^∗ = 0, 95) komt overeen met de ˆBh = 101, die berekend is in deze paragraaf [13, 33].

4.3 Analyse u

= u

(1 + u

sin ωt)

Om het model nog iets dichter bij de realiteit te brengen introduceren we een sei- zoenscyclus. Dit omdat gebleken is dat een El Ni˜no samenhangt met de tijd van het jaar. De kans op een El Ni˜no is namelijk veel groter wanneer de passaatwinden zwak zijn of oostelijk in tegenstelling tot de gemiddelde westelijke wind. Om deze seizoensgebondenheid door te voeren in ons systeem laten we u^∗ vari¨eren met t:

u^∗ = u₀(1 + u_ssin ωt).

De nieuwe parameters in dit model zijn u_s, die staat voor de amplitude en de periode ω, die we op 2π per jaar stellen. Hoe we deze u_s moeten kiezen is nog discutabel.

Enerzijds moeten we deze niet te groot kiezen, omdat de relatieve invloed van u₀ dan te verwaarlozen is en daarmee de oost-west symmetrie terugkeert, net als bij u^∗ = 0.

Anderzijds moet hij groot genoeg zijn om de passaatwinden om te keren, dus u_s > 1.

We gaan nu kijken wat er met het model gebeurt voor verschillende waardes voor B en xˆ _s (u₀ houden we op 0, 95). Om nog even te herhalen ˆB geeft aan hoe sterk de winden u en het temperatuurverschil y elkaar be¨ınvloeden. Dit zien we ook terug in figuur 4.5. In de onderste twee plots ( ˆB = 191 in plaats van ˆB = 127), zien we namelijk meer pieken en dus El Ni˜no’s.

Figuur 4.5: ut-plots, met a) ˆB = 127, u_s = 0, b) ˆB = 127, u_s= 3, c) ˆB = 191, u_s = 0 en d) ˆB = 191, u_s= 3.

Waar we nu juist ge¨ıntereseerd in zijn is de samenhang tussen de seizoensgebonden cycli die we hebben ge¨ıntroduceerd en het begin van een El Ni˜no. In figuur 4.6 zie je dat er inderdaad een samenhang is tussen de waarde van u^∗ en het begin van een El Ni˜no. Daaruit blijkt meteen dat het introduceren van goniometrische functie in u^∗ een goede keuze was om het model realistischer te maken. Wel in opvallen

dat door de seizoensgebondencyclus ook bij dubbel zoveel El Ni˜no ontstaan. Hieruit kan je concluderen dat dit toch nog niet een optimale keuze is of dit model zelf niet realistische genoeg is [27, 35, 41, 42].

Figuur 4.6: Op de x-as zie je het aantal El Ni˜no’s en op de y-as de maanden van het jaar. Verder is de onderbroken lijn de ontwikkeling van u^∗ gedurende het jaar. Voor beide figuren geldt ˆB = 191 en u₀ = 0, 95. Verder is voor a) u_s = 0 en voor b) u_s= 3 [35].

Hoofdstuk 5

Conclusie en discussie

Met de analyse van het model van Vallis in ons achterhoofd kunnen we zeggen dat we nog niet zijn beland bij een realistische voorspelling. Zo verdubbelt bijvoorbeeld het introduceren van een seizoensgebondencyclus het aantal El Ni˜no’s. Dit is natuurlijk niet de bedoeling. Een beter keuze zou misschien zijn om T^∗ met de tijd te laten vari¨eren, door T^∗ te vervangen door T^∗(1 + sin(ωt)). Toch zal ook dit waarschijnlijk niet een goede voorspelling geven.

Zoals in hoofdstuk 3 ook al aangegeven, zijn er 3 types modellen. Vallis’ model is simpel model waarbij we interessante eigenschappen van El Ni˜no hebben bekeken.

Willen we optimale voorspelling doen, zullen we moeten kijken naar de zogenaamde CGCM’s.

De huidige instuten die op korte termijn goede voorspelling maken zijn European Centre for Medium-Range Weather Forecasts (ECMWF) en National Centre of Envi- ronmental Prediction (NCEP). In figuur 5.1 zien we dat de gemiddelde temperatuur van het oppervlakte water in juni al ruim 0, 6^◦C boven de normale termperatuur ligt.

Dit geeft een sterke aanleiding tot een relatief sterke El Ni˜no. Of deze voorspelling daadwerkelijk uit zal komen is nog de vraag. Wat wel zeker is, is dat de modellen steeds meer duidelijkheid kunnen geven. Laten we het meemaken [9].

Figuur 5.1: Voorspelling van de SST in regio NINO3.4 van ECMWF (1 juni 2014) [9].

Bibliografie

[1] Art and tels diary of New Zealand (NZ) - Walker Circulation, El Ni˜noand La Ni˜na.

[2] National Weather Service - Monthly Atmospheric & SST Indices.

[3] Nation of change, URL: http://www.nationofchange.org/el-ni-o-could- make- 2014-hottest-year-record-1402502208, 2014.

[4] K. AchutaRao and K. Sperber. Simulation of the El Ni˜no Southern Oscillationn:

Results from the coupled model intercomparison project. Climate Dynamics, 19(3-4), 2002.

[5] K. AchutaRao and K. Sperber. Enso simulation in coupled ocean-atmosphere models: are the current models better? Climate Dynamics, 27, 2006.

[6] D.L.T. Anderson and J.P. McCreary. On the role of the indian ocean in a coupled ocean-atmosphere model of El Ni˜no and the Southern Oscillation. Journal of the Atmospheric Sciences, 1985.

[7] D.L.T. Anderson and J.P. McCreary. Slowly propagating disturbances in a cou- pled ocean-atmosphere model. Journal of the Atmospheric Sciences, 1985.

[8] Several authors. El Ni˜no, URL: http://en.wikipedia.org/.

[9] Several authors. Website van ecmwf, http://www.ecmwf.int.

[10] J. Bjerknes. A possible response of the atmospheric Hadley circulation to equa- torial anomalies of ocean temperature. Tellus, 1966.

[11] L. Bunge and A.J. Clarke. A verified estimation of the El Ni˜no index Nino-3.4 since 1877. Department of Oceanography, The Florida State University, Talla- hassee, Florida.

[12] G.J.H. Burgers and G.J. van Oldenborgh. El Ni˜no en La Ni˜na. www.knmi.nl.

[13] M. Cross. Introduction to Chaos - Chapter 7 Lyapunov exponents. internet.

[14] en A. von der Heydt Dijkstra, H.A. El Ni˜no: processes, dynamics and predicta- bility. Nieuw Archief voor Wiskunde, 2013.

[15] H.A. Dijkstra. Nonlinear Physical Oceanography. Springer Berlin Heidelberg, 2005.

[16] H.A. Dijkstra. Dynamical Oceanography. Springer Berlin Heidelberg, 2008.

[17] H.A. Dijkstra. Nonlinear Climate Dynamics. Cambridge, 2013.

[18] H.A. Dijkstra and G. Burgers. Fluid dynamics of El Ni˜no variability. Annual Review of Fluid Mechanics, 34, 2002.

[19] R.D. Euz´ebio and J. Llibre. Periodic solutions of El Ni˜no model through the vallis differential system. Discrete and continuous dynamical systems, 34(9), 2014. Only introduction.

[20] V. Govorukhin. Calculation lyapunov exponents for ode. internet.

[21] N.E. Graham and W.B. White. The El Ni˜no cycle: A natural oscillator of the pacific ocean-atmosphere. Science, 1998.

[22] A.C. Hirst. Unstable and damped equatorial modes in simple coupled ocean- atmosphere models. Journal of the Atmospheric Sciences, 1986.

[23] et al Hirst, A.C. Coupled Ocean-Atmosphere Models - Free equatorial instabilites in Simple coupled atmosphere-ocean models. Elsevier, 1985.

[24] F.F. Jin. Tropical ocean-atmosphere interaction, the pacific cold tongue, and the El Ni˜no-Southern Oscillation. Science, New Series, 274, 1996.

[25] M. Latif. ENSIP: the El Ni˜no simulation intercomparison project. Springer, 2001.

[26] K.M. Lau. Oscillation in a simple equatorial climate system. Journal of the Atmospheric Sciences, 1981.

[27] J. Llibre. Periodic solutions via averaging theory., Januari 2014.

[28] Google Maps. Google maps. www.maps.google.nl.

[29] J.P. McCreary. A model of tropical ocean-atmosphere interaction. Monthly Wea- ther Review, 1983.

[30] J.P. McCreary and D.L.T. Anderson. A simple model of El Ni˜no and the Southern Oscillation. Monthly Weather Review, 1984.

[31] J.P.Jr. McCreary. An overview of coupled ocean-atmosphere model of El Ni˜no and the Southern Oscillation. Journal of Geophysical Research: Oceans, 96, 1991.

[32] S.G.H. Philander, N.C. Lau, R.C. Pacanowski, and M.J. Nath. Two different si- mulations of the Southern Oscillation and El Ni˜no with coupled ocean-atmosphere general circulation models. Phil. Trans. Roy. Soc. (submitted), 1989.

[33] M.S. Pilant. Calculating the entire lyapunov spectra of the lorenz attractor.

internet.

[34] P.S. Schopf and M.J. Suarez. Vacillations in a coupled ocean-atmosphere model.

Journal of the Atmospheric Sciences, 1988.

[35] D. Strozzi. On the origin of interannual and irregular behavior in the El Ni˜no.

Master’s thesis, Princetion University, 1999.

[36] M.J. Suarez and P.S. Schopf. A delayed action oscillatior for enso. Journal of the Atmospheric Sciences, 1988.

[37] Taichiatduke. El Ni˜no 1. introduction, 2. beginning of El Ni˜no, 3. peru during El Ni˜no, 4. recovery of El Ni˜no., 2009. www.youtube.nl.

[38] A. Timmermann. Decadal enso amplitude modulations: a nonlinear paradigm.

Global and Planetary Change, 37, 2013.

[39] A. Timmermann, F.F. Jin, and J. Abshagen. A nonlinear theory for El Ni˜no bursting. J. Atmos. Sci., 60(1), 2003.

[40] F.F. Jin Timmermann, A. A nonlinear mechanism for decadal El Ni˜no amplitude changes. Geophysical Research Letters, 29(1), 2002.

[41] G.K. Vallis. El Ni˜no: A chaotic dynamical system? Science, 243(232), 1986.

[42] G.K. Vallis. Conceptual models of El Ni˜no and the southern oscillation. Journal of Geophysical Research, 93(13), 1988.

[43] T. Yamagata. Coupled Ocean-Atmosphere Models - Stability of a simple air-sea coupled model in the tropics. Elsevier, 1985.

[44] S.E. Zebiak and M.A. Cane. A model of El Ni˜no-Southern Oscillation. Monthly Weather Review, 1987.

Bijlage A

Routh-Hurwitz criterium

Het karakteristieke polynoom wat we gaan bekijken is de volgende:

λ³+ ( ˆC + 2)λ²+ ( ˆC + ˆB/ ˆC)λ + 2( ˆB − ˆC) = 0.

Hiervoor is het niet handig om eerst de nulpunten te gaan zoeken. Vandaar dat we gebruik maken van het zogenaamde Routh-Hurwitz criterium. Hierbij gebruiken we een transferfunctie waarvan ons karakteristieke polynoom de noemer is.

Theorem A.0.1 Routh-Hurwitz criterium Alle nulpunten van het karakteristieke po- lynoom hebben een negatief re¨eel deel dan en slechts dan als een bepaalde set van algebra¨ısche combinaties van co¨effici¨enten allemaal hetzelfde teken hebben.

Hoe deze bepaalde set gedefinieerd is, zullen we nader toelichten. In het algemeen krijgen we voor de polynoom a₃λ³+ a₂λ²+ a₁λ + a₀ de volgende tabel:

λ³ a3 a1

λ² a₂ a₀

λ¹ (a₁a₂− a₀a₃)/a₂ 0

λ⁰ a₀ 0

De bepaalde set die in het theorema genoemd is, bestaat uit alle waarden in de tweede kolom van de bovenstaande tabel. Wanneer deze waarden dus allemaal het- zelfde teken hebben, hebben alle nulpunten van het karakteristieke polynoom een negatief re¨eel deel. Wanneer dit het geval is, is het systeem dus stabiel.

Wij zijn ge¨ınteresseerd in het omslagpunt van stabiel naar onstabiel. Daarvoor moeten we dus opzoek naar een waarde van ˆB waarbij het teken omslaat. De tabel die bij ons polynoom hoort is de volgende:

λ³ 1 C +ˆ ^Bˆ_ˆ

C

λ² C + 2ˆ 2( ˆB − ˆC)

λ¹ ? 0

λ⁰ 2( ˆB − ˆC) 0

waarbij ? = ^{( ˆC+2)( ˆC+}

Bˆ

Cˆ)−2( ˆB− ˆC)

C+2ˆ .

Van de eerste twee waarden in de tweede kolom, namelijk 1 en ˆC + 2 weten we dat ze positief zijn. Ook van 2( ˆB − ˆC) weten we dat deze groter is dan nul, aangezien B > ˆˆ C. Om stabiliteit te cre¨eren is dus nodig dat ook ? groter is dan nul. Het omslagpunt ligt dus bij ? = 0.

( ˆC + 2)( ˆC + ^Bˆ_ˆ

C) − 2( ˆB − ˆC)

C + 2ˆ = 0,

( ˆC + 2)( ˆC + Bˆ

Cˆ) − 2( ˆB − ˆC) = 0, Cˆ²+ 4 ˆC − ˆB + 2Bˆ

Cˆ = 0, Cˆ²+ 4 ˆC = ˆB − 2Bˆ

Cˆ Cˆ²( ˆC + 4) = ˆB( ˆC − 2),

B =ˆ

Cˆ²( ˆC + 4) C − 2ˆ . Dit omslagpunt noemen we B_h = ^{Cˆ2( ˆ}_ˆ^C+4)

C−2 .

Bijlage B

Mathematica notebook codes

B.1 ut-plot, uy-plot en 3D plot

Clear[x, y, z]

x0 = 1; y0 = 1; z0 = 1;

Tend = 50;

initialConditions = {x[0] == x0, y[0] == y0, z[0] == z0};

LorenzEquations = {x’[t] ==

127*y[t] - 3*(x[t] + 0.95*(1 + 0*Sin[t*2 Pi])),

y’[t] == x[t]* z[t] - y[t], z’[t] == -x[t]*y[t] - z[t] + 1, initialConditions}; 12

s1 = NDSolve[LorenzEquations, {x[t], y[t], z[t]}, {t, 0, Tend}, MaxSteps -> \[Infinity]];

Plot[Evaluate[x[t] /. s1], {t, 0, Tend}, PlotRange -> All]

Plot[Evaluate[y[t] /. s1], {t, 0, Tend}, PlotRange -> All]

ParametricPlot[{x[t], y[t]} /. s1, {t, 0, Tend}, AspectRatio -> 1]

ParametricPlot3D[{x[t], y[t], z[t]} /. s1, {t, 0, Tend}, AspectRatio -> 1]

B.2 Bifurcatie plot

Clear[f];

Needs["Graphics‘Colors"]

f[x_, B_] := (x + 0.021) (x^2 + 1) - B/3 x Bmin = 0;

Bmax = 50;

xmin = -10;

xmax = 10;

sols = Solve[f[x, B] == 0, x];

p2 = Plot[ Evaluate[Table[y /. sol[[i]], {i, 1, Length[sol]}]], {B, Bmin, Bmax}, PlotStyle -> {{AbsoluteThickness[3.5], White}}, DisplayFunction -> Identity];

p1 = ContourPlot[f[x, B], {B, Bmin, Bmax}, {x, xmin, xmax},

Contours -> {0}, ContourShading -> True, FrameLabel -> {B, x}, PlotPoints -> 50,

ContourStyle -> {{White, AbsoluteThickness[1]}},

ColorFunction -> (RGBColor[Sign[#1], 0, 1 - Sign[#1]] &), DisplayFunction -> Identity];

Show[p1, p2, PlotRange -> {xmin, xmax}, DisplayFunction -> $DisplayFunction]

Bijlage C

Matlab codes

C.1 Vallis system

function f = vallis_ext(t,X)

% Vallis system

% Parameters B = 127;

C = 3;

x = X(1); y = X(2); z = X(3);

Y = [X(4), X(7), X(10);

X(5), X(8), X(11);

X(6), X(9), X(12)];

f = zeros(9,1);

% Vallis vergelijkingen f(1)=B*y-C*(x+0.95);

f(2)=x*z-y;

f(3)=-x*y-z+1;

% Gelineariseerd systeem Jac = [-C, B, 0;

z, -1, x;

-y, -x, -1];

f(4:12) = Jac*Y;

bron: [20].

C.2 Lyapunov exponenten

function [Texp,Lexp] = lyapunov(n,rhs_ext_fcn,fcn_integrator,tstart, stept,tend,ystart,ioutp);

% Input parameters:

% n - aantal vergelijkingen

% rhs_ext_fcn - handle of functions met de rechterkant van de de ODE’s.

% fcn_integrator - handle of ODE integrator function, wij gebruiken:

@ode45

% tstart - beginwaarde van t.

% stept - tijdstappen

% tend - eindwaarde van t

% ystart - beginwaarden in ODE.

% ioutp - step of print to MATLAB main window. ioutp==0 - no print,

% if ioutp>0 then each ioutp-th point will be print.

% Output parameters:

% Texp - tijdstip

% Lexp - Lyapunov exponents op dat tijdstip

% n= aantal niet-lineaire ODE’s

% n2=n*(n+1)= totaal aantal ODE’s n1=n; n2=n1*(n1+1);

% Aantal stappen

nit = round((tend-tstart)/stept);

y=zeros(n2,1); cum=zeros(n1,1); y0=y;

gsc=cum; znorm=cum;

% Beginwaarden y(1:n)=ystart(:);

for

i=1:n1 y((n1+1)*i)=1.0; end;

t=tstart;

% Hoofd loop for

ITERLYAP=1:nit

[T,Y] = feval(fcn_integrator,rhs_ext_fcn,[t t+stept],y);

t=t+stept;

y=Y(size(Y,1),:);

for i=1:n1

for j=1:n1 y0(n1*i+j)=y(n1*j+i); end;

end;

%

% Construeer een nieuwe basis mbv Gram-Schmidt znorm(1)=0.0;

for j=1:n1 znorm(1)=znorm(1)+y0(n1*j+1)^2; end;

znorm(1)=sqrt(znorm(1));

for j=1:n1 y0(n1*j+1)=y0(n1*j+1)/znorm(1); end;

for j=2:n1 for k=1:(j-1) gsc(k)=0.0;

for l=1:n1 gsc(k)=gsc(k)+y0(n1*l+j)*y0(n1*l+k); end;

end;

for k=1:n1 for l=1:(j-1)

y0(n1*k+j)=y0(n1*k+j)-gsc(l)*y0(n1*k+l);

end;

end;

znorm(j)=0.0;

for k=1:n1 znorm(j)=znorm(j)+y0(n1*k+j)^2; end;

znorm(j)=sqrt(znorm(j));

for k=1:n1 y0(n1*k+j)=y0(n1*k+j)/znorm(j); end;

end;

% Update grote van de vector

for k=1:n1 cum(k)=cum(k)+log(znorm(k));

end;

% normalize exponent for k=1:n1

lp(k)=cum(k)/(t-tstart);

end;

% Aanpassen output if ITERLYAP==1 Lexp=lp;

Texp=t;

else

Lexp=[Lexp; lp];

Texp=[Texp; t];

end;

if (mod(ITERLYAP,ioutp)==0) fprintf(

’t=%6.4f’,t);

for k=1:n1 fprintf(’ %10.6f’,lp(k)); end;

fprintf(

’\n’);

end;

for i=1:n1 for j=1:n1

y(n1*j+i)=y0(n1*i+j);

end;

end;

end;

bron: [20].

C.3 Run script

[T,Res] = lyapunov(3,@vallis_ext,@0de45,0,1,100000, [1,1,1], 1);

plot(T,Res);

title(’Dynamica van de Lyapunov exponenten) xlabel(’tijd’);

ylabel(’Lyapunov exponent’);

bron: [20].

Bijlage D

Data NINO3.4 en SOI

D.1 NINO 3.4

Bron: [2].

Bron: [2].

D.2 SOI

Bron: [2].

Bron: [2].