Studie av jonosfären och troposfärens påverkan av vulkanutbrott med GNSS Undersökning av hur GNSS-data kan användas för att studera vulkanut- brottet av Hunga Tonga Hunga Ha’apais påverkan på jonosfären och tro- posfären Kandidatarbete MATILDA HELLSTRÖM, MOA KALLÉN EDVIN KARLSSON, AXEL SKOOG INSTUTIONEN FÖR RYMD-, GEO- OCH MILJÖVETENSKAP CHALMERS TEKNISKA HÖGSKOLA Göteborg, Sverige 2024 www.chalmers.se www.chalmers.se Sammanfattning Den här rapporten utreder möjligheterna att kartlägga utbredningen av en tryckvåg genom att använda GNSS för att betrakta dess påverkan på troposfärens signal- fördröjning och elektrontätheten i jonosfären. Metoden har huvudsakligen utgjorts av bearbetning av GNSS- samt tryckdata från stationer geografiskt utspridda kring vulkanen Hunga Tonga Hunga Ha’apai. Vid analys av troposfärsfördröjningen har data bearbetats med mjukvaran GipsyX och till elektrontätheten har data förädlats med Matlabkod utvecklad av Long Tang [1]. Genom att framställa erhållna data för stationerna grafiskt har signalfördröjningen, Zenith Total Delay (ZTD), i troposfären samt den vinklade totala elektrontätheten (STEC) i jonosfären kunnat jämföras och analyseras. Eftersom tryckvariationerna orsakade av tryckvågen långt från utbrottsplatsen inte var stora nog för att ge upp- hov till en märkbar förändring av ZTD fanns svårigheter att använda detta mått vid analys av stationer mer avlägsna från vulkanen. Vidare förekom avvikelser i ZTD- variationen beroende på huruvida landmassa eller enbart vatten var beläget mellan utbrottsplatsen och mätstationen. Eftersom tryckvariationen föreföll minska i om- fattning då vågen propagerade över landmassor. Däremot kunde måttet bekräftas vara användbart för att erhålla hög precision för hastighetsbestämningen i närhe- ten av utbrottsplatsen. STEC visade sig vara mindre avståndberoende jämfört med ZTD men, varför denna bedömdes fungera väl för att bestämma en approximativ tryckvågshastighet, även på längre avstånd. Däremot skedde detta på bekostnad av precsionen. Med en estimerad tidsstämpel för när vågen bör ha passerat respektive station samt med överensstämmande resultat för både elektrontäthet och signalfördröjning kunde information om tryckvågens propagering genom atmosfären, så som utbredningshas- tighet, kartläggas tämligen väl. Sammanfattningsvis har studien bekräftat använd- barheten av att kombinera analys av både ZTD- samt STEC-data såvida tryckvågs- utbredningen orsakad av ett vulkanutbrott ska studeras. Därav bedöms rapportens syfte samt frågeställningar ha besvarats. i Abstract This report investigates the possibilities of using GNSS to map a volcanic pressure wave propagation in the troposphere and its impact on electron density in the io- nosphere. The method consisted of processing GNSS data and pressure data from stations geographically spread around the Hunga Tonga Hunga Ha’apai volcano using the GipsyX software and Matlab code developed by Long Tang [1]. By graphically presenting the obtained data for the stations, the signal delay, Ze- nith Total Delay (ZTD), in the troposphere, and the slanted total electron content (STEC) in the ionosphere could be compared and analyzed. Since the pressure varia- tions far from the eruption site caused by the pressure wave were not large enough to cause a noticeable change in ZTD, there were difficulties in using this measurement when analyzing stations further away from the volcano. Furthermore, the deviations in the ZTD depended on whether landmasses or only water was located between the eruption site and the measurement station due to smaller pressure variations when the wave propagated over landmasses. However, the measure could be confirmed to be useful whrn determining a high precision velocity near the eruption site. STEC was found to be less distance dependent compared to ZTD and thus, it was con- sidered to work better for determining an approximate pressure wave velocity over longer distances compared to the ZTD. The precision of the result was, however, lower for STEC when compared to ZTD. With an approximate timestamp for when the wave should have passed each station, and with consistent results for both STEC and ZTD, information about the propa- gation of the pressure wave through the atmosphere, such as propagation velocity, could be mapped. In summary, the study has confirmed the usefulness of combi- ning both ZTD and STEC data when studying pressure wave propagation caused by a volcanic eruption. Hence, the purpose and research questions of the report are considered to have been answered. Tillkännagivanden Tack till vår handledare Rüdiger Haas för fantastiskt stöd och hjälp under arbetets gång [2]. Tack till vår examinator Jan Johansson för all hjälp och skriptet till Vidar [3]. Tack till Onsala Rymdobservatorium för tillgång till datorn Vidar. Tack till Long Tang för utvecklingen av den Matlabmjukvara som använts i projektet för beräkning av den totala elektrontätheten [1]. Tack till NASA och JPL för mjukvaran GipsyX. Tack till Fackspråk och de studenter som bidragit med respons som underlättat skrivprocessen. Tack till GMT Steering Committee för The Generic Mapping Tools. - Matilda Hellström, Moa Kallén, Edvin Karlsson, Axel Skoog, Göteborg, maj 2024 iii Lista med Akronymer Nedan följer en lista med de akronymer som används i rapporten. dSTEC Förändring i vinklad total elektrontäthet, delta slant Total Electron Content GLONASS Rysslands Globala navigationssatellitsystem, GLObalnaya NAvi- gatsionnaya Sputnikovaya Sistema GMT Generiskt kartverktyg, Generic Mapping Tool GNSS Globalt navigations satellit system, Global Navigation Satellite Systems GPS Globalt Positioneringssystem, Global Positioning System HTHH Hunga Tonga–Hunga Ha’apai IGS Internationella GNSS tjänsten, International GNSS Service LP-filter Lågpass filter, Lowpass Filter PRN Pseudoslumpmässigt brus, Pseudo Random Noise STEC Vinklad total elektrontäthet, Slant Total Electron Content STECU Vinklad total elektrontäthet enhet, Slant Total Electron Content Unit TEC Total elektrontäthet, Total Electron Content TECU Total elektrontäthet enhet, Total Electron Content Unit ZHD Hydrostatisk zenit fördröjning, Zenith Hydrostatic Delay ZTD Total zenit fördröjning, Zenith Total Delay ZWD Zenit fuktighetsfördröjning, Zenith Wet Delay iv Innehåll Lista med Akronymer iv 1 Introduktion 1 1.1 Syfte . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2 Frågeställningar . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.3 Avgränsningar . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 2 Teori 4 2.1 Global Navigation Satellite System . . . . . . . . . . . . . . . . . . . 4 2.1.1 Uppbyggnad och funktion . . . . . . . . . . . . . . . . . . . . 4 2.1.2 Global Positioning System . . . . . . . . . . . . . . . . . . . . 5 2.1.3 Globalnaya Navigatsionnaya Sputnikovaya Sistem . . . . . . . 6 2.1.4 Galileo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.2 Vulkanism . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.3 Hunga Tonga–Hunga Ha’apai . . . . . . . . . . . . . . . . . . . . . . 8 2.4 Atmosfär . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.4.1 Troposfär . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.4.2 Jonosfär . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.5 Mjukvara . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.5.1 RINEX . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.5.2 GipsyX . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.5.3 IonKit-NH . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.6 Beräkningsmetoder och formler . . . . . . . . . . . . . . . . . . . . . 10 2.6.1 Zenith Total Delay, Zenith Wet Delay och Zenith Hydrostatic Delay . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.6.2 Saastamoinens modell . . . . . . . . . . . . . . . . . . . . . . 10 2.6.3 Savitzky-Golay filter . . . . . . . . . . . . . . . . . . . . . . . 11 2.6.4 Förväntad tryckvåg och avståndsbestämning . . . . . . . . . . 11 2.6.5 Total Electron Content . . . . . . . . . . . . . . . . . . . . . . 12 3 Metod 13 3.1 Datainsamling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 3.1.1 Val av tidsupplösning . . . . . . . . . . . . . . . . . . . . . . . 16 3.1.2 Val av elevationsvinkel . . . . . . . . . . . . . . . . . . . . . . 16 3.2 Användning av GipsyX . . . . . . . . . . . . . . . . . . . . . . . . . . 17 3.3 Grafisk framställning av Zenith Total Delay . . . . . . . . . . . . . . 17 v Innehåll 3.4 Beräkning av Total Electron Content . . . . . . . . . . . . . . . . . . 18 3.5 Analys av tryckvågens utbredning . . . . . . . . . . . . . . . . . . . . 18 3.6 Förväntade värden . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 4 Resultat 20 4.1 Uppmätt Zenith Total Delay . . . . . . . . . . . . . . . . . . . . . . . 20 4.2 Uppmätt Zenith Total Delay samt lufttryck . . . . . . . . . . . . . . 23 4.3 Eliminering av Zenith Wet Delay . . . . . . . . . . . . . . . . . . . . 24 4.4 Total Electron Content för 12 stationer . . . . . . . . . . . . . . . . . 25 4.5 Total Electron Content för enskilda satelliter . . . . . . . . . . . . . . 27 5 Diskussion 28 5.1 Zenith Total Delay . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 5.1.1 Zenith Total Delay i kombination med tryckdata . . . . . . . . 28 5.2 12 mätstationers uppmätta elektrontäthet . . . . . . . . . . . . . . . 29 5.3 Vinklade totala elektrontäthet för respektive station . . . . . . . . . . 29 5.4 Datainsamling och felkällor . . . . . . . . . . . . . . . . . . . . . . . 30 6 Slutsats 32 Referenser 33 A Python skript I A.1 GMT karta . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . I A.2 Distans . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . II A.3 ZTD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . V A.4 ZTD och tryck . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . XII B Matlabkod XV B.1 Tid- och distanskarta . . . . . . . . . . . . . . . . . . . . . . . . . . . XV B.2 dTEC-filer, antalet unika satelliter . . . . . . . . . . . . . . . . . . . XVI B.3 dTEC filter, alla stationer . . . . . . . . . . . . . . . . . . . . . . . . XVI B.4 dTEC, tre satelliter . . . . . . . . . . . . . . . . . . . . . . . . . . . . XIX vi 1 Introduktion Att kunna förutse naturkatastrofer möjliggör proaktivt arbete som kan bidra till att minska omfattningen av de problem som drabbade samhällen tvingas konfrontera. I de fall där katastrofer inte lyckas förutsägas i tid riskerar samhällen att bemöta kriser helt utan förberedelser, vilket kan resultera i ansenliga katastrofer, decimerad levnadsstandard och ihärdiga ekonomiska problem. Metoder med vilka katastrofer kan förutsägas är följaktligen av yttersta intresse för att garantera säkerhet på likväl individ- som samhällsnivå. Ett exempel på ett system med vilket exempelvis aktivitet i atmosfären kan detek- teras är Global Navigation Satellite System, GNSS [4]. Redan år 1964 detekterades förändringar i jonosfären på grund av en jordbävning i Alaska [5]. I dryga femtio år har naturfenomens inverkan på jonosfären därav kunnat mätas och genom att data insamlas både innan, under och efter katastroferna går det till viss mån att förutsäga uppkomsten av likande skeenden. Sedan 90-talet har stora datamängder samlats in från mätstationer och satelliter jorden runt med hjälp av GNSS, där majoriteten av alla data publicerats åtkomliga för allmänheten på internet [4]. Insamlade data innehåller information som möj- liggör kartläggning av fenomen som kontinentaldrift liksom studier av atmosfären och dess sammansättning [6]. Metoden utnyttjar att det genom analys av uppmätt tidsfördröjning hos en utsänd signal går att dra slutsatser om det medium i vilket signalen propagerar. Genom att använda satellitdata från ett flertal mätstationer kan information så som utbredningsriktning och -hastighet av vindar eller moln erhållas ur upptagna data, vilket endast är ett av GNSS många användningsområden [4]. I vissa avseenden krävs mer avancerade metoder för analys, varför kraftfull mjukvara för nästintill automatiserad analys av GNSS-data har utvecklats. GipsyX är ett exempel på en mjukvara som framställts av Jet Propulsion Labo- ratory (JPL) för underlättad dataanalys. Användningsområdena är många och då mjukvaran använts i studier i över 35 år bedöms tillförlitligheten vara stor [7]. I det här arbetet kommer 2000-talets mest omfattande vulkanutbrott studeras med GNSS-data som finns att tillgå från ett flertal mätstationer belägna på olika avstånd från vulkanen Hunga Tonga-Hunga Ha’apai (HTHH) [8]. Med en kombination av Matlabkod samt GipsyX kan både jonosfär och troposfär studeras och jämföras 1 1. Introduktion med tryckdata. Det ger, till skillnad från tidigare studier som vanligen fokuserar på endera atmosfärdelen, en mer omfattande förståelse för skeendena i fråga. 1.1 Syfte Arbetets syfte är att undersöka huruvida det är möjligt att studera vulkanutbrott med hjälp av GNSS-data. Mer specifikt kommer utbrottet av vulkanen Hunga Tonga Hunga Ha’apai kartläggas med data från ett antal mätstationer geografiskt spridda kring utbrottsplatsen. Studierna baseras på analys av jonosfärens elektrontäthet samt signalfördröjningen i troposfären, vilket möjliggör undersökning av tryckvågens utbredning. 1.2 Frågeställningar Rapporten kommer att undersöka utbredningen av den tryckvåg som vulkanutbrot- tet orsakade genom att kartlägga variationer i elektrontäthet samt signalfördröjning. Med utgångspunkt i erhållna resultat kommer följande frågeställningar besvaras. 1. Är det möjligt att med hjälp av den totala elektrontätheten i jonosfären kart- lägga när utbrottets tryckvåg passerade respektive station? 2. Vilka faktorer avgör huruvida ZTD eller STEC är lämpligast att använda för analys av tryckvågens utbredning? 3. Är det möjligt att använda GNSS-data för att bestämma tryckvågens utbred- ningshastighet? 1.3 Avgränsningar I projektet har data från ett begränsat antal stationer studerats eftersom det med- för stora, men hanterbara datamängder för projektets omfattning. De data som analyserats har primärt omfattat den dag då vulkanutbrottet skedde, men i viss utsträckning även dagen innan samt efter utbrottet. Upplösningen av GNSS-data, som innebär hur ofta markstationen hämtar data från satelliterna, har inte understigit 30 sekunder. Anledningen är att data upptagna var- je sekund, inte fanns tillgängliga för samtliga stationer. Vidare har en cut-off vinkel på 7° använts. Lägre elevationvinklar resulterar i att data från en längre sträcka inkluderas, vilket resulterar i ökad mätosäkerhet [9]. Endast data från Galileo, GLONASS och GPS har använts. Ytterligare satellitsy- stem hade eventuellt komplicerat datainsamlingen och hade sannolikt inte fungerat felfritt vid användningen av GipsyX. Dessutom har datamängderna från de tre sy- stemen varit tillräckliga för ändamålet i fråga. 2 1. Introduktion Troposfären har här valts att starta vid havsnivå och sluta vid 20 km över havet, medan motsvarande gränser för jonosfären har satts från 60 km till 1000 km ovanför havsytan. Eftersom övergången mellan atmosfärens delar är förhållandevis diffusa är gränserna inte allmängiltiga, men av praktiska skäl kommer övergången här be- handlas som en skarp gräns. Vidare har rapporten inte behandlat miljöaspekter så som utsläpp, skador på infra- struktur och byggnader samt eventuella dödsfall. Orsaken till att det uteslutits är dels att det inte kan mätas med GNSS-data, dels att det inte bedöms nödvändigt för att besvara de givna frågeställningarna. Vid analys av GNSS-data har mätosäkerheter gällande eventuella klockfel från sa- telliter och markstationer ej behandlats. Det beror på att erhållna satellitdata antas vara tillräckligt korrekt för att uppfylla syftet. 3 2 Teori I följande avsnitt presenteras den teoretiska bakgrund som krävs för att tillgodogöra sig resterande del av rapportens innehåll. Inledningsvis ges bakgrund till GNSS, både vad gäller uppbyggnad och funktion, vilket följs av information om vulkanism och vulkanen HTHH. Vidare beskrivs atmosfärens uppbyggnad samt den mjukvara som använts. Avslutningsvis presenteras de beräkningsmetoder och formler som ingår i projektets beräkningar. 2.1 Global Navigation Satellite System Global Navigation Satellite System, förkortat GNSS, är ett samlingsnamn för satel- litbaserade navigationssystem som möjliggör positionsbestämning med hög precision över hela världen [10]. GNSS utgörs av USA:s GPS, Rysslands GLONASS, Europas Galileo och Kinas Beidou. Till GNSS kan också ett fåtal regionala system inräknas. Samtliga av dessa system består av ett flertal satelliter som kretsar runt jorden och skickar ut signaler [11]. I det här avsnittet beskrivs GNSS uppbyggnad och funktion samt systemen GPS, GLONASS och Galileo i närmare detalj. 2.1.1 Uppbyggnad och funktion Satelliter tillhörandes GNSS kretsar runt jorden i omloppsbana med en höjd på runt 20 000 km över havsytan. Höjden och formen på omloppsbanan varierar visserligen något från system till system [12]. Satelliterna sänder ut signaler i form av elektro- magnetiska vågor i frekvensintervallet 1,2 till 1,6 GHz, vilket är frekvenser som utgör en del av det så kallade L-bandet [11]. Frekvenserna i fråga motsvarar en våglängd mellan ungefär 19 cm och 25 cm och används då våglängder i denna storleksordning kan tränga igenom moln och regn utan större påverkan [13]. För att kompensera för jonosfärens fördröjning skickar respektive satellit dessutom ut minst två olika frekvenser. 4 2. Teori Figur 2.1: Skiss över principen för positionsbestämning, där tre satellitsignaler används för att bestämma tidsskillnader, vilka i sin tur ger det så kallade pseudo- avståendet. Större delen av den positionsbestämning som genomförs baseras på så kallad trilate- ration, där avstånd kan fastställas givet satellitkoordinater [14]. Principen framgår av Figur 2.1 där tre satelliters signaler når en mätstation vars position kan bestäm- mas. Eftersom avståndsbestämningen är given av produkten av transmissionstiden och ljushastigheten krävs hög stabilitet och noggrann synkronisering av tid [11]. Samtliga satelliter är därav utrustade med atomklockor, vars precision är mycket hög [15]. Genom att använda överlagrade vågor, signalerna med positionsinformation och tids- stämpeln för signalutsändningen, kan så kallad pseudoslumpmässig bruskod (PRN) skapas. Denna binära signal innehåller nödvändig information för att avgöra identi- teten på satelliten samt information om satellitens position och omloppsbana [11]. Genom att mottagare på jorden skapar egna kopior av PRN kan de mottagna sig- nalerna jämföras med motsvarigheten som skapats på jorden. De tidsskillnader som observeras ger, såvida de multipliceras med ljushastigheten, det så kallade pseudoav- ståndet till satelliten i fråga [16]. Mer information om hur beräkningarna realiseras finns att tillgå i Avsnitt 2.6 Beräkningsmetoder och formler. 2.1.2 Global Positioning System Global Positioning System (GPS) är ett amerikanskt system som hanteras av den amerikanska regeringen. GPS utvecklades i första hand för den amerikanska mili- tären under 1970-talet men gjordes senare tillgängligt för civilt bruk [17]. Systemet består av 31 satelliter varav minst 24 ska vara fungerande och tillgängliga 95 procent av tiden [18]. Samtliga satelliter sänder ut modulerade vågor i samma frekvens där endast PRN varierar mellan satelliterna, vilket möjliggör identifiering av dessa [11]. GPS har en konstellation av 24 satelliter i sex olika cirkulära omloppsbanor lutade 5 2. Teori 55° mot ekvatorplanet [11]. Resterande sju satelliter fungerar som reserv. Omlopps- banorna har en höjd på 20 180 kilometer över havet, som resulterar i en omloppstid på 11 timmar och 58 sekunder [11]. Denna konstellation resulterar i att fyra till åtta satelliter alltid är i siktlinje 15° ovanför horisonten oavsett position på jordklotet [19]. 2.1.3 Globalnaya Navigatsionnaya Sputnikovaya Sistem Globalnaya Nsvigatsionnaya Sputnikovaya Sistem (GLONASS) är Rysslands mot- svarighet till GPS. Den första GLONASS-satelliten kallades Uragan och sköts upp i oktober 1982. Inledningsvis utvecklades systemet för att användas inom militären, men med tiden kom det att göras tillgängligt för delar allmänheten även om det fortfarande styrs under militär kontroll [20]. Efter ytterligare uppskjutningar var sy- stemet fullt funktionellt år 1995 [11]. Till skillnad från GPS sänder alla GLONASS- satelliter ut samma PRN-modulering men på något olika frekvenser. Satelliterna kan därav identifieras givet frekvensen, vilket kan jämföras med GPS-satelliternas PRN [11]. GLONASS består av 24 satelliter jämnt fördelade i tre olika omloppsbanor med en vinkel på ungefär 65° från ekvatorialplanet [11]. Denna konstellation ger, jämfört med andra system, en ökad synlighet över Ryssland [11]. Satelliterna kretsar vidare på en lite lägre höjd, 19 100 kilometer över havet, vilket ger en kortare omloppstid på ungefär 11 timmar och 16 minuter jämfört med GPS [11]. 2.1.4 Galileo Galileo är det europeiska navigationssystemet namngett efter den välkände veten- skapsmannen Galileo Galilei [19]. Det ägs och underhålls av European Space Agency (ESA) och är, till skillnad från GPS och GLONASS, främst avsett för civilt bruk. Den första satelliten tillhörandes Galileo sköts upp 2011 och därefter har uppskjut- ningarna fortsatt kontinuerligt. Den nuvarande uppsättningen av Galileo består av 28 satelliter fördelade i tre olika omloppsbanor [21]. Alla förutom två av dessa ligger i en omloppsbana 23 222 km från jordens yta med en vinkel på 56° från ekvatorial- planet, vilket resulterar i en omloppstid på ungefär 14 timmar och 5 minuter [11]. Resterande satelliter var feluppskjutna och används enbart till Search and Rescue, vilket de andra 26 också gör men då utöver sina primära användningsområden. Se- arch and Rescue innebär att satelliterna kan vidarebefordra en nödsignal till lokala räddnings- och operationscenter [21]. Galileo använder sig av två frekvenser som standard och kan därmed ge en posi- tionsbestämning i realtid med en säkerhet ned till en meter. Liksom GPS används PRN för att identifiera och skicka signalerna [21]. 6 2. Teori 2.2 Vulkanism En vulkan är en plats med en öppning som möjliggör utströmning av gaser samt fasta och smälta bergarter från jordens inre [22]. Magma kallas den blandning av gaser och smälta mineraler som finns i manteln, delen under jordskorpan. Magma som emanerat från manteln till ovanför jordytan klassificeras som lava [23]. En rad upptäckter runt år 1965 gjorde att den geologiska teorin om att jordens yta består av flera tektoniska plattor i rörelse blev etablerad fakta. Vulkanerna är oftast belägna längs med dessa tektoniska plattors kanter, där plattorna kan verka mot varandra, en konvergent gräns, eller ifrån varandra, en divergent gräns. Av Figur 2.2 framgår hur en gräns kan se ut. Förekomst av vulkaner som är geografiskt placerade innanför kanterna kallas för hot-spots [22]. Figur 2.2: Skiss som redogör för hur plattgränsernas verkan mot varandra kan resultera i rörelse av magma, vilket i sin tur kan mynna ut i ett vulkanutbrott. Beroende på vad som orsakar ett utbrott kan vulkanerna formas på huvudsakligen två olika sätt. Stratovulkan förkommer oftast vid konvergenta plattgränser och har en spetsigare topp samt ett explosionsartat utbrott. Sköldvulkan uppstår ofta vid hot-spots och vid divergenta plattgränser, där lavan i stället tenderar att flyta ut snabbt och breda ut sig över större ytor. Från en vulkanisk eruption härstammar magma i form av lava, följt av askmoln och gasutsläpp [22]. Det är vanligt att de tektoniska plattgränserna är belägna under havsytan. Vulkanis- ka processer i form av submarina utbrott, där lava, pyroklastiska flöden och kollaps av jordytan kan orsaka uppkomsten av tsunamier [24]. Pyroklastiska material är fasta partiklar och flytande material som kastas upp ur en vulkankrater i samband med explosionsartade utbrott [25]. 7 2. Teori 2.3 Hunga Tonga–Hunga Ha’apai Hunga Tonga Hunga Ha’apai är en stratovulkan belägen vid den konvergenta platt- gränsen mellan Stilla havets och den australienska plattan. Aktivitet hos vulkanen detekterades redan 2009, varefter utbrotten kom att öka i omfattning för att nå sin kulmen den 15 januari 2022 mellan klockan 04 och 05 då flera mycket kraftiga explo- sioner inträffade [8]. Vulkanutbrottet, vilket kom att bli det mest omfattande som inträffat under 2000-talet, omfattade totalt fem explosioner som gav upphov till en våldsam spridning av stora vattenmängder och en tryckvåg, varav den sistnämnda uppmättes nå jorden runt två gånger [8]. Tack vare GNSS finns stora datamängder från den relevanta perioden tillgängliga, varför utbrottet kan analyseras i detalj. I kombination med utbrottets omfattning är fallet beaktansvärt i en mer generell diskussion, vilket är ytterligare en aspekt som motiverar ämnets relevans. 2.4 Atmosfär Atmosfären delas in i ett flertal delar beroende på temperatur. I temperaturord- ning, från högst till lägst, kallas de olika lagern för troposfär, stratosfär, mesosfär, termosfär och allt ovanför termosfär kallas för exosfär [26]. Beaktas andra storheter än temperatur, så som elektrontäthet, brukar distinktionen mellan de kyligare av tidigare nämnda lagren inte vara lika explicit varför dessa går under benämningen jonosfär. I Figur 2.3 framgår de relevanta atmosfärdelarna för projektet. Figur 2.3: Av figuren framgår hur atmosfären vanligen uppdelas, varav större delen av de tre yttre lagren, exosfären, termosfären samt mesosfären, motsvarar jonosfären. Vilken uppdelning som är att föredra beror på vilka storheter som undersöks, varav den senare uppdelningen är av intresse i detta arbete. 8 2. Teori 2.4.1 Troposfär Troposfären är den del av atmosfären som är belägen närmst jordytan, varför luft- trycket är som högst här. Då drygt 75 % av all planetens luft liksom majoriteten av all vattenånga och vattenmolekyler finns i troposfären innesluter denna del av atmosfären majoriteten av allt väder. Det finns inga finita gränser som skiljer at- mosfärdelarna åt, i stället fasas de ut och därav är den övre gränsen av troposfären ungefärlig. Troposfären sträcker sig exempelvis ungefär 7 km från jordytan vid po- lerna respektive 18 km vid ekvatorn [26]. 2.4.2 Jonosfär Den del av atmosfären som utgör jonosfären är belägen mellan drygt 48 km och 985 km över havet [26]. Detta höjdintervall innefattar den så kallade mesosfären, termosfären och exosfären, men såvida det är av intresse att studera elektroner och joner är indelningen i jonosfär att föredra. Namnet jonosfär kommer av att solens radioaktiva strålning joniserar atomer och molekyler, vilket resulterar i uppkomsten av joner och fria elektroner. Dessa rekom- binerar under nattetid, varför processen är cyklisk [26]. Vidare gäller att jonosfären på drygt 1000 km höjd kommer vara fullt joniserad, varför denna del är fullkomligt reflekterande för radiovågor under 40 MHz [27]. 2.5 Mjukvara I det här avsnittet återfinns information om de olika mjukvaror som använts i pro- jektets olika delar. Detta inkluderar både GipsyX, RINEX samt IonKit-NH. 2.5.1 RINEX Receiver INdependent EXchange (RINEX) är ett standardiserat filformat som ofta används vid GNSS observationer [28]. Standarden föreslogs år 1989 av det Astrono- miska Institutet vid universitetet i Bern, Schweiz [29]. Målet var att standardisera datan som samlades in av de olika typer av satellitsystem samt mottagare. Tidiga- re hade varje tillverkare och system sina egna filformat vilket gjorde det svårt att analysera och jämföra olika data från olika stationer och satelliter [29]. Filformatet används idag vid nästan alla GNSS observationer. RINEX finns i olika versioner och lagrar data i textformat. Filformatet finns i dagsläget i fyra versioner varav version två och tre är vanligast förekommande [29]. 2.5.2 GipsyX GipsyX är en mjukvara utvecklad av JPL, vilket utgör en del av NASA, och grun- dar sig på över 35 år av utveckling och förfining [7]. Från bland annat GNSS-data kan information förädlas och ge en uppskattning på parametrar som satellitbanor, 9 2. Teori klockor, jordens orientering, hastigheten och fördröjningar i troposfären samt jono- sfären [30]. Mjukvaran underlättar analys av data från GNSS så som exempelvis GPS, GLONASS och Galileo. Vidare används GipsyX till geofysisk forskning om jordens förändringar, isens utbredning, klimatstudier av troposfär samt jonosfär och GNSS-nätverk. Därtill är den även ett verktyg vid omloppsbestämning av satelliter och flygplanspositionering [7]. 2.5.3 IonKit-NH Den mjukvara som använts för att bearbeta den totala elektrontätheten i jonosfären heter IonKit-NH och är utvecklad av Long Tang. Den består av ett bibliotek med flertalet Matlabkoder, som finns åtkomlig på GitHub [31]. Tidigare har koden har använts i två publikationer, varav en avhandlar TEC-analys av jordbävningen i Tohoku år 2011, medan den andra innehåller motsvarande för HTHHs utbrott 2022 [32]. För att göra koden kompatibel med filhanteringen som används på Mac gjordes viss modifikation av koden, men i övrigt är mjukvaran som använts utvecklad av Tang. 2.6 Beräkningsmetoder och formler I det här avsnittet kommer formler och beräkningsmetoder som använts under ar- betets gång presenteras och förklaras. Det som diskuteras är ZTD, olika modeller och filter samt avståndsbestämning och TEC. 2.6.1 Zenith Total Delay, Zenith Wet Delay och Zenith Hyd- rostatic Delay I nedre delen av atmosfären kan signalfördröjningen betraktas som en summa av en hydrostatisk fördröjning samt en ’wet delay’ som påverkas av lufttrycket respektive fuktigheten i troposfären [33]. Summan ger den totala fördröjningen i troposfären, vilket ger information om atmosfärsinnehåll. Genom att projicera satellitsignalen på zenit i förhållande till den station som mot- tager den beaktade signalen ges ZTD. Det här måttet kan därefter delas upp i zenit hydrostatisk fördröjning (ZHD) och zenit fuktighetsfördröjning (ZWD) [33]. ZTD plottas med fördel mot tid, varefter avvikelser i grafen ger information om vad sig- nalen propagerar genom. Är propagationsmediet svårgenomträngligt ges uppmätbar fördröjning av signalen. I fallet med en tryckvåg kan avvikelser i ZTD ge en upp- fattning om tryckvågens passage. 2.6.2 Saastamoinens modell Saastamoinens formel bygger på den totala zenit fördröjningen och den tolkas som det förväntade värdet på ZTD utifrån lufttryck, höjd på markstationen och fuk- tighet. Uttryckt i andra termer beror ZTD på den hydrostatiska fördröjningen och fuktighetsfördröjningen. Den ges av följande. 10 2. Teori ZTD = ZHD + ZWD (2.1) Saastamoinens formel kan skrivas mer utförligt enligt Ekvation 2.2. ZTD = 0.0022768 1 − 0.00266 cos (2φ)(P + 1255 TC + 273.15 · 6.1078· exp ( 17.269TC TC + 237.3) · rh 100) (2.2) För att beräkna ZHD måste även ZWD bestämmas. ZWD ges enligt följande. ZWD = 0.0022768 · ( 1255 TC + 273.15) · 6.1078 · exp ( 17.269TC TC + 237.3) · rh 100 (2.3) Här φ betecknar stationens latitud i radianer, P är det atmosfäriska trycket (hPa) vid den specifika markstationen, TC är temperaturen (°C) vid den specifika mark- stationen och rh är den relativa luftfuktigheten (%) [34]. 2.6.3 Savitzky-Golay filter Savitzky-Golay filter är ett approximationsfilter som gör en signal med mycket stör- ningar till en slätare signal. Den utvecklades utav Marcel E. Golay och Abraham Savitzky på 1960-talet inom analytisk kemi och används idag flitigt inom signal- behandlning och kan tolkas som ett lågpassfilter. Filtreringen baseras på två para- metrar p samt q. Variabeln p motsvarar här polynomets ordning och q är storleken på det tidsintervall signalen ska filtreras inom. Värdet på q måste vara större än p. Därefter sätts i(m) till den brusiga signalen och f(m) är Savitzky-Golay filtret av ordning av p. Vidare ansätts polynominalkoefficienterna α0, α1, ..., αp, varvid två ekvationer erhålls. Den första ger ett mått på standardavvikelsen δ, vilket definie- rar hur filtreringen sker. I samma ekvation sker även ett skift till höger där nästa utslätning görs med respektive q tills hela signalen är uppfylld [35]. δ = q−1 2∑ m= 1−q 2 [f(m) − i(m)]2 (2.4) Nedan framgår hur f(m) definieras [35]. f(m) = p∑ j=0 ajm j (2.5) 2.6.4 Förväntad tryckvåg och avståndsbestämning Haversine-formeln är en matematisk formel som används till att beräkna avståndet mellan två punkter på en sfär, vanligtvis jordens yta. Det görs med hjälp av longitud och latitud [36]. Formeln är given av Ekvation 2.6. d = 2R · arcsin  √√√√sin2 ( la2 − la1 2 ) + cos(la1) cos(la2) sin2 ( lo2 − lo1 2 ) (2.6) 11 2. Teori där d är avståndet mellan koordinaterna i kilometer, R är klotets radie i kilometer, la1, la2 är latitudkoordinater och lo1, lo2 är longitudkoordinater. 2.6.5 Total Electron Content Den totala elektrontätheten i jonosfären, uppmätt från satellit till markstation, är given av följande formel. TEC = 1 A f1 2f2 2 f1 2 − f2 2 (L1 − L2 + N + ε) (2.7) Här är A en konstant med värdet 40,3 m3/s2, f1, f2 är de olika frekvenserna för respektive satellit, ε är mätstörningar och L1, L2 är bärvågsobservationer från re- spektive satellit. Vidare är N bärvågens osäkerhet givet i meter. Den är konstant och approximeras med hjälp av följande formel. N = ⟨L1 − L2 − (P1 − P2)⟩, (2.8) där L1, L2 är bärvågsobservationer, P1, P2 är motsvarande för PRN, vilken genereras av både satellit och markstation för att tidsskillnaden som kan utläsas då signalerna sedan jämförs kan användas för distansbestämning. Beteckningen ⟨·⟩ används för att beskriva att den genomsnittliga operationen [1]. 12 3 Metod I det här avsnittet presenteras de metoder som använts för att analysera data och erhålla de resultat som presenteras i kommande avsnitt. Här behandlas metoden för hur data införskaffades, en redogörelse för hur mjukvaran GipsyX användes samt hur det grafiska resultatet framställdes med hjälp av programmering. Slutligen be- handlas framtagandet av förväntade värden baserat på föregående rapporter samt formler. 3.1 Datainsamling För att bestämma vilka GNSS-mätstationers data som är relevanta att analysera och använda i projektet användes det grafiska verktyget Network på den av Inter- national GNSS Service ägda hemsidan www.igs.org [4]. Med verktyget lokaliserades mätstationer i närområdet kring HTHH samt ytterligare stationer med längre av- stånd till HTHH. Ett krav på samtliga stationer, bortsett från TONG00TON, var att det skulle finnas tillgängliga data från stationerna under de tre dagar projektet valt att behandla. Stationerna som behandlas i rapporten presenteras i Tabell 3.1, där även avstånden från HTHH till respektive station finns med. 13 3. Metod Tabell 3.1: Samtliga stationer vars data har använts samt respektive avstånd till HTHH. Sorterat efter avstånd till HTHH. Station Avstånd till HTHH (km) TONG00TON 70,31 USP100FJI 699,06 FTNA00WLF 751,40 LAUT00FJI 820,06 SAMO00WSM 839,51 TUVA00TUV 1 457,34 CKIS00COK 1 620,83 PTVL00VUT 1 738,51 CHTI00NZL 2 580,98 THTI00PYF 2 727,13 OUS200NZL 3 095,33 TOW200AUS 3 919,84 GAMB00PYF 4 169,45 KOKB00USA 5 042,37 SCTB00ATA 6 443,17 MIZU00JPN 8 028,64 UCLU00CAN 9 162,79 CHUR00CAN 11 463,46 TASH00UZB 13 583,69 De fullständiga stationsnamnen består av fyra tecken för namnet av stationen i fråga, två siffror och tre tecken som representerar vilket land stationen finns i. Hädanefter kommer endast de fyra första tecknen för det fullständiga stationsnamnet användas i rapporten. Av Figur 3.1 återfinns en grafisk överblick av alla GNSS-stationer som använts i projektet. I Figur 3.2 framgår en inzoomad bild på de använda stationerna som är belägna i närheten av HTHH. Den grafiska illustrationen är genererad med hjälp av Generic Mapping Tool (GMT) [37] 14 3. Metod Figur 3.1: Överblick över samtliga GNSS-markstationer som använts i rapporten. Figur 3.2: De GNSS-stationer i närheten av HTHH som använts i projektet. Fastän TONG inte har tillgängliga data från hela den undersökta perioden är statio- 15 3. Metod nen inkluderad. Det beror på dess geografiska läge i förhållande till vulkanutbrottet förser analysen med värdefull information om tryckvågens utbredning samt jonosfä- rens elektrontäthet i närheten av utbrottet. I fallet med STEC behandlades RINEX-filerna inte i GipsyX. Istället laddades de ned för att behandlas med Matlabkod. Vidare gäller att data från endast GPS, GLONASS och Galileo användes för samtli- ga av de listade stationerna. Det beror på att övriga system bedömdes vara antingen otillräckligt testade gentemot GipsyX eller hade ett format som försvårade databe- handlingen. 3.1.1 Val av tidsupplösning Data som upptagits varje sekund fanns inte att tillgå för samtliga av de valda sta- tionerna under hela den relevanta perioden. Det resulterade, i kombination med att data upptagna var 30:e sekund gav tillräckliga resultat utan att skilja sig avseende- värt med sekundsdata, i att grafer genererades enligt resultatdelen med data från var 30:e sekund. Vidare fanns enbart tillgängliga tryckdata med en upplösning av vart 30:e sekund. För att erhålla samma upplösning motiverades valet av data från var 30:e sekund ytterligare. Även för STEC-graferna användes 30 sekunders data då det möjliggjorde användning av samma RINEX-filer i fallet med ZTD samt STEC. För tester och framställning av Pythonkod användes data från stationen ONSA belägen i Onsala, Sverige, från juli 2021. Därefter kördes koden med data upptagna varje 300:e sekund från de utvalda stationerna för att säkerställa kompatibiliteten. Efter framgångsrika tester användes data tagna av stationerna var 30:e sekund, vars resultat framställdes grafiskt. Därpå användes data upptagna varje sekund under de aktuella dagarna, varvid beslutet om att enbart använda 30 sekunders data kunde fattas efter jämförelse av de erhållna resultaten. 3.1.2 Val av elevationsvinkel På grund av att atmosfären påverkar de data mätstationen insamlar gäller att den vinkel med vilken satelliten står i förhållande till mätstationen måste beaktas. Re- dan vid 20° är den väg vilken signalen färdas genom atmosfären så lång att det uppkommer ett märkbart fel i förhållande till den modell av troposfären som mjuk- varan använder [9]. I detta projekt har elevationsvinkeln 7◦ valts som gräns för att minska felet i fråga, varför vinklar som underskrider 7◦ exkluderats. Innan valet av vid vilken vinkel som exklusionen skulle börja genomfördes tester för att undersöka den skillnad valet av 7◦ kontra 10◦ resulterade i. 16 3. Metod 3.2 Användning av GipsyX Mjukvaran GipsyX användes för databearbetning. Då mjukvaran fanns nedladdad på datorn Vidar tillhörandes Onsala Rymdobservatorium skapades en remotean- slutning till denna och genom att använda ett skript utvecklat av Jan Johansson underlättades användningen av GipsyX ytterligare [3]. Beroende på de nedladdade datafilernas tidsintervall mellan upptagna data tog re- spektive kommandosekvens olika lång tid att slutföra, där 90 minuter var snittet på tidsåtgången för de filerna med data upptagna varje sekund. Såvida körningar- na skedde utan anmärkning lade programskriptet resultaten i en mapp som sedan laddades ned för vidare analys. 3.3 Grafisk framställning av Zenith Total Delay Med den kod som finns att tillgå i Avsnitt A i Appendix framställdes grafer över de erhållna resultaten för ZTD som funktion av tid. För att underlätta identifiering av vilka avvikelser som sannolikt orsakats av tryckvågen adderades markörer vid de tider explosionernas tryckvågor estimerades passera respektive station. Det gjordes utefter tidigare rapporters resultat på 310-315 m/s för tryckvågshastigheten [38], [8]. För att jämföra de erhållna resultaten med förutsägelser av ZTD plottades me- teorologiska prognoser för ZTD i de skapade graferna. Därefter kunde grafernas principiella utseende avläsas och jämföras. Eftersom det enbart fanns förutsägelser av ZTD för var sjätte timme för ett fåtal av de stationer som undersökts var jämfö- relserna ytterst grova, varför meteorologiska data inte kom att användas i projektet. Förhoppningen med de framställda graferna där ZTD var plottad som funktion av tid var att urskilja tydliga avvikelser i ZTD. Avvikelserna skulle ge en indikation om hur tryckvågen utbreddes. Eftersom samtliga dagar innehöll avvikelser av annat ursprung än tryckvågen kunde en tydlig passage av tryckvågen inte urskiljas för många stationer, utan primärt för de stationer i närheten av HTHH varför vidare analys genomfördes. Genom att använda tryckdata för att undersöka och jämföra ZTD i relation till tryck kunde även mindre ZTD-avvikelser urskiljas. Detta på grund av att tryckskill- nader och ZTD-variationer korrelerar, varför en tydlig tryckförändring indikerar att ZTD-förändras. Vidare användes Saastamoinens formel för att beräkna vilken avvikelse i ZTD som kunde förväntas. Beräkningarna, vars resultat plottades som komplement till be- handlade data, gjorde det således möjligt att avgöra huruvida det överhuvudtaget är möjligt att urskilja ZTD-förändringen bland brus eller ej, vilket underlättade analysen. 17 3. Metod Då tryckvågen inte påverkar ZWD subtraherades den från ZTD för att ge ZHD. Förhoppningen var att kunna urskilja en större skillnad och därav erhålla en tydli- gare indikation för när tryckvågen passerade. Data filtrerades även med Savitzky-Golay-filter för att ta bort högfrekventa varia- tioner och brus. Förhoppningen var att erhålla en tydligare och mindre brusig graf. 3.4 Beräkning av Total Electron Content Användingen av Tangs mjukvara IonKit-NH implementerades genom att tillföra den samma obehandlade RINEX-filer som tillförts GipsyX. Utdata bestod av en fil inne- hållandes avvikelsen i STEC för stationen, delta Slant TEC (dSTEC) [39]. Filerna användes sedan vid framställning av en grafisk illustration, där dSTEC, angivet i STECU (Slant Total Electron Content Unit), plottades som funktion av tid. Markstationerna som ingår i analysen är totalt 12 stycken, som vidare genererar 12 dSTEC filer. I respektive dSTEC-fil återfanns data för flertalet satelliter, mellan 18-75 stycken, i relation till markstationen vilket framgår av Avsnitt B.2 i Appendix. Genom data från samtliga stationer kunde en omfattande graf framställas, där avvi- kelser i STEC angavs med hjälp av en färgskala. Matlabkoden för framställningen av grafen finns att tillgå i Avsnitt B.1 i Appendix. Ur grafen framgick hastigheten hos vågen, och genom att testa olika lutning på infogade linjer kunde en approximativ hastighet härledas. Därefter analyserades alla samtliga dSTEC värden för alla stationer i ytterligare Matlabkod som återfinns i Avsnitt B.3 i Appendix. Det gjordes för att hitta da- ta vid rätt avstånd och tid, som tydligt kan visa en skillnad i dSTEC. Eftersom satelliterna är i ständig rörelse finns ett intervall så att data endast tas från 100- 200 km från HTHH. Data valdes här från ett representativt urval satelliter, vars signalskärningspunkter med jonosfären var belägna på olika avstånd från HTHH. Grafisk framställning av dess värden i en gemensam figur gav tydligare och mer lät- tolkade resultat. Koden som användes för framställningen finns att tillgå i Avsnitt B.4 i Appendix. Ur figuren från graferna kunde även en approximativ hastighet för tryckvågen bestämmas, eftersom tiden för passage samt avståndet till HTHH var givna. 3.5 Analys av tryckvågens utbredning För att bestämma en utbredning av tryckvågen med hjälp av STEC har data för samtliga stationer framställts grafiskt. Avståndet från HTHH är här givet som funk- tion av tid. Genom att låta tryckvågens avvikelse vara rödmarkerad, medan ingen avvikelse framgår av gröna markeringar, har en rät linje vars lutning svarar mot utbredningshastigheten kunnat framställas. Av grafen framgår flera räta linjer, där respektive linje representerar de olika explosionerna. 18 3. Metod 3.6 Förväntade värden Koordinaterna för varje station har hämtats från International GNSS Service, vil- ket även är källan satellitdatafilerna hämtats från [4]. Med hjälp av Python-kod, se Avsnitt A.2 i Appendix, beräknades avståndet mellan varje station och HTHH vulkanen med Haversine-formeln. Avståndet användes vidare för att skapa en förväntad tid för när tryckvågen bör ha nått respektive station. Den estimerade utbredningshastigheten hos tryckvågen är 310-315 m/s, vilket baserat på tidigare studier av HTHH-utbrottet [38], [8]. Med utgångspunkt i att vulkanutbrottet bestod av fem explosioner har fem markerade vertikala intervall lagts till i graferna, vilket ger en estimerad tid för när tryckvågen borde nått respektive mätstation [40]. Saastamoinens formel har använts för att beräkna ZTD med hjälp av tryckdata, varefter resultatet har förts in i samma graf som uppmätt ZTD och använts som en riktlinje för vilken avvikelse i ZTD bör förväntas. Genom att beräkna tryckavvikelsen för respektive estimerad tryckvåg erhålls ett preliminärt värde på avvikelsen i ZTD. Då tryckavvikelserna minskar med ökande avstånd från HTHH är det med denna metod möjligt att avgöra huruvida förändringen av ZTD är urskiljbar bland brus. Detta kan i sin tur användas för att avgöra om det är nödvändigt att inkludera tryckdata för att finna vilka trender i grafen som motsvarar tryckvågens påverkan på signalfördröjningen. 19 4 Resultat Erhållet resultat presenteras i denna del fördelat på fyra delar; resulterande ZTD för stationerna belägna närmast utbrottsplatsen, jämförelse mellan tryckdata samt ZTD för avlägsna stationer, eliminering av ZWD och jonosfärens totala elektrontäthet i vertikalled under den undersökta tidperioden. 4.1 Uppmätt Zenith Total Delay De stationer vars ZTD undersökts är de sju stationer som ligger närmast HTHH. Stationerna har delats in i tre större avståndsspann, ca 70 km, 699 − 850 km och 1450 − 1650 km. Explosionerna är markerade med orange streck i figurerna nedan. I det första intervallet återfinns stationen TONG, vilken ligger 70 km från HTHH. I grafen som återger ZTD, se Figur 4.1, för stationen i fråga framgår att ZTD bå- de ökar och minskar under de estimerade intervallen för tryckvågen. Inledningsvis minskar den, för att sedan öka kraftigt efter andra utbrottet. Därefter en minskning, varpå en markant ökning igen. LP-filtrerad kurvan, är den ZTD-data som inhämtas där en Savitzky-Golay filtrering implementeras där alla mindre avvikelser filtrerar bort och kurvan blir slätare. Figur 4.1: ZTD [m] vid stationen TONG, ca 70 km från utbrottet. 20 4. Resultat I andra intervallet, 699-850 km, finns stationerna USP1, FTNA, LAUT och SAMO, presenterat i Figur 4.2. För både USP1 och FTNA sker en minskning varpå en ök- ning vid de första två ankomsterna av estimerad tryckvåg, samt ökning av ZTD vid resterade intervall. I tidsintervallet 04:30 till 07:00 uppmäts ett maximalt värde på ZTD vid 07-tiden för USP1 och vid 06.30-tiden för FTNA. För stationerna LAUT och SAMO sker det en knappt märkbar ökning vid de två första utbrotten. För stationerna längre bort än 800 km är det svårt att urskilja en förändring i ZTD på grund av tryckvågen, då ändringen inte är större än brus och mindre variationer i ZTD. I det tredje avståndsspannet 1455-1620 km från utbrottet finns stationerna TUVA och CKIS, se Figur 4.2. Gemensamt för dessa är att det inte förekommer några större variationer jämfört med brus. 21 4. Resultat (a) USP1, 699 km (b) FTNA, 751 km (c) LAUT, 820 km (d) SAMO, 839 km (e) TUVA, 1457 km (f) CKIS, 1620 km Figur 4.2: ZTD [m] för närliggande stationer. Det givna avståndet under varje figur är avståndet från vulkanutbrottet till stationen. 22 4. Resultat 4.2 Uppmätt Zenith Total Delay samt lufttryck I det här avsnittet presenteras ZTD tillsammans med lufttryck för stationerna TONG, MIZU, OUS2 samt TASH. Den ljusblå grafen motsvarar barometriskt tryck i hPa, medan den mörkgröna svarar mot ZTD. I graferna har även orangea, strec- kade linjer i vertikalled infogats vid de tidpunkter tryckvågorna estimerats anlända. Av Figur 4.3a framgår ZTD samt tryck för stationen TONG. Notera hur båda stor- heter ökar kraftigt med start vid tiden för den första anländande tryckvågen, där någon märkbar fördröjning inte tycks föreligga. För stationerna på längre avstånd från utbrottsplatsen, se Figur 4.3b, 4.3c, 4.3d, kan ingen tydlig trend i förändring av ZTD observeras. 03:00 04:00 05:00 06:00 07:00 08:00 09:00 10:00 Tid 985 990 995 1000 1005 1010 Ba ro m et ris kt tr yc k P (h Pa ) 2.54 2.56 2.58 2.60 2.62 2.64 2.66 ZT D TONG (a) TONG, 70 km 06:00 07:00 08:00 09:00 10:00 11:00 Tid 1015.0 1015.5 1016.0 1016.5 1017.0 1017.5 1018.0 1018.5 1019.0 Ba ro m et ris kt tr yc k P (h Pa ) 2.405 2.410 2.415 2.420 2.425 2.430 2.435 ZT D OUS2 (b) OUS2, 3 095 km 10:00 11:00 12:00 13:00 14:00 15:00 16:00 Tid 1006 1007 1008 1009 1010 1011 1012 Ba ro m et ris kt tr yc k P (h Pa ) 2.32 2.33 2.34 2.35 2.36 ZT D MIZU (c) MIZU, 8 028 km 16:00 16:30 17:00 17:30 18:00 Tid 947.5 948.0 948.5 949.0 949.5 950.0 Ba ro m et ris kt tr yc k P (h Pa ) 2.2400 2.2425 2.2450 2.2475 2.2500 2.2525 2.2550 2.2575 ZT D TASH (d) TASH, 13 584 km Figur 4.3: ZTD [m] och lufttryck [hPa] för stationer längre från vulkanutbrottet. Avståndet under varje figur är distansen från HTHH till stationen. 23 4. Resultat 4.3 Eliminering av Zenith Wet Delay Detta avsnitt innehåller de grafer som erhållits då ZWD subtraheras från ZTD, var- efter ZHD erhållits. Beräkningarna är gjorda på stationerna TONG, OUS2, MIZU samt TASH och resultatet är sammanställt i Figur 4.4. Av graferna framgår även resultat beräknade med Saastamoinens formel, de lågpassfiltrerade resultaten samt de tider tryckvågorna för respektive explosion förväntas nå stationen. Av graferna för TONG, se Figur 4.4a, framgår att trenderna för ZTD och ZHD är relativt lika. Både ZTD och ZHD varierar därav på liknande, förväntade sätt vid de estimerade tiderna för när tryckvågorna för respektive station ska anlända. Grafen för OUS2, Figur 4.4b, redogör för hur variationen i ZTD samt ZHD båda varierar mycket litet. Notera att OUS2 är belägen 3095 km från HTHH, men att det däremellan förekommer både land samt högre berg. Motsvarande grafer för MIZU, se Figur 4.4c, är desto mer varierande. Värt att notera är att MIZU är belägen 8028∼km från HTHH men att ingen landmassa är belägen mellan utbrottsplatsen och mätstationen. Av Figur 4.4d framgår hur avvikelserna för både ZTD och ZHD inte tycks vara stora nog för att enstaka avvikelser ska gå att urskilja. 24 4. Resultat (a) TONG, 70 km (b) OUS2, 3 095 km (c) MIZU, 8 028 km (d) TASH, 13 584 km Figur 4.4: ZTD [m] och ZHD [m] för stationer längre från vulkanutbrottet. Av- ståndet under varje figur är distansen från HTHH till stationen. 4.4 Total Electron Content för 12 stationer Den totala elektrontätheten som redovisas har bestämts med data från 12 GNSS- markstationer, vardera med tre stationer i olika riktningar. Efter bearbetning av RINEX filer i mjukvaran IonKit-NH, erhålls 12 filer med information om varje sa- tellits avvikelse i STEC, dSTEC. Notera att avvikelsen inte är projicerad i vertikalled utan i riktning mot satelliten, varför dSTEC erhållits. De stationer som använts och hur många unika satelliter som förekommer i respektive dSTEC-fil presenteras i Tabell 4.1. 25 4. Resultat Tabell 4.1: GNSS markstationer som använts, samt hur många satelliter respektive station har dSTEC-data på. Station Antal satelliter TONG 18 CKIS 21 THTI 34 CHUR 50 GAMB 60 PTVL 61 CHTI 67 TOW2 69 UCLU 70 KOKB 70 LAUT 71 SCTB 75 I Figur 4.5 framgår satelliternas banor i förhållande till distans och med tillhörande avvikelser i STEC. Färgvariationerna tycks bilda ett mönster. Notera att tre linjer svarandes mot hastigheterna 180 m/s, 250 m/s respektive 310 m/s har infogats i figuren. Samtliga linjer startar vid kl. 04:25, vilket antas vara tiden för när den största explosionen tog vid [40]. Figur 4.5: Karta över satellitbanornas avvikelser illustrerat i dSTECU. Notera linjerna och de mönster som framgår av färgskalan. 26 4. Resultat 4.5 Total Electron Content för enskilda satelliter Här visas dSTEC, för ett representativt urval av stationer på olika avstånd från HTHH. Avstånden har beräknats utifrån den longitud och latitud för signalens skär- ningspunkt med jonosfären. Av Figur 4.6 framgår tre grafer för uppmätt avvikelse i STEC, vilket är STEC som inte projicerats i vertikal-led utan taget direkt i satellitens riktning. Avvikelsen, dSTECU, varierar kring ett normaliserat värde som här satts till noll. Av figuren framgår även medelvärdet av avståndet från skärningspunkt till utbrottsplats under perioden. För den första grafen, hörandes till station-satellit kombinationen LAUT-81, där 81 motsvarar R21 (GLONASS satellit 21), ges som tydligast avvikelser efter kl. 05 och under en dryg halvtimme därefter, där den största avvikelsen erhålls vid ungefär 05:45. Detta kan svara mot hastigheter så låga som 180 m/s upp mot 500 m/s. Motsvarande värden för grafen CHTI-123, där 123 motsvarar E03 (Galileo satellit tre), respektive KOKB-145, där 145 motsvarar E25 (Galileo satellit 25), har beräknats uppgå till drygt 300 m/s samt 280 m/s då avvikelserna ges vid kl. 06:47 respektive 09:31. Figur 4.6: dSTEC för tre satelliter belägna vid olika avstånd från HTHH. Avståndet är angivet i kilometer och dSTEC uppmätt i dSTECU. 27 5 Diskussion I det här avsnittet tolkas och diskuteras det resultat som presenterats i föregående avsnitt. Exempelvis kommer ZTD att diskuteras med hjälp av tryckdata och Saas- tamoinens formel, varefter STEC analyseras och diskuteras. Vidare avhandlas de felkällor som kan ha påverkat det resultat som erhållits. 5.1 Zenith Total Delay Resultatet talar för att ZTD går att använda för att undersöka utbrottet av HTHH, men att det primärt fungerar för de mätstationer som är belägna nära utbrottsplat- sen. Av observationen framgår hur förändringen av ZTD vid den estimerade tiden för tryckvågens ankomst tycks variera som kraftigast för TONG, medan motsvaran- de resultat för övriga stationer är mindre explicit. Som exempel varierar ZTD för CKIS, se Figur 4.2f, mellan 2,43 m och 2,45 m medan motsvarande värden är cirka 2,56 m och 2,66 m för TONG. Att variationen är mindre för mer avlägsna stationer är ett förväntat resultat då tryckvariationerna i atmosfären blir mindre med ökande avstånd från utbrottsplat- sen, varför detsamma gäller ZTD-variationen. Det resultatet erhölls tämligen tidigt under projektets gång, varpå metoden kom att utökas med analys av tryckdata. 5.1.1 Zenith Total Delay i kombination med tryckdata De avståndsbegränsningar som präglar ZTD framgår även av figurerna i Avsnitt 4.2. Figurerna redogör för hur trenderna i framtagna tryckdata i varierande grad efterföljs av ZTD, tämligen omedelbart. Vid jämförelse av stationerna vars tryck- data och ZTD plottats går det att observera att ju längre en station är belägen från utbrottsplatsen, desto mindre uttalade är likheterna mellan trenderna för tryck respektive ZTD. Tryckdata och ZTD för stationen TONG, vilket framgår i Figur 4.3a, uppvisar en mycket likartad förändring strax innan, under och efter de estimerade tiderna för vågens ankomst. Motsvarande mönster för OUS2 som är belägen 3095 km från vulka- nen är avsevärt mindre utpräglat. Värt att notera i det här fallet är att motsvarande graf för MIZU visar på ett resultat där trenderna återigen är tydliga. Genom att betrakta mätstationernas geografiska lägen framgår att OUS2 först nås av tryck- vågen efter att den passerat större landmassor, däribland bergskedjor som skärmar 28 5. Diskussion mätstationen från vågen. Att avvikelsen i ZTD är tydligare för MIZU är således väntat, eftersom stationen nås av tryckvågen efter att den passerat öppet hav. Vid beräkning av förändringen i ZTD med Saastamoinens formel framgår att resul- tatet är förväntat. Det beror på den tryckförändring den annalkande vågen kommer orsaka resulterar i en ZTD-förändring som inte går att särskilja från de variationer som förekommer på grund av vädervariationer i troposfären, ZWD-delen av ZTD, liksom brus. För att eliminera påverkan av ZWD inhämtades meteorologiska data. Att subtra- hera dessa data från ZTD-data tycktes däremot inte göra trenderna tydligare. Som exempel kvarstod de många brusliknande förändringarna, varför arbetet med att urskilja tryckvågens passage inte förenklades. Därav undersöktes möjligheterna att använda meteodata inte ytterligare. 5.2 12 mätstationers uppmätta elektrontäthet Med färgskala har 12 mätstationers uppmätta dSTEC framställts grafiskt som funk- tion av både tid och avstånd från utbrottets epicentrum. Av resultaten framgår att ett antal linjelikannde former kan observeras ur färgskalan. Tre linjer med olika lut- ning 310 m/s, 250 m/s och 180 m/s har lagts till grafen, vilket var olika hastigheter som genom prövning visades överensstämma bäst med de erhållna STEC-linjerna. Då tidigare studier erhållit ett motsvarande värde på drygt 300-320 m/s bedöm- des resultatet rimligt. Däremot skulle även hastigheter så låga som 180 m/s kunna bedömas giltiga, varför resultaten är svårtolkade. 5.3 Vinklade totala elektrontäthet för respektive station Genom att bestämma den vinklade totala elektrontätheten för ett flertal satelliter med olika avstånd för skärningsplatsen med jonosfären kunde STEC användas för analys av tryckvågen. Av Figur 4.6 framgår ett representativt urval av de grafer som erhölls. Med ett estimerat värde för tidsintervallet då vågen passerade respekti- ve punkt kunde ett approximativt värde för dess hastighet bestämmas. Genom att iterera hastighetsbestämningen för ett flertal punkter gavs en approximativ hastig- het inom intervallen cirka 180-500 m/s. Det stämmer delvis överens med de linjer som approximerades efter de mönster som kunde observeras i dSTEC-avvikelsen. Troligen är hastigheterna över 300 m/s feluppmätta på grund av att det var syn- nerligen svårt att finna ett tillvägagångssätt där avvikelsen för rätt exposion kunde urskiljas. Vidare framgår av Figur 4.6 hur kartläggning av tryckvågen även är möjligt på stör- re avstånd jämfört med ZTD. dSTEC förefaller således vara att föredra för analys av tryckvågen vid avlägsna stationer. Däremot framgår inte avvikelser i dSTEC or- sakade av passage över landmassor eller skärmning av bergskedjor lika explicit som 29 5. Diskussion för ZTD. Att undersöka dSTEC och hur satellitdata skiljer sig för de fyra väderstrecken i för- hållande till en specifik mätstation hade med fördel kunnat göras i framtida studier. Resultaten i fråga hade då kunnat ge ytterligare information om hur tryckvågens utbredning varierar över avstånd och i olika riktningar. 5.4 Datainsamling och felkällor Enbart data från den 15:e januari är inkluderade i de slutgiltiga resultaten. Att ha exluderat ZTD för den 14:e och 16:e januari bedömdes berättigat då den tryckför- änding den andra passagen av tryckvågen orsakade inte gav en förändring av ZTD stor nog att urskilja bland brus eller lokala förändringar av ZWD. Värt att notera är att mer väl upplösta data troligen inte hade resulterat i ett tyd- ligare resultat. I stället erhölls endast fler variationer, varför analysen försvårades. Då tidigare studier valt att använda data upptagen var 30:e sekund samt eftersom tillgängliga tryckdata hade en upplösning svarandes mot var 30:e sekund valdes den upplösningen. Valet kan visserligen bidragit till att snabba variationer inte togs i beaktning. Emellertid bedöms felkällan i fråga inte ha kunnat decimera erhållet re- sultat avsevärt eftersom det primärt är ihållande trender som är av intresse. Det faktum att tryckvågens hastighet estimerades uppgå till 310 m/s är ytterligare en felkälla att ta i beaktning. Även om tidigare studier, likt omnämnt i teorin, funnit att storleken av hastigheten uppgått till approximativt 310 m/s är det sannolikt att den varierar mellan de fyra väderstrecken, likväl som på olika avstånd från utbrotts- platsen. Att 310 m/s må vara en inkorrekt approximation av utbredningshastigheten i samtliga riktningar innebär att de vertikala linjer som infogats i samtliga grafer kan vara missvisande. Vidare kan hastigheten i troposfären och jonosfären vara olika, varför hastigheten, vilka studierna bestämt i jonosfären, eventuellt är missvisande. Att variationen i tryckvågens hastighet i jonosfären är förhållandevis stor tros be- ro på två faktorer. Dels har en elevationsvinkel på 7° använts, vilket tros ha givit upphov till att atmosfärsinnehållet genom vilket signalen passerat ökat mätosäker- heten. Faktorer likt luftfuktighet och lokala tryckvariationer skulle kunna resultera i en stor mätosäkerhet som kan ha decimerat de erhållna resultaten validitet. Att använda en större elevationsvinkel i framtida studier hade därav kunnat ge mindre mätosäkerhet och en mer precis hastighetsbestämning. Det hade visserligen skett på bekostnad av mindre datamängder, varför en avvägning måste göras. Därtill gäller att metoden som användes vid hastighetsbestämningen i jonosfären var approximativ. Att bedöma vilken avvikelse som korresponderade mot rätt explosion var synnerligen svårt och att fel uppstått under denna procedur är således ytterst sannolikt. Inför framtida studier bör en bättre metod utvecklas eftersom precisio- nen då förbättras avsevärt. I denna studie har hastighetsbestämningen därav främst genomförts med ZTD och att i framtida studier jämföra hastigheten i troposfären 30 5. Diskussion med motsvarande i jonosfären vore följaktligen av intresse. Vidare har mätfel för klockor samt mätstationers signalmottagning försummats. Att fel som antagits vara försumbara haft en signifikant påverkan på det erhållna resultatet är osannolikt, men bör tas i beaktning fastän eliminering av felen är svårgenomförbart. 31 6 Slutsats I den här rapporten har möjligheterna att undersöka utbredningen av en tryckvåg med GNSS-data avhandlats. Bearbetade GNSS-data har hämtats från mätstationer geografiskt utspridda kring HTHH, som är vulkanen vars utbrott har undersökts. För bearbetningen användes dels med mjukvaran GipsyX, dels Matlabkod utvecklad av Long Tang. Tillsammans gav dessa analyserbara data för jonosfärens elektron- täthet samt signalfördröjningen i troposfären. De resultat som erhållits påvisar användbarheten av att kombinera både ZTD- samt STEC-data vid analys av tryckvågsutbredning. Erhållna resultat skulle således kun- na bidra till att underlätta valet av att använda ZTD eller STEC vid undersökning av vulkanutbrott, eftersom rapporten påvisat hur de är användbara i olika avseen- den. STEC är att föredra såvida ett vulkanutbrott ska studeras med data upptagna på ett stort avstånd från utbrottet. Fördelarna kommer av att de tryckvariationer tryckvågen från HTHH orsakade inte var stora nog att ge upphov till en märkbar förändring av ZTD på stora avstånd från utbrottsplatsen. Hade vågen dessutom passerat landmassa innan mätstationen nåddes var variationen ännu svårare att ur- skilja bland brus. Det var inte fallet för STEC eftersom även små tryckvariationer gav märkbara avvikelser, varför STEC var att föredra på större avstånd eller för stationer som skärmas av landmassa. Däremot bedömdes ZTD vara användbart för kartläggning av tryckvågens utbred- ning i närheten av utbrottsplatsen och för hastighetsbestämning med hög precision. Hastigheten hos den tryckvåg HTHH orsakade med ZTD-data kunde fastställas va- ra drygt 310 m/s. Vid analys av STEC-data erhölls ett tämligen stort intervall för vilken hastighet vågen propagerade med. Med kombinationen av ZTD- samt STEC- analys kunde tryckvågens utbredning kartläggas med hjälp av GNSS, varför syfte och frågeställningar bedöms vara uppfyllda och besvarade. 32 Referenser [1] L. Tang, ”IonKit-NH: a MATLAB-based toolkit for ionospheric detection of natural hazard IonKit-NH: a MATLAB-based toolkit for ionospheric detection of,” 2023, [Online]. doi: 10.21203/rs.3.rs-2308071/v2. [2] Chalmers Tekniska Högskola, ”Rüdiger Haas: Research,” [Online]. Tillgänglig: https://research.chalmers.se/person/rudiger (hämtad 2024-05-07). [3] Chalmers Tekniska Högskola, ”Jan Johansson: Research,” [Online]. Tillgäng- lig: https://research.chalmers.se/person/jmj (hämtad 2024-05-07). [4] International GNSS Service, ”Products,” [Online]. Tillgänglig: https://igs. org/products/ (hämtad 2024-02-09). [5] E. Astafyeva, ”Ionospheric Detection of Natural Hazards,” 4, dec. 2019, s. 1265– 1288. doi: 10.1029/2019RG000668. [6] S. Madry, ”Top Ten Things to Know About GNSS,” 2015, s. 99–100. doi: 10.1007/978-1-4939-2608-4_8. [7] National Aeronautics and Space Administration, ”GipsyX,” [Online]. Tillgäng- lig: https://gipsyx.jpl.nasa.gov/ (hämtad 2024-02-19). [8] E. Aa, S. R. Zhang, P. J. Erickson m. fl., ”Significant Ionospheric Hole and Equatorial Plasma Bubbles After the 2022 Tonga Volcano Eruption,” Space Weather, årg. 20, nr 7, juli 2022, [Online], issn: 15427390. doi: 10.1029/ 2022SW003101. (hämtad 2024-02-09). [9] H. Donovan, ”Reduction of the Minimum Elevation Angle for NASA Satel- lite Laser Ranging Tracking Operations,” Lanham, april 2001, [Online]. Till- gänglig: https://ilrs.cddis.eosdis.nasa.gov/lw12/docs/Donovan_ Reduction%20in%20the%20Minimum%20Elevation.pdf (hämtad 2024-02-19). [10] European Space Agency for the Space Program, ”Whats GNSS?” [Online]. Tillgänglig: https://www.euspa.europa.eu/european-space/eu-space- programme/what-gnss (hämtad 2024-02-09). [11] P. J. Teunissen och O. Montenbruck, ”Springer Handbook oƒ Global Naviga- tion Satellite Systems,” 2017. doi: 10.1007/978-3-319-42928-1. [12] National Aeronautics and Space Administration, ”CDDIS | | Techniques | GNSS | GNSS Overview,” [Online]. Tillgänglig: https://cddis.nasa.gov/ Techniques/GNSS/GNSS_Overview.html (hämtad 2024-02-17). 33 https://doi.org/10.21203/rs.3.rs-2308071/v2 https://research.chalmers.se/person/rudiger https://research.chalmers.se/person/jmj https://igs.org/products/ https://igs.org/products/ https://doi.org/10.1029/2019RG000668 https://doi.org/10.1007/978-1-4939-2608-4_8 https://gipsyx.jpl.nasa.gov/ https://doi.org/10.1029/2022SW003101 https://doi.org/10.1029/2022SW003101 https://ilrs.cddis.eosdis.nasa.gov/lw12/docs/Donovan_Reduction%20in%20the%20Minimum%20Elevation.pdf https://ilrs.cddis.eosdis.nasa.gov/lw12/docs/Donovan_Reduction%20in%20the%20Minimum%20Elevation.pdf https://www.euspa.europa.eu/european-space/eu-space-programme/what-gnss https://www.euspa.europa.eu/european-space/eu-space-programme/what-gnss https://doi.org/10.1007/978-3-319-42928-1 https://cddis.nasa.gov/Techniques/GNSS/GNSS_Overview.html https://cddis.nasa.gov/Techniques/GNSS/GNSS_Overview.html Referenser [13] R. Balasubramaniam och C. Ruf, ”Characterization of rain impact on L-Band GNSS-R ocean surface measurements,” mars 2020, s. 111 607. doi: 10.1016/ J.RSE.2019.111607. [14] Britannica Academic, ”Geoid,” [Online]. Tillgänglig: https://academic-eb- com.eu1.proxy.openathens.net/levels/collegiate/article/geoid/ 36465 (hämtad 2024-04-09). [15] European Space Agency, ”An intuitive approach to the GNSS positioning,” [Online]. Tillgänglig: https : / / gssc . esa . int / navipedia / index . php ? title = An _ intuitive _ approach _ to _ the _ GNSS _ positioning (hämtad 2024-02-17). [16] European Space Agency, ”GNSS signal,” [Online]. Tillgänglig: https://gssc. esa.int/navipedia/index.php?title=GNSS_signal (hämtad 2024-04-09). [17] Lantmäteriet, ”GPS,” [Online]. Tillgänglig: https://www.lantmateriet.se/ sv/geodata/gps-geodesi-och-swepos/GPS-och-satellitpositionering/ GPS-och-andra-GNSS/GPS/ (hämtad 2024-03-03). [18] GPS: The Global Positioning System, ”Space Segment,” [Online]. Tillgänglig: https://www.gps.gov/systems/gps/space/ (hämtad 2024-02-09). [19] J. Awange, ”GNSS Environmental Sensing,” 2018. doi: 10.1007/978-3-319- 58418-8. Tillgänglig: http://www.springer.com/series/3234. [20] T. S. Logsdon, ”GPS,” Britannica Academic, [Online]. Tillgänglig: https: //academic-eb-com.eu1.proxy.openathens.net/levels/collegiate/ article/GPS/396001 (hämtad 2024-04-15). [21] European Space Agency, ”ESA - What is Galileo?” [Online]. Tillgänglig: https: //www.esa.int/Applications/Navigation/Galileo/What_is_Galileo (hämtad 2024-03-24). [22] F. Dobran, Volcanic Processes: mechanisms in material transport, 1. utg. New York: Kluwer Academic/Plenum Publishers, 2001, vol. 590, s. 1–590. (hämtad 2024-02-12). [23] Britannica Academic, ”Magma,” [Online]. Tillgänglig: https://academic- eb-com.eu1.proxy.openathens.net/levels/collegiate/article/magma/ 50002 (hämtad 2024-04-26). [24] S. J. Day, ”Volcanic Tsunamis,” i The Encyclopedia of Volcanoes, [Online], Elsevier, 2015, s. 993–1009. doi: 10.1016/B978-0-12-385938-9.00058-4. (hämtad 2024-03-19). [25] Nationalencyklopedin, ”Pyroklastiskt material,” [Online]. Tillgänglig: http: //www.ne.se/uppslagsverk/encyklopedi/l%C3%A5ng/pyroklastiskt- material (hämtad 2024-03-25). [26] National Institute of Water and Atmospheric Research, ”Layers of the at- mosphere,” [Online]. Tillgänglig: https://niwa.co.nz/education- and- training/schools/students/layers (hämtad 2024-02-14). [27] K.G. Budden, The propagation of radio waves, 1. utg. Cambridge University Press, 1985, s. 11–76, isbn: 0-521-25461-2. 34 https://doi.org/10.1016/J.RSE.2019.111607 https://doi.org/10.1016/J.RSE.2019.111607 https://academic-eb-com.eu1.proxy.openathens.net/levels/collegiate/article/geoid/36465 https://academic-eb-com.eu1.proxy.openathens.net/levels/collegiate/article/geoid/36465 https://academic-eb-com.eu1.proxy.openathens.net/levels/collegiate/article/geoid/36465 https://gssc.esa.int/navipedia/index.php?title=An_intuitive_approach_to_the_GNSS_positioning https://gssc.esa.int/navipedia/index.php?title=An_intuitive_approach_to_the_GNSS_positioning https://gssc.esa.int/navipedia/index.php?title=GNSS_signal https://gssc.esa.int/navipedia/index.php?title=GNSS_signal https://www.lantmateriet.se/sv/geodata/gps-geodesi-och-swepos/GPS-och-satellitpositionering/GPS-och-andra-GNSS/GPS/ https://www.lantmateriet.se/sv/geodata/gps-geodesi-och-swepos/GPS-och-satellitpositionering/GPS-och-andra-GNSS/GPS/ https://www.lantmateriet.se/sv/geodata/gps-geodesi-och-swepos/GPS-och-satellitpositionering/GPS-och-andra-GNSS/GPS/ https://www.gps.gov/systems/gps/space/ https://doi.org/10.1007/978-3-319-58418-8 https://doi.org/10.1007/978-3-319-58418-8 http://www.springer.com/series/3234 https://academic-eb-com.eu1.proxy.openathens.net/levels/collegiate/article/GPS/396001 https://academic-eb-com.eu1.proxy.openathens.net/levels/collegiate/article/GPS/396001 https://academic-eb-com.eu1.proxy.openathens.net/levels/collegiate/article/GPS/396001 https://www.esa.int/Applications/Navigation/Galileo/What_is_Galileo https://www.esa.int/Applications/Navigation/Galileo/What_is_Galileo https://academic-eb-com.eu1.proxy.openathens.net/levels/collegiate/article/magma/50002 https://academic-eb-com.eu1.proxy.openathens.net/levels/collegiate/article/magma/50002 https://academic-eb-com.eu1.proxy.openathens.net/levels/collegiate/article/magma/50002 https://doi.org/10.1016/B978-0-12-385938-9.00058-4 http://www.ne.se/uppslagsverk/encyklopedi/l%C3%A5ng/pyroklastiskt-material http://www.ne.se/uppslagsverk/encyklopedi/l%C3%A5ng/pyroklastiskt-material http://www.ne.se/uppslagsverk/encyklopedi/l%C3%A5ng/pyroklastiskt-material https://niwa.co.nz/education-and-training/schools/students/layers https://niwa.co.nz/education-and-training/schools/students/layers Referenser [28] International GNSS Service, ”Formats and Standards,” [Online]. Tillgänglig: https://igs.org/formats-and-standards/ (hämtad 2024-05-07). [29] M. Chen, ”What is RINEX format and its latest version?” [Online]. Tillgänglig: https://igs.org/formats-and-standards/ (hämtad 2024-05-07). [30] W. Bertiger, Y. Bar-Sever, A. Dorsey m. fl., ”GipsyX/RTGx, a new tool set for space geodetic operations and research,” Advances in Space Research, årg. 66, s. 469–489, 2020, [Online]. doi: 10.1016/j.asr.2020.04.015. Tillgänglig: http://creativecommons.org/licenses/by/4.0/ (hämtad 2025-02-26). [31] L. Tang, ”IonKit-NH,” [Online]. Tillgänglig: https://github.com/tanggdut/ IonKit-NH/tree/main (hämtad 2024-04-15). [32] L. Tang, Z. Li och B. Zhou, ”Large-area tsunami signatures in ionosphere ob- served by GPS TEC after the 2011 Tohoku earthquake,” s. 93, 2018, [Online]. Tillgänglig: https://doi.org/10.1007/s10291-018-0759-1. [33] P. Benevides, J. Catalão, P. Miranda och M. J. Chinita, ”Analysis of the relation between GPS tropospheric delay and intense precipitation,” 2013, [Online]. doi: 10.1117/12.2028732. [34] A. Sam-Khaniani och R. Naeijian, ”Evaluation of modified Saastamoinen ZTD model using ground-based GPS observation over Iran,” Earth Science Infor- matics, årg. 16, nr 3, s. 2339–2353, sept. 2023, [Online], issn: 18650481. doi: 10.1007/S12145-023-01033-4/TABLES/4. [35] S. Kumar Sunori Mehul Manu Amit Mittal Pradeep Juneja, ”Savitzky-Golay Filter Design for Removing Noise from EEG Signal,” [Online]. doi: 10.1109/ ICIDCA56705.2023.10100215. [36] H. Alkan och H. Celebi, ”The Implementation of Positioning System with Trilateration of Haversine Distance,” 2019, [Online], s. 1–6. doi: 10.1109/ PIMRC.2019.8904289. [37] D. Tian, L. Uieda, W. J. Leong m. fl., PyGMT: A Python interface for the Generic Mapping Tools, version 0.12.0, [Online]. doi: 10 . 5281 / zenodo . 11062720. Tillgänglig: https://doi.org/10.5281/zenodo.11062720. [38] T. Tonegawa och Y. Fukao, ”Mesospheric pressure source from the 2022 Hunga, Tonga eruption excites 3.6-mHz air-sea coupled waves,” 2023, [Online]. doi: 10.1126/sciadv.adg8036. [39] L. Tang, ”Ionospheric disturbances of the January 15, 2022, Tonga volcanic eruption observed using the GNSS network in New Zealand,” GPS Solutions, årg. 27, s. 53, 2023. Tillgänglig: https://doi.org/10.1007/s10291-023- 01395-8 (hämtad 2024-04-15). [40] E. Astafyeva, B. Maletckii, T. D. Mikesell m. fl., ”The 15 January 2022 Hunga Tonga Eruption History as Inferred From Ionospheric Observations,” 10, maj 2022, [Online]. doi: 10.1029/2022GL098827. 35 https://igs.org/formats-and-standards/ https://igs.org/formats-and-standards/ https://doi.org/10.1016/j.asr.2020.04.015 http://creativecommons.org/licenses/by/4.0/ https://github.com/tanggdut/IonKit-NH/tree/main https://github.com/tanggdut/IonKit-NH/tree/main https://doi.org/10.1007/s10291-018-0759-1 https://doi.org/10.1117/12.2028732 https://doi.org/10.1007/S12145-023-01033-4/TABLES/4 https://doi.org/10.1109/ICIDCA56705.2023.10100215 https://doi.org/10.1109/ICIDCA56705.2023.10100215 https://doi.org/10.1109/PIMRC.2019.8904289 https://doi.org/10.1109/PIMRC.2019.8904289 https://doi.org/10.5281/zenodo.11062720 https://doi.org/10.5281/zenodo.11062720 https://doi.org/10.5281/zenodo.11062720 https://doi.org/10.1126/sciadv.adg8036 https://doi.org/10.1007/s10291-023-01395-8 https://doi.org/10.1007/s10291-023-01395-8 https://doi.org/10.1029/2022GL098827 A Python skript A.1 GMT karta import numpy as np import pygmt import os # Path to current directory current_directory = os.path. dirname (os.path. abspath ( __file__ )) # Find cordinates document and variables for latitude and longitude file_cord = os.path.join( current_directory , ’koordinater . txt ’) file_names = os.path.join( current_directory , ’ koordinater_pos .txt ’) data = np. loadtxt ( file_cord ) data_names = np. genfromtxt (file_names , dtype=’str ’) la , lo = data [:, 0], data [:, 1] names = data_names [:,] # Honga tunga ha apai cordinates la_hthh = -20.545 lo_hthh = -175.3925 fig = pygmt.Figure () fig.coast( # Set the x-range (East/West) and the y-range (North/ South) region="173/203/ -23/ -7", # Set projection to Mercator , and the figure size to 15 centimeters projection ="M15c", # Set the color of the land land=" forestgreen ", # Set the color of the water I A. Python skript water=" skyblue ", # Display the national borders and set the pen thickness to 0.5p borders ="1/0.5p", # Display the shorelines and set the pen thickness to 0.5p shorelines ="1/0.5p", # Set the frame to display annotations and gridlines frame="af", ) fig.plot(x=lo , y=la , style="s0.2c", fill="yellow", pen=" black", label = ’GNSS markstationer ’) fig.plot(x=lo_hthh , y=la_hthh , style="t0.2c", fill="red", pen="black", label = ’Hunga Tonga Hunga Ha ‘ apai ’) fig.legend () fig.show () A.2 Distans from numpy import sin , cos , arcsin , pi , round , sqrt import numpy as np import os # Path to current directory current_directory = os.path. dirname (os.path. abspath ( __file__ )) # Find cordinates document and variables for latitud and longitud file = os.path.join( current_directory , ’koordinater .txt ’) data = np. loadtxt (file) la2 , lo2 = data [:, 0], data [:, 1] # Honga tunga ha apai cordinates la1= -20.545 lo1= -175.3925 # Ground stations name pos= [’KOKB:’, ’UCLU ’, ’CHUR ’, ’CKIS ’, ’THTI ’, ’GAMB ’, II A. Python skript ’TONG ’, ’CHTI ’, ’SCTB ’, ’LAUT ’, ’USP1 ’, ’PTVL ’, ’TOW2 ’, ’MIZU ’, ’FTNA ’, ’SAMO ’, ’TUVA ’, ’TASH ’, ’OUS2 ’ ] # Degrees to radians def deg2rad ( degrees ): radians = degrees * pi / 180 return radians # Haversive formula def getDistance (latitude1 , longitude1 , latitude2 , longitude2 ): delta_lo = longitude2 - longitude1 delta_la = latitude2 - latitude1 r = 6371 # earth radius distance = (2 * r * arcsin(sqrt(sin( deg2rad ( delta_la /2))**2 + cos( deg2rad (la1)) * cos( deg2rad (la2)) * sin( deg2rad ( delta_lo /2))**2))) return round(distance ,2) # kilometers # Stores estimated time it takes for a wave (from HTHH to several ground stations ) est_time = [] # Calculating estimated time for i in range(data.shape [0]): # Check the number of rows in data file distance = getDistance (la1 ,lo1 ,la2 ,lo2) # Input of the cordinates p_wave = ((( distance [i ]*1000) /310) /60) # Velocity 310 m/s or 340 m/s. /60 because we want minutes . est_time .append(p_wave) print(pos[i] ,’Avstand :’, distance [i], ’km , ljudvag : ’, p_wave , ’min ’) III A. Python skript # Make a dictonary , to easy copy into other python scripts distance_dict = dict(zip(pos , est_time )) print( ’\n’, distance_dict , ’\n’) print(data.shape [0]) IV A. Python skript A.3 ZTD import matplotlib .pyplot as plt import numpy as np import math from scipy.signal import savgol_filter from scipy.fft import fft from datetime import datetime , timedelta from operator import add import os # ZTD Stations : # TONG , OUS2 , MIZU , TASH , FTNA , LAUT , SAMO , TUVA , USP1 , CKIS # Choose Station : station = "TONG" # Settings for plot: # Enable ZTD Plot: ztd = True # Enable Saastamoinen plot. Only available for: TONG , OUS2 , MIZU , TASH saastamoinen = False # Enable Savitzky Golay filter: savitzky = True # Enable Errorbars for ZTD: error_ZTD = False # Enable explosions estimate : estimate = True # Enable png save: save = False # Enable ZHD. Only available for: TONG , OUS2 , MIZU , TASH enable_zhd = False # Dictionary of minutes until pressurewave at 310 m/s: est_time_dict = {’TONG ’: 3.7801075268817206 , ’LAUT ’: 44.08924731182796 , ’USP1 ’: 37.58387096774194 , ’MIZU ’: 431.647311827957 , ’FTNA ’: 40.39784946236559 , ’SAMO ’: 45.13494623655914 , ’TUVA ’: 78.3516129032258 , ’TASH ’: 730.3059139784946 , ’OUS2 ’: 166.41559139784945 , ’CKIS ’: 87.14139784946236} # Convert time tag V A. Python skript sec = ’30s’ def sec2date (sec): ref_time = datetime (2000 , 1, 1, 12, 0, 0) return ref_time + timedelta ( seconds =sec) # Get path current_directory = os.path. dirname (os.path. abspath ( __file__ )) subfolder_path = os.path.join( current_directory , ’ZTD_ ’+ sec+"_ny") # Save location for plots autosave_path = os.path.join( current_directory , ’ Autosave_plots ’ ) # Read contents file_path = os.path.join( subfolder_path , ’2022 -01 -15 _7’+ station +’.txt ’) data = np. loadtxt ( file_path ) # Extract time , second column and fourth column. Delta is the error time_raw , column2 , delta , column4 = data [:, 0], data [:, 1], data [:, 2], data [:, 3] time = list(map(sec2date , time_raw )) # Sum the columns sum= column2 + column4 # A variable for estimated wave arrival minutes =0 minutes_est = est_time_dict [ station ] # Plot: plt.figure (1) if ztd: plt.plot(time , sum , label="ZTD ", color=’c’) if savitzky : # Low pass filter of ZTD , Saastamoinen ZTD and ZHD savitzky_sum = savgol_filter (sum , 40, 2) plt.plot(time , savitzky_sum , label="LP filtrerad ZTD", color=’mediumblue ’) VI A. Python skript # Only if Saastamoinen plot is enabled : if saastamoinen : # Gets location of Meteo -Rinex file subfolder_path_t = os.path.join( current_directory , ’ Meteo_data ’) file_path_t = os.path.join( subfolder_path_t , station + ’_meteo_20220115 .txt ’) data_t = np. loadtxt ( file_path_t ) # Reads rinex file # Creates lists of the different parameters in Rinex file year , month , day , hour , minutes , pressure , temp , rhum = data_t [:, 0], data_t [:, 1], data_t [:, 2], data_t [:, 3], data_t [:, 4], data_t [:, 6], data_t [:, 7], data_t [:, 8] # Extracts and creates correct date -time format time_Meteo = [ datetime (year =2022 , month =1, day=int(d) , hour=int(h), minute=int(m)) for d, h, m in zip( day , hour , minutes )] # Different coordinates for the different stations if station == "TONG": phi = -0.369 height = 0.00563 elif station == "OUS2": phi = -0.8005 height = 0.00261 elif station == "MIZU": phi = 0.683 height = 0.0117 elif station == "TASH": phi = 0.721 height = 0.04397 i = 0 # e is a part of the complete Saastamoinen equation to simplify the code e = [] for t in temp: e.append ((6.1078 * math.exp ((17.269 * float(t))/( float(t) + 273.3))) * (float(rhum[i]) /100)) i = i + 1 # Calculate the Zenit wet delay j = 0 ZWD_saas = [] for z in e: ZWD_saas .append (0.0022768 * ((0.05 + (1255/( temp[ VII A. Python skript j] + 273.15) )) * z) / (1 -0.00266 * math.cos (2 * phi - 0.00028 * height))) j = j + 1 # Calculate the Zenit hydrostatic delay ZHD_saas = [] for p in pressure : ZHD_saas .append (0.0022768 * p / (1 -0.00266 * math .cos (2 * phi) - 0.00028 * height)) # Add Wet delay and Hydrostatic delay to get Total delay ZTD_saas = list(map(add , ZWD_saas , ZHD_saas )) # Plot Saastamoinen calculated total delay plt.plot(time_Meteo , ZTD_saas , label=" Saastamoinen ZTD", color=’sandybrown ’) if savitzky : # Low pass filter of Saastamoinen ZTD savitzky_saas = savgol_filter (ZTD_saas , 40, 2) plt.plot(time_Meteo , savitzky_saas , label="LP filtrerad Saastamoinen ZTD", color=’tab:red ’) if error_ZTD : plt. fill_between (time , sum - delta , sum + delta , alpha =0.3 , label=’OsÃďkerhetsintervall ’) # Save PNG file of plot if save: plt. savefig ( autosave_path +’/’+ station +’.png ’, format= ’png ’) if enable_zhd == True: subfolder_path_t = os.path.join( current_directory , ’ Meteo_data ’) # Meteodata from stations is 5 min interval exept for TONG (30 sec) if station == "TONG": subfolder_path_ZTD = os.path.join( current_directory , ’ZTD_30s ’) else: subfolder_path_ZTD = os.path.join( current_directory , ’ZTD_300s ’) file_path_t = os.path.join( subfolder_path_t , station + ’_meteo_20220115 .txt ’) file_path_ZTD = os.path.join( subfolder_path_ZTD , ’ 2022 -01 -15 _7’+ station +’.txt ’) data_t = np. loadtxt ( file_path_t ) VIII A. Python skript data_ZTD = np. loadtxt ( file_path_ZTD ) # Data from Meteo year , month , day , hour , minutes , pressure , temp , rhum = data_t [:, 0], data_t [:, 1], data_t [:, 2], data_t [:, 3], data_t [:, 4], data_t [:, 6], data_t [:, 7], data_t [:, 8] # Data from ZTD time_raw , column2 , delta1 , column4 = data_ZTD [:, 0], data_ZTD [:, 1], data_ZTD [:, 2], data_ZTD [:, 3] # Create the ZTD from the columns ZTD = column2 + column4 # Time series for ZTD def sec2date (sec): # Converts time tag ref_time = datetime (2000 , 1, 1, 12, 0, 0) return ref_time + timedelta ( seconds =sec) time_ZTD = list(map(sec2date , time_raw )) # Time series for Meteo time_Meteo = [ datetime (year =2022 , month =1, day=int(d) , hour=int(h), minute=int(m)) for d, h, m in zip( day , hour , minutes )] # Different parameters for the stations if station == "TONG": phi = -0.369 height = 0.00563 elif station == "OUS2": phi = -0.8005 height = 0.00261 elif station == "MIZU": phi = 0.683 height = 0.0117 elif station == "TASH": phi = 0.721 height = 0.04397 i = 0 # e is a part of the complete Saastamoinen equation e = [] for t in temp: e.append ((6.1078 * math.exp ((17.269 * float(t))/( float(t) + 273.3))) * (float(rhum[i]) /100)) i = i + 1 j = 0 # Calculate wet dealy ZWD = [] IX A. Python skript for z in e: ZWD.append (0.0022768 * ((0.05 + (1255/( temp[j] + 273.15) )) * z) / (1 -0.00266 * math.cos (2 * phi - 0.00028 * height))) j = j + 1 # Remove the Wet delay from Total delay to get Hydrostatic delay ZHD = [] for k in range(len(ZWD)): ZHD.append(ZTD[k] - ZWD[k]) # Plot Hydrostatic delay plt.plot(time_Meteo , ZHD , label="ZHD", color=’ palegreen ’) # Plot Savitzky Golay: if savitzky : # Low pass filter of ZHD savitzky_ZHD = savgol_filter (ZHD , 40, 2) plt.plot(time_Meteo , savitzky_ZHD , label="LP filtrerad ZHD", color=’g’) # Estimated time of explosions est1_start = datetime (2022 , 1, 15, 4, 20, 0) + timedelta ( minutes = minutes_est ) est1_end = datetime (2022 , 1, 15, 4, 22, 30) + timedelta ( minutes = minutes_est ) est2_start = datetime (2022 , 1, 15, 4, 25, 0) + timedelta ( minutes = minutes_est ) est2_end = datetime (2022 , 1, 15, 4, 28, 0) + timedelta ( minutes = minutes_est ) est3_start = datetime (2022 , 1, 15, 4, 51, 30) + timedelta ( minutes = minutes_est ) est3_end = datetime (2022 , 1, 15, 4, 53, 0) + timedelta ( minutes = minutes_est ) est4_start = datetime (2022 , 1, 15, 5, 5, 30) + timedelta ( minutes = minutes_est ) est4_end = datetime (2022 , 1, 15, 5, 7, 30) + timedelta ( minutes = minutes_est ) est5_start = datetime (2022 , 1, 15, 5, 8, 30) + timedelta ( minutes = minutes_est ) est5_end = datetime (2022 , 1, 15, 5, 15, 30) + timedelta ( minutes = minutes_est ) # Plot estimated bars if estimate : # Plot vertical bars where the change is expected X A. Python skript col = ’sandybrown ’ plt. fill_betweenx (y=[ min(sum) -0.5, max(sum)+0.2] , x1= est1_start , x2=est1_end , color=col , alpha =0.3 , label=’ Explosioner ’) plt. fill_betweenx (y=[ min(sum) -0.5, max(sum)+0.2] , x1= est2_start , x2=est2_end , color=col , alpha =0.3 ) plt. fill_betweenx (y=[ min(sum) -0.5, max(sum)+0.2] , x1= est3_start , x2=est3_end , color=col , alpha =0.3 ) plt. fill_betweenx (y=[ min(sum) -0.5, max(sum)+0.2] , x1= est4_start , x2=est4_end , color=col , alpha =0.3 ) plt. fill_betweenx (y=[ min(sum) -0.5, max(sum)+0.2] , x1= est5_start , x2=est5_end , color=col , alpha =0.3 ) # Window config plt.ylabel("ZTD") plt.xlabel("Tid") plt.title("ZTD , " + station + ’, ’ + sec + ’, ’ + " 2022 -01 -15") plt.xticks( rotation =45) plt.gca ().xaxis. set_major_formatter (plt. matplotlib .dates. DateFormatter (’%H:%M ’)) plt.legend () plt.show () XI A. Python skript A.4 ZTD och tryck import matplotlib .pyplot as plt import numpy as np from datetime import datetime , timedelta import os # Choose station station = ’TONG ’ # Dictionary of minutes until pressurewave at 310 m/s: est_tid = {’TONG ’: 3.7801075268817206 , ’MIZU ’: 431.647311827957 , ’TASH ’: 730.3059139784946 , ’OUS2 ’: 166.41559139784945} # Get estimated minutes for the choosen station minutes_est = est_tid [ station ] # Get path current_directory = os.path. dirname (os.path. abspath ( __file__ )) subfolder_path_t = os.path.join( current_directory , ’ Tryck_data ’) subfolder_path_ZTD = os.path.join( current_directory , ’ ZTD_data ’) file_path_t = os.path.join( subfolder_path_t , station +’ _tryck.txt ’) file_path_ZTD = os.path.join( subfolder_path_ZTD , ’ 2022 -01 -15 _7’+ station +’.txt ’) # Read content data_p = np. loadtxt ( file_path_t ) data_ZTD = np. loadtxt ( file_path_ZTD ) # Data from file year , month , day , hour , minutes , pressure = data_p [:, 0], data_p [:, 1], data_p [:, 2], data_p [:, 3], data_p [:, 4], data_p [:, 6] # Timeline tidtagg = [ datetime (year =2022 , month =1, day=int(d), hour= int(h), minute=int(m)) for d, h, m in zip(day , hour , minutes )] tidtagg_est = datetime (2022 , 1, 15, 5, 0, 0) + timedelta ( minutes = est_tid [ station ]) XII A. Python skript # Data from ZTD time1_raw , column2_1 , delta1 , column4_1 = data_ZTD [:, 0], data_ZTD [:, 1], data_ZTD [:, 2], data_ZTD [:, 3] # Create new array , the sum of second and fourth column sum = column2_1 + column4_1 # Timeline for ZTD def sec2date (sec): # Transform the time ref_time = datetime (2000 , 1, 1, 12, 0, 0) return ref_time + timedelta ( seconds =sec) time_ZTD = list(map(sec2date , time1_raw )) # Volcanic explosions intervals est1_start = datetime (2022 , 1, 15, 4, 20, 0) + timedelta ( minutes = minutes_est ) est1_end = datetime (2022 , 1, 15, 4, 22, 30) + timedelta ( minutes = minutes_est ) est2_start = datetime (2022 , 1, 15, 4, 25, 0) + timedelta ( minutes = minutes_est ) est2_end = datetime (2022 , 1, 15, 4, 28, 0) + timedelta ( minutes = minutes_est ) est3_start = datetime (2022 , 1, 15, 4, 51, 30) + timedelta ( minutes = minutes_est ) est3_end = datetime (2022 , 1, 15, 4, 53, 0) + timedelta ( minutes = minutes_est ) est4_start = datetime (2022 , 1, 15, 5, 5, 30) + timedelta ( minutes = minutes_est ) est4_end = datetime (2022 , 1, 15, 5, 7, 30) + timedelta ( minutes = minutes_est ) est5_start = datetime (2022 , 1, 15, 5, 8, 30) + timedelta ( minutes = minutes_est ) est5_end = datetime (2022 , 1, 15, 5, 15, 30) + timedelta ( minutes = minutes_est ) style=’dashed ’ # Plot figure fig , ax1 = plt. subplots () tids_est = datetime (2022 , 1, 15, 4, 20, 0) + timedelta ( minutes = minutes_est ) plt. axvline (x=tids_est , color=’sandybrown ’, linestyle = style , label=’Explosioner ’) tids_est = datetime (2022 , 1, 15, 4, 27, 0) + timedelta ( XIII A. Python skript minutes = minutes_est ) plt. axvline (x=tids_est , color=’sandybrown ’, linestyle = style ) tids_est = datetime (2022 , 1, 15, 4, 51, 0) + timedelta ( minutes = minutes_est ) plt. axvline (x=tids_est , color=’sandybrown ’, linestyle = style ) tids_est = datetime (2022 , 1, 15, 5, 5, 0) + timedelta ( minutes = minutes_est ) plt. axvline (x=tids_est , color=’sandybrown ’, linestyle = style ) tids_est = datetime (2022 , 1, 15, 5, 8, 0) + timedelta ( minutes = minutes_est ) plt. axvline (x=tids_est , color=’sandybrown ’, linestyle = style ) # Plot barometric pressure on the first y-axis ax1.plot(tidtagg , pressure , label=’Barometriskt tryck ’, color=’skyblue ’) ax1. set_xlabel (’Tid ’) ax1. set_ylabel (’Barometriskt tryck P (hPa)’, color=’ skyblue ’) ax1. tick_params (axis=’y’, labelcolor =’skyblue ’) # Choose color on y-axis # Create second y-axis based on the first one ax2 = ax1.twinx () ax2.plot(time_ZTD , sum , label="ZTD", color=’forestgreen ’) ax2. set_ylabel (’ZTD ’, color=’forestgreen ’) ax2. tick_params (axis=’y’, labelcolor =’forestgreen ’) plt.title( station ) plt.xticks( rotation =45) plt.gca ().xaxis. set_major_formatter (plt. matplotlib .dates. DateFormatter (’%H:%M ’)) plt. tight_layout () plt.show () XIV B Matlabkod B.1 Tid- och distanskarta clear all; close all; % ------------------------------------------------- % time - distance plot % ------------------------------------------------- dtecfile =’dTEC_12stationer / allstations_12 .dte ’; data=load( dtecfile ); x=data (: ,2) /3600; % s-->h y=data (: ,9); % km z=data (: ,5); % TECU % -------------------------------------------------- epi =4.25; % Eruption , time 4.25 xt1 =[ epi :0.5:16]; % timeline x-axis % Linear lines representing velocity yt1 =310*( xt1 -epi) *3600/1000; % distance in km yt2 =250*( xt1 -epi) *3600/1000; % distance in km yt3 =150*( xt1 -epi) *3600/1000; % distance in km figure hold on scatter (x,y,1,z,’.’, ’DisplayName ’, ’Satellitdata ’); plot(xt1 ,yt1 ,’k--’, ’DisplayName ’, ’310 m/s’) plot(xt1 ,yt2 ,’k-’, ’DisplayName ’, ’250 m/s’) plot(xt1 ,yt3 ,’k.-’, ’DisplayName ’, ’180 m/s’) colormap (jet); colorbar xlabel(’UT (h)’); ylabel(’Epicenter distance (km)’); legend(’show ’); XV B. Matlabkod set(gca ,’FontSize ’ ,15) set(gca ,’YLim ’ ,[0 10000]) B.2 dTEC-filer, antalet unika satelliter clear all; close all; % ---------------------------------------------- % function : Check how many diffent satellites % appearing in a dTEC files. % input: dte files % ---------------------------------------------- % Name of GNSS stations stations = { ’samo ’, ’uclu ’, ’chur ’, ... ’ckis ’, ’thti ’, ’gamb ’, ... ’tong ’, ’chti ’, ’sctb ’, ... ’laut ’, ’ptvl ’, ’tow2 ’ }; for i = 1: length( stations ); % Load datafiles filename = [’dTEC_12stationer /’ stations {i} ’0150.22 o .dte ’]; data = load( filename ); % Unique values [ unique_sat ] = unique(data (: ,3)); disp ([ stations {i}, ’ ’, num2str (length( unique_sat ))] ) end B.3 dTEC filter, alla stationer clear all; close all; % ---------------------------------------------- % function : Make a figure containing satellite dTEC % graphs , based on three distances from % satellite to HTHH , from larger dte file from % several stations with several satellites . % input: dte file % ---------------------------------------------- % Load file XVI B. Matlabkod filename = [’dTEC_12stationer / allstations_12 .dte ’]; data = load( filename ); % Colors to the graphs colors = {’b’,’b’, ’r’, ’g’, ’k’, ’c’, ’m’, ’b’, ’r’, ’g’ , ’k’, ’c’, ’g’, ’k’, ’c’}; c = 1; % Set three minimum and maximum distance min1 = 800; max1 = 1300; min2 = 2500; max2 = 2800; min3 = 4400; max3 = 4900; % Open the plot window figure; hold on; % Three forloops which iterate satellites 1 -180 and collect data from three % different distances due to min and max variable . for i = 1:180; % Find all row index to i-number satellite , data before 12.05 pm and % distance interval . [ii] = find(data (: ,3) == i & data (: ,2) < 43300 & round (data (: ,9)) > min1 & round(data (: ,9)) < max1 ); [unique_stn , ~, freq_stn ] = unique(data(ii ,4)); % Unique GNSS stations that appear [unique_sat , ~, freq_sat ] = unique(data(ii ,3)); % Unique satellites that appear % Want to have longer data series to plot if length(ii) > 150; date = datetime ( num2str (data(ii , 1)), ’ InputFormat ’, ’yyyyMMdd ’); timetag = date + seconds (data(ii ,2)); if ~ isempty (ii); % For different colors in the plot c = c + 1; XVII B. Matlabkod end dis = mean(data(ii ,9)); plot( datetime (timetag , ’format ’, ’hh:mm:ss’), data(ii ,5) ,[colors{c} ’.-’], ’DisplayName ’, [ ’St ’ num2str ( unique_stn (1)) ’ ,Sa ’ num2str ( unique_sat (1)) ’ Dist ’ num2str (dis) ]); % disp(’Unika stationer :’); % disp( unique_stn ); end end % ----------------------------- 2nd distance % Satellites from 1 to 180 for i = 1:180; % Find all row index to i-number satellite , data before 12.05 pm and % distance interval . [ii] = find(data (: ,3) == i & data (: ,2) < 30000 & round(data (: ,9)) > min2 & round(data (: ,9)) < max2 ); % Find all row index to i-number satellite % Want to have longer data series to plot if length(ii) > 250 & length(ii) < 350; date = datetime ( num2str (data(ii , 1)), ’ InputFormat ’, ’yyyyMMdd ’); timetag = date + seconds (data(ii ,2)); if ~ isempty (ii); % For different colors in the plot c = c + 1; end [unique_stn , ~, freq_stn ] = unique(data(ii ,4)); % Unique GNSS stations that appear [unique_sat , ~, freq_sat ] = unique(data(ii ,3)); % Unique satellites that appear dis = mean(data(ii ,9)); plot( datetime (timetag , ’format ’, ’hh:mm:ss’), data(ii ,5) ,[colors{c} ’.-’], ’DisplayName ’, [ ’St ’ num2str ( unique_stn (1)) ’ ,Sa ’ num2str ( unique_sat (1)) ’ Dist ’ num2str (dis) ]); end end XVIII B. Matlabkod % ----------------------------------- 3rd distance % Satellites from 1 to 180 for i = 1:180; % Find all row index to i-number satellite , data before 12.05 pm and % distance interval . [ii] = find(data (: ,3) == i & data (: ,2) < 40000 & data (: ,2) > 21500 & round(data (: ,9)) > min3 & round( data (: ,9)) < max3 ); % Find all row index to i- number satellite % Want to have longer data series to plot if length(ii) > 350 ; date = datetime ( num2str (data(ii , 1)), ’ InputFormat ’, ’yyyyMMdd ’); timetag = date + seconds (data(ii ,2)); if ~ isempty (ii); % For different colors in the plot c = c + 1; end [unique_stn , ~, freq_stn ] = unique(data(ii ,4)); % Unique GNSS stations that appear [unique_sat , ~, freq_sat ] = unique(data(ii ,3)); % Unique satellites that appear dis = mean(data(ii ,9)); plot( datetime (timetag , ’format ’, ’hh:mm:ss’), data(ii ,5) ,[colors{c} ’.-’], ’DisplayName ’, [ ’St ’ num2str ( unique_stn (1)) ’ ,Sa ’ num2str ( unique_sat (1)) ’ Dist ’ num2str (dis) ]); end end xlabel(’Tid (h)’); ylabel(’dTEC ’); title ([’Delta total electron content , ’ ]); legend(’show ’); B.4 dTEC, tre satelliter clear all; close all; XIX B. Matlabkod % ---------------------------------------------- % function : Make a figure containing 3 satellite % dTEC graphs , based on three distances from % satellite point to HTHH. % input: dte files % ---------------------------------------------- % Name of GNSS stations station1 = ’laut ’ ; station2 = ’chti ’ ; station3 = ’kokb ’ ; % Load datafiles filename1 = [’dTEC_12stationer /’ station1 ’0150.22 o.dte ’ ]; data1 = load( filename1 ); filename2 = [’dTEC_12stationer /’ station2 ’0150.22 o.dte ’ ]; data2 = load( filename2 ); filename3 = [’dTEC_12stationer /’ station3 ’0150.22 o.dte ’ ]; data3 = load( filename3 ); % Set the minimum and maximum distance for the satellites min1 = 900; max1 = 1100; min2 = 2500; max2 = 2700; min3 = 4300; max3 = 4900; % Colors in the plot colors = {’b’, ’r’, ’g’, ’k’}; % Bring up plot window figure; hold on; % These following rows below are collecting data index vectors and plot the data. % Its same structure for all satellites so the comments for the code for satellite 3 is also relevant to the XX B. Matlabkod other 2. % ---------------- Satellite 3 ------------------ % Satellite number i_3 = 145; % ii_3 find index from data3 file which matches the satellite number , time before % 12.05 pm , distance from hthh >min3 and min3 & round(data3 (: ,9)) < max3 ); [ unique_sat ] = unique(data3(ii_3 ,3)); % Unique satellites that appear date = datetime ( num2str (data3(ii_3 , 1)), ’InputFormat ’, ’ yyyyMMdd ’); % Get the date timetag = date + seconds (data3(ii_3 , 2)); % Date + get the time dis = mean(data3(ii_3 ,9)); % Mean distance km plot( datetime (timetag , ’format ’, ’hh:mm:ss’), data3(ii_3 ,5) ,[colors {3} ’.-’], ’DisplayName ’, [ station3 ’-’ num2str ( unique_sat (1)) ’ , ’ num2str (dis) ’ km’]); % ---------------- Satellite 2 ------------------ i_2 = 123; [ii_2] = find(data2 (: ,3) == i_2 & data2 (: ,2) < 43300 & round(data2 (: ,9)) > min2 & round(data2 (: ,9)) < max2 ); [ unique_sat ] = unique(data2(ii_2 ,3)); date = datetime ( num2str (data2(ii_2 , 1)), ’InputFormat ’, ’ yyyyMMdd ’); timetag = date + seconds (data2(ii_2 ,2)); dis = mean(data2(ii_2 ,9)); plot( datetime (timetag , ’format ’, ’hh:mm:ss’), data2(ii_2 ,5) ,[colors {2} ’.-’], ’DisplayName ’, [ station2 ’-’ num2str ( unique_sat (1)) ’ , ’ num2str (dis) ’ km’ ]); % ---------------- Satellite 1 ------------------ i_1 = 81; [ii_1] = find(data1 (: ,3) == i_1 & data1 (: ,2) < 43300 & round(data1 (: ,9)) > min1 & round(data1 (: ,