Monitoring Erebus Volcano’s
Active Lava Lake; Tools,
Techniques and Observations.
Nial John Peters
St Catharine’s College
April 2014
A dissertation presented for the degree of Doctor of Philosophy
Aaron, this is all your fault...
Cover Image: A sequence of false-colour images captured with a thermal infra-red camera in
December 2010 showing a large gas bubble arriving in the Erebus lava lake.
i
Declaration
This dissertation is the result of my own work and includes nothing which is the
outcome of work done in collaboration except where specifically indicated in the text.
I confirm that it does not exceed the word limit required by the degree committee of
Earth Sciences and Geography and does not contain work that has been submitted
for a degree, diploma or other qualification at any other university.
Nial Peters
ii
Foreword
No one ever stays for the ending credits, so I felt impelled to write something mean-
ingful at the beginning instead. Having had whatever little literary flair I might
once have possessed crushed from me by too much scientific-style writing (or was it
too little? I’m not sure any more), it will be neither fluent nor eloquent, but anyway
here goes...
I cannot say that I have researched this thoroughly, but it seems a safe assumption
to make that not many caving accidents result in PhDs being awarded∗. However,
had it not been for that unfortunate event in the Austrian Alps several years ago,
I suspect that I would never have been invited to Antarctica as “caving support”,
would therefore never have met my supervisor, and probably would never have em-
barked on a PhD at all. A series of unfortunate and regrettable events perhaps?
Well, I leave that for my examiners to determine. Just a few months into my PhD,
I found myself playing on the swings in the Christchurch botanic gardens with my
supervisor. “You know, this is exactly what I imagined doing a PhD would be like”
I said†. Admittedly, my initial ideas about what was involved were somewhat wide
of the mark, the remainder of the PhD certainly involved far fewer toys designed
for young children. However, toys or no toys, the whole experience (with just a few
notable exceptions) has been hugely enjoyable and, truthfully, I have absolutely no
regrets. For fear of sounding vindictive, I will refrain from wishing you as much
fun reading this dissertation as I had writing it, but you should at least glance-over
the acknowledgements section at the end. To me at least, it is the most meaningful
part of this thesis and accredits the greatest things that I take away from this – the
experiences shared with my friends and colleagues in some of the most spectacular
places on Earth.
∗The “awarded” part is obviously somewhat unclear at the time of writing.
†His reply was: “Yeah, it’s pretty much how I imagined supervising would be too!”
iii
Summary
Active lava lakes present a rare opportunity to observe directly the complex pro-
cesses occurring within a magma body. Situated on Ross Island, Antarctica, the
3794-m-high crater of Erebus volcano has hosted a phonolite lava lake for decades.
Previous studies have shown that many of the lake’s characteristics, such as surface
velocity, gas flux and gas composition, exhibit a pronounced pulsatory behaviour
on a time-scale of ∼10 min. Focusing primarily on the analysis of infra-red (IR)
imagery acquired from the crater rim, this dissertation considers how the periodic
behaviour of the Erebus lava lake evolves over decadal time periods, how the cyclic
fluctuations of the different properties are interrelated and what can be inferred
about the mechanisms occurring beneath the surface of the lake from these obser-
vations. Creation of new hardware, software and methodologies to facilitate these
types of observations is a strong focus of this work.
Chapter 1 introduces the nature of active lava lakes, reviews previous studies of
Erebus and presents in detail the research objectives that are addressed by the
subsequent chapters. In Chapter 2, a new thermal camera system that was developed
as part of this study is described. Designed to run autonomously at the crater-rim
of Erebus, this system was installed in December 2012 and has enabled, for the first
time, extended time-series of images to be acquired. Chapter 3 briefly describes
some of the other hardware and software that was developed as part of this study
and outlines how it has been utilised for volcano monitoring. In Chapter 4, a dataset
of IR images collected between 2004–2011 is used to assess inter-annual variability
in the pulsatory behaviour of the surface motion of the Erebus lava lake. The cyclic
behaviour is found to be a sustained feature of the lake, and no obvious changes are
observed across the time period analysed. Data collected with the camera system
described in Chapter 2 are analysed in Chapter 5 and combined with measurements
from other instruments to assess the correlation between the cyclic behaviours of
different lake properties. Cycles in surface speed, surface elevation, gas flux and gas
composition are found to be highly correlated with each other. In Chapter 6, the
surface velocities calculated in the preceding chapters are revisited, and the two-
dimensional structure of the flow field is analysed. Chapter 7 demonstrates how the
motion tracking methodologies developed for studying the Erebus lava lake can be
used to improve high time resolution sulphur dioxide flux estimates - a significant
challenge faced in the study presented in Chapter 5. Finally, Chapter 8 presents a
synthesis of the key findings and conclusions from the preceding chapters.
iv
Contents
1 Introduction 1
1.1 Open-Vent Volcanism . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.2 Active Lava Lakes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.3 Erebus Volcano . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.4 Degassing Mechanisms . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.4.1 Gas Pistoning . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.4.2 Pressurisation Changes . . . . . . . . . . . . . . . . . . . . . . 7
1.4.3 Bi-Directional Flow Instability . . . . . . . . . . . . . . . . . . 8
1.4.4 Bubble Wakes . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.5 Fieldwork . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.6 Summary of Activity . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
1.7 Research Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.7.1 Objective 1 – Improved Monitoring . . . . . . . . . . . . . . . 12
1.7.2 Objective 2 – Assessment of Annual Variability of the Lava
Lake Cycles . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
1.7.3 Objective 3 – Multi-parameter Observations of Cyclic Behaviour 13
1.7.4 Thesis Overview . . . . . . . . . . . . . . . . . . . . . . . . . 13
1.8 Published and Submitted Materials . . . . . . . . . . . . . . . . . . . 14
2 Development of an Autonomous Thermal Camera System for Mon-
itoring the Lava Lake at Erebus Volcano 16
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.2 Hardware . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.2.1 Camera . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.2.2 Control system . . . . . . . . . . . . . . . . . . . . . . . . . . 22
2.2.3 Power supply . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
2.2.4 Time synchronisation . . . . . . . . . . . . . . . . . . . . . . . 23
2.2.5 Sensors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
2.2.6 Data storage . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
v
2.2.7 Telemetry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
2.2.8 Enclosure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
2.3 Software . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
2.3.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
2.3.2 Implementation and architecture . . . . . . . . . . . . . . . . 27
2.3.3 Reuse potential . . . . . . . . . . . . . . . . . . . . . . . . . . 29
2.4 Results and discussion . . . . . . . . . . . . . . . . . . . . . . . . . . 29
2.4.1 The camera system . . . . . . . . . . . . . . . . . . . . . . . . 29
2.4.2 The dataset . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
2.4.3 Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
3 Other Volcanic Monitoring Hardware and Software 33
3.1 Data Visualisation Software - AvoPlot . . . . . . . . . . . . . . . . . 34
3.1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
3.1.2 Implementation and architecture . . . . . . . . . . . . . . . . 35
3.1.3 Reuse potential . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.1.4 Volcanological Applications . . . . . . . . . . . . . . . . . . . 37
3.2 Portable Spectrometer Scanners . . . . . . . . . . . . . . . . . . . . . 39
3.2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
3.2.2 Hardware and Software . . . . . . . . . . . . . . . . . . . . . . 40
3.2.3 Measuring SO2 Flux . . . . . . . . . . . . . . . . . . . . . . . 43
3.2.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
4 Decadal Persistance of Cycles in Lava Lake Motion 46
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
4.2 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
4.2.1 Camera Hardware . . . . . . . . . . . . . . . . . . . . . . . . . 48
4.2.2 Data Selection . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
4.2.3 Motion Tracking . . . . . . . . . . . . . . . . . . . . . . . . . 51
4.2.4 Time-Series Analysis . . . . . . . . . . . . . . . . . . . . . . . 53
4.2.5 Bubbles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
4.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
4.3.1 Spectrograms . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
4.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
4.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
5 Multi-parameter Observations of Cyclic Lava Lake Behaviour 75
5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
vi
5.2 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
5.2.1 Lava Lake Surface Velocity . . . . . . . . . . . . . . . . . . . . 77
5.2.2 SO2 Flux . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
5.2.3 Gas Composition . . . . . . . . . . . . . . . . . . . . . . . . . 81
5.2.4 Lava Lake Surface Elevation . . . . . . . . . . . . . . . . . . . 82
5.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
5.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
5.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93
6 Analysis of the 2D Surface Velocity of the Erebus Lava Lake 95
6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
6.2 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
6.2.1 Data Selection . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
6.2.2 Divergence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
6.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
6.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
6.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
7 Use of Motion Estimation Algorithms for Improved SO2 Measure-
ments Using UV Cameras 104
7.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
7.2 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107
7.2.1 Motion Estimation Algorithms . . . . . . . . . . . . . . . . . . 107
7.2.2 Algorithm Parameters . . . . . . . . . . . . . . . . . . . . . . 109
7.2.3 Using ATHAM as a Validation Tool . . . . . . . . . . . . . . . 110
7.2.4 Villarrica Dataset . . . . . . . . . . . . . . . . . . . . . . . . . 113
7.2.5 Image Registration . . . . . . . . . . . . . . . . . . . . . . . . 114
7.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115
7.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122
7.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125
8 Conclusions 127
8.1 Improved Monitoring . . . . . . . . . . . . . . . . . . . . . . . . . . . 128
8.2 Assessment of Annual Variability of the Lava Lake Cycles . . . . . . . 129
8.3 Multi-parameter Observations of Cyclic Behaviour . . . . . . . . . . . 129
8.4 The Future . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130
Acknowledgments 131
vii
Appendix A Data Availability 132
Bibliography 135
viii
Chapter 1
Introduction
“I can’t see any possible reason why this won’t work.”
– Bill McIntosh, Erebus, December 2010.
1.1 Open-Vent Volcanism
Broadly speaking, “open-vent” volcanoes are those that persistently expose magma
at the Earth’s surface. In this respect, they encompass many of the popular stereo-
types of volcanic activity, exhibiting sustained passive degassing, which is often ac-
companied sporadically by varying degrees of low-level explosive, active-degassing,
activity [Tamburello et al., 2012]. At many open-vent volcanoes, the magma-
atmosphere interface is very shallow, often with visible incandescence at the surface.
Well-known examples of such systems include Villarrica in Chile, and Stromboli in
the Aeolian Islands, both of which attract many thousands of tourists per year due to
their impressive, but relatively benign, activity. The persistent activity also attracts
volcanologists, in particular those interested in degassing processes [e.g. Witter et al.,
2004, Shinohara and Witter, 2005, Mather et al., 2004a,b, Sawyer et al., 2011] and
also those wishing to test new monitoring technologies [e.g. McGonigle et al., 2009,
Mori and Burton, 2009, Bluth et al., 2007].
1.2 Active Lava Lakes
Persistently active lava lakes are a spectacular but rare form of open-vent volcanism
found at only a handful of volcanoes around the world. Notable examples include
Erta 'Ale in Ethiopia, Kı¯lauea in Hawaii (since 2008), Nyiragongo in the Democratic
Republic of the Congo and, the focus of this study, Erebus in Antarctica. An active
lava lake is the exposed top of a volcano’s magmatic plumbing system. It is impor-
tant to draw the distinction between active lava lakes, which maintain a link with
the deeper magmatic system, and passive lava lakes which do not [Tilling, 1987].
The latter are formed due to ponding of erupted lava in a topographical depression,
which then rapidly cools and solidifies giving the lake a lifespan of hours to days.
Active lava lakes however, may have lifespans of many decades.
Longevity of a lava lake has been argued to reflect either effective transfer of magma
between the lake and the deeper system [e.g. Oppenheimer et al., 2004, Francis et al.,
1993], or a supply of gas bubbles from depth [Witham and Llewellin, 2006, Bouche
et al., 2010]. It can be shown experimentally that processes occurring at depth will
manifest themselves at the surface as changes in the lake’s behaviour, for example
its surface level [Witham et al., 2006] or gas flux [Divoux et al., 2009]. It follows
therefore, that observations of lake properties can yield valuable insights into the
processes occurring at depth in the magmatic system, where direct measurements
2
are not possible. This link with the deeper magmatic system makes the study of
active lava lakes of particular importance.
1.3 Erebus Volcano
Erebus is a 3794-m-high stratovolcano located on Ross Island, Antarctica (Fig. 1.1).
It is often claimed to be the southernmost active volcano in the world and is known
to have hosted an active phonolite lava lake since at least 1972 [Giggenbach et al.,
1973]. Although other small lakes have appeared intermittently over this period, the
main, “Ray” lake (Fig. 1.2), has been a permanent feature of the crater throughout
(with the notable exception of 1984–1985 when it was buried following sustained
explosive eruptions) [Kyle et al., 1990].
Often described as an “ideal laboratory volcano”, the persistent, low-level activity
of Erebus, does indeed make it particularly amenable to detailed close-range study.
Ground-based measurements can be made from the rim of its main crater, which
commands a direct view of the active lava lake. As one of only a handful of active
lava lakes in the world [Tazieff, 1994], and being unique in its phonolitic composi-
tion [Kelly et al., 2008], it has received significant attention from the volcanological
community, and a wide range of studies have been conducted to understand the
magmatic processes responsible for its decadal persistence.
A high degree of thermal structure in the crust, which covers the lake, has enabled
feature tracking algorithms to be applied to time series of infrared images so as to
investigate the lake’s surface motion [Oppenheimer et al., 2009]. The same ther-
mal images have also been used to study the thermal radiance of the lake [Calkins
et al., 2008]. The continuity of thermal radiation from the lake also enables it to be
used as a source for open-path Fourier Transform Infra-Red (FTIR) spectroscopy,
allowing detailed analysis of the plume composition at a high temporal resolution
[Oppenheimer and Kyle, 2008]. The extremely dry atmosphere of Antarctica is
highly advantageous for FTIR measurements and has also enabled terrestrial laser
scanning (TLS) observations to be made of the lake surface [Jones, 2013, Jones et al.,
2014] during periods of particularly low humidity inside the crater.
As with all open-vent volcanoes, Erebus’s volcanic plume is a sustained feature.
During favourable weather conditions it is often coherent and vertically rising over
3
Figure 1.1: Map showing the location of Ross Island and Erebus volcano. Graphic
created by the cartography section of the Cambridge University Geography Depart-
ment.
4
Figure 1.2: The active lava lake at Erebus volcano as it appeared in December 2010.
The lake is ∼40 m across its long axis. Photo: Clive Oppenheimer
a sufficient distance to allow UV spectroscopy to be used to ascertain SO2 flux
[Kyle et al., 1990, Sweeney et al., 2008]. The use of a Dual Wide Field of View
(DW-FOV) instrument has allowed rapid variations in flux to be studied [Boichu
et al., 2010]. During less favourable weather conditions, the plume is often grounded
across the northern side of the crater allowing direct measurements using “multi-
gas” instruments [Shinohara, 2005, Aiuppa et al., 2005] to be made. In addition to
complementing the FTIR measurements, these have also enabled the measurement
of hydrogen emissions [Moussallam et al., 2012].
A common finding of the aforementioned studies is that many of the lava lake’s
properties (SO2 flux, gas composition, radiant heat, surface velocity) demonstrate
pseudo-periodic fluctuations with a period of 5–18 min. In light of the observed
behaviour of lab-based flow experiments [Huppert and Hallworth, 2007], it was pro-
posed by Oppenheimer et al. [2009] that these fluctuations may be due to instability
in bi-directional flow in the lake’s feeder conduit causing a pulsatory flow of hot,
volatile-rich magma from the conduit into the lake. Pseudo-periodic behaviour has
also been observed at Kı¯lauea, where it has been attributed to gas pistoning [Patrick
et al., 2011, Orr and Rea, 2012], and at Erta 'Ale, where it has been explained by
the regular arrival of large bubbles [Bouche et al., 2010]. Experiments conducted by
Witham et al. [2006] have also suggested that bubble driven pressurisation changes
could account for pulsatory behaviour. These mechanisms are discussed in greater
detail and in the context of Erebus in Section 4.4.
The stable convective behaviour of the Erebus lava lake is punctuated by intermit-
tent [De Lauro et al., 2009] Strombolian eruptions associated with the rupture of
large (decametric) gas bubbles at the lake surface [Dibble et al., 2008, Gerst et al.,
5
2013]. Phases of increased and more intense Strombolian activity (several explosions
per day, with ejecta escaping the main crater) recur, lasting 1–10 months and are
followed by more extended intervals during which gas bubble bursts are less frequent
and of a smaller size (a few explosions per week, with all ejecta being confined to
the crater – see for example Jones et al. [2008]). The chemical and mineralogical
composition of erupted lavas has remained constant for approximately 17 ka and the
abundance of unusually large anorthoclase crystals is indicative of sustained shallow
magma convection throughout this period [Kelly et al., 2008]. Indeed, the presence
of such large crystals may be a significant influence on the behaviour of the shallow
convection at Erebus [Molina et al., 2012]. Other properties of the lake also demon-
strate remarkably consistent long-term behaviour, for example SO2 flux [Sweeney
et al., 2008] and radiant heat output [Wright and Pilger, 2008].
1.4 Degassing Mechanisms
As mentioned above, several different mechanisms have been proposed to explain
the degassing behaviours observed at different lava lakes around the world. These
are summarised in the following sections. Distinguishing between these mechanisms
as being the process responsible for the behaviour of the Erebus lava lake is a major
motivation behind the study of the surface velocity of the lake presented in this
thesis.
1.4.1 Gas Pistoning
Gas pistoning has been proposed as a mechanism to explain observed surface level
fluctuations both in the Pu'u 'O¯'o¯ lava lake at Kı¯lauea [Orr and Rea, 2012] and also
in a perched lava channel on the flanks of Kı¯lauea during 2007 [Patrick et al., 2011].
Figure 1.3 shows the different stages of gas pistoning activity which, collectively,
lead to cyclic variations in the surface level of the lava lake. The “filling” stage of
the cycle is gradual, as small gas bubbles rising through the lake are trapped either
beneath a solidified surface crust or a foam layer close to the surface. The increasing
gas volume in the lake causes a rise in surface level. At a critical thickness, the layer
of bubbles collapses, forming large gas slugs which escape from the lake causing a
very rapid reduction in surface level. As small bubbles continue to rise through the
lake, the process begins again.
6
(a) (b) (c) (d)
Figure 1.3: Cartoon of the gas pistoning mechanism (this is a modified version of
the figure in Patrick et al. [2011]). (a) gas bubbles are trapped near the surface of
the lake, either beneath a cooled crust or in a foam layer. (b) Gas accumulation in
the lake causes an increase in surface level. (c) The foam layer collapses into a few
large gas slugs which break the surface, sometimes explosively. (d) With the loss of
gas, the lake surface returns to its normal level, and the process begins again.
Figure 1.4: Cartoon of the bubble-driven pressurisation changes mechanism (figure
taken from Witham et al. [2006]). (a) A regular supply of bubbles from depth
drives upflow of magma into the lake. (b) The upflow of magma is impeded by the
hydrostatic pressure in the lake, allowing bubbles to coalesce into large slugs. (c)
The release of large gas slugs decreases the hydrostatic pressure at the base of the
conduit relative to the lake, triggering a downflow of magma.
1.4.2 Pressurisation Changes
The bubble-driven pressurisation changes mechanism was proposed by Witham et al.
[2006] based on laboratory based analogue models of lava lake systems. In their
experiments, gas was injected into a vertical pipe (the conduit) filled with water,
with a tank on the top (the lake) and connected to a large reservoir of water (such
that the volume of water in the lake-conduit system could change). When the gas
supply was first turned on, bubbles rising through the conduit drove an upflow of
liquid from the conduit into the lake. As the lake filled, the increase in hydrostatic
head stalled this upflow, allowing the bubbles to coalesce into large gas slugs. Release
of these slugs causes fluctuations in the pressure at base of the conduit, and when
these are sufficiently large the pressure in the lake exceeds that at the bottom of
the conduit resulting in downflow of liquid. Downflow continues until hydrostatic
equilibrium is regained at which point the cycle repeats.
7
Figure 1.5: Cartoon of the bi-directional flow mechanism (figure taken from Oppen-
heimer et al. [2009]). Flow instability associated with bi-directional flow of magma
in the conduit results in a pulsatory supply of bubble-rich magma from depth into
the lake. (a) In between arrival of pulses of magma, convection in the lake is slug-
gish and the gas emissions are dominated by gas exolved at depth which is able to
rise to the surface independent of the melt. (b) Pulses of magma result in increased
convection speeds and a change in gas signature due to increased shallow exsolution.
1.4.3 Bi-Directional Flow Instability
Huppert and Hallworth [2007] demonstrated that, at certain Reynolds numbers, two
fluids of different densities flowing in opposite directions in a vertical pipe develop
a flow instability that leads to a highly periodic pulsed flow. Oppenheimer et al.
[2009] proposed that such a mechanism may be applicable to volcanic conduit-lake
systems. In this case, the density disparity between upflow and downflow is due
to gas loss from magma in the lake. Hence, the upflow consists of bubble rich,
low density magma from depth whereas the downflow consists of denser, degassed
magma from the lake. The pulsatory arrival of fresh batches of magma into the
lake gives rise to cyclic variations in the chemistry of the gas emissions, the surface
velocity of the lake and its surface level.
8
Figure 1.6: Sketch of the bubble wake mechanism (figure taken from Bouche et al.
[2010]). Large bubbles rising through the conduit entrain small bubbles in their
turbulent wake. The wake rises with the large bubble but can become unstable
and detach. In this case the wake rises more slowly and independently of the large
bubble, behaving as a bubble cluster.
1.4.4 Bubble Wakes
Although not directly related to cyclic behaviour of lava lakes, the degassing mech-
anism proposed by Bouche et al. [2010] is also of interest. Based on observations of
the Erta 'Ale lava lake, the principal idea is that large bubbles whose buoyancy is
sufficient to allow them to rise through the conduit independently of the melt create
a turbulent wake behind them. Small bubbles which are largely static with respect
to the melt, become entrained in the wake of the large bubbles and rise with them. If
the wake grows beyond a critical size (proportional to the size of the leading bubble)
it will detach from the bubble and rise independently through the melt, behaving as
a bubble cluster. At Erta 'Ale, such bubble clusters are attributed with the sporadic
fire-fountaining activity that is sometimes observed.
1.5 Fieldwork
All data relating to Erebus presented in this thesis have been collected as part of the
Mount Erebus Volcano Observatory’s field campaigns. Fieldwork on Erebus volcano
is limited to the austral summer, and typically takes place between late-November
9
Figure 1.7: The Lower Erebus Hut (just visible in the bottom right) with the main
Erebus crater clearly visible in the background.
and early-January. Where I refer to a field season by year, I am referring to the year
in which it began.
Field operations are based at LEH, the Lower Erebus Hut, (Fig. 1.7), which sits
∼2 km from the crater rim on the north-west side of the caldera. Logistical support
for the field campaigns is provided by the United States Antarctic Program, whose
McMurdo station is a short helicopter flight from LEH. From LEH, the crater rim is
reachable by first driving by snowmobile to Nausea Knob (NKB) and then walking
the final 0.5 km.
I have participated in five consecutive field seasons at LEH, from 2009–2013. How-
ever, the 2009 field season was prior to my involvement in the study being presented
here, and I worked primarily in the ice caves during this season.
1.6 Summary of Activity
In the following analyses, data from field campaigns between 2004 and 2013 have
been used. Although the general behaviour of the lava lake at Erebus is fairly
consistent from year to year, there are some observable variations. It is therefore
important to set the results presented within the context of the state of activity of
the lake during each of the respective field campaigns.
During the 2004 field season, there were two separate lava lakes present in the crater
10
Figure 1.8: A “large” bubble burst in the Erebus lava lake during the 2012 field cam-
paign. The lake is ∼20 m across and there is a ∼2 s time delay between consecutive
frames. Photos: Tim Burton
[Calkins et al., 2008]. By the 2005 field season, only the “Ray” lake remained [Davies
et al., 2008], and no additional lakes have since been observed. All data presented
here are from the “Ray” lake, and henceforth, I refer to it simply as the lava lake.
Since 2004, the visible surface area of the lake has reduced by a factor of approxi-
mately four. Possible reasons for this change are discussed in detail in Section 4.3.
Despite the reduction in visible surface area, there have been no observable changes
in the behaviour of the lava lake. Stable convective behaviour has been maintained
throughout. The convective behaviour is characterised by the lateral migration of
cracks in the lake’s crust across surface. These typically move radially outwards
from the approximate centre of the lake, but more complex flow patterns with sev-
eral “convection cells” and abrupt reversals of motion direction are also common.
Lobes of fresh lava are occasionally observed to flow across the surface of the lake
from upwelling regions, and there is visible subduction of the surface crust at the
downwelling regions.
Bubbles of a variety of sizes are observed to surface in the lake. I describe bubbles as
“large” if they result in significant ejection of material from the lake. Such bubbles
(e.g. Fig. 1.8) are typically 10–30 m in diameter and cause a visible emptying of the
lake. I classify such events as being distinct from the far more frequently occurring
“small”, metre and sub-metre scale bubbles (e.g. panel (d) in Fig. 4.9) which arrive
at the surface of the lake, but do not result in explosive activity.
A study of explosive events (due to large bubbles) between 2003–2011 using seismic
data [Knox, 2012] shows that, with the notable exception of the period from late
2005 to early 2007, the frequency of explosive events has remained fairly constant
11
at a few per week, with ejecta being entirely confined to the crater. During 2005–
2006, however, there were several explosions per day, often of sufficient magnitude
to propel ejecta out of the crater, the frequency of these events then gradually de-
clined and by the 2007 field season the lake had returned to its usual level of activity.
In 2013, Erebus once again entered a particularly active phase, with two or three
explosions each day that were large enough for ejecta to exit the crater completely.
Bombs several metres in length were observed to land up to ∼400 m down-slope
from the crater rim. At the time of writing in 2014, data from the seismic network
deployed around Erebus [Aster et al., 2004] shows that large explosions continue
to occur several times a day (although with a brief period in early March where
activity subsided).
1.7 Research Objectives
The following sections outline the three key research objectives which this study
strives to address. Section 1.7.4 then summarises how these fit together and how
they have been structured within the context of this thesis. It will be seen that
this work has involved the development of an imaging system and of data analysis
techniques, coupled with detailed spatial and temporal analysis of image data and
their interpretation.
1.7.1 Objective 1 – Improved Monitoring
Although field campaigns on Erebus typically last up to two months, the complex
logistics involved in reaching LEH mean that typically data are only collected for a
few weeks. The crater rim of Erebus is a formidable environment in which to conduct
scientific measurements. High winds, corrosive gases and extremely low tempera-
tures (<-50 °C during the winter) pose serious technical challenges for personnel and
instruments alike, compounding the problem of short time series as frequent equip-
ment failures lead to fragmented datasets. Gaps in the data make extended time
series analysis problematic at best, and often impossible.
The first objective therefore, is to improve the data collection capabilities
of the Mount Erebus Volcano Observatory.
12
1.7.2 Objective 2 – Assessment of Annual Variability of the
Lava Lake Cycles
The first ground-based infra-red images of the Erebus lava lake were collected in
2004 [Calkins et al., 2008] and the dataset has been extended during each subse-
quent field campaign (except for the 2005 season, for which the data have not been
made available to the author). While the data are far from continuous, they do
provide suitable resolution and range for a multi-year analysis of the lava lake’s
cyclic behaviour. Despite this, previous studies have only analysed very short time
series (a few hours), and little is known about the persistence and/or evolution of
the pulsatory behaviour on longer timescales.
The second objective therefore, is to characterise and evaluate the be-
haviour of lake on a decadal time scale.
1.7.3 Objective 3 – Multi-parameter Observations of Cyclic
Behaviour
As noted in Section 1.3, cyclic behaviour has been observed in many different prop-
erties of the lake. However, due to the many problems associated with obtaining
continuous datasets outlined above, the relationships between the cycles in differ-
ent properties is largely unknown. These relationships are key in understanding
the mechanism(s) responsible for the pulsatory behaviour of the lake. Furthermore,
since the observed behaviour at the surface is likely controlled by processes occurring
at depth, this presents a valuable opportunity to better understand the dynamics
of the subsurface magmatic system at both Erebus and other open-vent volcanoes
around the world.
The final objective therefore, is to definitively link the cycles in surface
motion, gas composition, gas flux and lake surface level.
1.7.4 Thesis Overview
Each of the research objectives detailed above is addressed in a separate chapter
of this thesis. The hardware and software that I have developed to facilitate im-
proved data collection at Erebus (and other volcanoes around the world) is detailed
in Chapters 2 and 3. In Chapter 4, I use the dataset of IR images collected between
2004–2011 to evaluate the pulsatory behaviour of the lava lake on an inter-annual
13
timescale. More recent data, collected using the camera system described in Chap-
ter 2, is analysed in Chapter 5 and combined with data from other instruments to
show the remarkable correlation between cycles in different properties of the lake
(surface motion, gas composition, gas flux, and surface elevation). Chapter 7 demon-
strates how some of the methods developed for analysing the IR images from Erebus
may be used to address problems encountered in other fields of volcanic monitoring,
specifically sulphur dioxide measurements using UV cameras. Finally, Chapter 8
gives an overview of the key findings in each chapter and discusses what constraints
can be placed on the subsurface magmatic system of Erebus given the results pre-
sented in Chapters 4 and 5.
1.8 Published and Submitted Materials
A significant portion of this thesis is based on articles that have been published in
peer-reviewed journals. As is the nature of modern science, parts of these studies
have involved collaborative work and contributions from co-authors. A statement
of co-author involvement is included at the beginning of each chapter, detailing the
specifics of any collaborations involved. A summary of the publications arising from
this study is given below, along with references to their positions within the body
of this thesis.
Chapter 2 has been published as:
Peters, N., C. Oppenheimer, and P. Kyle. “Autonomous Thermal Camera
System for Monitoring the Active Lava Lake at Erebus Volcano, Antarctica”.
Geoscientific Instrumentation, Methods and Data Systems 3, no. 1 (February
5, 2014): 13–20.
Section 3.1 has been published as:
Peters, Nial. “AvoPlot: An Extensible Scientific Plotting Tool Based on Mat-
plotlib”. Journal of Open Research Software 2, no. 1 (February 10, 2014).
Chapter 4 has been published as:
Peters, Nial, Clive Oppenheimer, Philip Kyle, and Nick Kingsbury. “Decadal
Persistence of Cycles in Lava Lake Motion at Erebus Volcano, Antarctica”.
Earth and Planetary Science Letters 395 (June 1, 2014): 1–12.
Chapter 5 has been published as:
Peters, N., C. Oppenheimer, D. R. Killingsworth, J. Frechette, and P. Kyle,
14
“Correlation of cycles in Lava Lake motion and degassing at Erebus Volcano,
Antarctica”, Geochemistry, Geophysics, Geosystems, 2014
Chapter 7 has been published as:
Peters, N., A. Hoffmann, T. Barnie, M. Herzog, C. Oppenheimer, “Use of
Motion Estimation Algorithms for Improved Flux Measurements Using SO2
Cameras”, Journal of Volcanology and Geothermal Research, in press 2014.
These articles have not been reproduced verbatim within the thesis. In addition
to the necessary changes to combine them into a single cohesive document, where
appropriate they have also been extended to cover greater detail than was possible
within the word- and figure-counts of the journals.
15
Chapter 2
Development of an Autonomous
Thermal Camera System for
Monitoring the Lava Lake at
Erebus Volcano
The work presented in this chapter was undertaken in order to meet the “Improved
Monitoring” research objective of this study (Section 1.7.1). In order to obtain longer
continuous time series of IR images of the Erebus lava lake, a new thermal cam-
era system has been developed and installed at the crater rim. The new system is
designed to be autonomous, and capable of capturing images of the lava lake contin-
uously throughout the year. This represents a significant improvement over previous
systems which required the frequent attention of observatory researchers and could
therefore only be operated during a few weeks of the annual field campaigns. The
extreme environmental conditions at the summit of Erebus pose significant chal-
lenges for continuous monitoring equipment, and a custom made system was the
only viable solution. In this chapter the hardware and software of the new system
is described in detail.
The majority of content from this chapter has been published as: Peters, N., C. Op-
penheimer, and P. Kyle. “Autonomous Thermal Camera System for Monitoring the
Active Lava Lake at Erebus Volcano, Antarctica”. Geoscientific Instrumentation,
Methods and Data Systems 3, no. 1 (February 5, 2014): 13–20.
All text and figures presented in this chapter were created by me, and the role of
the co-authors did not extend beyond that of normal supervision. All design, con-
struction, assembly and on-site installation of the camera system, telemetry system
and power system was undertaken by me, with the notable exceptions of: (i) the
most recent load-controller installed at the Nausea Knob “power station”, which
was designed and installed by Nick Salava from the US Antarctic Program, (ii) the
tripod mount for the camera, which was constructed by Kostas Milonidis. All the
computer codes used for the camera system were written by me.
Figure 2.1: Simultaneous images from a digital SLR camera (note that the visibility
of the lake in this image is exceptional) and the IR camera. The high degree of
thermal structure in the lake’s crust is clearly visible in the IR image and can be
used in time-series of such data to estimate the surface motion of the lake (see
Chapter 4). The ability of the IR camera to image through the volcanic plume is
also clearly demonstrated. The lake is ∼ 40 m across in these images.
2.1 Introduction
Thermal infrared (IR) cameras are an important tool in the study of lava lakes and
have been used extensively in a number of lake studies around the world, ranging
from quantification of radiative heat output [e.g. Calkins et al., 2008, Oppenheimer
et al., 2004] to studying surface velocity [Oppenheimer et al., 2009]. They are also
used more widely as an operational tool in volcano monitoring [e.g. Spampinato
et al., 2011]. Although absorption of infrared radiation by volcanic gases makes
accurate temperature measurements difficult to achieve [Sawyer and Burton, 2006],
the ability to image through optically opaque volcanic plumes (Fig. 2.1) gives IR
cameras a significant advantage over conventional cameras for monitoring lava lakes.
17
The Mount Erebus Volcano Observatory (MEVO) has been operating a thermal
camera on Erebus during annual field campaigns (typically from late November to
early January) since 2004. However, due to the many challenges of working in such
an extreme environment, the time series of IR images tend to be very short and frag-
mented. In 2010, the decision was made to upgrade the thermal camera to a system
that could potentially (contingent on a reliable power supply) run autonomously
and continuously year-round. After two years of development and testing, the new
system was installed during the 2012 field season.
Autonomous instruments have become a popular choice for researchers from a va-
riety of disciplines working in Antarctica. Although it might seem an obvious so-
lution to obtaining long time series measurements without the need for extended
field campaigns, the technical challenges of creating such systems are numerous [see
for example Bauguitte et al., 2011, Lawrence et al., 2004]. MEVO already operates
several year-round seismic stations [Aster et al., 2004], however, the location of the
thermal camera system and its high data rate meant that it presented many new
challenges.
Autonomous thermal camera systems developed by other volcano observatories [e.g.
Patrick et al., 2014] often rely on year-round telemetry to transmit data back to a
centralised server in real-time. However, the power requirements of telemetry sys-
tems presently makes this approach unsuitable at Erebus. The system developed
for Erebus also had to withstand far more extreme environmental conditions than
are found at most other volcanoes that have year-round monitoring systems in place.
2.2 Hardware
The camera system is mostly an aggregation of off-the-shelf parts, with very little
custom-made hardware involved. Although this has meant that some compromises
have had to be made in terms of specifications, the reduction in cost and develop-
ment time more than outweighs the disadvantages. Figures 2.2 and 2.3 show the
main components used, which are discussed in detail in the following sections. The
approximate cost of the components are outlined in Table 2.1. It should be noted
however, that the control system could be used with any (possibly cheaper) camera
unit, provided that it supports the GenICam interface.
18
Component Approximate Cost
Thermal Camera £20,000
Single Board Computer £300
GPS £80
Solid-state Hard-disc £300
Sensors, SSR and Interface £200
Enclosure (inc. IR window) £400
Table 2.1: Approximate cost breakdown of thermal camera system.
Figure 2.2: The main components of the camera system: (A) FLIR SC645 IR-
camera, (B) BCT-RE2 single board computer, (C) Garmin 18 LVC GPS, (D) Intel
256 Gb solid-state hard-disc, (E) Phidgets solid-state relay.
19
BCT-RE2
U
S
B
R
S
2
3
2
GPIO
Ethernet
SC645
GPS SSD
SSR
LEDs
Phidget
Interface
Sensors
12 V
USB->Ethernet
To telemetry 
system
Figure 2.3: Block diagram showing the main components of the camera system.
20
2.2.1 Camera
The thermal camera unit itself is a FLIR SC645, with an uncooled 640×480 element
detector, spectral range of 7.5–13µm, and a spatial resolution (instantaneous field
of view) of 0.69 mrad. This unit was chosen because, at the time of purchase, it
was one of the few thermal cameras that allowed digital output of images in real-
time (rather than storing them on an on-board storage medium such as a compact
flash card). Furthermore, it supports the non-proprietary GenICam interface, which
allows capture parameters to be set on the camera and images to be received via
Ethernet.
The camera has three possible temperature ranges that it may be set to. In general
we use the range 273-923 K, however, we have occasionally used the higher temper-
ature range of 573-2273 K. The temperature range used for each image is stored in
the PNG header data of the image file.
The camera has a built-in calibration source and surface temperatures could, in
principle, be obtained per pixel. However, this would require parameterisation of
in-band emissivity of the lava lake surface and atmospheric transmittance. While a
reasonable estimate might be made for the former, the atmospheric transmittance
is highly variable in space and time due to the fluctuating column amounts of vol-
canic gases in the optical path between the camera and the lake surface (∼ 350 m
range; see Fig. 4.1). Several magmatic gas species emitted from the lake absorb in
the detector’s waveband. In an instantaneous image, the transmittance will likely
differ widely for each pixel due to the inhomogeneity of the volcanic plume. And
from frame to frame, the plume structure may vary considerably. Thus, no unique
value of transmittance would be appropriate to correct all the images to surface
temperatures. Sawyer and Burton [2006] have shown how use of open-path Fourier
transform spectroscopy to measure absorbing species in the optical path can pro-
vide some traction on this problem but only an imaging system would be able to
constrain both the spatial and temporal variability in transmittance (and it would
require a vast amount of data processing). Consequently, we have recorded images
in the “radiometric” setting on the camera - in which the recorded signal is linearly
proportional to in-band radiance reaching the detector. While these data are not
well suited to estimating surface temperatures on the lake or radiative power out-
puts, they are valuable for studying lava lake motion, as discussed in Section 2.4.3.
21
2.2.2 Control system
The control system comprises of an ARM-based single board computer (SBC) which
is connected to the camera via Ethernet. We used a Blue Chip Technology (BCT)
RE2 SBC. Although there is a huge range of SBCs on the market that cover
most combinations of interfaces, performance and power consumption, the choice
of boards with an extended operating temperature range is considerably more lim-
ited. Due to the high computational expense of image compression, we found in
testing that we needed a processor of at least 500 MHz. We also required multi-
ple USB host ports and several GPIO (general purpose input/output) pins. With
its 600 MHz ARM OMAP Cortex A8 processor, −40 °C operating temperature and
large range of interfaces, the BCT-RE2 met all our requirements at a power con-
sumption of 2 W. The system runs Ubuntu Linux (version 10.10). Additional kernel
modules had to be compiled to enable user-space access to the GPIO pins and also
to facilitate use of a USB to Ethernet converter.
2.2.3 Power supply
Between 2011–2013 the entire power system which serves the crater rim was replaced.
The old system was approaching ten years of service, and had become overly com-
plex, and unreliable. A combination of extreme winds and corrosive gases makes the
crater rim of Erebus an unsuitable site for solar panels and wind generators. Instead,
power is generated 0.5 km down-slope at the Nausea Knob (NKB) seismic station
site, where a 1000 Ah battery bank is charged using a ∼ 0.5 kW array of photovoltaic
panels, and a 100 W wind turbine. A Schaefer AEP-1500 inverter (chosen for its
extended temperature range) is then used to transmit 230 V through a heavy-duty
cable to the crater rim, where it is converted back to 12 V DC using an REL-185-1004
power supply (also chosen for its extended temperature range). The 230 V cable is
buried beneath the surface to help protect it from UV radiation and volcanic bombs.
Sufficient slack has been left in the cable that deformation of the ground where it is
buried (e.g. due to a bomb impact) will not pull it taught and potentially damage it.
Based on past experience of unfamiliar researchers causing damage to equipment
by connecting it to the wrong power supply, or with reversed polarity, the camera
system incorporates a custom made protection circuit. Reverse polarity protection
is provided by a MOSFET, and over-voltage protection is provided by a “crowbar”
circuit, which short-circuits a resettable fuse (PTC) at supply voltages above 15 V.
22
It is our belief that all electronic field equipment should be protected in a similar
way. Protection circuits are straightforward to make, cost almost nothing, and pre-
vent simple mistakes from totally destroying equipment. We cannot stress this point
enough!
Power to the camera itself is controlled via one of the SBC’s GPIO pins using a Phid-
gets 3053 solid-state relay (SSR). This allows the control software to “hard-restart”
the camera in case of problems.
Overall, the system requires ∼ 11 W of power. Between 6–7 W is used by the camera
itself, ∼ 2 W by the SBC and the remainder by the hard disc.
2.2.4 Time synchronisation
The majority of volcanological research requires the inter-comparison of data from
many different instruments, making accurate time stamping of the utmost impor-
tance. To ensure that the system time on the SBC did not drift, we linked it via
its RS232 serial port to a Garmin GPS 18 LVC receiver. The system time was then
synchronised with the GPS using standard Linux utilities. Initially we tried the
gpsd program, however, we found it to consume excessive system resources and so
we switched to the built-in GPS support in ntpd. We had hoped to use the Lin-
uxPPS kernel module to provide very accurate time synchronisation using the PPS
(pulse per second) signal output from the GPS. Unfortunately, the RS232 port on
the BCT-RE2 is implemented over the I2C bus, and the driver for this does not cur-
rently support LinuxPPS. Despite this, we found that synchronisation was possible
to within a few milliseconds just using the NMEA data from the GPS. The control
software for the camera has been designed to allow a sufficient delay, following a re-
boot, for GPS time synchronisation to occur before image capture resumes.
2.2.5 Sensors
The power consumption of the camera system, together with the temperature in-
side and outside of its enclosure are monitored using Phidgets sensors. These are
connected to the SBC via USB using a Phidgets 1018 2 interface board. Phidgets
sensors were chosen for their ease of use and simple software interface.
23
Dec 25 2012 Jan 08 2013 Jan 22 2013 Feb 05 2013 Feb 19 2013 Mar 05 2013 Mar 19 2013 Apr 02 2013 Apr 16 2013
Date (UTC)
40
30
20
10
0
10
20
30
40
T
e
m
p
e
ra
tu
re
 (
°C
)
Figure 2.4: Time series of temperature readings both inside (upper series) and
outside (lower series) of the camera enclosure. It should be noted that the outside
temperatures reported appear to be ∼ 10 °C higher than the actual air temperature.
The data collected by the sensors between December 2012 and April 2013 show a
stable power consumption of ∼11 W throughout. Figure 2.4 shows the recorded tem-
peratures. The sharp drops in internal temperature are caused by power outages,
and the subsequent rapid heating clearly demonstrates the ability of the system to
keep itself warm simply using waste heat from the SBC and hard disc. It should be
noted that the external temperature sensor seems to over-read by ∼10 °C compared
to measurements made with a handheld thermocouple. This is likely caused by heat
leakage from the camera enclosure into the external temperature sensor enclosure.
2.2.6 Data storage
The low temperature performance of hard disc drives (HDDs) seems to be poorly re-
searched and/or documented. When choosing a data storage solution for the camera
system we were unable to find any high-capacity discs with an extended temperature
range. In the end we chose a 250 Gb Intel solid state drive (SSD). While it was not
rated for low operating temperatures, it did have a lower power consumption than
conventional HDDs. Some consideration was given to the idea of storing batches
of images on the SBC’s internal flash memory and then periodically transferring
them to the SSD. Although this would have reduced the power consumption of the
24
ABC
D
E
F
Figure 2.5: The main components of the new network switch installed at the crater
rim in 2011: (A) Power-over-Ethernet (POE) port for TL-45, (B) 12 V power input,
(C) weather-proof Ethernet ports, (D) 12 V to 24 V DC/DC converter to power TL-
45, (E) Parvus PRV-1059 Ethernet switch, (F) POE power injector.
system, since the SSD would only need to be powered up occasionally, we decided
that the risk of data loss due to the increased complexity outweighed the potential
benefits.
2.2.7 Telemetry
During the annual MEVO field campaign, a Trango TL-45 5 GHz wireless Ethernet
link is used to telemeter data from instruments at the crater rim back to servers at
our field camp. Following numerous failures of the system due to cold-damage of ca-
bles and connectors, the whole network system surrounding the TL-45 was replaced
in 2011. In the new system, instruments are connected to the TL-45 via a low-
power, low-temperature Parvus PRV-1059 Ethernet switch using low-temperature
ethernet cable and Amphenol—weather-proof connectors (Fig. 2.5). The high power
consumption of the telemetry system (∼ 12 W) means that it is only operated during
the field campaign (when there is ample solar charging). The thermal camera sys-
tem incorporates a USB to Ethernet converter to provide it with a second Ethernet
interface and allow connection to the telemetry system during this period. Images
can then be streamed to the field camp in addition to being stored on the SSD,
providing some redundancy in data storage and the ability to monitor the lava lake
in real-time.
25
Figure 2.6: The new thermal camera system installed at the crater rim of Erebus
volcano. Photo: Clive Oppenheimer
2.2.8 Enclosure
Previous thermal cameras operated by MEVO were simply wrapped in bubble-wrap
to protect them from the cold. Somewhat surprisingly this proved to be very effec-
tive. However, the corrosive gases emitted from the volcano rapidly took their toll,
and metal parts (including electrical contacts) soon degraded and became unusable.
An additional problem with previous systems was that standard tripod mounts were
not rigid enough to prevent considerable camera shake caused by the high winds
which frequent the summit of Erebus. An air-tight plastic Peli— case was used as
the enclosure for the new system, to which a rigid aluminium frame was bolted.
The frame was attached to a tripod comprised of scaffolding bars to form a solid
mounting which does not vibrate in the wind (Fig. 2.6). The IR window used in the
camera enclosure is 3 mm thick Germanium with a 8-12µm anti-reflective coating.
All external fixtures are either aluminium or stainless steel to reduce corrosion. The
camera casing has been insulated with expanded polystyrene and all internal fix-
ings are made of plastic to reduce conductive heat loss. Combined with the heat
generated by the enclosed electronics, this meant that the temperature inside the en-
closure remained above -5 °C throughout its operation in 2013 (except during power
outages), despite external temperatures of less than −50 °C (Fig. 2.4).
26
2.3 Software
2.3.1 Overview
The control software for the camera system is comprised of several small programs
which are managed by an initialisation script. All programs share a single configu-
ration file, which can be edited to change the operational parameters of the system.
The main capture program is responsible for setting the capture parameters of the
camera, and then entering an infinite loop of capturing images, performing some
basic pre-processing on them and then compressing them as PNG files. During the
field campaigns when the telemetry system is in place, two separate server programs
are used to send images and sensor data to a viewer program which runs on a com-
puter at the field camp. This allows real-time monitoring of the lava lake and of the
thermal camera system itself.
2.3.2 Implementation and architecture
The initialisation script is written in Bash. It is launched automatically on system
boot by upstart (a replacement for the System V init daemon) once the networking
and GPIO services have been started. The script performs a series of system checks
such as ensuring that the SSD is mounted correctly, initialises GPS time synchro-
nisation and then starts the capture and server programs. If errors are encountered
that would prevent the capture of images, then the whole system is restarted. An
error handling safety net is provided by the kernel watchdog module, which will
restart the system if the capture software crashes for some reason or if the system
“freezes” (runs out of memory, or has excessive CPU usage). As a final layer of fault
tolerance, a crontab entry is used to restart the whole system once a week.
The main capture program is written in C++ and uses libconfig++ for reading
the configuration file, libpng for compressing and storing the images, and libaravis
[Pacaud, 2011] for interfacing with the camera. The images output from the camera
have a bit depth of 16. To reduce the size of the PNG files produced, they are
scaled to a bit depth of 8 before being compressed. In order to preserve as much
dynamic range in the pixels representing the lava lake as possible, the images are
first thresholded to set all low valued pixels to zero. They are then linearly scaled
27
to 8 bit values as follows:
x8 = (x16 −min16)× 255
max16 −min16 (2.1)
where x8 is the 8 bit pixel value, x16 is the 16 bit pixel value and min16 and max16
are the minimum and maximum pixel values in the 16 bit image respectively.
Converting the images to 8 bit using this algorithm is sensitive to dead or hot pixels,
since these will skew the minimum and maximum pixel values. This is potentially
disastrous, as the dynamic range of pixels within the lake in the resulting 8 bit im-
ages could be severely impacted. To prevent this from happening, a simple dead/hot
pixel detection algorithm is applied to each image before conversion to 8 bit. This
works by comparing each pixel to its neighbours and looking for very strong con-
trasts in value. In order to implement this algorithm efficiently, it ignores the edge
pixels of the image, meaning that every pixel analysed has eight neighbours. Edge
pixels are simply set to zero. Image files are written to disc using incremental num-
bers as file names. This method was chosen as it will not result in images being
overwritten in the event of unexpected system restarts (e.g. due to power cuts) or
loss of system time.
The server programs are written in C. The server for the sensor data uses libphidgets
to interface with the sensors. Values read are recorded in a log file and also passed
to any connected clients via a socket. The server program for the images works by
using inotify to monitor the output directory for the images from the capture pro-
gram (which is read from the capture program’s configuration file). Each new image
detected is transmitted to any connected clients as a string of raw PNG data, and
the user-gpio kernel module is used to flash a status LED connected to one of the
GPIO pins of the SBC. The server programs were written as separate programs from
the capture program in order to keep the latter as simple, and therefore reliable, as
possible.
The viewer program is written in Python and connects to the server programs via
a standard socket. Multiple instances of the viewer can connect to the server pro-
grams simultaneously.
28
2.3.3 Reuse potential
All the software created for this project is freely available under the terms of the
Gnu Public License from http://dx.doi.org/10.6084/m9.figshare.784942. Al-
though the image acquisition software is specific to the SC645 camera, fairly minor
modifications would enable it to work with any camera that supports the GenICam
interface. It is not however, suitable for applications that need capture rates above
a few frames per second. The image server and viewer programs could be reused
“as is” for any image capture application that outputs PNG files, provided that the
server side system has inotify.
2.4 Results and discussion
2.4.1 The camera system
The system was installed at the crater rim of Erebus in November 2012 and, other
than a few minor issues with the camera being unable to focus to infinity when it
was first mounted in the enclosure, it operated without any problems up until the
end of the 2012 field season. It was then left running with the intention that it con-
tinue to acquire images continuously until the following field season. On returning
to the volcano in November 2013, we discovered that the battery bank which serves
the camera system had failed. Several of the batteries had cracked casings and all
registered ∼ 6 V. The exact cause for this failure was never determined, however,
we suspect that it was caused by a faulty wind turbine system which allowed the
batteries to over-discharge and subsequently freeze. As a result, power was never
restored to the camera after the winter, and therefore images were only acquired up
until 26 April 2013. It is worth noting however, that when power was restored, the
system resumed acquisition without any further intervention. The wind turbine and
associated controller have now been replaced with a new, redundant system which
we hope will provide a greater degree of protection for the battery bank.
A further failure of the system was that the stainless-steel mounting bolt which se-
cures the camera enclosure to the tripod had bent, allowing the enclosure to rotate
such that the camera was no longer pointing at the lake. Fortunately, this had hap-
pened after the power had failed, and so no additional data were lost. It is somewhat
surprising that such a small enclosure can provide enough wind-resistance to bend
a 5 mm diameter bolt. We have replaced the bolt with a much larger one, and plans
29
2013
bomb
c
a
b
le
Figure 2.7: Photograph showing the landing site of a large (∼1.5 m diameter) vol-
canic bomb which was erupted in 2013 and landed directly on top of the 230 V cable
which takes power from Nausea Knob to the crater rim. Luckily the cable seemed
to be undamaged.
are underway to replace the entire tripod setup with a more substantial mounting
system during the 2014 field campaign.
Despite the problems outlined above, the camera system worked well until late April
when the sun stopped rising and solar charging ended. Almost four million images
were acquired; the largest thermal infrared dataset ever recorded on Erebus. The
system was left running again after the 2013 field season and, with the improved
power system in place, we are hopeful of achieving year-round data acquisition dur-
ing 2014. However, the increased level of explosive activity which began in 2013
and continued until late-February 2014 (Section 1.6) has added significantly to the
risk of terminal failure of the system. During this period, large bombs, up to 4 m in
length, were being erupted from the crater two to three times a day. One of these
landed within a few metres of the camera system itself, and several had struck the
cable that carries power to the crater rim (Fig. 2.7). Fortunately, the cable appeared
undamaged by the impacts – rewarding the huge efforts that were taken to bury it!
We have no way to tell if the camera system is still operational until we return to
Erebus in November 2014.
2.4.2 The dataset
Figure 2.8 shows a short sequence of images captured by the camera system dur-
ing testing in 2010, which document a large gas slug arriving in the lake. These
30
Figure 2.8: A sequence of three IR images recorded ∼ 0.17 s apart, showing the
arrival of a large gas slug into the lava lake.
images were captured at 6 Hz, whereas the camera system ordinarily operates at
0.5 Hz. The dataset from the new camera system is freely available at http:
//www.usap-data.org/entry/NSF-ANT11-42083/. All images in the archive are
in PNG format, and contain important meta-data in their file headers. Information
about extracting these meta-data can be found in the README file which is part
of the archive. At the time of writing the archive contains all the images that were
captured by the camera system during the 2012 field season. This will be expanded
upon annually as new data are recovered from the camera. Older datasets from
previous camera systems will also be made available once the images have been
converted into a suitable archive format. It should be noted that the images are
obtained with an oblique viewing geometry, and require orthorectification before
they are used. This has not been applied to the images in the archive. Techniques
for correcting the images are discussed in Section 4.2.3.
2.4.3 Applications
As discussed in detail in Chapter 4, the Erebus lava lake is covered by a thin crust
of cooler material. As the lake convects, this crust cracks and the cracks migrate
across the surface driven by the motion of the underlying lava. Figure 2.1 shows
a comparison between the visible appearance of the lake and an IR image (it should
be noted that such a clear view of the lake is exceptionally rare). As can be seen
in the figure, the cracks are clearly visible in the IR image and, as demonstrated by
Oppenheimer et al. [2009], can be tracked using a wavelets based motion estimation
algorithm [Magarey and Kingsbury, 1998] to determine the surface velocity of the
lake. Figure 2.9 shows a short time series of mean lake surface speed calculated in
this way. The pulsatory behaviour of the surface motion is thought to be associ-
ated with a density-driven bi-directional flow of magma in the conduit feeding the
lake [Oppenheimer et al., 2009]. A primary goal of the new thermal camera sys-
tem is to provide long, uninterrupted observations of the pulsatory behaviour of the
31
23:36:00 23:46:00 23:56:00 00:06:00 00:16:00 00:26:00 00:36:00 00:46:00
Time (UT)
0.01
0.02
0.03
0.04
0.05
0.06
0.07
0.08
0.09
0.10
M
e
a
n
 l
a
ke
 s
p
e
e
d
 (
m
s−
1
)
Figure 2.9: Time series of mean lake surface speed calculated from IR images cap-
tured on 22 December 2010. The pulsatory behaviour of the lake is clearly visible
with a period of ∼ 12 min. The sharp spikes in speed (“shot noise”) are caused by
bubbles reaching the surface of the lake.
lake to facilitate a detailed comparison with numerical simulations [Molina et al.,
2012], analogue experiments [Witham et al., 2006, Huppert and Hallworth, 2007]
and datasets from other instruments deployed at Erebus. The latter is covered in
Chapter 5.
The episodic arrival of large gas bubbles into the lake cause spectacular Strombolian
eruptions (Fig. 2.8). Although detailed studies have been made of these events using
other techniques [Gerst et al., 2013, Jones et al., 2008], little has yet been learnt of
how they influence the pulsatory behaviour of the lake. Previous thermal camera
systems on Erebus did not record the majority of the bubble events since their cap-
ture rate was too slow. However, the new system will allow explosive episodes to be
recorded. Explosions can easily be identified in the data either as “shot noise” in
the calculated surface motion or, for large events, as a sharp (factor of ∼ 5) increase
in the compressed size of the image files. Although some analysis of bubble events
is presented in Chapter 4, sufficiently long datasets were not available to conduct an
in-depth study at the time of writing. This is therefore left as a promising avenue
of research for future studies.
32
Chapter 3
Other Volcanic Monitoring
Hardware and Software
The work presented in this chapter was undertaken as part of the “Improved Mon-
itoring” research objective of this study (Section 1.7.1). Although the outputs have
not featured significantly in other sections of this dissertation, they do constitute
a large body of work, which is relevant to the monitoring of open-vent volcanoes
more generally. As such, we feel that it is necessary to present them here. However,
this chapter should be viewed as a summary of the work, rather than an in-depth
analysis. We present new data visualisation software, which, in addition to provid-
ing generic plotting capabilities, can be customised using plug-ins to provide very
specific data processing capabilities such as real-time SO2 measurements. We also
present portable spectrometer scanning units for making spectroscopic measure-
ments at remote field sites.
Section 3.1 of this chapter has been published as: Peters, Nial. “AvoPlot: An Ex-
tensible Scientific Plotting Tool Based on Matplotlib”. Journal of Open Research
Software 2, no. 1 (February 10, 2014).
All the hardware and software described in the chapter was created independently,
by me.
3.1 Data Visualisation Software - AvoPlot
3.1.1 Introduction
AvoPlot is a simple-to-use graphical plotting program written in Python [van Rossum,
2013] and making extensive use of the matplotlib plotting library [Hunter, 2007].The
software was originally created as a real-time visualisation program for monitoring
sulphur dioxide emissions from volcanoes. However, with several different projects
under way in our research group, all requiring a different set of data manipulation
tools but similar user interfaces, it soon became clear that a more generic solution
would be extremely useful; AvoPlot was designed to fill this role. Essentially, Avo-
Plot aims to achieve two things: (i) To provide a simple and intuitive graphical user
interface for the matplotlib plotting library, with tools for importing common data
formats. (ii) To provide an extensible platform that researchers can easily customise
to suit their specific data visualisation and analysis requirements.
Although there are numerous open-source plotting programs already available, none
(that the author is aware of) provide the same functionality as AvoPlot. Spread-
sheet programs such as OpenOffice Calc [Apache, 2013] and Gnumeric [Gnome,
2013] cannot load arbitrary data types, and also lack good data navigation tools
(panning, zooming etc.) for their plots. The plot viewer built into matplotlib, pro-
vides good data navigation tools, but has very few graphical facilities for editing plot
parameters. The Veusz [Sanders, 2013] program incorporates most of the features of
AvoPlot. It has a comprehensive graphical interface for creating and editing plots,
as well as providing a plug-in interface to allow users to load arbitrary data types
and define new analysis functions. However, Veusz implements its own plotting
backend which is unlikely to be as familiar to developers as matplotlib and has a
smaller community providing support. Furthermore, the scope for writing complex
plug-ins is far more limited than in AvoPlot, which allows developers the freedom to
use any wxPython [Dunn, 2003, Rappin and Dunn, 2006] constructs when writing
their plug-in interfaces.
At the time of writing, AvoPlot is still at an early stage of development. Although
two stable versions have been released, the available functionality is somewhat lim-
ited. Despite this, it has already been used in several different volcanological projects
and was used extensively by the author during the preparation of this thesis. Plug-
ins for working with differential optical absorption spectroscopy (DOAS) data and
Fourier transform infrared (FTIR) spectra are currently under development and will
34
be released in due course.
3.1.2 Implementation and architecture
AvoPlot is written in pure Python. Graphical interface elements are implemented
using wxPython and plots are created and displayed using matplotlib and its wx-agg
backend.
The software consists of two major components; the user interface (UI) and the plug-
in system. The UI provides a generic graphical interface to much of matplotlib’s
functionality, allowing users to change plot parameters (such as line styles, colours
etc.) without needing to write any code. It also allows users to interact with any
plug-ins that have been installed. An overview of the UI is shown in Fig. 3.1.
The plug-in system is intended to give developers a simple way both to import their
specific data into the UI, and to provide tools for manipulating them. Details of
how to write plug-ins are given in the AvoPlot documentation which can be found in
the “doc” folder of the source distribution, or online as part of the software archive.
Furthermore, example plug-ins can be found in the “examples” folder of the source
distribution.
AvoPlot comes with a simple file import plug-in which will be installed automati-
cally. Currently, this can only load text files containing columns of numerical data.
However, it will be extended in future releases to support other file types such as
CSV, Excel spreadsheets etc. Sample data for creating a plot with this plug-in
are provided in the “spectrum data.txt” file of the “examples” folder of the source
distribution and step-by-step instructions for its use can be found in the AvoPlot
documentation.
3.1.3 Reuse potential
As a standalone program, AvoPlot provides an easy to use, generic plotting tool
capable of producing publication-quality figures. Although the current feature set
is relatively limited, future development will focus on exposing more of matplotlib’s
functionality through AvoPlot’s UI. This will give users access to matplotlib’s ex-
tensive capabilities without requiring them to write Python scripts. In this respect,
35
Figure 3.1: Screen-shot of the AvoPlot user interface showing its main components:
A) Toolbar providing the “standard” matplotlib pan, zoom and save controls. B)
Navigation panel allowing selection of different plots, axes and data series. C)
Control panel providing controls for editing plot properties such as axis titles, line
styles etc. Any tools provided by the plug-in for the currently selected plot are
also displayed here under a different tab. D) Tabbed plot window with a variable
number of panes allowing multiple plots to be viewed simultaneously.
36
the program has a high reuse potential “as is” for a wide spectrum of end-users
regardless of coding experience.
However, the main reuse potential for AvoPlot comes from its plug-in based ar-
chitecture. With just a few lines of relatively simple code, researchers can write a
plug-in that will plot their specific data set (which may be in a format difficult to
import with mainstream plotting programs) in the AvoPlot UI. They can then take
advantage of a feature-rich plot editing environment to customise their plot without
having to edit or run scripts. Furthermore, AvoPlot plug-ins are distributed using
Python’s distutils system [Ward and Baxter, 2013], making it simple for other re-
searchers (who may have no programming experience) to install plug-ins written by
their peers.
In addition to data import, AvoPlot plug-ins also allow researchers to easily inte-
grate their own data manipulation and analysis tools into the AvoPlot UI, enabling
them to create a graphical data processing software without having to write exten-
sive UI code of their own. The aim is to encourage researchers to make the analysis
tools which they create user-friendly and accessible for others, who may not be pro-
grammers.
The software has been successfully tested (functional testing) on Linux (Ubuntu
and OpenSUSE), Microsoft Windows (XP, 7 and 64bit 8) and MacOSX. There is
no reason however, why AvoPlot cannot be installed on any system that can run
Python.
3.1.4 Volcanological Applications
As already discussed, AvoPlot is designed to be a generic scientific plotting tool that
can be specialised to suit particular applications through the use of plug-ins. Within
a volcanological context, there are currently two such plug-ins under active devel-
opment. The first is for analysing FTIR spectra, and the second is for performing
real-time estimates of SO2 emissions using UV spectrometers.
The FTIR plug-in aims to automate the background subtraction tasks which are
currently performed manually. Very early versions of the plug-in were used in the
study conducted by Iacovino et al. [2014] and development work continues to create
a more advanced and user-friendly tool [Iacovino et al., 2013].
37
14:54:00 14:59:00 15:04:00 15:09:00
Time (UTC)
6
7
8
9
10
11
12
13
I 3
35
n
m
/I
3
10
n
m
Figure 3.2: Plot produced using the AvoPlot SO2 plug-in to visualise UV spectra
collected with a scanning instrument at Villarrica volcano in February 2012. The
peaks correspond to in-plume measurements as the spectrometer is scanned back
and forth in a transect of the plume.
The SO2 emission plug-in is designed to provide a faster and more advanced replace-
ment for the Vispect software [Wright, 2008], which is currently the “standard” tool
used by volcanologists for UV spectroscopy in the field. The principal of this soft-
ware is to provide real-time estimations of the column amount of SO2 in the field of
view of the spectrometer, thus enabling researchers to rapidly assess the data they
are collecting. Such a tool is vital for effective data collection, since it is otherwise
impossible to tell if the field of view of the spectrometer is intersecting the plume.
This is especially true for data collected during traverses beneath the plume [e.g.
McGonigle, 2002], and for scanning units such as those described in Section 3.2. Vis-
pect achieves this by plotting a ratio of measured intensities at wavelengths outside
and inside an absorption band of SO2 (typically 335 nm and 310 nm respectively), a
peak in the ratio represents increased amounts of SO2. Although Vispect is a useful
tool and has been widely adopted by researchers around the world, it has several
key limitations. The procedures it uses for loading spectral data are extremely slow,
making the assessment of large data sets impossible. Furthermore, the approach used
to detect newly written spectrum files is flawed, and causes the program to crash
once more than a few thousand spectra have been collected. The current implemen-
tation of the AvoPlot plug-in offers the same functionality as Vispect (Fig. 3.2), but
uses optimised spectrum loading routines which can run concurrently on multiple
CPUs, resulting in processing times that are over fifty times faster than those of
Vispect (even on a single-core CPU). System-level directory-watching utilities (in-
38
otify in Linux and ReadDirectoryChangesW from the native API in Windows) are
used to detect newly created spectrum files, ensuring fast, efficient detection that
scales to any number of spectra. The plug-in has been used extensively during field
campaigns on Erebus as well as during spectroscopic studies of volcanoes in Costa
Rica [Moussallam et al., 2014b], Chile (unpublished) and Mexico [Cassidy et al.,
2014].
In addition to the studies making use of these specific plug-ins, the main AvoPlot
program itself has also been used in several projects. Its fitting functions, have been
used for performing spectrometer calibrations [Section 5.2.2, Peters et al., 2014b]
and for analysing sulphur data [Iacovino, 2014], and its general plotting capabilities
have been used to produce many of the figures in this dissertation.
3.2 Portable Spectrometer Scanners
3.2.1 Introduction
SO2 flux is an important parameter in volcanic monitoring and is commonly mea-
sured using UV spectrometers. In order to calculate flux, spectra must be collected
across a full transect of the plume perpendicular to its transport direction. For
horizontally drifting plumes, this is often achieved by carrying the spectrometer
(with its field of view at zenith) in a traverse beneath the plume [e.g. McGonigle,
2002]. However, this technique requires there to be an unobstructed path beneath
the plume, and cannot be used for vertically rising plumes. In these situations scan-
ning units are often used instead [e.g. McGonigle, 2007]. Although different designs
for scanning units exist, their basic operation is the same; a stepper motor is used
to rotate the the optical axis of the spectrometer’s field of view such that spectra
can be collected across a full transect of the plume.
Recently, much emphasis has been put on developing scanning units for permanent
installation around volcanoes [e.g. Salerno et al., 2009, Edmonds et al., 2003]. As
such, the scanners themselves tend to be large, heavy and require an external power
supply. This makes them unsuitable for studies at inaccessible volcanoes where
equipment must be carried on foot. The scanning units presented here were specif-
ically designed with this in mind and as a result are light, cheap and low-power.
39
Figure 3.3: Screen-shot of the main interface window for the scanner control pro-
gram.
3.2.2 Hardware and Software
The software for the scanners is written in Python, using wxPython for the user-
interface. It is split into three main components: (i) the control program for the
scanners, (ii) a script for assigning angles to spectrum files and (iii) ica inspect
– a graphical program for calculating fluxes from processed spectral data. All
the software is available open-source under the terms of Gnu Public License from
http://code.google.com/p/avoscan/.
Figure 3.3 shows a screen-shot of the main interface window of the control program
for the scanners. Start- and end-angles for scanning are set relative to a user de-
fined zenith. For vertical scans, this should usually be straight-upwards, which can
be obtained accurately using a spirit-level placed on the top of the spectrometer’s
telescope. Two types of scan are possible, “continuous” in which the stepper motor
is always in motion and “stepped” in which the stepper pauses for a user-defined
amount of time at a certain number of increments between the start- and end-angles.
In both cases the stepper reverses direction at the end points of the scan. Although
this introduces a significant error into the recorded angles (due to the slack in the
gearbox), it does not affect the accuracy of SO2 flux calculations (see Section 3.2.3)
and is a necessary compromise to avoid having to use a mirror for scanning. In
addition to setting scan parameters, the control program allows field notes to be
40
recorded in the header data of the angle-versus-time file, ensuring that the notes do
not become separated from the data to which they refer.
The ica inspect program was the original basis for the AvoPlot software [Section 3.1,
Peters, 2014] and has also been used, in a slightly modified form, for processing
XANES (X-ray Absorption Near Edge Structure) data [Moussallam et al., 2014a].
With respect to the scanning units, it is used to calculate fluxes from the processed
spectral data. Note that it does not perform the actual SO2 column amount re-
trievals, that must be done separately (see Section 5.2.2 for a description of the
retrieval process). Rather, ica inspect takes the retrieved SO2 column amounts and
their corresponding angles and integrates them (see Eq. 3.2 below) to give a flux for
each complete scan. A key step in this process is the subtraction of the background
(i.e. the out-of-plume) SO2 column amount. The distinction between in-plume and
out-of-plume regions is estimated by fitting a Gaussian to the data from each indi-
vidual scan. Out-of-plume column amounts are taken to be any that are more than
three standard-deviations from the mean of the Gaussian fit. Although this gives
a reasonable estimate, some manual adjustment of the bounds is usually required.
The background is then approximated by a linear fit through the out-of-plume data.
However, it is common for these data to contain many outliers, for example caused
by the ground entering the field of view of the spectrometer at the extrema of the
scan, and these must be removed so as not to skew the background estimation. The
background is therefore calculated using an iterative approach; an initial linear fit
is made through all the data, points that could not be fitted within their associated
error (where the error is determined by the retrieval software) are identified and
the furthest outlier is removed from the fit, the process is then repeated until all
remaining points can be fitted within error. Despite the simplicity of this approach,
it has been found to produce very robust background estimates, even when the data
are exceptionally noisy.
The main hardware components of the scanning units are shown in Fig. 3.4. Even at
fast scanning speeds, power consumption is sufficiently low that the scanners can be
powered via their USB interface. Since the majority of spectrometers used for mea-
suring volcanic SO2 are also USB-powered, this is a huge advantage; removing the
need to carry heavy batteries into the field. The power supply board includes reverse-
polarity protection (as described in Section 2.2.3), with a standard fuse providing
over-current protection. The stepper-motor and controller are standard components
produced by Phidgets (part numbers 3319 1 and 1063 1 respectively). It should be
41
AB
C
D
E
F
G
Figure 3.4: Photograph showing the key hardware components of the scanning units:
(A) Phidgets 1063 1 stepper controller, (B) Phidgets 3319 1 stepper motor, (C)
Power supply (5–12 V DCDC converter, reverse polarity and over-current protec-
tion), (D) CMPS10 tilt-compensated compass, (E) FT-MOD-4232HUB USB-Serial
hub, (F) Tripod mount, (G) Waterproof USB ports.
noted that although these components do not have an extended operating temper-
ature range, the scanners have been used, without fault, for many hours at −30‰
during field campaigns on Erebus. The FT-MOD-4232HUB USB-Serial hub is used
to provide additional USB-host ports on the scanner which can be used, for example,
for connecting the spectrometer. It also provides four serial ports. One of these is
exposed on the back of the scanner to allow a handheld GPS unit to be attached
for accurate time-keeping and automatic location logging (although this is not yet
supported by the software). Another is exposed to allow connection of a web-camera
which could collect a visual record of each scan, allowing rapid assessment of plume
conditions in post-processing (also not yet supported by the software). A serial port
is also used to connect a three-axis, tilt-compensated CMPS10 electronic compass.
This allows the software to automatically record the angle from horizontal that the
scan was made at (for horizontal scans), and the bearing of the scan plane (for
vertical scans). All the electronics are sealed in a waterproof box, and waterproof
connectors are used to ensure that dust, ash or moisture cannot enter.
42
3.2.3 Measuring SO2 Flux
Scans are made by mounting a telescope (coupled to the spectrometer by an optical-
fibre) to the scanner. The spectrometer records continuously (the control software
for the spectrometer is independent of that of the scanner) as the scanner rotates
back and forth. Angles are assigned to each spectrum in post-processing using their
timestamps and an angle-versus-time log file created by the scanning software. The
angle (θτ ) associated with a spectrum captured at time τ is given by:
θτ = θ0 + (τ − τ0) · dθ
dt
(3.1)
where dθ
dt
is the angular velocity of the scanner, τ0 is the time record closest to τ in
the angle-versus-time log file and θ0 is the corresponding angle.
Using the ica inspect program, processed spectra are split into individual scans
based on their corresponding angles and each scan is manually divided into in-
plume and out-of-plume regions. A linear fit through the out-of-plume regions is
used as an estimate of the background and subtracted from the scan. SO2 flux (Φ) is
then calculated using the standard scanning DOAS integration method [e.g. Salerno
et al., 2009]. For a scan containing n spectra, the flux in tonnes per day is given by:
Φ = κυ
n−1∑
i=0
r (θi+1 − θi) · Ai + Ai+1
2
cosα (3.2)
where κ is the conversion factor from ppm m3 s−1 to tonnes per day (2.3× 10−4), υ
is the plume transport speed perpendicular to the scan direction, r is the distance
from the spectrometer to the plume, Ai is the SO2 column amount recorded at scan
angle θi and α is the inclination angle of the scan from horizontal.
The stepper motor in the scanners has a 99:1 planetary gearbox, giving it an accuracy
of ± 0.01 ◦. Note that errors due to crash-back in the gearbox are systematic for a
single scan, and therefore do not affect the integrated column amount. However, the
precision of the timestamps in the spectrum files (1 second) does have an effect on
the accuracy of the angles. Using Equation 3.1 and the standard error propagation
equation [e.g. Boas, 1983a] the error in the angle (δθτ ) of a spectrum recorded at
time τ can be calculated by:
δθτ = |θτ | ·
√
δ2θ +
(
(τ − τ0) · dθ
dt
)2
·
(
2δ2θ +
δ2τ
(τ − τ0)2
)
(3.3)
43
where δθ = 0.01
◦ is the uncertainty in the angle of the stepper motor and δτ = 1 s
is the uncertainty in the capture time of the spectrum.
Error in the retrieved SO2 column amount is usually estimated by the retrieval
software by analysing the residual from the fitting of the reference and Ring spectra.
Propagating these errors along with those in the angles through Equation 3.2 gives
the following as an estimate of the error in the flux (δΦ):
δ2Φ = Φ
2
[(
κυδICA
Φ
)2
+
(
δυ
υ
)2]
(3.4)
where δυ = 0.3 m s
−1 is the uncertainty in the plume rise speed and δICA is the error
in the integrated column amount given by:
δ2ICA =
n−1∑
i=0
A2i


√
δ2θi + δ
2
θi+1
θi+1 − θi
2 +

√
δ2Ai + δ
2
Ai+1
Ai+1 + Ai
2
 (3.5)
As can be seen in Figure 3.5, the resulting errors are rather large (in the order of
± 50%). The error is dominated by the uncertainty in the angle for each spectrum,
which is large due to the single second precision of the timestamps used in the spec-
trum files. Improvements in the software have since resolved this issue by using
the same directory watching code as AvoPlot (Section 3.1.4) to achieve millisecond
precision.
3.2.4 Results
Figure 3.5 shows a comparison of SO2 fluxes recorded at Villarrica volcano in Febru-
ary 2012 using the scanning units described above, a UV camera and pedestrian
traverses beneath the plume. The data shows that despite the very different time
resolutions of the three techniques, the calculated fluxes are comparable. These
data were collected as part of a study into using UV spectrometers to calibrate
images from UV cameras. Unfortunately, the same methodology was developed in-
dependently by Lu¨bcke et al. [2013] and so the study was abandoned. However, the
scanning units have since been used during several field campaigns on Erebus and
studies at Turrialba volcano in Costa Rica [Moussallam et al., 2014b] and Colima
volcano in Mexico [Cassidy et al., 2014].
44
14:00 15:00 16:00 17:00 18:00 19:00
Time (UTC)
0
50
100
150
200
250
300
350
S
O
2
 F
lu
x
 (
to
n
n
e
s 
p
e
r 
d
a
y
)
UV Camera
DOAS Traverses
DOAS Scans
Figure 3.5: SO2 flux data from Villarrica volcano, collected on 10 February 2012.
The data show a comparison of the SO2 fluxes calculated from UV camera images,
UV spectrometer scans (DOAS scans) and UV spectrometer pedestrian traverses
(DOAS traverses).
45
Chapter 4
Decadal Persistance of Cycles in
Lava Lake Motion
The work presented in this chapter was undertaken in order to meet the “Assess-
ment of Annual Variability of the Lava Lake Cycles” research objective of this study
(Section 1.7.2). Although previous studies of Erebus volcano’s active lava lake have
shown that many of its observable properties (gas composition, surface motion and
radiant heat output) exhibit cyclic behaviour with a period of ∼10 min, they have
typically focused on short data sets (a few hours). Subsequently, little is known
about the longer-term behaviour of these cycles. In this chapter we investigate the
multi-year progression of the cycles in surface motion of the lake using an extended
(but intermittent) dataset of thermal infrared images collected between 2004 and
2011. Cycles with a period of ∼5–18 min are found to be a persistent feature of
the lake’s behaviour and no obvious long-term change is observed despite variations
in lake level and surface area. The times at which gas bubbles arrive at the lake’s
surface are found to be random with respect to the phase of the motion cycles, sug-
gesting that the remarkable behaviour of the lake is governed by magma exchange
rather than an intermittent flux of gases from the underlying magma reservoir.
The majority of content from this chapter has been published as: Peters, Nial, Clive
Oppenheimer, Philip Kyle, and Nick Kingsbury. Decadal Persistence of Cycles in
Lava Lake Motion at Erebus Volcano, Antarctica. Earth and Planetary Science
Letters 395 (June 1, 2014): 1–12.
All content of this chapter was created by me, and the role of the co-authors did not
extend beyond that of normal supervision. Terrestrial laser scan data were collected
and pre-processed by Jed Frechette, Laura Jones and Drea Rae Killingsworth. How-
ever, subsequent calculation of lake surface areas were performed by me. Computer
codes implementing the motion estimation algorithm were provided by Prof. Nick
Kingsbury, but all other code was created by me. My implementation of the Morlet
wavelet transform was, however, closely based on code provided by Marie Boichu.
4.1 Introduction
As discussed in Section 1.3, many of the properties of the Erebus lava lake exhibit
a pronounced pulsatory behaviour. Oppenheimer et al. [2009] observed that the
radiative heat loss, surface velocity and certain magmatic gas ratios all oscillated
with a period of ∼10 min. The cycles appeared to be phase locked with each other,
suggesting a common mechanism was responsible for the oscillations in each prop-
erty. Evidence of similar cyclicity has been observed in the SO2 flux [Boichu et al.,
2010] and this is explored in more detail in Chapter 5. Fluctuations in the H2/SO2
ratio have also been measured [Moussallam et al., 2012], but these have yet to be
linked definitively to the cycles observed by Oppenheimer et al. [2009].
One possible explanation for the observed behaviour is pulsatory exchange flow of
hot, degassing magma into the lake from the subjacent conduit. It has been shown
experimentally that given two liquids flowing in opposite directions in a vertical pipe
(for example driven by a density difference between them), under certain flow con-
ditions an instability occurs which results in a pulsed flow [Huppert and Hallworth,
2007]. Oppenheimer et al. [2009] suggested that such behaviour may explain the cy-
cles at Erebus volcano, with bubbly and degassing, low density magma rising up the
conduit into the lake whilst degassed, denser magma sinks back down the conduit
again. The resulting pulsatory flow delivers packets of fresh magma into the lake
quasi-periodically, giving rise to the observed cycles in lake properties. The period
of the cycles would be expected to reflect the rheological properties and velocity of
the bubbly flow and geometry of the conduit.
The previous studies at Erebus have analysed only very short time-series of data
(circa 1 h), and no investigation of the long-term behaviour of the cycles has yet
been conducted. However, thermal infrared (IR) images of the Erebus lava lake
have been collected almost every year since 2004 during the Mount Erebus Volcano
Observatory’s annual austral summer field campaigns. Using a similar technique to
that of Oppenheimer et al. [2009], we have extracted mean surface speed estimates
from the usable portions of the now substantial IR dataset. Using the mean surface
speed as a proxy to assess the cyclicity of the lake motion, we present an overview
of its behaviour between 2004 and 2011 and compare this to visible changes in the
lake’s appearance. Using a dataset recorded at higher time resolution in 2010, we
identify times when bubbles arrive at the surface of the lake and compare this to
the phase of the cycles.
47
4.2 Methodology
Fieldwork on Erebus volcano is limited to the austral summer, and typically takes
place from late-November to early January. Where we refer to a field season by year,
we are referring to the year in which it began. The logistics involved in reaching
the crater rim, combined with frequent bad weather conspire to limit the interval
of IR image data acquisition to a few weeks each year. The intervals of useful data
are further reduced due to fluctuations in the IR transmission between camera and
lava lake. When the gas/aerosol plume is highly condensed (high relative humidity)
the IR transmission in the camera waveband is poor and the images of the lake are
of unusable quality. The latest IR camera system, which was deployed in December
2012, is capable of year-round operation (dependent on power) [Chapter 2, Peters
et al., 2014c]. The data from this fully automated system was not available at the
time of this investigation, but is analysed in Chapter 5.
4.2.1 Camera Hardware
All IR images of the Erebus lava lake used in this chapter were recorded by tripod-
mounted camera systems installed at the Shackleton’s cairn site on the northern
side of the Main Crater (Fig. 4.1). Three different IR camera systems have been
used on Erebus since 2004 (summarised in Table 4.1). The first of these was an
Agema Thermovision 550 mid-infrared (3–5µm) camera, as described by Oppen-
heimer et al. [2004], which acquired images with a 4 s time-step. The second was
a FLIR P25 camera equipped with a 72 mm IR lens. This camera has an uncooled
320×240 element detector with a waveband of 7.5–13µm and a spatial resolution
(instantaneous field of view) of 0.31 mrad. The low temperatures at the crater rim
of Erebus severely impacted the P25’s ability to write images to its compact flash
card, and as a result the interval between successive images varies between 8 and
20 s. The most recent IR camera model we have employed is a FLIR SC645, with
an uncooled 640×480 element detector, waveband of 7.5–13µm, and a spatial reso-
lution of 0.69 mrad. Although the SC645 is capable of frame rates of up to 25 Hz,
we have typically acquired images at 2 s intervals. It was, however, operated with a
0.5 s time-step in acquisitions for a part of the 2010 field season and we have used
these higher time resolution data for bubble event analysis (Section 4.2.5).
48
43°
240 m
224 m
Figure 4.1: Photograph showing the IR camera setup on the crater rim of Erebus
and cartoon showing the viewing geometry of the lake (calculated from TLS data).
Year IR Camera Waveband (µm) Capture Rate (Hz)
2004 AGEMA 3 - 5 0.25
2005 No data collected - -
2006 - 2009 FLIR P25 7.5 - 13 0.125 - 0.05
2010 - 2011 FLIR P25 and FLIR SC645 7.5 - 13 0.125 - 0.05 and 0.5
2011 - 2014 FLIR SC645 7.5 - 13 0.5
Table 4.1: IR camera systems used at Erebus volcano by year.
4.2.2 Data Selection
Interruptions to the recording of IR images on Erebus are common. The Agema
and P25 cameras both required their memory cards to be changed regularly and
equipment failure was frequent due to the harsh operating conditions. These fac-
tors have resulted in a segmented data set with many gaps. The first step in data
selection was to split the data into groups of continuous acquisition that contained
no two images more than 40 s apart. Groups spanning less than one hour of acquisi-
tion were discarded. Subsequent data processing was performed on a per group basis.
High winds at the summit of Erebus cause camera shake, potentially introducing
large errors into the velocity estimates calculated by the motion tracking algorithm.
This problem is particularly acute in data from the Agema and P25 cameras, which
did not have such stable tripod mounts as does the new SC645 system. Attempted
stabilisation of the images in post-processing failed due to the lack of distinctive sta-
tionary features in the images. Instead, a simpler approach was followed, in which
only periods of data with little or no camera shake were analysed. Due to the large
volume of available data, an automated routine for identifying such periods was
49
110 120 130 140 150 160 170
x coordinate
80
90
100
110
120
130
140
150
160
y
 c
o
o
rd
in
a
te
Figure 4.2: Positions (in pixels) of the lake bounding box within images from 2010
showing their assignment into clusters. The scatter (standard deviation) of points
within one cluster is used as an indicator of camera shake. Only images assigned
to clusters with a standard deviation of less than 1.0 were used for surface motion
estimation.
developed. This involved first defining the bounding box of the lake in each image
by thresholding the image at a predetermined level, and identifying the non-zero re-
gion. Images in which the bounding box could not be found, or was unusually small
were rejected, as these characteristics point to poor visibility of the lake (typically
caused by high relative humidity, blowing snow, or rime-ice accumulation on the
lens). The centre coordinates of the bounding boxes were then assigned to clusters
using SciPy’s fclusterdata function [Jones et al., 2001]. To reduce the run time of
the clustering algorithm, duplicate bounding box positions were discarded before
clusters were computed. Using the standard deviation of the bounding box coordi-
nates in each cluster as an indicator of camera shake, the best clusters for each year
(typically with a standard deviation of < 1.0 pixels) were selected (Fig. 4.2). As a
final check of data quality, the images in each cluster were compiled into a video,
which was then viewed to ensure good visibility of the lake and minimal camera
shake throughout.
50
4.2.3 Motion Tracking
The steps involved in obtaining a time-series of mean lake surface speed from a set
of IR images are shown in Fig. 4.3 and are described in detail below.
Since the focal plane of the thermal camera is not parallel to the surface of the lava
lake, perspective effects mean that each pixel in the image represents a different
distance in the plane of the lake. To correct for this distortion, each image was
rectified before the motion tracking was carried out. The required transformation
was calculated by matching points in the image to points in a terrestrial laser scan
of the lake. OpenCV’s [Bradski, 2000] cvFindHomography function was then used to
calculate the required transformation matrix, and the cvWarpPerspective function
used to apply it (see Bradski and Kaehler [2008]). Correcting the images in this
way also accounts for any lens distortion. Terrestrial laser scan data of the lava
lake were only available for 2008 onwards. For thermal images from earlier years,
the homography matrix was calculated from the viewing angle of the camera and
the size of the lake (which had been estimated with a handheld laser range finder).
Although this method neglects lens distortion, we expect the effects to have little
impact on the results obtained.
The significant temperature contrast between the lake and the surrounding crater
causes problems for the feature tracking algorithm. As the strongest feature in the
image, the lake boundary tends to dominate over the structure within the lake that
we are actually interested in. This issue can be overcome by masking the regions
outside of the lake with Gaussian-distributed white noise with a mean and variance
similar to that of the pixels within the lake. Random noise is used rather than a
fixed value to prevent the output of the bandpass filters used in the wavelet decom-
position from being exactly zero, as this causes the algorithm to fail.
The feature tracking algorithm itself is based on the Dual-Tree Complex Wavelet
Transform (DT-CWT) [Kingsbury, 2001]. Unlike the widely used Discrete Wavelet
Transform (DWT) (see for example Polikar [2010]), DT-CWT is approximately shift
invariant, meaning that shifts in input signal do not cause changes in energy dis-
tribution between wavelet coefficients at different scales. Instead, it is found that
shifts in the input signal manifest themselves as a phase shift between the wavelet
coefficients. By decomposing each frame in a series of images using the DT-CWT
and comparing the phase shifts between the subimages, it is possible to estimate
the displacement field that maps features in one frame to the next. The estimate
51
1700 1750 1800 1850 1900 1950
Image  Number
0.02
0.03
0.04
0.05
0.06
0.07
M
e
a
n
 s
u
rf
a
c
e
 v
e
lo
c
it
y
 (
m
s
−
1
)
A
B
C
D
Figure 4.3: Flow diagram showing the steps in extracting mean suface velocity
from IR images. A) Perspective correction to account for viewing angle of camera.
B) Random noise masking to reduce contrast between lake and surroundings. C)
DT-CWT feature tracking estimates motion vectors between consecutive frames by
calculating phase shift between complex wavelet coefficients. D) Mean velocity found
by averaging motion vectors in the central region of the lake.
52
is first made at a coarse level of decomposition, with subsequent estimates made at
finer levels so as to refine the result [Magarey and Kingsbury, 1998]. Oppenheimer
et al. [2009] tuned the method specifically for working with IR images from Erebus
volcano, and we adopted the same parameters for the motion estimation reported
here. The algorithm was verified by passing in sets of two identical images, one of
which was shifted by a known amount and checking that the shift was correctly de-
tected. Further verification was achieved by comparing velocity estimates obtained
from the algorithm with visual estimates found by counting pixels.
Finally, the mean surface speed of the lake was found by averaging the magnitudes
of the computed velocity vectors. To avoid possible edge effects, only velocity vec-
tors from the central region (at least 3 pixels inside the lake boundary) of the lake
were included in the averaging.
An animation showing results of the motion tracking methodology described above
(from rectified images to mean lake speed) is provided in the supplementary material
for the electronic version of Peters et al. [2014d].
4.2.4 Time-Series Analysis
The Fourier Transform is a popular tool for evaluating frequency components in time
series data. Although it has perfect frequency resolution, it does not provide any
time resolution and is therefore not suitable for evaluating non-stationary signals. It
is possible to achieve some time resolution by splitting the signal under analysis into
windows of fixed length and applying a Fourier Transform to each window. However,
the choice of window length is critical since although short windows provide better
time resolution, they also inhibit the frequency resolution of the Fourier Transform.
Wavelet transforms partially solve this resolution problem by using a variable length
window to achieve higher frequency resolution at low frequencies, and higher time
resolution at high frequencies. This is desirable, since high frequency components
tend to be short-lived in “real world” signals. Referred to as “Multi-Resolution
Analysis”, this approach makes wavelet transforms particularly suitable for assess-
ing the time-frequency structure of non-stationary signals. An excellent summary
of wavelet transforms is given by Polikar [2010].
As can be seen in Fig. 4.4e, the mean surface speed time-series obtained are highly
non-stationary. To evaluate the periodic components of the series with time, we
53
therefore use a Morlet wavelet transform to produce spectrograms of the data. Our
implementation of the Morlet transform is the same as that of Boichu et al. [2010].
The mean speed data were interpolated to a uniform 1 s time step prior to the Mor-
let transform using simple linear interpolation.
As illustrated by Figs. 4.4a and 4.4b, some of the ∼5–18 min cycles are of much
greater amplitude than others, and will result in a very high modulus in the Morlet
transform. Longer time-series tend to exacerbate this problem, since they often con-
tain at least a few very high amplitude oscillations, which then saturate the colour
scale and mask much of the other detail. In this way, the cyclicity of the lake may
not be apparent even if it exists. However, creating a spectrogram of just the data
from the “non-cyclic” time period, reveals that there are indeed still ∼5–18 min pe-
riod components present – they are simply of lower amplitude.
4.2.5 Bubbles
Bubbles breaking the surface of the lake manifest themselves as sharp peaks in the
mean surface speed time-series. The poor time resolution of the Agema and P25
datasets mean that most bubbles are not recorded. However, much of the SC645
data from 2010 was recorded at 2 Hz, which is more than sufficient to capture the
arrival of bubbles at the surface.
Bubble events were located in time by comparing the mean speed time-series to a
low-pass filtered copy of itself. Bubbles were classified as events where the speed
was greater than 1.2 standard deviations above the low-pass filtered value. The
value of 1.2 was chosen by comparing bubble events detected by the algorithm to
those located manually in a test set of data spanning three hours. The analysis was
conducted on a continuous time-series of good quality data from 24 December 2010,
spanning approximately 13 h. By visually inspecting the IR images corresponding
to each of the bubble events, we determined that all events were small (metre-scale,
with no ejection of material from the lake).
The bubble events detected are uniformly distributed in time. However, this tells
us nothing of how they are related to the pulsatory behaviour of the lake. What
is of real interest is how bubble events relate to the phase of the speed cycles; for
example, do more bubbles surface during periods of fast surface movement? In order
to evaluate a possible relationship between the cyclicity and bubble events we use
54
00:02 05:35 11:08 16:42 22:15 03:48
Time (UT)
100
600
1100
1600
2100
P
e
ri
o
d
 (
se
co
n
d
s)
0.00
0.15
0.30
0.45
0.60
0.75
0.90
1.05
1.20
×10 5
100
600
1100
1600
2100
P
e
ri
o
d
 (
se
co
n
d
s)
0.0
1.5
3.0
4.5
6.0
7.5
9.0
×10 5
100
600
1100
1600
2100
0.0
0.3
0.6
0.9
1.2
1.5
1.8
2.1
2.4
×10 5
0.02
0.04
0.06
0.08
0.02
0.04
0.06
0.08
0.02
0.04
0.06
0.08
0.10
M
e
a
n
 v
e
lo
ci
ty
 (
m
s−
1
)
(a) (b)
(c) (d)
(e)
(f)
Figure 4.4: (e) Selection of mean surface velocity data from December 2010 and
(f) the corresponding Morlet wavelet transform showing the periodicities present.
(c, d) Expanded sections of the time series showing how an increase in period is
accompanied by an increase in amplitude. (a, b) Morlet wavelet transforms of the
expanded sections showing evidence of periodicity even when it is not clearly present
in Morlet transform of the full time series (f). The colour scales represent the Morlet
transform modulus.
55
the method of delays (e.g. Kantz and Schreiber [2003]) to reconstruct our time-series
data into a phase space representation. If the bubble events are somehow correlated
to the phase of the speed cycles then we argue that their distribution in phase space
will differ from that of a random sample taken from the time-series. We can imagine
this as being due to a clustering of bubble events at certain positions in phase space.
A clear description of the method of delays for phase space reconstruction is given
by Kantz and Schreiber [2003] and also by Richter and Schreiber [1998]. We will
therefore not attempt to describe the technique in any detail here, instead we give
the specifics of how the parameters required for delay reconstruction were calculated
for our time-series. Essentially, phase space reconstruction of a time-series involves
mapping each sample in the series to a vector in phase space. For a scalar time-
series x1, x2, x3, . . . the method of delays can be used to calculate the corresponding
phase space vectors xn = (xn, xn−l, xn−2l, . . . , xn−(d−1)l) where l is known as the lag,
and d is the embedding dimension. Figure 4.5 shows an example of phase space
reconstruction using the method of delays with an embedding dimension of 2, and
demonstrates how points in the time-series are mapped into phase space.
Careful selection of both the lag and the embedding dimension is of paramount
importance to obtain an effective phase space reconstruction of the original time-
series. Estimation of a suitable lag for our data was performed using two different
techniques. The first involved calculation of the autocorrelation function of the
time-series. The time required for the autocorrelation function to decay by a factor
of e is stated as being a reasonable estimate for the lag [Kantz and Schreiber, 2003].
For our data, this gave a lag of 142 s. A second technique for lag estimation, also
described by Kantz and Schreiber [2003], is to find the minimum of the mutual infor-
mation of the time-series. We calculated the mutual information of our time-series
using the TISEAN software package [Hegger et al., 1998], and found the minimum
to be at 200 s. We repeated our bubble event analysis for several different values of
lag between these two estimates and found no appreciable difference in the results.
The results presented here are from the analysis using a lag of 150 s.
To determine a good embedding dimension, we followed the approach proposed by
Hegger and Kantz [1999], in which the false nearest neighbour method [Kennel et al.,
1992] is combined with surrogate data tests. The false nearest neighbour method
compares the ratio of distances between points and their nearest neighbours between
an embedding dimension of n and n+1. If the ratio is greater than a threshold (s),
56
t1−τ t1 t2−τ t2
t
y(t1−τ)
y(t1 )
y(t2−τ)
y(t2 )y
(t
)
Time series
y(t1 ) y(t2 )
y(t)
y(t1−τ)
y(t2−τ)
y
(t
-τ
)
Phase space reconstruction
Figure 4.5: Phase space reconstruction of a simple time-series using the method of
delays. The original time-series is shown in the left hand plot and the phase space
reconstruction using a lag of τ and an embedding dimension of 2 is shown in the
right hand plot. Samples from the time-series (e.g. t1 and t2) are mapped to phase
space by plotting them against the value of the time-series one lag prior to their
sample time.
the points are said to be “false neighbours”. A high percentage of false nearest neigh-
bours is indicative of too low a choice for the embedding dimension. The TISEAN
software package was used both for the calculation of false nearest neighbours and
for the creation of surrogate data [Schreiber and Schmitz, 1999]. Figure 4.6 shows
the calculated percentage of false nearest neighbours for different embedding di-
mensions and thresholds. There is a clear difference in behaviour between the real
data and the surrogate data. This indicates that the loss of false neighbours as we
move to higher embedding dimensions is not simply due to linear correlations in
the data. Furthermore, we can see that, for an embedding dimension of 4, the per-
centage of false nearest neighbours falls away very rapidly for low threshold values.
We therefore used 4 as our embedding dimension when performing the phase space
reconstruction. The data were low-pass filtered prior to phase space reconstruction
to remove noise and the spikes due to bubbles.
The time-series analysed contains 141 bubble events (Fig. 4.7). We compared the
cumulative distribution function (CDF) of the bubble events to a reference CDF in
each of the phase space dimensions. The reference CDF is the CDF of the time-series
itself. As an indicator of the expected variation in CDFs, the standard deviation of
10,000 CDFs, each constructed from 141 points randomly sampled from the time-
series, was computed. A significant variation of the bubble event CDF from that
of the reference in any of the dimensions, would indicate some correlation to the
57
0 2 4 6 8 10 12 14
s
0.0
0.2
0.4
0.6
0.8
1.0
FN
N
 f
ra
ct
io
n
Real data
0 2 4 6 8 10 12 14
s
0.0
0.2
0.4
0.6
0.8
1.0
FN
N
 f
ra
ct
io
n
Surrogate data
Figure 4.6: The fraction of false nearest neighbours (FNN) against the threshold
value for embedding dimensions from 1 to 6 (top to bottom) and a lag of 150 s, for
both the real data (left hand plot) and a surrogate data series (right hand plot). The
difference between the real and surrogate data shows that the reduction in FNN as
we increase the embedding dimension is not simply due to linear correlations in the
data. For an embedding dimension of 4, the FNN fraction reduces very quickly for
small threshold values and there is little difference in behaviour compared to larger
embedding dimensions. We therefore chose 4 as a suitable embedding dimension for
our analysis.
phase of the cycle. Differences between CDFs were quantified using the two-sample
Kolmogorov-Smirnov test (K-S test). The computed critical value for the K-S test
at 90% confidence (based on a reference sample size of 90901, and a bubble event
sample size of 141) is 0.102.
To verify the technique, we created a set of 95 fake bubble events located at the
peaks of the mean speed cycles (Fig. 4.7). These events were then subjected to the
same analysis as the real bubble events. The critical value for the K-S test at 90%
confidence is 0.125 for the fake bubble sample size of 95. As shown in Fig. 4.8, the
CDFs for the fake bubble events show a strong deviation from that of the random
samples in each of the phase space dimensions, with K-S test results of 0.50, 0.15,
0.16 and 0.13 respectively (i.e. all above the critical K-S value, suggesting the two
samples came from different distributions). Hence, the technique correctly identi-
fied the correlation between the fake bubble events and the phase of the speed cycles.
58
0.00
0.02
0.04
0.06
0.08
0.10
0.12
0.14
Raw time series
0.01
0.02
0.03
0.04
0.05
0.06
0.07
Low-pass filtered, showing bubble events
04:00 05:00 06:00 07:00 08:00 09:00 10:00 11:00 12:00 13:00 14:00 15:00 16:00
Time (UT)
0.01
0.02
0.03
0.04
0.05
0.06
0.07
Low-pass filtered, showing fake bubble events
M
e
a
n
 s
u
rf
a
ce
 s
p
e
e
d
 (
m
s−
1
)
Figure 4.7: Top panel: time-series of mean lake surface speeds from 24 December
2010. Middle panel: the same data, low-pass filtered to remove noise and spikes due
to bubbles and with bubble events marked. Bottom panel: low-pass filtered data
with the fake bubble events used for testing marked.
59
0.02 0.03 0.04 0.05 0.06
V(t )
0.0
0.2
0.4
0.6
0.8
1.0
C
u
m
u
la
ti
v
e
 p
ro
b
a
b
il
it
y
0.02 0.03 0.04 0.05 0.06
V(t  - τ)
0.0
0.2
0.4
0.6
0.8
1.0
C
u
m
u
la
ti
v
e
 p
ro
b
a
b
il
it
y
0.02 0.03 0.04 0.05 0.06
V(t  - 2τ)
0.0
0.2
0.4
0.6
0.8
1.0
C
u
m
u
la
ti
v
e
 p
ro
b
a
b
il
it
y
0.02 0.03 0.04 0.05 0.06
V(t  - 3τ)
0.0
0.2
0.4
0.6
0.8
1.0
C
u
m
u
la
ti
v
e
 p
ro
b
a
b
il
it
y
(a) (b)
(c) (d)
Figure 4.8: The cumulative distribution functions (CDFs) in each of the four phase
space dimensions for the bubble events (solid line) and fake bubble events (dashed
line). The x-axes represent the coordinates of the events in the corresponding phase
space dimension, from the zero-lags dimension, V(t), up to the three-lags dimension,
V(t - 3τ), where τ=150 s). The shaded region represents one standard deviation on
either side of the reference CDF (i.e. the CDF of the mean speed data). Large
deviations from the reference CDF are indicative of a correlation with the phase of
the speed cycles, as can be seen in the fake bubble data. The deviation from the
reference in the first dimension (a) of the bubbles CDF is attributed to imperfect
filtering of the signal rather than a phase dependence.
60
4.3 Results
The 2010 field season was characterised by exceptional visibility of the lava lake. In
addition to the IR images captured, several short time-series of visible images were
captured using a digital SLR camera equipped with a telephoto lens. Figure 4.9
shows a short time-series of mean surface speed and mean surface temperature data
calculated from IR images, with visible images corresponding to peaks and troughs in
the speed also shown. There are no consistent differences observed between the ap-
pearance of the lake surface during periods of high speeds and periods of low speeds.
Oppenheimer et al. [2009] found a strong correlation between the phase of cycles
in mean surface speed and radiative heat loss in their data set. This correlation
is further demonstrated by the time-series shown in Fig. 4.9 and Fig. 4.10. Note
however, that since we have not attempted an accurate temperature calibration of
the IR images, we present mean surface temperatures normalised to their maximum
value. What is not clear from these data alone, is whether the temperature varia-
tions observed are due to an increase in the surface temperature of the lava in the
lake, or an increase in the area of cracks in the surface crust of the lake caused
by the increased motion [Oppenheimer et al., 2009]. Additional cracks will expose
more of the underlying lava to the surface and will therefore cause an increase in the
mean temperature recorded by the IR camera. Increased cracking during periods
of higher surface speed is not obvious in the images shown in Fig. 4.9, suggesting
that the changes in recorded temperature are indeed due to variations in the mean
temperature of the lake surface. However, we feel that a qualitative argument such
as this is insufficient to rule out increased cracking as a cause.
In an attempt to identify more rigorously the reason for the temperature cycles,
we compared the histograms of the thermal images at the minima and maxima of
the cycles. If the cycles are caused by an increase in surface temperature, then we
would expect the histograms at the cycle maxima to be shifted relative to those at
the minima. If increased cracking is the cause, we would expect more high tem-
perature pixels, resulting in a skewing of the histograms at the maxima compared
to those at the minima. Unfortunately, the results obtained were ambiguous, with
greater differences between histograms from the same point in the cycles than found
by comparing those at maxima to those at minima. The cause of the measured tem-
perature fluctuations remains elusive, however, it seems likely that they are caused
by a combination of both increased surface cracking and increased surface temper-
ature.
61
09:30:00 10:00:00 10:30:00 11:00:00 11:30:00
Time (UTC)
0.02
0.03
0.04
0.05
0.06
0.07
0.08
S
p
e
e
d
 (
m
s
−
1
)
0.92
0.94
0.96
0.98
1.00
1.02
N
o
rm
a
li
s
e
d
 s
u
rf
a
c
e
 t
e
m
p
e
ra
tu
re
Temperature
Speed
(a) (b)
(a)
(c) (d)
Figure 4.9: Time series of mean lake surface speed and mean surface temperature
(arbitrary units) from 17 December 2010 calculated from images captured with the
P25 IR camera, with photographs showing the appearance of the lake surface at
periods of high surface speed (a,b) and low surface speed (c,d). There is no distinct
difference between the appearance of the lake surface during periods of high or low
surface motion. Note that in (d), a small gas bubble can be seen just reaching the
surface. Note that the bubble does not appear as “shot noise” in the speed data
because of the poor time resolution of the P25 camera. The lake is approximately
40 m across its long axis.
62
23:57 00:02 00:07 00:12 00:17 00:22 00:27 00:32 00:37 00:42
Time (UT)
0.01
0.02
0.03
0.04
0.05
0.06
0.07
0.08
0.09
0.10
M
e
a
n
 s
u
rf
a
ce
 s
p
e
e
d
 (
m
s−
1
)
0.88
0.90
0.92
0.94
0.96
0.98
1.00
N
o
rm
a
lis
e
d
 s
u
rf
a
ce
 t
e
m
p
e
ra
tu
re
Temperature
Speed
Figure 4.10: Short time-series of mean surface velocity and mean surface temper-
ature data from 21–22 December 2010 calculated from images acquired with the
SC645 camera. We have not attempted to retrieve accurate temperatures from the
images, and instead report unitless temperature values normalised by the maximum
value in the time-series. The pulsatory behaviour is particularly clear during this
period, and the symmetry of the peaks in time (i.e. about their horizontal centre)
is evident.
It is important to realise that even if increased surface cracking is insignificant, cy-
cles in the surface temperature do not necessarily reflect periodic variations in the
internal temperature of the lake itself. It is possible for example, that during periods
of increased surface motion the mean surface age of the lake is lower. Younger crust
is likely to be thinner, and hence the conductive heat flow through it will be larger,
resulting in higher measured surface temperatures despite the bulk temperature of
the lake remaining static.
Figure 4.10 shows a short time-series of mean surface speed and mean surface tem-
perature calculated from IR images captured in 2010. The pulsatory behaviour was
particularly pronounced during the period shown, and the waveform of the cycles is
clear. The peaks in speed and temperature are approximately Gaussian in shape,
with rising and falling edges that are symmetric about the centre of the peak. The
peaks tend to be shorter lived than the troughs, suggesting a system with a sta-
ble baseline state that is being perturbed, rather than a system that is oscillating
about a mid-point. There also appears to be a correlation between the magnitude
of the cycles and their period, with longer period cycles having a greater amplitude.
However, such a relationship is not always observed, and there are many instances
where the opposite is true.
Morlet spectrograms of the mean speed data from the 2007–2011 field seasons are
63
provided in Section 4.3.1. What is clear from the data is that the cycles in speed are
not strictly periodic. Instead, there tends to be a broad range of periodic components
present, centred at around 900 s. However, these components appear to be fairly
consistent across the dataset and have not changed appreciably during the period
of study. Figure 4.11 further illustrates this point, showing the time average (nor-
malised to have a mean of zero and standard deviation of one) of the modulus of all
the Morlet spectrograms from each field season. The general trend towards higher
modulus at longer periods is due to the fact that long period variations in mean
speed tend to be of greater amplitude than short period variations (as is typical for
most time-series data from natural systems). Despite this, the broad peak around
900 s is evident in the data from the 2007–2011 field seasons. The time-series from
the 2004 and 2006 field seasons were of insufficient duration to allow analysis for long
period behaviour, and as a result do not show the same behaviour as the other years.
It is unfortunate that during the 2005 and 2006 field seasons, which covered the
period of increased explosive activity, IR data were either not collected at all, or are
of insufficient length to compare to other years. However, as shown in Fig. 4.12, the
pulsatory behaviour of the lake appears to be unperturbed by large bubble bursts.
The figure shows a short time-series of mean surface velocity data from 29–30 De-
cember 2010, during which a large (∼30 m) bubble arrived at the surface of the lake.
Despite a significant ejection of material from the lake, the mean speed data show
that the pulsatory behaviour appears to be uninterrupted on timescales greater than
that of lake refill (several minutes). It is interesting to note that at the time of the
explosion, the Morlet spectrogram shows a particularly strong periodic component
at ∼1000 s. We believe that this may be caused by increased surface speeds in the
build-up to the explosion and also during the recovery phase as the lake refills. The
IR images show that the lake level rises rapidly immediately prior to a large bubble
reaching the surface, likely causing an increase in the recorded surface speed. Rapid
flow of lava into the lake during the refill phase of an explosive event is also likely
to cause elevated surface speeds.
In addition to the apparent stability of cycles in surface speed, the magnitude of
the surface speed has also remained approximately unchanged since 2004. Although
the mean surface speed can exhibit considerable variability (∼3–20 cm s−1) on a
timescale of days, no systematic change was observed over the period of study.
Whilst the behaviour of the mean surface speed has remained remarkably stable,
64
0 500 1000 1500 2000
Period (seconds)
2011
2010
2009
2008
2007
2006
2004
M
e
a
n
 M
o
rl
e
t 
tr
a
n
s
fo
rm
 m
o
d
u
lu
s
 (
n
o
rm
a
lis
e
d
)
Figure 4.11: Normalised time averages of the Morlet transform modulus of all avail-
able data from each field season. The data from each year have been vertically offset
from each other for clarity. A broad peak corresponding to the period of the lake
cycles (∼900 s) is evident in the 2007–2011 data confirming that the cyclic behaviour
is a persistent and largely unchanging feature of the lake. The time-series from the
2004 and 2006 field seasons were of insufficient length to be able to resolve long
period fluctuations.
65
Figure 4.12: A 5 h time-series of mean surface speed data from the 29–30 December
2010 calculated from images acquired with the SC645 camera and the corresponding
Morlet spectrogram. The IR images above show the state of the lake before, during
and after a large (∼30 m) bubble burst. Cycles of between ∼600–1100 s are visible
across the time-series both in the Morlet spectrogram and in the mean speed data.
The strong ∼1000 s signal in the spectrogram corresponding to the explosion itself
may be caused by elevated speeds due to the rising of the lake level prior to the
explosion and the subsequent refilling of the lake afterwards.
66
2004 2005 2006 2007 2008 2009 2010 2011
Year
500
1000
1500
2000
2500
La
ke
 s
u
rf
a
ce
 a
re
a
 (
m
2
)
IR image estimate
LiDAR estimate
Figure 4.13: Surface area of the Erebus lava lake by field season (i.e. December of the
year shown). The areas have been estimated from a combination of terrestrial laser
scan data (Jones [2013], Frechette and Killingsworth pers. comm.) and rectified IR
images.
the visual appearance of the lava lake has changed significantly. Figure 4.13 shows
how the surface area of the lake (calculated from IR images and terrestrial laser scan
data) has decreased since the first measurements in 2004. Overall the surface area
has reduced by a factor of approximately four. The terrestrial laser scan data also
show that since at least 2008 (when the first such data were recorded), the decrease
in area has been accompanied by a 3–4 m per year drop in the lake’s mean surface
elevation (Jones [2013], Frechette and Killingsworth pers. comm.). The dramatic
reduction in surface area cannot be accounted for by the drop in surface elevation
(i.e. due to the lake receding into a conical basin) since the lake walls are observed
(from terrestrial laser scan data and visual observations) to have a near-vertical
profile. It seems likely therefore, that the observed decrease in surface area is due
to a change in the geometry of the lake basin, either due to accretion of lava onto
its sides, or deformation of the surrounding rock. The apparent lack of influence
of lake geometry on the cyclic behaviour would tend to suggest that the cycles are
driven by processes occurring in the magma conduit or a connected reservoir rather
than in the lake itself.
Figure 4.8 shows the cumulative distribution functions (CDFs) of the bubble events,
and fake bubble events, in each of the four phase space dimensions. The shaded areas
67
delimit one standard deviation on either side of the reference CDFs. As discussed
in Section 4.2.5, the CDFs for the fake bubbles show a strong deviation from the
reference, correctly identifying the correlation between the phase of the speed data
and the fake bubble events. In contrast, the CDFs for the real bubble events are
very similar to the reference in all but the first dimension. The K-S test gives values
of 0.15, 0.05, 0.07 and 0.06 for the four dimensions, respectively. Apart from the
first dimension, these are all below the critical K-S value at 90% confidence (0.102),
indicating that the bubble events are from the same distribution as the speed data
itself and that there is, therefore, no correlation between the phase of the speed
cycles and the bubbles.
In the first dimension, the CDF of the bubble events appears to be the same shape
as that of the mean speed data, but shifted slightly to the right (higher mean speed).
We believe that this is caused by the failure of our low-pass filtering to remove fully
the spikes caused by bubble events in the mean speed data, rather than any cor-
relation with the phase of the cycles. As a result, bubble events appear to occur
at slightly higher speeds than they actually do, shifting the CDF to the right. We
tested this hypothesis by plotting the CDF of 141 randomly selected points from
the speed data with a linear offset of 0.002 ms−1 added. The results showed a CDF
that matched that of the mean speed data in all dimensions except the first, which
showed a linear offset to the right as expected. We therefore conclude that the
bubble events are not correlated to the phase of the velocity cycles and that the
deviation we observe in the first dimension is due to the low-pass filter’s inability
to remove completely the effects of bubble events from the underlying mean speed
signal.
4.3.1 Spectrograms
Shown below are representative Morlet spectrograms of the mean surface speed of
the lava lake from 2007-2011. The spectrograms shown here are from the longest
continuous time series of IR images from the P25 camera from each year. Data from
2004 and 2006 have not been included, since the time series are too short to evaluate
periods of above ∼600 s.
The spectrograms demonstrate that the cycles are not perfectly periodic, nor is there
any discernible pattern to the variation in their period. Instead, strong wavelet re-
sponses are apparent, ranging in period between ∼500–1100 s, and often with sharp
68
transitions between them. The complexity of the spectrograms (in addition to the
issues discussed in Fig. 4.4) demonstrates the difficulty of interpreting extended time
series of speed data due to the high degree of non-stationarity.
03:01 04:24 05:47 07:11 08:34 09:57 11:21 12:44 14:07
Time (UT)
100.0
300.0
500.0
700.0
900.0
1100.0
1300.0
1500.0
P
e
ri
o
d
 (
se
co
n
d
s)
2007
0.0
0.4
0.8
1.2
1.6
2.0
2.4
2.8
3.2
×10 5
11:43 14:30 17:16 20:03 22:50
Time (UT)
100.0
300.0
500.0
700.0
900.0
1100.0
1300.0
1500.0
P
e
ri
o
d
 (
se
co
n
d
s)
2008
0.0
0.3
0.6
0.9
1.2
1.5
1.8
2.1
2.4
2.7
×10 5
10:02 11:25 12:49 14:12 15:35 16:59
Time (UT)
100.0
300.0
500.0
700.0
900.0
1100.0
1300.0
1500.0
P
e
ri
o
d
 (
se
co
n
d
s)
2009
0.0
0.8
1.6
2.4
3.2
4.0
4.8
5.6
6.4
×10 5
69
00:02 05:35 11:08 16:42 22:15 03:48
Time (UT)
100.0
300.0
500.0
700.0
900.0
1100.0
1300.0
1500.0
P
e
ri
o
d
 (
se
co
n
d
s)
2010
0.00
0.15
0.30
0.45
0.60
0.75
0.90
1.05
1.20
×10 5
06:40 09:27 12:13 15:00 17:47 20:33
Time (UT)
100.0
300.0
500.0
700.0
900.0
1100.0
1300.0
1500.0
P
e
ri
o
d
 (
se
co
n
d
s)
2011
0.0
0.2
0.4
0.6
0.8
1.0
1.2
1.4
1.6
×10 5
4.4 Discussion
A common conclusion of multi-year studies conducted at Erebus volcano is that
its behaviour is remarkably stable. Observations of radiant heat output [Wright
and Pilger, 2008], SO2 flux [Sweeney et al., 2008] and seismicity [Aster et al., 2003]
have all found very little variation during the past decade. Our findings that the
pulsatory behaviour of the lava lake has been a persistent and unchanging feature
(both on a daily and a yearly time-scale) since at least 2004 fit well with these pre-
vious findings and further emphasise the remarkable stability of Erebus’s magmatic
system. The preservation of cycles in surface speed despite large perturbations to
the system (decametric bubble bursts) is indicative not only of the stability of the
responsible mechanism, but also that it is likely sourced at a deeper level than the
lake itself. This argument is further supported by the consistency of the motion
cycles despite a reduction in the area of the lava lake, exposed at the surface, over
70
the period of observation.
The broad width of the peak in the spectrograms that we observe is consistent with
the findings of Oppenheimer et al. [2009] who found the period of the fluctuations in
mean lake speed to vary between ∼5–18 min for the 2004 dataset. Although short
sections of our data appear to contain several discrete frequency bands within this
range, such fine scale structure is never observed consistently over time periods of
more than a few hours. No clear pattern to the variation of the period of fluctuations
measured is evident from the spectrograms. However, it is important to consider
how well the mean speed of the lake surface represents the underlying process re-
sponsible for the pulsatory behaviour. Even if this process contains well defined,
discrete periods, complex flow dynamics and other forcing mechanisms (such as
thermal convection within the lake) may result in a highly non-linear response in
the surface speed. It is possible that the broad distribution in period of the cycles in
speed observed is due to the complex coupling between a periodic driving mechanism
and the dynamics of the bubbly flow within the lake. Given the correlation between
surface motion and gas composition ratios (which must have different couplings to
the driving mechanism) reported by Oppenheimer et al. [2009], we believe that the
variability in period stems primarily from the variability in the underlying driving
mechanism.
Current theories on driving mechanisms for lava lake fluctuations can be grouped
into three main categories; (i) instability in density driven, bi-directional flow of
magma in the conduit feeding the lake [Oppenheimer et al., 2009], (ii) “gas piston-
ing” caused by gas accumulation either beneath a solidified crust on the surface of
the lake [Patrick et al., 2011] or as a foam layer at the top of the lava column [Orr
and Rea, 2012], and (iii) gas bubble-driven pressurisation changes [Witham et al.,
2006]. In the latter mechanism, the (uni-directional) upflow of bubbly magma in
the conduit is interrupted by excess hydrostatic pressure in the lake. Stagnation of
the flow allows bubbles in the conduit to coalesce into large gas slugs that rise to
the surface independently of the melt. The loss of large gas slugs at the surface of
the lake causes an increase in pressure at the base of the conduit. If this exceeds
the pressure in the magma chamber then downflow occurs, suppressing the ascent
of bubbles in the conduit. As the lake drains, the downflow reduces until it can no
longer suppress the ascent of bubbles, and the cycle repeats. Witham et al. [2006]
were able to demonstrate this mechanism by bubbling air through a water column
with a basin attached to the top to represent the lake. They observed cyclic varia-
71
tions in the depth of water in the basin, consisting of a logarithmic increase in depth
followed by a rapid, linear decrease.
As shown by Orr and Rea [2012], gas pistoning, as observed at Kı¯lauea, is also an
asymmetric process, consisting of a relatively slow, cumulative deviation from the
baseline state of the system as bubbles are trapped in the foam layer or beneath
the solidified crust, followed by a sudden release of the accumulated gas and rapid
return to the baseline state. The temporal symmetry of the perturbations in the
Erebus lava lake is not consistent with either of these mechanism. It may be argued
that the complex geometry of the upper magmatic system of Erebus could lead to
a more symmetric variation than observed by Witham et al. [2006] and Orr and
Rea [2012]. However, our finding that the arrival of small (metre-scale) bubbles at
the surface of the lake is uncorrelated with the phase of the speed cycles is only
consistent with the bi-directional flow mechanism. Both bubble-driven mechanisms
require a periodic release of bubbles prior to lake draining and, in the case of the
Witham et al. [2006] mechanism, a significant decrease in the number of bubbles
during lake draining. Large (decametric) bubbles, typically occur at Erebus only a
few times per week, and cannot therefore be responsible for the ∼5–18 min cycles.
Since no periodic release of small bubbles is observed either (visually e.g. Fig. 4.9,
or statistically Fig. 4.8), we argue that the pulsatory behaviour of the lava lake at
Erebus volcano is driven by magma exchange between a shallow magma chamber
[Zandomeneghi et al., 2013] and the lake via bi-directional flow in the connecting
conduit.
It is interesting to note that, on average, bubble events in the data presented in
Fig. 4.7 occur every 5.5 min. This is comparable to the cycles in surface speed,
which range from ∼5–18 min in period. However, given that some cycles occur
without any bubbles surfacing (e.g. 09:00–09:15 in Fig. 4.7) and given the random
distribution of bubbles with respect to the phase of the cycles (Fig. 4.8), we believe
the similarity in timescales to be coincidental.
Pulsatory behaviour deriving from bi-directional flow in a conduit has been demon-
strated for single-phase systems using two fluids of different densities [Huppert and
Hallworth, 2007]. However, any exchange of magma occurring at the Erebus lava
lake will clearly be multi-phase, and its dynamics will be influenced not only by the
presence of gas bubbles but also by the large anorthoclase crystals which constitute
30–40% of the melt volume [Kelly et al., 2008]. Indeed, numerical simulations of
72
the Erebus magmatic system indicate that the inclusion of crystals has a very sig-
nificant effect on the flow dynamics [Molina et al., 2012]. While it is likely that gas
bubbles play an even more significant role than the crystals, a complete multi-phase
flow model of the Erebus system is not yet available. Whilst it is possible that the
dynamics observed by Huppert and Hallworth [2007] may not be applicable to a
complex multi-phase system such as that at Erebus, the lack of compelling evidence
for an alternative mechanism leads us to conclude that instability associated with
density-driven bi-directional flow is the most likely explanation for the observed
cyclic behaviour. As noted by Oppenheimer et al. [2009], the density contrast driv-
ing the flow is likely to result primarily from magma degassing.
Bouche et al. [2010] observed that bubbles in the lava lake at Erta ‘Ale volcano
may be trapped beneath the cooled crust at the surface of the lake and forced to
travel laterally until they encounter a crack in the crust before they can surface.
If such a process were also occurring in the Erebus lake, then it would invalidate
our comparison of the bubble events to the cycles in surface speed. The variable
duration of lateral migration of bubbles would prevent any direct comparison of the
timings of the bubble events and the phase of the cycles, since it would tend to
randomise their arrival at the surface. However, it can be observed in the IR images
that even small bubbles (1 m in diameter) break the surface of the Erebus lake
in areas with no visible cracks. We do not therefore believe that the crust on the
Erebus lake inhibits bubble ascent, nor that it causes significant lateral displacement
of bubbles. These differences likely reflect the contrasting rheologies of the magmas
involved, which in turn reflect the differences in composition, temperature, viscosity
and crystal content.
In our analysis of the correlation of bubble events to lake cycles (Section 4.2.5), we
have only looked at small bubbles in detail, since the dataset did not contain any
large events. Small bubbles may be sourced within the lake itself, whereas large
(decametric, causing ejection of material out of the lake) bubbles are thought to
have originated at greater depths [Oppenheimer et al., 2011, Burgisser et al., 2012].
It is possible that the passage of large bubbles through the conduit may perturb
the bi-directional flow of magma, causing variations in the period of lake surface
speed fluctuations. Although no such variation was observed in Fig. 4.12, we do not
believe this to be sufficient evidence to discount such a possibility. Since the arrival
of large bubbles is relatively infrequent, a time-series spanning several months would
need to be analysed to achieve a statistically significant sample size with which to
73
investigate possible effects of large bubbles on the lake’s motion.
4.5 Conclusions
In this chapter, we have reported an analysis of thermal infrared image data of the
active lava lake at Erebus volcano that spans seven field campaigns from 2004–2011.
In total 370,000 useful images were acquired representing 42 “field days” of obser-
vations and spanning contiguous observations of up to 44 h duration. The images
were analysed using a feature-tracking algorithm to determine surface motion vec-
tors from which the mean speed of the surface of the lake was derived, and this
parameter used to monitor the lake’s pulsatory behaviour. Shot noise in the mean-
speed data was found to indicate bubbles arriving at the surface of the lake, allowing
an analysis of the timings of bubble events with respect to the phase of the surface
speed cycles.
Since 2004, the apparent size (surface area) of the Erebus lava lake has decreased
by a factor of four. Despite these changes in the lake’s appearance, its pulsatory
behaviour has remained constant over the period of study, exhibiting cycles in mean
surface speed with periods in the range ∼5–18 min. Mean surface speeds are typ-
ically between 3 and 20 cm s−1. No obvious long-term progression of the cycles
was observed. Surface speed time-series are not symmetrical about their mean
(the troughs in speed are broader than the peaks), suggesting that the pulsatory
behaviour is due to intermittent perturbations of the system, rather than an oscil-
latory mechanism.
Bubbles arriving at the surface of the lake show no correlation to the phase of the
surface speed cycles. We therefore conclude that the pulsatory behaviour of the
lake is associated primarily with the flow dynamics of magma exchange within the
shallow plumbing system rather than by a flux of bubbles.
While we have analysed a substantially larger dataset than Oppenheimer et al.
[2009], we have still been limited by the intermittent coverage. We hope that our
recently-installed autonomous thermal camera system will yield much more extended
time-series, facilitating investigations into the effect of large (decametric) bubbles
on the pulsatory behaviour of the lake.
74
Chapter 5
Multi-parameter Observations of
Cyclic Lava Lake Behaviour
The work presented in this chapter was undertaken in order to meet the “Multi-
parameter Observations of Cyclic Behaviour” research objective (Section 1.7.3). As
discussed in chapters 1 and 4, numerous studies at Erebus volcano have recorded
pulsatory behaviour in many of the observable properties of its active lava lake.
However, although the periods of the cycles in the different properties recorded are
similar, suggesting they are linked, a lack of overlap between the different mea-
surements has prevented direct comparisons from being made. Using high time-
resolution measurements of surface elevation, surface speed, gas composition and
SO2 flux we demonstrate for the first time a definitive link between the cyclic be-
haviour in each of these properties. The cycles are found to be in-phase with each
other, with a small but consistent lag of 1–3 min between the peaks in surface ele-
vation and surface speed. Explosive events are found to have no observable effect
on the pulsatory behaviour beyond the ∼5 min period required for lake refill.
The majority of content from this chapter has been published as: Peters, N., C.
Oppenheimer, D. R. Killingsworth, J. Frechette, and P. Kyle, “Correlation of cycles
in Lava Lake motion and degassing at Erebus Volcano, Antarctica”, Geochemistry,
Geophysics, Geosystems, 2014
The collection and processing of the FTIR data presented in this chapter was per-
formed by Clive Oppenheimer. Collection and analysis of the TLS data was un-
dertaken by Drea Rae Killingsworth and Jed Frechette. The section describing the
FTIR methodology (Section 5.2.3) was written by Clive Oppenheimer, and is repro-
duced here with his consent.
5.1 Introduction
As discussed in Section 1.3, the pulsatory behaviour of the Erebus lava lake manifests
itself in many of its observable properties. In Chapter 4, we used the surface motion
of the lake as a metric to asses the inter-annual variability of the cyclic behaviour.
However, the choice of surface speed as a parameter was made purely on the basis of
an extended set of IR images being available. Observing how the cycles in different
parameters relate to each other is key in understanding the underlying mechanism
responsible for them. Since the cycles are believed to derive from processes occur-
ring at depth, a detailed knowledge of their behaviour will undoubtedly lead to an
improved understanding of how Erebus’s magmatic system as a whole functions.
Unfortunately, coordinating multi-instrument studies of volcanoes, such that the
time series of different data overlap, is a non-trivial exercise. As a result, previous
field campaigns at Erebus volcano have largely failed to concurrently record fluc-
tuations in more than one parameter. Similarities between surface motion and gas
composition were observed by Oppenheimer et al. [2009] in a short time series from
2004. However, due to the limited available data and a lack of high precision time-
stamping, a detailed analysis of the correlation was not possible. Thus, although it
has been widely assumed that the cyclic behaviour of the different parameters must
be linked, to date there has been no significant body of data to confirm this.
As the infrastructure installed by the Mount Erebus Volcano Observatory (MEVO),
such as power and data telemetry, continues to be improved, longer continuous
time series are becoming possible [Chapter 2, Peters et al., 2014c] and with them,
a greater degree of overlap between different observations. Here we present data
from December 2012 and December 2013 which show for the first time a definitive
link between the pseudo-periodic behaviour of the lava lake’s surface speed, surface
elevation, SO2 flux and plume composition.
In December 2013, the lake was approximately 30 m across (Fig. 5.1), with a sur-
face area of ∼580 m2 (measured using TLS). This is slightly larger than its surface
area in 2012 (∼400 m2), but is just a quarter of its 2004 size (Fig. 4.13). As dis-
cussed in Section 1.6, in 2013, Erebus was in a particularly active phase, with several
large explosive events occurring per day. Despite the increased activity, the lake’s
behaviour in periods of quiescence between explosions was indistinguishable from
other years. We believe therefore that the data presented here are representative of
Erebus’s behaviour in general rather than only during its periods of elevated activity.
76
30 m
Figure 5.1: Visible and infrared images of the lava lake on 15 December 2013. Freshly
ejected material surrounding the lake is visible in both images.
The findings of this chapter provide an important cornerstone for future studies,
and will be of particular relevance to future efforts to model pulsatory behaviour
either through experiments or computer simulations. They also provide a greater
body of evidence in support of the conclusions reached in Chapter 4.
5.2 Methodology
5.2.1 Lava Lake Surface Velocity
Using a thermal infrared (IR) camera situated on the northern rim of Erebus’s main
crater (∼330 m from the lake), images of the lava lake were recorded at a rate of
∼0.5 Hz. The thermal camera system used in this chapter is described in detail in
Chapter 2. Despite the small visible surface of the lake compared to many previous
years, estimates of the surface velocity using the same wavelet-transform feature
tracking algorithm as used in Chapter 4 were still possible. TLS data collected dur-
ing the 2013 field season were used to correct the IR images for perspective effects
and the mean lake surface speed was then calculated using the same methodology
as in Chapter 4.
Unlike the study presented in Chapter 4, there was no need to search the images for
periods of data with little camera shake. The new tripod system installed in 2012
[Chapter 2, Peters et al., 2014c] has solved this issue. Instead, periods of data with
good visibility of the lake (low humidity in the crater) which coincided with data
from the FTIR and the DW-FOV DOAS were chosen. This selection process was
77
done by manually inspecting the images.
Subsequent processing of the mean speed time series presented has been avoided.
The data shown have not been enhanced or filtered in any way. This is deliberate
so as to make the signal to noise ratio clear to see. The correlations between data
sets are still clearly visible “by eye” and we feel that this is sufficient to support our
conclusions.
5.2.2 SO2 Flux
Ultraviolet Differential Optical Absorption Spectroscopy (UV DOAS) is a powerful
tool for monitoring volcanic sulphur dioxide emissions [Galle et al., 2003]. Typically,
UV spectra are collected across a transect of the volcanic plume perpendicular to
its transport direction, either through the use of a scanning unit or by carrying
the spectrometer in a traverse [e.g. McGonigle, 2002, Salerno et al., 2009]. Sulphur
dioxide flux is calculated by integrating across the transect and multiplying by rise
speed (for a vertically rising plume) or by wind speed (for a horizontally drifting
plume). Although this technique is well established and commonly used within the
volcanological community, the necessity to record spectra at many positions across
the plume greatly restricts the time resolution of the flux measurements made. It is
therefore not possible to use the results in studies of dynamic phenomenon such as
Erebus’s pulsatory gas flux. Furthermore, the inclusion of plume transport speed,
which is often poorly constrained, in the calculation introduces large errors which
are difficult to quantify.
Increasingly, ultraviolet cameras [Mori and Burton, 2006, Bluth et al., 2007] are
being used to quantify SO2 emissions associated with short lived events such as
Strombolian eruptions [Tamburello et al., 2012, Dalton et al., 2009, Mori and Bur-
ton, 2009] and for comparison with other high time resolution monitoring data, such
as observations of seismic tremor [Nadeau et al., 2011]. Plume transport speeds can
be determined from the UV images either by a simple correlation of pixel values
between successive frames, or by more advanced optical flow methods [Chapter 7,
Kern et al., 2012, Peters et al., 2014a]. However, despite their promise as a tool for
high time resolution SO2 measurements, there remain many significant calibration
issues which must be addressed before accurate fluxes can be obtained [Kern et al.,
2010]. Although some of these issues have begun to be addressed by using a co-
located UV spectrometer to calibrate the images [Lu¨bcke et al., 2013], the technique
78
is still very much in its infancy. This, combined with the rather low SO2 flux of the
Erebus plume (∼0.5 kg s−1), makes measuring the periodic fluctuations using a UV
camera problematic. Instead, we took a similar approach to that of Boichu et al.
[2010], the so-called Dual Wide Field of View DOAS (DW-FOV DOAS) technique.
DW-FOV DOAS uses a pair of UV spectrometers connected by optical fibres to
two telescopes. Each telescope houses cylindrical lenses to give its spectrometer a
field of view that is very wide in one dimension (typically the horizontal) and very
narrow in the other. The fields of view of the two spectrometers are aligned parallel
along their long axes, with a small offset in the narrow field of view direction. The
width of the fields of view is such that they encompass the entire plume (the width
can be adjusted according to plume dimensions and viewing geometry) allowing the
column amount of sulphur dioxide to be measured instantaneously at two different
heights in the plume. Column amount measurements are made repeatedly with both
spectrometers and by cross-correlation of the two time series the plume transport
speed, and therefore flux, can be measured.
The spectrometers used in this study were OceanOptics USB2000+ instruments with
a spectral range of 283–427 nm. Mercury lamp calibration spectra were collected at
both the beginning and the end of the field campaign. We measured the instrument
response function (assuming a Gaussian shape) from the calibration spectra using
the AvoPlot software [Section 3.1, Peters, 2014] and found it to have a consistent
full width at half maximum of 0.6 nm for both spectrometers across the wavelength
range used for the SO2 column amount retrievals (∼310–324 nm).
SO2 column amounts were estimated using the standard DOAS methodology [Platt
and Stutz, 2008] whereby out-of-plume background spectra are averaged and used
to calculate the Ring spectrum (this is used to correct for the distortion of the so-
lar Fraunhofer lines due to rotational Raman scattering in the atmosphere). The
averaged background is then fitted to both a solar spectrum [Chance and Ku-
rucz, 2010] and the Ring spectrum calculated from the solar spectrum to esti-
mate the wavelength calibration errors in the spectrometer, the so-called “shift
and squeeze” correction. Each in-plume spectrum is then corrected for shift and
fitted to the Ring spectrum calculated from the averaged background and scaled
reference spectra for SO2 [Vandaele et al., 1994] and O3 [Burrows et al., 1999]
(the dominant absorbers in the wavelength range of the fitting window). The op-
timum scaling factor for the SO2 reference spectrum is directly proportional to
79
Figure 5.2: Photograph showing the plume conditions on 14 December 2013, with
annotations describing the viewing geometry of the DW-FOV DOAS instrument.
Although the plume is rising approximately vertically, there is a certain amount of
grounded plume around the crater which will cause the instrument to overestimate
the SO2 flux.
the amount of SO2 along the optical path of the spectrometer allowing column
amounts to be extracted. All DOAS retrievals were implemented using a com-
bination of the DOASIS software [Kraus, 2006] and scripts written by Tsanev
(http://www.geog.cam.ac.uk/research/projects/doasretrieval).
Plume rise speeds were determined using the same windowed cross-correlation ap-
proach as Boichu et al. [2010], whereby a fixed length window of column amount data
from the lower spectrometer is extracted at regular time intervals, normalised with
respect to its mean and standard deviation, and cross-correlated with the normalised
column amount data from the upper spectrometer (up to a certain maximum lag).
In this study we used a window length of 600 s and a maximum lag of 120 s. The
peak in the cross-correlation of the two signals is taken to represent the time taken
for the plume to rise the separation distance between the two fields of view, in our
case a distance of 36.5 m (Fig. 5.2). The rise speeds calculated were smoothed using
a low-pass filter before being used in the flux calculation.
The column amounts recorded by each spectrometer in the DW-FOV DOAS system
represent the mean column amount across its field of view at the distance of the
80
plume [Boichu et al., 2010]. The SO2 flux, φ, can therefore be expressed as:
φ =
104M
N
Av
A
lower
Dvθ
WFOV
(5.1)
where M is the molar mass of SO2 in kg mol
−1, N
Av
is Avogadro’s number, A
lower
is the column amount recorded by the lower spectrometer in molec cm−2, D is the
horizontal distance to the plume in metres (see Fig. 5.2), v is the rise speed in m s−1
and θ
WFOV
is the field of view of the spectrometer in the horizontal (wide) direction
in radians.
Note that our equation for flux differs slightly from that of Boichu et al. [2010],
since we have accounted for the increased path length through the plume due to the
inclined viewing angle. This results in an additional cos(α) term, which cancels with
that in the expression used by Boichu et al. [2010]. We use the column amounts
recorded by the lower spectrometer in our flux calculation to minimise the distance
above the lake that we are measuring.
The weather conditions at the summit of Erebus on 14 December 2013 were clear,
with no clouds within the fields of view of the spectrometers and light winds of be-
tween 1–3 ms−1. The plume was rising approximately vertically for a few hundred
metres above the crater rim, before drifting sideways with the wind. However, as
can be seen in Fig. 5.2, the plume was somewhat “messy” with a certain amount
of stagnated gas within the fields of view of the spectrometers. Whilst this will
not adversely affect the rapid variations in SO2 flux that we are interested in, it is
important to realise that the fluxes reported here may be overestimates compared
to the true values.
5.2.3 Gas Composition
Open-path absorption spectra of the gas emissions were collected with a MIDAC
Corporation M-4402-1 FTIR spectrometer sited on the crater rim, with the lava lake,
∼330 m line-of-sight distant, acting as infrared source [Oppenheimer and Kyle, 2008,
Oppenheimer et al., 2011]. All interferograms were recorded with a time step of 1.7 s
and a nominal 0.5 cm−1 spectral resolution. The measurements presented here were
made in December 2013. Retrievals of the species recognised in the absorption spec-
tra were made from transformed interferograms using a code that simulates and fits
atmospheric transmittance in discrete wavebands [Burton et al., 2000], specifically
81
2020–2100 cm−1 for H2O, CO2, CO and OCS. Uncertainties on these measurements
are typically about 5% [e.g. Horrocks et al., 2001]. Gas ratios are calculated from
the raw column amounts after corrections are made for contributions of atmospheric
water and CO2 [Oppenheimer and Kyle, 2008]. This procedure works particularly
well in this setting owing to the hyperarid background atmosphere.
We inspected the ratios of all possible pairs of the seven magmatic gas species de-
tected and found that the CO2/OCS ratio most clearly revealed the cyclic variation
sought. We have therefore focused on this quantity for the purposes of the multi-
parameter inter-comparison sought here.
5.2.4 Lava Lake Surface Elevation
An Optec ILRIS 3D terrestrial laser scanner was used to survey the lava lake sur-
face from the same location on the crater rim as the thermal camera and FTIR
spectrometer. The scanner uses a near-infrared laser (1535 nm wavelength), with a
scanning frequency of 2.5 kHz. At a distance of 300 m (the approximate distance
to the lake’s surface), the spot diameter of the laser beam is ∼7 cm. The scanner
was connected to a Trimble GPS pulse-per-second signal to provide accurate time
stamps for each data point.
The data used in this study consists of a series of repeated 30-s-duration scans,
focused on the lava lake surface. Scan data from within a circular region of 4 m
diameter, approximately centred in the middle of the lake, is averaged and used as a
measure of lake surface elevation. The density of points within the averaged region
is variable due to changing plume conditions, however, it typically contains around
90 measurements. The time taken to record these measurements (< 2 s) is small
compared to the timescale of change in lake surface elevation.
5.3 Results
Although cyclic behaviour is observed in many of the gas ratios measured with
the FTIR [Ilanko et al., 2014], in the 2013 dataset they are most apparent in the
CO2/OCS ratio time series and we therefore use this ratio to characterise the chang-
ing gas chemistry. Figure 5.3 shows four three-hour time series of mean lake surface
speed and CO2/OCS ratio from 6, 12 and 13 December 2013. Each three-hour
82
sequence contains a large (lake-emptying) explosion. These are highlighted, and
inset sequences of thermal images show the progression of each event with a ∼2 s
time-step. The correlation between the gas ratio and speed data is striking, and
clearly demonstrates that the cyclic behaviour of the two properties must be inher-
ently linked. Periods where there is no obvious correlation (e.g. around 16:00 on 12
December) correspond to times where there are no appreciable cycles in the time
series. Although there are periods where the peak in mean speed clearly precedes
the peak in CO2/OCS ratio (e.g. 18:30 on 6 December 2013), there are also periods
where the converse is true (e.g. 11:50 on 6 December 2013). A windowed correlation
analysis confirmed that there was no consistent lag between the signals. Given the
potentially complex coupling between the underlying mechanism driving the pul-
satory behaviour and the properties being measured, this is perhaps to be expected.
The large explosions shown in Fig. 5.3 do not appear to have any significant impact
on the cyclic behaviour of the lake. Following a short period (∼5 min, correspond-
ing to the time-scale of lake refill) of elevated speed and CO2/OCS ratio, both time
series are seen to return to similar values as prior to the explosion and resume their
pulsatory behaviour. No obvious changes in the period or amplitude of the cycles
following an explosion are apparent.
Figure 5.4 shows a little over two hours of SO2 flux, CO2/OCS ratio and mean lake
surface speed data collected on 14 December 2013. The SO2 flux data have been
offset by -3 min to account for the time taken for gas emitted at the lake’s surface to
rise to within the DW-FOV DOAS instrument’s field of view, a distance of ∼500 m
(the mean rise speed measured with the DW-FOV DOAS was 2.7 ms−1, giving a rise
time of ∼3 min from the lake’s surface to the field of view of the DW-FOV DOAS
instrument). However, it is important to realise that this is an approximation, and
the variability of the true delay time will result not only in a shift of the measured
flux compared to the real flux, but also in compression or rarefaction. The corre-
lation of the SO2 flux with the other signals is contingent on inhomogeneities in
gas flux at the lake’s surface rising coherently until they can be measured with the
DW-FOV DOAS. Turbulence and mixing in the lower portions of the plume are
likely to disrupt any structure in the emissions from the lake, and we believe that
this may account for the rather poor correlation of the data series at some times.
Despite this, the general agreement between the three time series is compelling, and
provides very strong evidence for the fluctuations in SO2 flux being driven by the
same underlying mechanism as those in the gas composition and lake surface mo-
83
Time(UTC)
Figure 5.3: Comparison between the CO2/OCS ratio measured with the FTIR and
the mean speed of the lake surface. IR images show the progression of explosive
events at a ∼2 s timestep.
84
Figure 5.4: Data from the 14 December 2013 showing the correlation between SO2
flux, CO2/OCS ratio and the mean speed of the lava lake surface. Cycles with
particularly good correlation are marked.
tion. Peaks in SO2 flux are observed to coincide with peaks in lake surface speed.
However, given the uncertainty in the time delay between gas being emitted from
the lake to when it is measured, particularly since this will vary with time, it is
not possible to investigate the alignment of peaks in any greater detail and there
may well be short lags between the signals which are not resolvable. There were no
explosions during the acquisition period of the SO2 flux data, so we cannot comment
on how they might influence the flux.
A ∼7 h time-series of mean lake surface speed and lake surface elevation recorded on
17 December 2012 is shown in Fig. 5.5. Again, the correlation between the cyclicity
in both properties is clearly visible. The sequence of IR images in the figure show
the progression of the large bubble burst that occurred at 17:24 UTC. Immediately
following the explosion, there is a ∼2 min gap in the TLS acquisition due to the
lake being obscured by the fumes from the explosion. The lake is then seen to refill
rapidly over the next ∼3 min and subsequently resume a pulsatory behaviour which
is indistinguishable from that exhibited prior to the explosion. This recovery time
is consistent with that of the large bubble events highlighted in Fig. 5.3 as is the
fact that the pulsatory behaviour does not appear to be influenced by the explosions
beyond the timescale of lake refill.
Also evident in Fig. 5.5 and shown more clearly in Fig. 5.6, is that the peak in surface
speed occurs after the peak in surface elevation. Figure 5.7a shows the results of
a windowed cross-correlation of the two series (as described in Section 5.2.2) using
85
15:00 16:00 17:00 18:00 19:00 20:00 21:00
Time (UTC on 17 December 2012)
6
5
4
3
2
1
0
1
2
R
e
la
ti
v
e
 s
u
rf
a
c
e
 e
le
v
a
ti
o
n
 (
m
)
0.00
0.02
0.04
0.06
0.08
0.10
0.12
M
e
a
n
 s
p
e
e
d
 (
m
s
−
1
)
Mean speed
Surface  e levation
Figure 5.5: Data from 17 December 2012, showing the surface elevation of the lake
(measured using TLS) and the mean surface speed (estimated from IR images). IR
images show the progression of a large explosion that occurs at 17:24 with a timestep
of ∼2 s.
18:33 18:43 18:53 19:03 19:13 19:23
Time (UTC on 17 December 2012)
3
2
1
0
1
2
R
e
la
ti
v
e
 s
u
rf
a
ce
 e
le
v
a
ti
o
n
 (
m
)
0.00
0.02
0.04
0.06
0.08
0.10
0.12
M
e
a
n
 s
p
e
e
d
 (
m
s−
1
)
Mean speed
Surface elevation
Figure 5.6: Expanded section of the data presented in Fig. 5.5, showing the lag
between the surface elevation and surface speed of the lake.
86
15:00 16:00 17:00 18:00 19:00 20:00
Time (UTC on 17 December 2012)
300
200
100
0
100
200
300
La
g 
(s
)
08:41 08:56 09:11 09:26 09:41 09:56 10:11 10:26 10:41 10:56
Time (UTC on 13 December 2013)
300
200
100
0
100
200
300
La
g 
(s
)
(a)
(b)
Figure 5.7: Results of a windowed cross-correlation using a window length of 1200 s
of: (a) the lake surface speed and the surface elevation on 17 December 2012; (b)
lake surface speed and CO2/OCS molar ratio on 13 December 2013. The plotted
data are the lags corresponding to the maximum in the cross-correlation at each
time step and therefore represent the lag between the two data series.
87
a window length of 1200 s and shows that the surface speed typically lags the the
surface elevation by 1–3 min. This contrasts with the relationship between the sur-
face speed and CO2/OCS ratio (Fig. 5.7b), which are generally in phase with each
other. Periods when the time series were very dissimilar and could therefore not
be cross-correlated successfully (e.g. due to the explosion at 17:24 on 17 December
2012) have been omitted from the figure for clarity.
5.4 Discussion
Given the apparent complexity of Erebus’s shallow magmatic system [e.g. Zan-
domeneghi et al., 2013] and the sensitivity of multi-phase flows to small changes in
system parameters [e.g. Molina et al., 2012], the regularity and longevity of the cyclic
behaviour of the lava lake is somewhat surprising. That such disparate methodolo-
gies as FTIR spectroscopy, motion estimation and DW-FOV DOAS are not only
all capable of independently detecting this behaviour but also demonstrate such a
striking correlation with each other, is remarkable, particularly given the challeng-
ing conditions under which the data were collected. The correlation between the
different data provides persuasive evidence that the cycles are not an artefact of the
measurement techniques.
Fig. 5.7 shows that the cycles in CO2/OCS ratio are in-phase with those in the lake’s
surface speed. In other words, the CO2/OCS ratio reaches its peak when the lake
motion is most vigorous, and its trough when the lake is most sluggish. This raises
the important question as to why the CO2/OCS ratio is varying. A comprehensive
analysis of the gas composition cycles is forthcoming in Ilanko et al. [2014] but we
note here that the two species are related through the following redox equilibrium:
3CO + SO2 
 2CO2 + OCS (5.2)
Whose equilibrium constant, K1 is given by:
K1 =
x
OCS
x
CO2
2
x
CO
3x
SO2
· P−1 (5.3)
where xa is the mole fraction of species a and P is the pressure.
We can also consider the equilibrium between CO2 and CO as follows:
88
CO +
1
2
O2 
 CO2 (5.4)
For which the equilibrium constant is given by:
K2 =
x
CO2
x
CO
√
f
O2
(5.5)
Rearranging, we can write the CO2/OCS molar ratio as:
x
CO2
x
OCS
=
K2
3
K1
· fO2
3
2
x
SO2
· P−1 (5.6)
There are thus several factors that can influence the CO2/OCS ratio, including
temperature (the equilibrium constants are temperature-dependent), oxygen fugac-
ity and pressure, in addition to the mole fraction of SO2 in the gas phase (itself
dependent on redox conditions, temperature and pressure). The system is under-
determined but can be solved if we assume (i) the measured gases represent an
equilibrium composition, and (ii) one of the unknown parameters. In this case, we
make the assumption that the gases are equilibrated at atmospheric pressure. This
enables us to compute, spectrum by spectrum, the equilibrium temperature and oxy-
gen fugacity. Under these assumptions, we may conclude that the cycles represent
variation of approximately 0.1 log units relative to the rock buffer Quartz-Fayalite-
Magnetite (QFM; Fig. 5.8). This does not explain the origin of this variation but
provides compelling evidence for a close association between magma degassing, re-
dox conditions and the fluid dynamics of the lava lake.
Measurements made by Ripepe et al. [2002] on Stromboli, revealed a correlation
between increased explosive activity and intensified gas puffing. Periods where de-
gassing consisted of vigorous puffing were found to coincide with a high number of
Strombolian explosions, whereas periods with a low number of Strombolian events
were characterised by fewer, lower-velocity gas puffs. Stromboli alternates between
these two different degassing regimes on timescales of 5–40 min, not dissimilar to
the ∼5–18 min cycles in SO2 flux observed on Erebus. Ripepe et al. [2002] proposed
that the observed behaviour at Stromboli was due to changes in the gas supply rate
into a shallow foam layer. On Erebus however, no correlation between explosive
activity and passive degassing behaviour is observed. Large (decametric) bubbles
are too infrequent to correlate with the cyclicity in the gas flux. For example, cyclic
variations in SO2 flux were observed in the data presented in Fig. 5.4, during a pe-
riod when no large bubbles were recorded. Although smaller (metre-scale) bubbles
89
08:33 09:03 09:33 10:03 10:33 11:03
Time (UTC on 13 December 2013)
0
2000
4000
6000
8000
10000
CO
2
/O
CS
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0.0
∆
QF
M
 (l
og
 u
ni
ts
)
CO2 /OCS
∆QFM
Figure 5.8: Comparison of the CO2/OCS molar ratio on 13 December 2013 and
the computed oxygen fugacity at which the gas equilibrated (relative to the QFM
buffer).
occur on a similar time-scale to that of the cycles in SO2 flux, these have been found
to be uncorrelated with the fluctuations in surface speed [Chapter 4, Peters et al.,
2014d] and by inference therefore, also with variations in the gas flux. Explosive
activity on Erebus appears to be entirely decoupled from its passive degassing be-
haviour. This is consistent with the findings of Burgisser et al. [2012], who propose
that while explosive degassing at Erebus is due to rupture of gas bubbles that rise
from depth independently of the melt, its passive degassing is dominated by ex-
solved and exsolving volatiles from magma in the lake itself. This would suggest
that the fluctuations in SO2 flux observed at Erebus (Fig. 5.4) are due to changes
in supply of gas-rich magma into the shallow system. However, in contrast to the
mechanism proposed by Ripepe et al. [2002] for Stromboli, we argue that significant
gas accumulation at a shallow level does not occur at Erebus. This prevents large
bubbles from forming in the shallow system due to coalescence, allowing significant
changes in SO2 flux without any changes in explosive activity.
A high degree of gas permeability of the Erebus lava lake was also proposed by
Peters et al. [2014d] (see also Chapter 4). This argument can be strengthened by
considering the emissions that occur during a single cycle. From the TLS data (e.g.
Fig. 5.5) we can see that a change in surface elevation of ∼2 m is typical during a
single cycle. Treating the upper part of the lake basin as cylindrical, with a surface
area of ∼400 m2 (its 2012 size), this represents a volume change of 800 m3. The
additional amount of SO2 released during a single cycle (above the background SO2
emission rate) is ∼5300 mol (found by integrating the flux data in Fig. 5.4). FTIR
90
0 200 400 600 800 1000 1200 1400
Time (seconds)
0
2000
4000
6000
8000
10000
12000
14000
V
o
lu
m
e
 (
m
3
)
Change in lake volume
Cumulative gas volume
Figure 5.9: Comparison of the volume change of the lava lake and the additional
volume (at atmospheric pressure) of gas released (above the background level) during
a single cycle.
measurements indicate a typical SO2/total-gas molar ratio of ∼113, giving a total
gas output of 1.4×104 m3 in a single cycle (assuming standard temperature and
pressure). These results are summarised in Fig. 5.9, which shows a comparison be-
tween the volume of gas released and the volume change of the lake during a single
cycle. It should be noted that overlapping SO2 flux data and TLS data are not
yet available, and so whilst the lake volume data are from 17 December 2012, the
gas volume data are calculated from SO2 flux data recorded on 14 December 2013.
However, since inter-annual variability in surface elevation variations and SO2 flux
is small at Erebus [Sweeney et al., 2008, Jones, 2013], it is reasonable to compare
these data. The data show that the total volume of gas released during a single cycle
is an order of magnitude greater than the volume change of the lake, supporting the
argument that the lake must be highly permeable to gas.
The 1–3 min lag between the peak in surface motion (and therefore also in gas
composition) and the peak in surface elevation may be indicative of the dwell-time
of gas in the lake. The lag may be explained by the time taken for fresh pulses of
magma injected into the lake to rise to the surface, but this would suggest that gas
escaping at the surface has travelled with the melt rather than rising independently
through the lake. It is important to remember however, that even at the low-point
of the cycles there is a net gas flux. If bubbles are not rising independently of the
melt, then there must be a sustained transfer of magma to the surface of the lake.
Since the surface motion is not seen to stagnate, this seems feasible. However, the
definition of “surface” in this argument is rather unclear, and it is not yet possible
to constrain either the depth at which magma motion might influence the surface
motion measured with the thermal camera, or the depth at which the gas might
91
separate from the magma.
Witham and Llewellin [2006] showed that stable lava lake systems would rapidly re-
turn to equilibrium following a perturbation. This is indeed what is observed in the
Erebus lake, with a rapid return to “normal” behaviour following explosive events.
It is interesting to note that the refill following an explosion does not appear to
“overshoot” the maximum lake level reached during periods of quiescence. Recov-
ery from an explosion is rapid, and shows that the response of the lake system to
a perturbation is much faster than the time-scale of the pulsatory behaviour. This
would suggest that the recorded waveform of the cycles in surface elevation (which
must be a convolution of the lake’s response function and the driving mechanism’s
waveform) are representative of the underlying process responsible for the cycles
and has not been overly smoothed due to a sluggish response of the lake system.
Although the waveforms of the surface motion and gas composition indicate that
the observed cyclicity is caused by periodic perturbations to a system that is other-
wise at equilibrium [Chapter 4, Peters et al., 2014d], the surface elevation appears to
be approximately symmetric about its mean, indicating a system that is oscillating
about its equilibrium state. It seems likely that the surface elevation of the lake
is a more “direct” measure of the underlying process responsible for the pulsatory
behaviour and we therefore believe it to be an oscillatory mechanism rather than an
impulsive one. One possible explanation is that the base of the lake is maintained
at magmastatic equilibrium. Fluctuations in surface elevation balance changes in
the bulk density of the lake caused by injection of bubbly (lower density) magma
from the conduit, and gas loss at the surface.
In Chapter 4, we noted that whilst the cycles in surface speed were not always clear,
there did always appear to be some periodicity present, concluding that the pul-
satory behaviour was a sustained feature of the lake. These findings are supported
by the data presented here. In particular, Fig. 5.5 (e.g. between 20:10 and 21:00)
demonstrates that even when pulsatory behaviour in the lake surface level is not
evident, there are still clear cycles in the surface speed. A similar phenomenon is
shown in Fig. 5.4 (e.g. between 04:08 and 05:08) where cycles are apparent in the
surface speed, but not in the CO2/OCS ratio. It is curious that at different times,
cyclic behaviour may be present in some properties but not others. We attribute
this behaviour to the fact that each property being measured is a different surface
manifestation of a common underlying periodic (or pseudo-periodic) process. The
coupling between the underlying process and the measured quantity could be ex-
92
tremely complex, and also susceptible to other forcing mechanisms. The clearest
example of this is the SO2 flux measurements. Even if the flux of SO2 is perfectly
periodic at the lake’s surface, changing atmospheric conditions within the crater will
alter how visible this periodicity is at the point of measurement and under certain
circumstances plume transport dynamics may dominate the flux, masking the peri-
odic signal entirely. A similar argument may be made for the FTIR measurements,
where the measured composition may be influenced by an accumulation of gases
within the crater over the time-scale of one lake cycle. In the case of the surface
elevation and surface speed measurements, it is possible that other processes within
the lake (e.g. thermal convection or foam formation) can overshadow the mechanism
responsible for the pulsatory behaviour, concealing its influence.
The discovery of periodic degassing at Mt. Etna in Sicily [Tamburello et al., 2013,
Pering et al., 2014], raises the possibility that, despite the differences in magma
viscosities, the models of magma flow being developed to explain the pulsatory
behaviour of the Erebus lava lake may be more widely applicable to open-vent vol-
canoes around the world. However, some care is needed when assessing periodic
degassing since it is possible that the development of turbulence within the plume
between the vent and the point of measurement (often a few hundred meters) could
cause a periodicity in the measured flux, even if the source is constant. In the data
presented here, we believe that this is not the case, as it would seem an unlikely
coincidence for turbulence effects to be so well correlated with the cycles occurring
in plume chemistry and lake surface motion.
5.5 Conclusions
In this chapter, we have presented measurements made during December 2012 and
December 2013 of the SO2 flux, CO2/OCS ratio, surface elevation and mean sur-
face speed of the Erebus lava lake. Well established methodologies (UV and FTIR
spectroscopy, TLS and motion estimation from thermal imagery), which have been
used during numerous previous field campaigns on Erebus, were employed. How-
ever, the unprecedented overlap in the time series presented has enabled us to show
for the first time a definitive link between the pseudo-periodic fluctuations observed
in each of these properties. Cycles (with an approximate period of 10 min) in SO2
flux, CO2/OCS ratio, oxygen fugacity and mean surface speed were shown to be
correlated, and in-phase with each other, providing strong evidence for a common
93
underlying mechanism. The cycles in surface speed were found to lag cycles in sur-
face elevation by 1–3 min. The volume of gas released during a single cycle (above
the background gas flux) is an order of magnitude greater than the associated change
in lake volume. This suggests that gas is effectively lost from the lake as soon as
it has entered. Large explosions appeared to have no appreciable effect on the pul-
satory behaviour observed, at least beyond the time taken for the lake to refill to its
normal level (∼5 min). However, insufficient events were analysed to rule out the
possibility of subtle changes to the behaviour caused by large gas bubbles, leaving
this as an interesting avenue for future investigation.
94
Chapter 6
Analysis of the 2D Surface
Velocity of the Erebus Lava Lake
This chapter looks in more detail at the structure of the velocity field over the
surface of the Erebus lava lake. Whereas the previous chapters have focused on
the mean velocity, here we examine the 2D velocity field using a combination of
vector plots and the divergence of the flow. Time periods from both 2004 and 2013
are analysed and found to have similar flow patterns despite the dramatic change
in lake surface area between the years. The flow typically consists of flow from a
single principal source region to a single principal sink region. No variation in flow
structure is observed during the cycles in the lake’s behaviour, nor does it appear
different during periods of little cyclic activity. The observations fit well with the bi-
directional flow mechanism proposed by Oppenheimer et al. [2009], however, they
also imply that there is a continuous flow of magma into the lake and that the
observed pulsatory behaviour is due to periodic intensifications of this flow.
6.1 Introduction
The previous chapters of this thesis have focused primarily on the study of spatially-
averaged metrics of the behaviour of the Erebus lava lake. In Chapter 4, we used
mean surface speed to analyse inter-annual variability in the lake’s behaviour and
in Chapter 5 we combined the mean speed data with other “bulk property” mea-
surements such as SO2 flux and mean surface elevation to demonstrate how the
behaviours of these properties are linked. However, inherent in the method of com-
puting the mean surface speed of the lake (see Section 4.2.3), is the calculation of
the two-dimensional velocity field over the whole of the lake’s surface. The analysis
presented previously has therefore discarded much of the information that is avail-
able from the IR dataset.
In this chapter, we revisit the velocity fields calculated in the previous chapters and
look in more detail at the two-dimensional structure of flow at the lake’s surface.
The flow structure is evaluated across complete cycles in mean lake surface speed
for data collected in both 2004 and 2013. As discussed in Chapter 4, despite the
significant change in lake surface area between 2004 and 2013, the overall behaviour
of the lake remained similar, with ∼5–20 min period cyclicity evident. We also com-
pare the flow structure between a time period in 2013 with no evident cyclicity and
a ensuing period demonstrating strongly cyclic behaviour. The structure of the flow
during these periods is analysed using both vector plots and plots of the divergence
of the flow field.
6.2 Methodology
The IR images analysed in this chapter were collected using the same thermal cam-
era systems as described in Chapter 4 [see also Peters et al., 2014d]. The calculation
of the lake surface velocities followed the same methodology as used in the previous
chapters (see Section 4.2.3).
6.2.1 Data Selection
Time periods for analysis were chosen by considering plots of the mean surface speed
of the lake (e.g. Fig. 5.3). Periods where the cyclic behaviour was particularly evi-
dent were used in order to demonstrate how the velocity field changed during the
96
progression of a cycle. Data from both December 2004 and December 2013 were
used in order to compare the velocity fields for different sizes of lake (∼2300 m2 in
2004 and ∼580 m2 in 2013).
A period of data from 15:45–16:45 UTC on 12 December 2013 was also chosen for
analysis as it shows a transition between a period of no evident cyclicity and well
defined cycles (see Fig. 5.3).
Visibility of the lake was good during all time periods used (low humidity in the
crater) and there was no apparent camera shake.
6.2.2 Divergence
Evaluating temporal changes in a 2D velocity field is challenging due to the dif-
ficulty of visualising such data. In previous chapters, this issue was addressed by
summarising the velocity field using its mean. However, this gives little information
about the structure of the velocity field. Here, we use the divergence of the velocity
field as a means of assessing the progression of the surface motion of the lake. Di-
vergence was chosen since it not only provides information about the 2D structure
of the velocity field, but also has a simple physical interpretation.
Given a 2D vector function of the form V(x, y), the divergence of the function is
defined [e.g. Boas, 1983b] as:
div V = ∇·V = ∂Vx
∂x
+
∂Vy
∂y
(6.1)
The physical interpretation of the divergence at a particular point is the extent to
which that point behaves as a source or sink. A point in a vector field with positive
divergence is a source, whereas negative divergence indicates a sink region. In the
context of lake surface motion (where the vector field under consideration is the
displacement field), source regions correspond to spreading centres associated with
lava up-welling and sink regions correspond to “subduction zones” associated with
lava down-welling.
97
6.3 Results
Figure 6.1 shows a short time series of mean lake surface speed calculated from IR
images recorded on 17 December 2004. On this date, the lake was ∼2300 m2 in size.
The time series spans a single cycle in mean speed and plots of the 2D flow field
and its associated divergence are shown at different points in the cycle. Although
the flow field is rather complex, with many source and sink regions, the majority of
these are transient and the flow structure is dominated by flow from a source region
in the upper right of the lake to a sink region in the lower left. This flow structure
is consistent across cycles in the mean surface speed of the lake.
Figure 6.2 shows how the 2D velocity field and its divergence changes across a single
cycle in mean speed on 13 December 2013. At this time, the lake was ∼580 m2 in
size, just a quarter of its size in 2004. The data are generally smoother than those
shown in Fig. 6.1. This is likely due to the superior temporal and spatial resolution
of the thermal camera system used in 2013 compared to that of 2004 (0.5 Hz and
640×480 pixels compared to 0.25 Hz and 320×240 pixels). The figure shows that the
lake surface has one pronounced source region (lower left) and one pronounced sink
region (upper right). The configuration of these regions does not change throughout
the duration of the cycle in mean speed. There is no reversal of the flow between the
leading and trailing edges of the peak in mean speed. The source and sink regions
are observed to intensify at the peak of the cycle. However, this is to be expected
given the increase in magnitude of the surface speeds.
Figure 6.3 shows the 2D velocity field and associated divergence at different times
during a ∼1.25 hr period of mean surface speed data on 12 December 2013. During
the first ∼30 min of the timeseries there is no apparent cyclic behaviour in the lake’s
mean surface speed, whereas, for the latter part of the timeseries there are clearly
defined cycles. Despite this dramatic change, the structure of the velocity field varies
little across the transition in behaviour of the lake, with flow being predominantly
from a source region in the lower left to a sink region in the upper right. As can
be seen from the plot of mean speed in the figure, during the period with no cyclic
behaviour, the mean speed of the lake’s surface is comparable to the troughs in mean
speed during the period with well defined cycles.
98
86
4
2
0
2
4
6
8
D
iv
e
rg
e
n
ce
×10 2
8
6
4
2
0
2
4
6
8
D
iv
e
rg
e
n
ce
×10 2
8
6
4
2
0
2
4
6
8
D
iv
e
rg
e
n
ce
×10 2
8
6
4
2
0
2
4
6
8
D
iv
e
rg
e
n
ce
×10 2
04:19 04:24 04:29
Time (UTC on 17 Dec. 2004)
0.00
0.05
0.10
0.15
0.20
M
e
a
n
 s
p
e
e
d
 (
m
s−
1
)
Figure 6.1: Comparison of the 2D velocity field (black arrows) and associated diver-
gence (colour scale) across a single cycle in the mean surface speed of the lake com-
puted from IR images recorded on 17 December 2004. A positive divergence (red)
is indicative of a spreading/up-welling region while a negative divergence (blue) in-
dicates a converging/down-welling region. Arrow lengths are proportional to speed.
99
6.0
4.5
3.0
1.5
0.0
1.5
3.0
4.5
6.0
D
iv
e
rg
e
n
ce
×10 2
6.0
4.5
3.0
1.5
0.0
1.5
3.0
4.5
6.0
D
iv
e
rg
e
n
ce
×10 2
6.0
4.5
3.0
1.5
0.0
1.5
3.0
4.5
6.0
D
iv
e
rg
e
n
ce
×10 2
6.0
4.5
3.0
1.5
0.0
1.5
3.0
4.5
6.0
D
iv
e
rg
e
n
ce
×10 2
09:08 09:13 09:18
Time (UTC on 13 Dec. 2013)
0.00
0.01
0.02
0.03
0.04
0.05
0.06
0.07
0.08
M
e
a
n
 s
p
e
e
d
 (
m
s−
1
)
Figure 6.2: Comparison of the 2D velocity field (black arrows) and associated diver-
gence (colour scale) across a single cycle in the mean surface speed of the lake com-
puted from IR images recorded on 13 December 2013. A positive divergence (red)
is indicative of a spreading/up-welling region while a negative divergence (blue) in-
dicates a converging/down-welling region. Arrow lengths are proportional to speed.
100
15:50 16:00 16:10 16:20 16:30 16:40 16:50
Time (UTC on 12 Dec. 2013)
0.00
0.01
0.02
0.03
0.04
0.05
0.06
0.07
0.08
M
ea
n 
sp
ee
d 
(m
s−
1
)
6.0
4.5
3.0
1.5
0.0
1.5
3.0
4.5
6.0
Di
ve
rg
en
ce
×10 2
6.0
4.5
3.0
1.5
0.0
1.5
3.0
4.5
6.0
Di
ve
rg
en
ce
×10 2
6.0
4.5
3.0
1.5
0.0
1.5
3.0
4.5
6.0
Di
ve
rg
en
ce
×10 2
6.0
4.5
3.0
1.5
0.0
1.5
3.0
4.5
6.0
Di
ve
rg
en
ce
×10 2
Figure 6.3: Comparison of the 2D velocity field (black arrows) and associated diver-
gence (colour scale) between a time period with no apparent cyclic behaviour (15:45
– 16:25 UTC) and a time period with well defined cycles (16:25 – 17:00 UTC) on 12
December 2013.
101
6.4 Discussion
When comparing the divergence fields calculated from images from 2004 to those of
2013 it is tempting to conclude that the principal up-welling region, and therefore
the location of the conduit, has moved from the upper-right to the lower-left of the
lake. However, it is important to remember that the fields of view of the camera
systems recording the images are different between the two years, and there is no
way to assess how the lake boundaries in each year relate to each other. Given the
fairly rapid variations in the topography of Erebus’s inner crater that have been ob-
served [Jones, 2013], it seems plausible that the positioning of the conduit could vary
between years. However, there is insufficient data to rule out the possibility that
the conduit remains fixed, and that it is the relative position of the lake that changes.
Peters et al. [2014d] (see also Section 4.3) found that the CDF of bubble events
is shifted towards higher mean speeds than a reference of random samples of lake
speeds, suggesting that bubbles tend to surface when the lake surface is moving
faster than average. This was attributed to the failure of the low-pass filter used in
the analysis to completely remove the spikes in speed caused by the bubble events.
However, animations of the 2D velocity plots and associated divergences (available
on request) show that, following a bubble, there is rapid flow into the void space
left by the bubble. Although very localised, this region of high velocity will skew
the mean speed value, and may account partially for the offset of the CDF.
In Section 1.4, proposed mechanisms for lava lake degassing were discussed. In Chap-
ter 4, we discounted the gas pistoning and bubble-driven pressurisation changes as
candidates for the Erebus system based on the absence of periodic gas releases
from the lake and also because of the symmetry of the Erebus cycles in time. The
structure of the velocity field supports these conclusions, since both aforementioned
mechanisms require uni-directional upflow followed by uni-directional downflow of
magma in the conduit. No reversal of the surface flow pattern of the lake is observed,
suggesting that these mechanisms are not active at Erebus.
The observed flow structure at Erebus is indicative of a stable and continuous ex-
change flow between the lava lake and the deeper magmatic system. The cycles that
are observed are due to an intensification of this flow rather than a change in flow
regime. Such observations are not at odds with the bi-directional flow mechanism
proposed by Oppenheimer et al. [2009], however, what is not made clear by Op-
penheimer et al. [2009] is that underpinning the pulsatory flow of magma into the
102
lake is a relatively stable baseline flow which is present even when cyclic behaviour
is not. Furthermore, it is possible, given the frequency of bubble events presented
in Chapter 4, that the perturbations to the baseline flow are caused by detached
bubble wakes rather than flow instability. However, this hypothesis seems less likely
since cycles are occasionally observed without any evidence of bubbles arriving at
the surface (e.g. Figure 4.7).
6.5 Conclusions
In this chapter, we have presented estimates of the surface velocity field of the Ere-
bus lava lake in December 2004 and December 2013. The same motion estimation
methodology as employed in Chapter 4 and Chapter 5 was used. However, here, the
structure of the flow field was analysed rather than just the mean surface velocity.
In addition to vector plots of the surface velocities, the divergence of the flow field
was used as a tool to assess changes in its structure. It was found that the flow field
did not change in structure during the course of a cycle. Specifically, there was no
reversal of flow or changes in the number of source and sink regions. The surface
flow of the lake in 2004 was broadly similar to that in 2013 with one predominant
source region and one predominant sink region, despite the large change in lake
surface area between the years. The structure of the surface flow field was found
to be similar between periods of strongly cyclic activity and periods of little cyclicity.
The continuity of the upwelling region in the lake, and lack of evidence for any flow
reversal precludes gas pistoning and the Witham et al. [2006] pressurisation change
mechanism as being responsible for the observed cyclicity of the lake. Although not
directly contradictory to the bi-directional flow model proposed by Oppenheimer
et al. [2009] (see also Section 1.4), the results suggest that the pulsatory delivery of
magma into the lake is superimposed on top of a relatively stable baseline flow. It
appears that the cyclic behaviour of the lake is caused by periodic intensifications
of this baseline flow.
103
Chapter 7
Use of Motion Estimation
Algorithms for Improved SO2
Measurements Using UV Cameras
UV cameras are rapidly gaining popularity as a tool for monitoring SO2 emissions
from volcanoes. A variety of different SO2 camera systems have been developed
with varying patterns of image acquisition in space, time and wavelength. Despite
this diversity, there are two steps common to the workflows of most of these sys-
tems; aligning images of different wavelengths to calculate apparent absorbance and
estimating plume transport speeds, both of which can be achieved using motion
estimation algorithms. In this chapter we present two such algorithms, the Dual
Tree Complex Wavelet Transform-based algorithm used in the previous chapters
and the Farneba¨ck Optical Flow algorithm. We assess their accuracy using a syn-
thetic dataset created using the numeric cloud-resolving model ATHAM, and then
apply them to real world data from Villarrica volcano. Both algorithms are found
to perform well and the ATHAM simulations offer useful datasets for benchmarking
and validating future algorithms.
The majority of content from this chapter has been published as: Peters, N., A.
Hoffmann, T. Barnie, M. Herzog, C. Oppenheimer, “Use of Motion Estimation
Algorithms for Improved Flux Measurements Using SO2 Cameras”, Journal of Vol-
canology and Geothermal Research, in press 2014
The development and execution of the ATHAM model to produce synthetic UV
images was conducted entirely by Alex Hoffmann (with supervision from Michael
Herzog), and other than providing some advice on initialisation parameters, I did
not contribute to this part of the work. The text describing the model (Section 7.2.3)
was also written by Alex Hoffmann, and is reproduced here with his consent. All
post-processing of the model data (including development of the processing codes),
and all other content of this chapter, is my own work.
7.1 Introduction
Ultraviolet (UV) cameras are the latest addition to the family of UV spectroscopy
techniques used to measure sulphur dioxide (SO2) emissions. Since the first vol-
canological demonstrations [Mori and Burton, 2006, Bluth et al., 2007] they have
been used in numerous studies of SO2 fluxes from volcanoes [e.g. Holland et al., 2011,
Tamburello et al., 2011, Campion et al., 2012, Smekens et al., 2013, La Spina et al.,
2013]. Using UV cameras, estimates of SO2 flux at ∼1 s temporal resolution are
possible, enabling studies of highly dynamic degassing events such as Strombolian
eruptions [Tamburello et al., 2012, Dalton et al., 2009, Mori and Burton, 2009] and
pulsatory degassing [Tamburello et al., 2013, Pering et al., 2014].
In order to measure SO2, UV cameras typically use bandpass filters to image at two
narrow wavebands, one within and one outside-of an absorption band of SO2 (typi-
cally 310 and 330 nm respectively). The cameras operate on the principle that the
apparent absorbance between the two bands will vary only as a function of atmo-
spheric species that absorb unequally across them [Mori and Burton, 2006]. How-
ever, the calibration of apparent absorbance to SO2 column amount is non-trivial
[Kern et al., 2010], and calibration curves are typically estimated from imaging gas
cells of known concentration, or UV spectrometer retrievals coincident with a region
in the absorption image [Lu¨bcke et al., 2013, Kern et al., 2013].
Images at the two wavebands are either captured simultaneously, necessitating two
separate cameras, or sequentially, using a filter wheel to change the recorded wave-
band for each image. The latter method permits a single camera unit to be used, but
at the expense of having a short time delay between the in-band and out-of-band im-
ages. Since volcanic plumes are seldom stationary in the atmosphere, this is clearly
undesirable. Heterogeneities may have shifted by several pixels between consecu-
tive images, compromising the SO2 retrieval method described above. Simultaneous
capture of images at two wavebands is also problematic, since the retrieval of SO2
column amounts requires the field of view of each pixel in the in-band image to
match that of the corresponding pixel in the out-of-band image. Any misalignment
of the optical planes of the two camera units and any differences in lens distortion
between them will mean that this is not the case and will lead to errors in the
SO2 column amounts calculated. UV cameras are typically used “in the field” and,
as such, are often subjected to rapid temperature changes, vibration and physical
shocks. Alignment of optical components to the degree of single-pixel correspon-
dence is simply not practical.
105
A further problem, also applicable to other gas flux measurement techniques [e.g.
McGonigle et al., 2005], is that of estimating plume transport speed. To calculate
SO2 flux from UV camera images, the SO2 column amounts must be integrated
across a transect of the plume and multiplied by the transport speed perpendicular
to the transect. Typically, the plume speed is either estimated manually [Bluth et al.,
2007], which is unfeasible for long datasets or, more commonly, by cross-correlating
the integrated column amount values from two complete transects of the plume at
different distances from the vent [Kantzas et al., 2010]. The latter method assumes
a uniform velocity field across the entire plume and does not account for complex
flow structures or vorticity. One has only to observe a volcanic plume, to realise that
this is a poor assumption to make. Boichu et al. [2010] discuss in more detail the
compromises that must be made between time- and motion-resolution when using
the cross-correlation method. An additional, and perhaps more serious issue, is that
the lack of reliable “ground truth” data means that the errors introduced by using
this method are largely unconstrained.
Estimating plume speed and accurate registration of images from two separate cam-
eras are essentially the same problem; that of calculating motion vectors which map
pixels in one image to their corresponding positions in another. This is a common
computer-vision problem, which is addressed by a large number of well-established
motion estimation algorithms of varying complexity and trade-offs [e.g. Geiger et al.,
2014]. To date, use of these techniques in volcanology has mostly been constrained to
estimating lava motion from infrared images [James et al., 2007, Oppenheimer et al.,
2009, Lev et al., 2012, Peters et al., 2014d], and their application to UV camera im-
ages has so far been limited to the work of Kern et al. [2012] who demonstrated how
optical-flow methods may be used to estimate SO2 flux. We extend this work here
by assessing the performance of two different motion estimation algorithms when
applied to UV images, both for flux estimation and for registering images from sep-
arate cameras. We use both a Dual-Tree Complex Wavelet Transform (DT-CWT)
algorithm [Magarey and Kingsbury, 1998], and the Farneba¨ck optical-flow algorithm
[Farneba¨ck, 2003].
As already discussed, a major problem with estimating plume transport speeds is
that of validation. With no possibility of accurately measuring plume speed by an
alternative method, it is not possible to assess the performance of the techniques
being developed to estimate plume speeds, nor is it possible to constrain the errors
106
which they may introduce. In an attempt to address this issue, we use the Active
Tracer High-Resolution Atmospheric Model (ATHAM) run in a two-dimensional
configuration to create a set of synthetic images of a volcanic plume from a pas-
sively degassing volcano. Algorithms for estimating plume speed can be applied to
the synthetic images and the results compared to the “real” velocity field computed
by the model, providing a simple method of assessing their performance in the con-
text of volcanic plume observations. We demonstrate the technique with both the
DT-CWT and Farneba¨ck motion estimation algorithms.
It is not our intent to conduct an in-depth comparison of the many different motion
estimation algorithms available, and we acknowledge that more suitable algorithms
may exist for use with UV camera images than the two which we present here.
Instead, we hope to demonstrate the potential applications and benefits of motion
estimation algorithms in volcanic-plume UV imaging, and establish a framework for
validating and comparing them within this context.
Our specific aims are threefold: (i) to introduce the ATHAM model as a tool for
validating plume transport speed estimation methods, (ii) to present two different
motion estimation algorithms which are suitable for use with images from UV cam-
eras and show how they perform on both real and synthetic data, (iii) to demonstrate
how these algorithms can be employed to provide accurate registration of images
from separate camera units.
7.2 Methodology
7.2.1 Motion Estimation Algorithms
The DT-CWT algorithm works by decomposing images using the Dual-Tree Com-
plex Wavelet Transform [Kingsbury, 2001]. The phase shift between the decomposed
subimages of two frames is related to the motion of features between those frames.
The algorithm employs an iterative approach, using coarse features to create an ini-
tial motion estimate and then refining this using progressively finer-scale features.
Full details of the algorithm are given by Magarey and Kingsbury [1998].
The DT-CWT algorithm has been used in previous volcanological studies to estimate
the surface motion of the lava lake at Erebus volcano from time series of infrared
images [Chapter 4 Oppenheimer et al., 2009, Peters et al., 2014d]. Whilst this is
107
clearly a very different scenario to that of estimating plume transport speeds from
UV images, many of the fundamental problems are the same; identifiable features are
often diffuse, transient, and may change significantly between frames. The algorithm
copes well under these conditions, and is also robust against the noise introduced
by varying amounts of volcanic plume, which degrades the camera’s view of the lake.
The Farneba¨ck algorithm works by expressing the local neighbourhoods of each pixel
in the images as polynomial expansions. The translations between the polynomial
expansions for each successive frame can then be calculated and, by enforcing a cer-
tain model of motion (e.g. the affine motion model) the two-dimensional displace-
ment field can be estimated. The Farneba¨ck algorithm also employs an iterative
approach, with a starting estimate for motion being derived from coarse features
and then refined using progressively finer features. The algorithm is described in
detail by Farneba¨ck [2003].
In this investigation we used the same implementation of the DT-CWT algorithm
as in Chapters 4 and 5. MatLab code for this was kindly provided by Prof. Nick
Kingsbury at the Engineering Department of the University of Cambridge. We
used the implementation of the Farneba¨ck algorithm provided by the open-source
computer-vision library OpenCV [Bradski, 2000].
Both of these implementations allow several parameters to be set by the user to
control the operation of the algorithm. A description of the parameters along with
the specific values used for this study are given in Section. 7.2.2. Broadly speaking
the parameters set both the averaging level and the scale of features used by the
algorithms. Selecting higher averaging and larger features will usually result in a
more robust (less prone to noise, more tolerant of transient features) motion esti-
mate, but at the expense of a less detailed motion field. Given a set of images for
which the motion vectors are known (e.g. the data presented in Figs. 7.2 and 7.4), it
would be possible to optimise an algorithm’s performance by systematically varying
its input parameters and minimising the mean-squared error with respect to the
known velocities (e.g. using Levenberg-Marquardt optimisation). However, for this
study we have not employed such a rigorous approach and instead selected suitable
parameters by “trial and error” using a small subset of our data. We believe this
approach to be sufficient to support our conclusions, however, it should be noted
that the velocity estimates presented could likely be improved upon by using a more
thorough parameter selection process.
108
7.2.2 Algorithm Parameters
Below, we describe the input parameters for the two motion estimation algorithms,
and list the values that were assigned to them for this investigation. We refer to the
parameters by the same names as are used in the source-code and documentation for
the specific implementations that we used (see Section 7.2.1). Optional parameters
that have not been defined below were left set to their default values. Different values
were used for the image registration (Fig. 7.10) and the plume speed estimations
(Figs. 7.2, 7.4 and 7.8). Both sets of values are listed in the tables below.
DT-CWT Algorithm
Parameter Value for Flux Estimation Value for Registration
w [−pi/2.15,−3pi/2.15] [−pi/2.15,−3pi/2.15]
nit 6 6
levelsel [[6, 4], [6, 3], [5, 3], [5, 3]] [[4, 4], [4, 3], [4, 3], [4, 3], [4, 2]]
sizeqfilt [4, 2, 2, 1] [48, 32, 32, 16]
avlevel 3 4
Table 7.1: Parameters used with DT-CWT algorithm for both image registration
and plume speed estimation.
w Expected phase rotation per sample of the high and low band filters used in the
DT-CWT decomposition.
nit Number of iterations of the algorithm. A higher value may lead to a more
refined motion estimate, but at the expense of increased computation time.
levelsel Defines the coarsest and finest CWT levels to be used during each iteration.
In general it is best to use only the coarse levels for the early iterations, and
then gradually expand the range to progressively finer levels. This parameter
effectively selects the “size” of the features that are being tracked. Use of larger
features often gives estimates that are less prone to noise, but also results in
a less detailed velocity field.
sizeqfilt Defines the size of the smoothing filter used at each iteration. A larger
filter gives a smoother velocity field, but also causes loss of detail.
avlevel Defines the CWT level which has the same resolution as the velocity esti-
mates.
109
Farneba¨ck Algorithm
Parameter Value for Flux Estimation Value for Registration
pyr scale 0.5 0.5
levels 4 3
winsize 20 15
iterations 5 10
poly n 7 5
poly sigma 1.5 1.2
flags OPTFLOW FARNEBACK GAUSSIAN OPTFLOW FARNEBACK GAUSSIAN
Table 7.2: Parameters used with Farneba¨ck algorithm for both image registration
and plume speed estimation.
pyr scale Scale factor for each pyramid level compared to the previous.
levels Number of pyramid levels to create. A higher number results in coarser
features being used for the initial motion estimate.
winsize Defines the size of the averaging filter used. A larger filter gives a smoother
velocity field, but also causes loss of detail.
iterations Defines the number of iterations of the algorithm at each pyramid level.
A higher value may lead to a more refined motion estimate, but at the expense
of increased computation time.
poly n Size of the neighbourhood used in the polynomial expansion for each pixel.
A higher value results in more robust, but less detailed motion estimates.
poly sigma The standard deviation of the Gaussian smoothing filter used during
the polynomial expansion.
flags Optional parameters. We used the OPTFLOW FARNEBACK GAUSSIAN option so
that the algorithm uses a Gaussian filter rather than a simple “box” filter. This
gives more robust motion estimates at the expense of increased computation
time.
7.2.3 Using ATHAM as a Validation Tool
We used the non-hydrostatic, fully compressible, multiphase, mesoscale Active Tracer
High-Resolution Atmospheric Model (ATHAM) to simulate a synthetic SO2 plume
from passive degassing of an idealised volcano, loosely based on Villarrica. Built
110
within the Large Eddy Simulation framework, ATHAM was originally designed for
simulating plinian eruption dynamics [Oberhuber et al., 1998, Herzog et al., 1998,
2003]. Within a volcanic context, it has been used for research into volcanic gas
scavenging, particle aggregation and interactions with microphysics [Textor et al.,
2006a,b], and more recently, into the dynamics of large co-ignimbrite [Herzog and
Graf, 2010] and phreatomagmatic [Van Eaton et al., 2012] eruption clouds.
ATHAM is based on a modular structure that facilitates the coupling of specific
process modules to the dynamical core. The core solves, within an Eulerian frame-
work, the Navier-Stokes equations of momentum, and the pressure, temperature
and tracer transport equations for a gas-particle mixture [Oberhuber et al., 1998].
Active tracers, i.e. solid, liquid or gaseous components that have their own heat ca-
pacities and densities, influence the mixture’s flow by changing its thermodynamic
and dynamic properties. Momentum and tracer equations are described in flux form
to guarantee conservation of momentum and mass; the heat transport equation is
in advective form.
Short- and longwave radiative transfer, as well as surface fluxes of sensible and latent
heat, have not been applied in our model runs, to reduce complexity and computa-
tional expense. Resolved and unresolved turbulent entrainment of ambient air into
the plume modifies its buoyancy and disperses tracers and temperature anomalies.
Subgrid-scale (SGS) turbulence is treated as eddy diffusivity and parametrised on
the basis of a 1.5-order prognostic SGS anisotropic turbulent kinetic energy (TKE)
closure scheme, described by Herzog et al. [2003]. The mixed-phase cloud micro-
physics included in the model correspond to a simple prognostic bulk Kessler-type
approach, predicting cloud droplets, rain, ice crystals and graupel as specific con-
centrations [Herzog et al., 1998]. The current model configuration does not include
aerosol microphysics and nucleation, cloud nuclei activation, or ash aggregation.
For our validation exercise, we used the same volcano model configuration as in the
previously cited studies. The volcanic topography was parametrised as a simple
Gaussian-shaped mountain with a height of 2850 m and a width parameter chosen
such that the height-to-base ratio (500/2000) captured within our model domain
roughly matched that of the Villarrica edifice. The crater was also chosen to be
Gaussian-shaped (inverted), with a depth of 100 m, a width parameter of 50 m, and
a flat bottom. Emissions from Villarrica typically emanate from large (5–20 m)
vents in the cooled-crust covering the ∼40 m-diameter lava lake that occupies the
111
base of its crater [Witter et al., 2004, Shinohara and Witter, 2005]. To emulate heat
and volcanic gas fluxes from the lake, we applied a vertical velocity forcing over
the cross-section of a lake, assumed to have a 50 m diameter. The baseline “exit”
velocity was prescribed at 1 ms−1, ramped up over 30 s, then modulated in time by
a sine function of the same amplitude and a period of 2 min, mimicking periodicity
in lake activity. Additionally, we superimposed random heterogeneity at each of
the lake’s 25 grid-points of ±50% of the “exit” velocity. The potential temperature
of the erupting gas-particle mixture above the lake was set to 1000 K, and specific
gas and aerosol concentrations were aligned with the values measured by Sawyer
et al. [2011]. From the tabled mass ratios of emitted gas species, we set the specific
concentrations of the two gaseous tracers included in our model runs, magmatic
H2O and SO2, to 782 and 79 g kg
−1, respectively. From the ratio of the measured
mass fluxes of SO2 (3.7 kg s
−1) to aerosol (1 kg s−1), we inferred a total aerosol active
tracer specific concentration of 21 g kg−1. The remainder was assumed to be dry air
from unresolved near-surface turbulent entrainment. All particulate aerosol was as-
signed to a single size bin with a radius of 0.5µm and in equal amounts of 7 g kg−1
to 3 size bins with radii of 0.1, 0.5 and 5µm in a second simulation as described
in Section 7.3. All aerosol was given a specific heat capacity of 1200 J kg−1K−1 and
a density of 2000 kg m−3. Note that the actual values chosen to characterise vol-
canic degassing are not overly important, as the primary objective of the modelling
exercise is to provide validation for motion estimation algorithms, rather than to
constrain accurate plume dynamics and composition of a specific volcano. The tur-
bulent length scale, and the horizontal and vertical TKE components over the lake
were set to 2 m and 2 m2s−2, respectively.
To account for wind-driven advection, without incurring excessive computational
expense, we opted for a two-dimensional Cartesian grid, acknowledging that this
approach precludes three-dimensional turbulent entrainment and therefore limits
effective plume dispersion. An equal and homogeneous grid-spacing of 2 m in both
the horizontal and vertical dimensions was used, roughly matching the 1 m2 pixel
size of our Villarrica dataset (Section 7.2.4). We used 1010×751 grid-points to cover
a domain 2018 m wide and 1500 m high, starting at an altitude of 2350 m. A no-
slip condition was applied at the ground-level lower boundary, characterised by a
roughness length of 10 cm; the lateral boundaries were open. The top boundary was
treated as a rigid lid. Open boundaries defy mass conservation within the domain,
requiring us to limit our simulations to 15 min. The minimum and maximum time
steps were 0.01 and 1 s, respectively, governed by a Courant-Friederichs-Lewy cri-
112
terion limited to 0.8 for advection. Complete model data-fields were output every
minute, and the full velocity, density and SO2 fields were output roughly every sec-
ond. Typical runtime on two computational cores was approximately one week.
Atmospheric initial conditions are specified in the model via a standard single verti-
cal sounding profile of temperature, relative humidity and horizontal winds. These
characterise atmospheric stability and determine potential advective transport. We
used a sounding from Puerto Montt (station SCTE/85799, 41.43° S 73.10°W;∼250 km
south of Villarrica) collected on 8 February 2012, at 12:00 UTC (-4 h time zone),
which sets the simulation start to 08:00:00 LT (local time). Since the crater extends
out of the atmospheric boundary layer and into the free troposphere, we estimate
that the use of an early morning profile is inconsequential to the evolution of the
plume. Only the west-east (u) component of the horizontal wind was retained. This
was aligned with the model’s orientation, and scaled down to 20% of the measured
values, in order to avoid excessively fast advection of the plume (and of turbulence)
out of our small model domain, and strong up- and down-slope flow.
The model output was used to create a set of synthetic calibrated UV images, where
pixel values represent column amounts of SO2 in kg m
−1. The motion estimation
algorithms were then applied to these images, and their performance evaluated by
comparing the calculated SO2 flux across a particular line of pixels to its real value
(calculated using the motion field from the model output). The fluxes were calcu-
lated by integrating the product of the SO2 column amounts across a transect of
the image with the components of the velocity perpendicular to that transect.
Prior to applying the motion estimation algorithms, background and foreground re-
gions (i.e. clear sky and the ground) were masked with Gaussian-distributed white
noise, with a mean value approximating that of the pixels within the plume. This
was found to greatly improve the performance of both algorithms, preventing spu-
rious motion vectors from being computed in regions of very little texture and at
high contrast boundaries.
7.2.4 Villarrica Dataset
A ∼20 min sequence of UV images recorded at Villarrica volcano, Chile, on 8 Febru-
ary 2012 was used to test both algorithms’ performance on “real” data. The im-
ages were captured using an Apogee Alta U47-UV camera with an E2V CCD47-10,
113
1024×1024 pixel CCD and fitted with a 50 mm lens, giving a single pixel “footprint”
of 1.02×1.02 m2 (distance to the plume was ∼4 km). The camera used a filter-wheel
to alternate between in-band and out-of-band imaging, with an in-band image be-
ing recorded approximately every 4 s. Only the in-band images (310 nm) have been
used in this study. We have not attempted to calibrate the images into SO2 column
amounts and the colour-scales on the images presented are in arbitrary units.
As in the case of the ATHAM data, random noise was used to mask background
and foreground regions prior to motion estimation.
7.2.5 Image Registration
To demonstrate how motion estimation algorithms may be used to ensure an accu-
rate mapping between pixels from different cameras (i.e. between the in-band and
out-of-band images), we used a single UV camera image of Villarrica and applied a
known distortion to it using OpenCV’s cvWarpImage function. As shown in Fig. 7.1,
the distortion was made up of a perspective shift (the expected effect if the optical
planes of the two cameras are not co-planar), and a rotation of 0.5° clockwise about
the top left corner of the image. The two motion tracking algorithms were then used
to calculate the shift vectors that mapped features in the distorted image to their
respective positions in the original image.
In order to distribute errors in the shift estimates across the whole of the corrected
image, we used the shift vectors to calculate a single corrective transform which we
then applied to the whole image. In this case we assumed a projective distortion,
and were therefore able to calculate the optimal (least-squares error) correction using
OpenCV’s cvFindHomography function [described in detail by Bradski and Kaehler,
2008]. However, more complex corrections (e.g. including lens distortion) could be
calculated in a similar manner using the correspondences found by the motion track-
ing algorithms.
To avoid errors caused by edge effects, we only included shifts from within a set
region of interest (ROI) when calculating the corrective transformation. The ROI
was selected so as to exclude the edges of the image, the crater region (which has
very high intensity contrasts) and the clear sky regions (which have almost no in-
tensity contrast) since motion estimates in these regions are often inaccurate. This
selection procedure could potentially be automated by thresholding the image based
114
5 pixels
5 pixels
θ
θ = 0.5°
A B C
Figure 7.1: Cartoon of the distortion used to create a test image for the image
registration technique. (A) The original (square) image, (B) perspective shift, (C)
rotation (0.5° clockwise about top left corner).
on texture i.e. excluding regions where variation in pixel value is very small.
7.3 Results
Figure 7.2 shows the results of applying the motion estimation algorithms to the
synthetic images produced using ATHAM. An animation showing the velocity fields
produced for the full time series of images is provided in the supplementary material
to the online version of this article. The SO2 fluxes through a horizontal line of pix-
els 170 m above the crater, calculated using the velocities from ATHAM and the two
motion estimation algorithms, are plotted against time. As shown in Fig. 7.3, the
agreement of the estimated fluxes with the fluxes computed from ATHAM is striking
and provides a strong argument for the use of these algorithms in computing plume
transport speeds. Both the DT-CWT and Farneba¨ck methods show a good corre-
lation across the full range of flux values, with R2 values of 0.95 and 0.97 respectively.
The model used to produce the images for Fig. 7.2 had a closed top (rigid lid) bound-
ary (as described in Section 7.2.3), acting essentially in a similar way as a very strong
temperature inversion layer. When the buoyant plume reached the top boundary of
the domain, it was forced to spread laterally and large down-welling regions quickly
led to a negative flux across the integration line (as can be seen beyond 08:04:31 in
Fig. 7.2). To overcome this issue, an additional 49 vertical grid-points were added,
extending the vertical model domain to a total elevation of more than 6 km on a
stretched grid. In addition, we increased vertical wind shear by applying a linear
multiplier function to the wind profile, varying between 1 and 1.5 over the lower
1.5 km vertical extent of the model domain. Furthermore, we applied ±20% random
white noise onto the initial horizontally homogeneous wind field, and initialised the
115
ATHAM DT-CWT Farneback
08:00:01 08:00:31 08:01:01 08:01:31 08:02:01 08:02:31 08:03:01 08:03:31 08:04:01 08:04:31
Simulation Time
3
2
1
0
1
2
3
4
5
6
S
O
2
 f
lu
x
 (
kg
 s
−1
)
Flux
ATHAM
DT-CWT
Farneback
Figure 7.2: Summary of results from the “closed-top” ATHAM simulation. SO2
fluxes calculated using the velocity fields from ATHAM, the DT-CWT algorithm,
and the Farneba¨ck algorithm are plotted against time. Snapshots of the velocity
fields (corresponding to the time denoted by the vertical grey line) superimposed
over the SO2 amounts are also shown. The line of pixels over which the fluxes were
calculated is highlighted in white.
2 1 0 1 2 3 4 5 6
ATHAM SO2  flux (kg s
−1 )
3
2
1
0
1
2
3
4
5
6
Fa
rn
e
b
a
ck
 S
O
2
 f
lu
x
 (
kg
 s
−1
)
y = 0.95x - 0.07
R2  = 0.968
Farneback
2 1 0 1 2 3 4 5 6
ATHAM SO2  flux (kg s
−1 )
2
1
0
1
2
3
4
5
6
D
T
-C
W
T
 S
O
2
 f
lu
x
 (
kg
 s
−1
)
y = 0.88x - 0.03
R2  = 0.953
DT-CWT
Figure 7.3: Scatter plots showing the relationship between the SO2 fluxes computed
using velocity estimates from the Farneba¨ck (left) and DT-CWT (right) algorithms
and the real fluxes from the “closed-top” ATHAM model. A linear regression line
and the associated R2 value are also shown.
116
ATHAM DT-CWT Farneback
08:03:00 08:04:00 08:05:00 08:06:00 08:07:00
Simulation Time
2
0
2
4
6
8
10
12
S
O
2
 f
lu
x
 (
kg
 s
−1
)
Flux
ATHAM
DT-CWT
Farneback
Figure 7.4: Summary of results from the “open-top” ATHAM simulation. SO2
fluxes calculated using the velocity fields from ATHAM, the DT-CWT algorithm,
and the Farneba¨ck algorithm are plotted against time. Snapshots of the velocity
fields (corresponding to the time denoted by the vertical grey line) superimposed
over the SO2 amounts are also shown. The line of pixels over which the fluxes were
calculated is highlighted in white.
SGS TKE with isotropic turbulence on the order of 1 m2 s−2, to promote resolved
turbulent effects. However, the increased randomisation of the initial conditions did
not produce a significant effect. The motion of the plume produced by this setup
was predominantly horizontal (and often grounded in the vicinity of the crater),
necessitating a vertical integration line. The results from this simulation are shown
in Fig. 7.4. Note however, that the extended vertical extent of the model domain
has been cropped from the plotted data such that the plotted extent matches that of
Fig. 7.2. Once again, the effectiveness of the motion estimation algorithms is clearly
demonstrated (Fig. 7.5). As shown in the animated sequence (see supplementary
materials), the periods in which the motion estimation algorithms significantly un-
derestimate the SO2 flux (e.g. 08:06) correspond to periods when large portions of
the plume are grounded at the integration line and lack any significant structure.
Using the same images as for Figs. 7.2 and 7.4, we made flux estimates using the
conventional cross-correlation technique often used in UV imaging [Kantzas et al.,
117
0 2 4 6 8 10 12
ATHAM SO2  flux (kg s
−1 )
2
0
2
4
6
8
10
12
Fa
rn
e
b
a
ck
 S
O
2
 f
lu
x
 (
kg
 s
−1
)
y = 0.97x - 0.12
R2  = 0.973
Farneback
0 2 4 6 8 10 12
ATHAM SO2  flux (kg s
−1 )
2
0
2
4
6
8
10
12
D
T
-C
W
T
 S
O
2
 f
lu
x
 (
kg
 s
−1
)
y = 1.01x - 0.10
R2  = 0.965
DT-CWT
Figure 7.5: Scatter plots showing the relationship between the SO2 fluxes computed
using velocity estimates from the Farneba¨ck (left) and DT-CWT (right) algorithms
and the real fluxes from the “open-top” ATHAM model. A linear regression line
and the associated R2 value are also shown.
2010]. The integrated SO2 column amounts for the cross-correlation were taken from
20 m either side of the integration lines used in Figs. 7.2 and 7.4. The correlation
window used spanned 20 s. Short time series of the resulting fluxes are shown in
Figs. 7.6 and 7.7. Figure 7.6 shows a period when the cross-correlation method per-
forms well, and Fig. 7.7 shows a period when it does not. Fluxes calculated using
velocity estimates from the Farneba¨ck algorithm are also shown for comparison. It
can be seen that although the cross-correlation method produces reasonable flux es-
timates when the plume consists of discrete puffs (Fig. 7.6), when the plume is less
structured it can produce errors of up to 300 % (Fig. 7.7). The Farneba¨ck algorithm
performs consistently better than the cross-correlation method, and appears to be
less sensitive to plume conditions.
As already discussed, it is not possible to quantitatively assess the performance
of motion estimation algorithms on real UV images. However, some insight can
be gained through a visual inspection of the motion vectors produced by the al-
gorithms. Figure 7.8 shows a snapshot of the motion vectors calculated by both
the DT-CWT and Farneba¨ck algorithms when applied to a sequence of UV images
of Villarrica (see Section 7.2.4). Animations of the motion estimates from the full
sequence of images are provided in the supplementary material to the online version
of this article. As shown by the histograms in Fig. 7.8, both algorithms produce a
similar distribution of velocity vectors over the plume region. In general the DT-
CWT technique produces a smoother velocity field. However, this is possibly due to
our choice of algorithm parameters rather than the algorithm itself and, as discussed
118
08:01:00 08:01:15 08:01:30 08:01:45 08:02:00 08:02:15 08:02:30 08:02:45 08:03:00
Simulation time
1
0
1
2
3
4
5
S
O
2
 f
lu
x
 (
kg
 s
−1
)
ATHAM
Cross-correlation
Farneback
Figure 7.6: Short time series comparing SO2 fluxes calculated from synthetic images
from the “closed-top” ATHAM simulation using a simple cross-correlation method
and the Farneba¨ck algorithm. The true fluxes (calculated using the ATHAM velocity
field) are also shown. This time period was chosen deliberately as an example of
good performance of the cross-correlation technique.
08:04:53 08:05:08 08:05:23 08:05:38 08:05:53 08:06:08 08:06:23 08:06:38 08:06:53
Simulation time
0
2
4
6
8
10
12
S
O
2
 f
lu
x
 (
kg
 s
−1
)
ATHAM
Cross-correlation
Farneback
Figure 7.7: Short time series comparing SO2 fluxes calculated from synthetic images
from the “open-top” ATHAM simulation using a simple cross-correlation method
and the Farneba¨ck algorithm. The true fluxes (calculated using the ATHAM velocity
field) are also shown. This time period was chosen deliberately as an example of
poor performance of the cross-correlation technique.
119
DT-CWT
0 1 2 3 4 5 6 7 8
Speed (m s−1 )
0.00
0.05
0.10
0.15
0.20
0.25
0.30
0.35
N
o
rm
a
lis
e
d
 c
o
u
n
t
Farneback
0 1 2 3 4 5 6 7 8
Speed (m s−1 )
0.00
0.05
0.10
0.15
0.20
0.25
0.30
Figure 7.8: A UV image of Villarrica captured on 8 February 2012, with motion
vectors estimated using the DT-CWT algorithm (left image) and the Farneba¨ck
algorithm (right image). Normalised histograms showing the distribution of mag-
nitudes of the velocity vectors produced by each method are also plotted. Motion
vectors from the outside 10 pixels of the image have been excluded from the his-
tograms, as have the zero velocity vectors.
in Section 7.2.1, a more rigorous methodology for selecting parameter values would
be needed to enable direct comparison of the methods. Both algorithms produce
a mean plume velocity of around 3.5 m s−1, which seems plausible for Villarrica.
Closer inspection of the estimated motion shows, as demonstrated by Fig. 7.9, that
the magnitudes and directions of the velocity vectors from both algorithms match
the motion of easily observable features.
The results of using the DT-CWT and Farneba¨ck algorithms to calculate the cor-
rective transform for a distorted image (as described in Section 7.2.5) are shown in
Fig. 7.10. To assess the performance of the registration we have plotted the percent-
age error of each pixel value in the corrected image with respect to its value in the
original image, i.e. the error associated with the pixel at coordinates (i, j), Ei,j, is
120
A B
DC
Figure 7.9: (A, B) Two successive UV images of Villarrica with a ∼4 s delay be-
tween them, the region displayed in C and D is denoted by the white rectangles.
(C, D) Motion vectors mapping pixels in A to their corresponding positions in B;
computed using the DT-CWT algorithm (C) and the Farneba¨ck algorithm (D). No-
table features in the plume from image A (blue) and B (red) have been enhanced
using Sobel edge detection and are displayed. Correct motion vectors are expected
to point from features shown in blue to the corresponding feature in red.
121
given by:
Ei,j =
|Oi,j − Ci,j|
Oi,j
× 100 (7.1)
where Oi,j is the pixel value in the original image, and Ci,j is its value in the cor-
rected image.
The correction to the image was applied using OpenCV’s cvWarpImage function (as
was the original distortion). It should be noted that distortion and subsequent cor-
rection, even when the transformation is known exactly, leads to errors of 30–40 %
in the out-of-plume region of the corrected image (where the pixel values are close
to zero). This is due to the interpolation required when performing the transforma-
tions. The striped pattern observed in the error plots is caused by a small-amplitude
periodic fluctuation in the intensities of pixels in the original image, and is not an
artifact of the registration technique. We are unclear on the exact origin of this
phenomenon but suggest that it may be caused by an instrumental effect due to the
particular hardware configuration of the camera. As can be seen in Fig. 7.10, both
algorithms produce similar results, with resultant errors being within the interpo-
lation error of the cvWarpImage function. In the plume region, the errors in the
corrected images are <3 % for both algorithms.
7.4 Discussion
In general both algorithms performed similarly, with equal resultant errors when
used for image registration, and comparable correlations with real flux values when
applied to our synthetic data (Figs. 7.3 and 7.5). However, we found the Farneba¨ck
algorithm to be less sensitive to input parameters, and therefore a little easier to
use than the DT-CWT algorithm. As previously mentioned, a wide range of motion
estimation algorithms exist and it is not our intent to promote the DT-CWT or
Farneba¨ck algorithms in particular as being superior for use with UV camera im-
ages. That said, the results produced by both when applied to the Villarrica and
ATHAM datasets were encouraging, and suggest that, even if not optimal, either
algorithm would be a reasonable choice for this application.
More specialised optical-flow algorithms exist, which are specifically designed for
tracking fluid flows. As such, they incorporate constraints from fluid mechanics,
such as continuity, to prevent un-physical flow estimations. The findings of Corpetti
122
Raw image showing ROI Error in uncorrected image ROI
0
10
20
30
40
50
60
70
80
90
E
rr
o
r 
(%
)
Error in ROI corrected
using DT-CWT
0
5
10
15
20
25
30
35
40
E
rr
o
r 
(%
)
Error in ROI corrected
using Farnebäck
0
5
10
15
20
25
30
35
40
E
rr
o
r 
(%
)
Figure 7.10: Results of using the DT-CWT and Farneba¨ck algorithms to estimate
a projective distortion to an image. Top left: the original, undistorted image with
the region of interest (ROI) over which motion vectors were estimated highlighted.
Top right: percentage pixel errors in the ROI of the distorted image. Lower left:
percentage pixel errors in the ROI of the image corrected using the DT-CWT algo-
rithm. Lower right: percentage pixel errors in the ROI of the image corrected using
the Farneba¨ck algorithm.
123
et al. [2002] show that inclusion of these constraints can lead to improved accuracy
of the motion estimates produced for fluid flows. Such algorithms may be more
suitable for estimating the plume transport speed than the DT-CWT or Farneba¨ck
algorithms. However, successful application of the continuity constraint relies on
accurate knowledge of the SO2 amounts in each pixel. Whilst this is the case for
the synthetic images produced by ATHAM, the calibration of real images is more
questionable and it is not clear how these errors might propagate into the motion
estimate. In the future it may be possible to combine the motion estimates and the
apparent absorbance data into a single coherent plume model which could then be
inverted for SO2 flux, minimising the errors from both the motion and the calibra-
tion. However, this is pure speculation, and is beyond the scope of this investigation.
Synthetic plume images, such as those created using ATHAM, are clearly a useful
validation tool for SO2 flux calculation techniques. However, the extent to which
they are representative of images of a real plume is questionable. In short, it is not
necessarily true that a motion estimation algorithm that can perfectly reproduce
the velocity field of the synthetic images will work well when applied to real data.
However, it seems fair to say that the synthetic images provide a useful starting
point for algorithm development and tuning.
The configuration and setup of the model used in this study was limited by comput-
ing time, and it is likely that with more resources its capabilities could be exploited
to a greater extent. In particular, our simulations fail to reproduce the very rapid
dilution of the plume that is seen in the Villarrica images. This may be due in part
to the fact that the model is run only in two-dimensions, greatly reducing the scope
for entrainment of the surrounding atmosphere into the plume. A three-dimensional
setup using ATHAM, whilst possible, would be very computationally expensive. It
would, however, facilitate very useful studies into the effects of treating a real vol-
canic plume as a two-dimensional object when calculating transport speeds. It is
important to remember that each pixel in a UV image represents an integrated
amount of SO2 in the plume along the optical path of that pixel. Hence, estimating
motion vectors from such an image can only ever yield a coarse approximation to the
true plume velocity. Three-dimensional plume model data would allow the errors
introduced by estimating transport speeds from a two-dimensional projection of the
plume to be investigated.
Although both algorithms were able to accurately compute the transformation re-
124
quired to correct our distorted image, in reality the problem is more complex. The
image pairs to be registered from UV camera systems employing two camera units
are captured at different wavebands and features visible in the in-band image may
not be present in the out-of-band image. This will likely degrade the performance
of the algorithms and may compromise their ability to compute the motion vectors
altogether. The extent of this problem will depend on the composition of the specific
plume being imaged, and could be alleviated by careful selection of regions of the
images over which to perform the motion estimation. An alternative, and possibly
superior approach, would be to periodically capture a set of images at the same
waveband from both cameras and use these to compute the registration parameters
required.
It is theoretically possible to use the motion estimation techniques to correct for
more complex image distortions than the projective transformation which we demon-
strated (e.g. radial lens distortions). However, allowing more degrees of freedom
when calculating the corrective transformation will lead to a less robust result,
which is more susceptible to errors in the motion estimates. Since the lens distor-
tion should be relatively constant, we believe that this correction would be better
calculated prior to field deployment of the camera, perhaps using the well-established
“chess board” technique [described in detail by Bradski and Kaehler, 2008].
7.5 Conclusions
We have presented two different motion estimation algorithms, a DT-CWT based
algorithm and the Farneba¨ck algorithm, and showed how they can be applied to
UV camera images to improve the calculation of SO2 flux. The ATHAM model was
used to create a set of synthetic UV images for which the SO2 amounts and velocity
field were known, and these were used to validate the estimates produced by the
two algorithms.
Motion estimation algorithms are clearly a useful tool when applied to SO2 mea-
surements using UV cameras. It has been found in this study that they can provide
a means to accurately register images from two separate camera units, in addition
to providing more robust estimates of plume velocity for use in flux calculations
than the currently-used cross-correlation method. The two algorithms presented
here were similar in performance and would both be suitable for UV camera appli-
125
cations, however, we found the Farneba¨ck algorithm a little easier to use. A more
thorough comparison of different algorithms is needed to determine the optimal al-
gorithm for UV camera applications and this will require a more rigorous approach
to input parameter selection than was used in this study.
Using plume models, such as ATHAM, to produce test datasets for motion estima-
tion algorithms provides a simple way to evaluate and compare their performance.
Further research is required to ensure that the data produced are an accurate rep-
resentation of real volcanic plumes, and extension of the model setup into three-
dimensions promises to provide valuable insights into the errors associated with flux
measurements.
126
Chapter 8
Conclusions
“Well if everything else wasn’t broken, we’d be sorted.”
– Bill McIntosh, Erebus, December 2010.
I have developed new volcanic monitoring technologies and techniques and also anal-
ysed the results of applying these to the Erebus lava lake. A new thermal camera
system was created and deployed at the crater rim of Erebus and new hardware and
software for recording SO2 fluxes was also developed. Additionally, a new method-
ology for estimating plume transport velocity from UV camera images was investi-
gated. Overall, IR images of the Erebus lava lake from 2004 up to 2013 have been
analysed to obtain the mean speed of the lake surface. The most recent IR images
were collected at concurrently with measurements of SO2 flux, gas composition and
surface elevation of the lake allowing, for the first time, a detailed inter-comparison
of the cyclic behaviour observed in each. The conclusions from these studies are
summarised below, organised in terms of the research objectives which were laid out
in Section 1.7.
8.1 Improved Monitoring
The new thermal camera system developed as part of this study and deployed on
Erebus in December 2012 (Chapter 2) is a great success. Despite the failure of the
power system during the winter months, the dataset collected surpasses any previ-
ously collected, both in terms of length and continuity. The camera has not only
greatly enhanced the monitoring capacity of the Mount Erebus Volcano Observa-
tory, but has also provided a valuable pilot project to test the feasibility of running
instruments at the crater rim outside of the short time window of the annual field
campaigns. Future instruments are already under consideration and plans are cur-
rently being made to install year-round telemetry to the crater rim during the 2014
field season, allowing “live” imagery of the lake to be collected. As well as the camera
system itself, the infrastructure that was installed to support it (power, Ethernet)
has been hugely beneficial to the field campaigns in general. Data loss due to power
and network failures is now largely a thing of the past, and time-series data contin-
ues to improve in terms of continuity, facilitating ever more detailed analyses.
High temporal resolution measurements of SO2 flux from volcanoes are challenging
to achieve for a variety of reasons. A major contribution to the problem is the diffi-
culty of accurately measuring of plume transport velocity. I have demonstrated that
using the motion tracking methodologies used during the IR image analysis part of
this study, more accurate estimates of plume transport velocity can be obtained
than is possible with the currently-used cross-correlation method. These findings
will not only be useful for future studies at Erebus, but also for researchers and
128
volcano observatories around the world who use UV cameras for making SO2 flux
measurements.
8.2 Assessment of Annual Variability of the Lava
Lake Cycles
Pulsatory behaviour in the mean surface speed of the lake has been sustained since
at least 2004, when the first measurements were made. Cycles, with a period of
5–18 min, are always present although with varying degrees of visibility in the data.
Despite significant changes in the visible size of the lake, no change in the cyclic be-
haviour was detectable across the period studied (2004–2013). This would suggest
not only that the cyclic behaviour is a long-term stable feature, but also that the
mechanism responsible occurs at a greater depth than the lake itself.
Large (decametric) bubbles are too infrequent to be considered as a possible driving
mechanism for the pulsatory behaviour and the timings of smaller bubbles are found
to be uncorrelated with the phase of the cycles in surface speed. These findings point
towards a driving mechanism based on the flow dynamics of magma exchange within
the shallow magmatic system and largely rule out gas-driven mechanisms.
8.3 Multi-parameter Observations of Cyclic Be-
haviour
In addition to the cyclic behaviour in the mean surface speed, similar behaviour is
found in the SO2 flux, plume composition and surface elevation of the lake. The
cycles in each of these parameters are well correlated with each other, suggesting
that they all derive from a common underlying process. It can be shown by con-
sidering the magma flux required to explain the mass of H2O released per cycle,
that the observed variations in surface elevation of the lake can be accounted for by
exchange of magma in and out of the lake. This supports the hypothesis that the
pulsatory behaviour of the Erebus lake is driven by magma exchange rather than gas
flux. Furthermore, it suggests that the lake is highly permeable to gas, and does not
present a significant barrier to its escape. However, this is somewhat contradicted
by the finding that the peak in surface motion and gas composition lag the peak
129
in surface elevation by 1–3 min. This may be explained by the time taken for fresh
magma injected into the lake to rise to the surface, suggesting that the gas rises
with the melt through the lake.
The behaviour of the lake following explosive events caused by large bubble bursts
is indistinguishable from its behaviour prior to the events, with the exception of
the ∼5 min required for the lake to refill. No change in amplitude or period of the
pulsatory behaviour is observed, suggesting a high degree of decoupling between the
explosive and passive degassing mechanisms of the system.
8.4 The Future
This study has done much to progress our understanding of the Erebus lava lake;
identifying cycles in the surface elevation, surface motion, SO2 flux and gas chem-
istry as a persistent and stable feature of its behaviour, which is decoupled from
the explosive activity. However, there are clearly a lot of questions left unanswered.
In particular, the exact mechanism responsible for the pulsatory behaviour remains
elusive. Whilst the observations presented here have provided some insight into
feasible mechanisms, it is simply not possible to test these hypotheses further based
on observations alone. It is our belief that the next breakthroughs in understanding
the Erebus magmatic system will come from combining the observations presented
here with lab-based flow models and computer simulations. This will allow the exact
conditions required to reproduce the observed behaviour to be deduced. The work
presented here has brought us several steps closer to understanding Erebus, but it
is a long road and the end is not yet in sight.
130
Acknowledgements
None of this would have been possible without the help and support of a huge
number of people. First and foremost, I thank my supervisor Clive Oppenheimer.
This PhD has played out almost exactly as I had hoped, and I owe that to you.
Aaron Curtis, thanks for getting me into this mess in the first place and I really
hope we manage to continue to find excuses to go into the field together. Phil Kyle,
I’m still suspicious that the only reason you ever took me to Erebus was because
you were embarrassed about nearly burning Clive’s house down. But whatever the
reason, thank you! Vitchko Tsanev, Talfan Barnie, Amy Donovan, Lex Hoffmann
and Nick Kingsbury, it has been a great pleasure to work with you all. I would
also like to thank my parents for their continued support and encouragement. A
huge and deeply heartfelt thanks to Andrew and Bela, my wonderful house-mates,
and to all the rest of the inhabitants and associates of the Long Road commune.
Cycle-cake Sundays, box-of-wine Wednesdays, slugs and electric taps; it really is a
“special” place to live and I hope we can keep it going for a bit longer at least. Last
but not least, I must thank The Volcanofiles (Kelby Hicks, Tehnuka Ilanko, Kayla
Iacovino, Yves Moussallam). It has certainly been a rough road guys, but the good
times were some of the best in my life, and it is those that I will remember. Kelby,
your passion and enthusiasm was truly inspirational; we all miss you dearly.
The Volcanofiles at Villarrica, Chile, 2012. Photo: Kayla Iacovino.
131
Appendix A
Data Availability
This appendix contains charts showing when thermal images of the Erebus lava lake
have been recorded. Images from both the P25 and the newer SC645 camera (see
Table 4.1) are included. It will be seen that the data are far from continuous, even
with the new camera system [Chapter 2, Peters et al., 2014c]. However, many of the
gaps in the the data are caused by a build-up of rime-ice on the lens of the camera,
obscuring its view of the lake entirely.
It should be noted that the time-spans of the data included here have not been
“cleaned” of unusable images (i.e. those in which the lake is not clearly visible).
Therefore, even if images are available for a particular time period, it may not be
possible to obtain useful results from them.
2
0
0
6
–
2
0
0
7
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
D
ec
2
0
0
7
–
2
0
0
8
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
D
ec
2
0
0
8
–
2
0
0
9
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
D
ec
2
0
0
9
–
2
0
1
0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
D
ec
2
0
1
0
–
2
0
1
1
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
D
ec
2
0
1
1
–
2
0
1
2
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
D
ec
F
ig
u
re
A
.1
:
D
at
a
av
ai
la
b
il
it
y
fo
r
th
e
P
25
th
er
m
al
ca
m
er
a.
2
0
1
0
–
2
0
1
1
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
D
ec
2
0
1
1
–
2
0
1
2
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
D
ec
2
0
1
2
–
2
0
1
3
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
N
ov
D
ec
J
a
n
F
eb
M
ar
A
p
r
F
ig
u
re
A
.2
:
D
at
a
av
ai
la
b
il
it
y
fo
r
th
e
S
C
64
5
th
er
m
al
ca
m
er
a.
134
Bibliography
A. Aiuppa, C. Federico, G. Giudice, and S. Gurrieri. Chemical mapping of a fu-
marolic field: La Fossa Crater, Vulcano Island (Aeolian Islands, Italy). Geophysi-
cal Research Letters, 32(13), 2005. ISSN 1944-8007. doi: 10.1029/2005GL023207.
Apache. Apache OpenOffice Calc, 2013. URL http://www.openoffice.org/
product/calc.html.
R. Aster, S. Mah, P. Kyle, W. McIntosh, N. Dunbar, J. Johnson, M. Ruiz, and
S. McNamara. Very long period oscillations of Mount Erebus volcano. Journal of
Geophysical Research, 108:22 PP., Nov. 2003. doi: 200310.1029/2002JB002101.
R. Aster, W. McIntosh, P. Kyle, R. Esser, B. A. Bartel, N. Dunbar, B. Johns, J. B.
Johnson, R. Karstens, C. Kurnik, M. McGowan, S. McNamara, C. Meertens,
B. Pauley, M. Richmond, and M. Ruiz. Real-time data received from Mount
Erebus volcano, Antarctica. Eos, 85(10):97–104, Mar. 2004. doi: 200410.1029/
2004EO100001.
S. J.-B. Bauguitte, N. Brough, M. M. Frey, A. E. Jones, D. J. Maxfield, H. K.
Roscoe, M. C. Rose, and E. W. Wolff. A network of autonomous surface
ozone monitors in Antarctica: technical description and first results. Atmo-
spheric Measurement Techniques, 4(4):645–658, Apr. 2011. ISSN 1867-8548. doi:
10.5194/amt-4-645-2011.
G. Bluth, J. Shannon, I. Watson, A. Prata, and V. Realmuto. Development of an
ultra-violet digital camera for volcanic SO2 imaging. Journal of Volcanology and
Geothermal Research, 161(1–2):47–56, Mar. 2007. ISSN 0377-0273. doi: 10.1016/
j.jvolgeores.2006.11.004.
M. L. Boas. Error Propagation Equation. In Mathematical Methods in the Physical
Sciences, pages 734–735. Wiley, 2nd edition, Apr. 1983a. ISBN 0471044091.
M. L. Boas. Mathematical Methods in the Physical Sciences. Wiley, 2nd edition,
Apr. 1983b. ISBN 0471044091.
135
M. Boichu, C. Oppenheimer, V. Tsanev, and P. R. Kyle. High temporal resolution
SO2 flux measurements at Erebus volcano, Antarctica. Journal of Volcanology and
Geothermal Research, 190(3-4):325–336, 2010. ISSN 0377-0273. doi: 10.1016/j.
jvolgeores.2009.11.020.
E. Bouche, S. Vergniolle, T. Staudacher, A. Nercessian, J.-C. Delmont, M. Frogneux,
F. Cartault, and A. Le Pichon. The role of large bubbles detected from acoustic
measurements on the dynamics of Erta ’Ale lava lake (Ethiopia). Earth and
Planetary Science Letters, 295(1–2):37–48, June 2010. ISSN 0012-821X. doi:
10.1016/j.epsl.2010.03.020.
G. Bradski. The OpenCV Library. Dr. Dobb’s Journal of Software Tools, 25(11):
120, 122–125, Nov. 2000. ISSN 1044-789X.
G. Bradski and A. Kaehler. Learning OpenCV: Computer Vision with the OpenCV
Library. O’Reilly Media, 1st edition, Sept. 2008. ISBN 0596516134.
A. Burgisser, C. Oppenheimer, M. Alletti, P. R. Kyle, B. Scaillet, and M. R.
Carroll. Backward tracking of gas chemistry measurements at Erebus volcano.
Geochemistry, Geophysics, Geosystems, 13(11), 2012. ISSN 1525-2027. doi:
10.1029/2012GC004243.
J. Burrows, A. Richter, A. Dehn, B. Deters, S. Himmelmann, S. Voigt, and J. Or-
phal. Atmospheric Remote-Sensing Reference Data from GOME: 2. Temperature-
Dependent Absorption Cross Sections of O3 in the 231-794 nm Range. Journal
of Quantitative Spectroscopy and Radiative Transfer, 61(4):509–517, Mar. 1999.
ISSN 0022-4073. doi: 10.1016/S0022-4073(98)00037-5.
M. R. Burton, C. Oppenheimer, L. A. Horrocks, and P. W. Francis. Remote sensing
of CO2 and H2O emission rates from Masaya volcano, Nicaragua. Geology, 28
(10):915, 2000. ISSN 0091-7613. doi: 10.1130/0091-7613(2000)28〈915:RSOCAH〉
2.0.CO;2.
J. Calkins, C. Oppenheimer, and P. Kyle. Ground-based thermal imaging of lava
lakes at Erebus volcano, Antarctica. Journal of Volcanology and Geothermal Re-
search, 177(3):695–704, 2008. ISSN 0377-0273. doi: DOI:10.1016/j.jvolgeores.
2008.02.002. Volcanology of Erebus volcano, Antarctica.
R. Campion, M. Martinez-Cruz, T. Lecocq, C. Caudron, J. Pacheco, G. Pinardi,
C. Hermans, S. Carn, and A. Bernard. Space- and ground-based measurements
of sulphur dioxide emissions from Turrialba Volcano (Costa Rica). Bulletin of
136
Volcanology, 74(7):1757–1770, June 2012. ISSN 0258-8900, 1432-0819. doi: 10.
1007/s00445-012-0631-z.
M. Cassidy, P. D. Cole, K. Hicks, N. Varley, N. Peters, and A. Lerner. Rapid and
slow: Varying magma ascent rates provide the mechanism for small Vulcanian
eruptions. Earth and Planetary Science Letters (in review), 2014.
K. Chance and R. Kurucz. An improved high-resolution solar reference spectrum
for earth’s atmosphere measurements in the ultraviolet, visible, and near infrared.
Journal of Quantitative Spectroscopy and Radiative Transfer, 111(9):1289–1295,
June 2010. ISSN 0022-4073. doi: 10.1016/j.jqsrt.2010.01.036.
T. Corpetti, E. Memin, and P. Perez. Dense estimation of fluid flows. IEEE Trans-
actions on Pattern Analysis and Machine Intelligence, 24(3):365–380, Mar. 2002.
ISSN 0162-8828. doi: 10.1109/34.990137.
M. P. Dalton, I. M. Watson, P. A. Nadeau, C. Werner, W. Morrow, and J. M.
Shannon. Assessment of the UV camera sulfur dioxide retrieval for point source
plumes. Journal of Volcanology and Geothermal Research, 188(4):358–366, Dec.
2009. ISSN 0377-0273. doi: 10.1016/j.jvolgeores.2009.09.013.
A. G. Davies, J. Calkins, L. Scharenbroich, R. G. Vaughan, R. Wright, P. Kyle,
R. Castano, S. Chien, and D. Tran. Multi-instrument remote and in situ ob-
servations of the Erebus Volcano (Antarctica) lava lake in 2005: A compari-
son with the Pele lava lake on the jovian moon Io. Journal of Volcanology
and Geothermal Research, 177(3):705–724, Nov. 2008. ISSN 0377-0273. doi:
16/j.jvolgeores.2008.02.010.
E. De Lauro, S. De Martino, M. Falanga, and M. Palo. Modelling the macroscopic
behavior of Strombolian explosions at Erebus volcano. Physics of the Earth and
Planetary Interiors, pages 174–186, 2009. ISSN 0031-9201. doi: 10.1016/j.pepi.
2009.05.003.
R. Dibble, P. Kyle, and C. Rowe. Video and seismic observations of Strombolian
eruptions at Erebus volcano, Antarctica. Journal of Volcanology and Geothermal
Research, 177(3):619–634, Nov. 2008. ISSN 0377-0273. doi: 10.1016/j.jvolgeores.
2008.07.020.
T. Divoux, E. Bertin, V. Vidal, and J. C. Ge´minard. Intermittent outgassing through
a non-Newtonian fluid. Physical Review E, 79(5):56204, 2009.
R. Dunn. wxpython: wxwindows for python. http://www.wxpython.org. 2003.
137
M. Edmonds, R. A. Herd, B. Galle, and C. M. Oppenheimer. Automated, high
time-resolution measurements of SO2 flux at Soufrie`re Hills Volcano, Montserrat.
Bulletin of Volcanology, 65(8):578–586, Nov. 2003. ISSN 0258-8900, 1432-0819.
doi: 10.1007/s00445-003-0286-x.
G. Farneba¨ck. Two-Frame Motion Estimation Based on Polynomial Expansion.
In J. Bigun and T. Gustavsson, editors, Image Analysis, number 2749 in Lecture
Notes in Computer Science, pages 363–370. Springer Berlin Heidelberg, Jan. 2003.
ISBN 978-3-540-40601-3, 978-3-540-45103-7. URL http://link.springer.com/
chapter/10.1007/3-540-45103-X_50.
P. Francis, C. Oppenheimer, and D. Stevenson. Endogenous growth of persistently
active volcanoes. Nature, 366(6455):554–557, Dec. 1993. doi: 10.1038/366554a0.
B. Galle, C. Oppenheimer, A. Geyer, A. J. McGonigle, M. Edmonds, and L. Hor-
rocks. A miniaturised ultraviolet spectrometer for remote sensing of SO2 fluxes:
a new tool for volcano surveillance. Journal of Volcanology and Geother-
mal Research, 119(1–4):241–254, Jan. 2003. ISSN 0377-0273. doi: 10.1016/
S0377-0273(02)00356-6.
A. Geiger, P. Lenz, C. Stiller, and R. Urtasun. The KITTI Vision Benchmark
Suite, 2014. URL http://www.cvlibs.net/datasets/kitti/eval_stereo_
flow.php?benchmark=flow.
A. Gerst, M. Hort, R. C. Aster, J. B. Johnson, and P. R. Kyle. The first second
of volcanic eruptions from the Erebus volcano lava lake, Antarctica - energies,
pressures, seismology, and infrasound. Journal of Geophysical Research: Solid
Earth, 118(7):3318–3340, 2013. ISSN 2169-9356. doi: 10.1002/jgrb.50234.
W. F. Giggenbach, P. R. Kyle, and G. L. Lyon. Present volcanic activity on Mount
Erebus, Ross Island, Antarctica. Geology, 1(3):135–136, Nov. 1973.
Gnome. GNOME Office / Gnumeric - Welcome to Gnumeric!, 2013. URL https:
//projects.gnome.org/gnumeric/.
R. Hegger and H. Kantz. Improved false nearest neighbor method to detect deter-
minism in time series data. Physical Review E, 60(4):4970–4973, Oct. 1999. doi:
10.1103/PhysRevE.60.4970.
R. Hegger, H. Kantz, and T. Schreiber. Practical implementation of nonlinear time
series methods: The TISEAN package. arXiv e-print chao-dyn/9810005, Sept.
1998. URL http://arxiv.org/abs/chao-dyn/9810005. CHAOS 9 (1999) 413.
138
M. Herzog and H.-F. Graf. Applying the three-dimensional model ATHAM to vol-
canic plumes: Dynamic of large co-ignimbrite eruptions and associated injection
heights for volcanic gases. Geophysical Research Letters, 37(19):L19807, Oct. 2010.
ISSN 1944-8007. doi: 10.1029/2010GL044986.
M. Herzog, H.-F. Graf, C. Textor, and J. M. Oberhuber. The effect of phase changes
of water on the development of volcanic plumes. Journal of Volcanology and
Geothermal Research, 87(1–4):55–74, Dec. 1998. ISSN 0377-0273. doi: 10.1016/
S0377-0273(98)00100-0.
M. Herzog, J. M. Oberhuber, and H.-F. Graf. A Prognostic Turbulence Scheme
for the Nonhydrostatic Plume Model ATHAM. Journal of the Atmospheric Sci-
ences, 60(22):2783–2796, Nov. 2003. ISSN 0022-4928, 1520-0469. doi: 10.1175/
1520-0469(2003)060〈2783:APTSFT〉2.0.CO;2.
A. P. Holland, I. M. Watson, J. C. Phillips, L. Caricchi, and M. P. Dalton. De-
gassing processes during lava dome growth: Insights from Santiaguito lava dome,
Guatemala. Journal of Volcanology and Geothermal Research, 202(1–2):153–166,
Apr. 2011. ISSN 0377-0273. doi: 10.1016/j.jvolgeores.2011.02.004.
L. A. Horrocks, C. Oppenheimer, M. R. Burton, H. J. Duffell, N. M. Davies, N. A.
Martin, and W. Bell. Open-path Fourier transform infrared spectroscopy of
SO2: An empirical error budget analysis, with implications for volcano moni-
toring. Journal of Geophysical Research: Atmospheres, 106(D21):27647–27659,
Nov. 2001. ISSN 2156-2202. doi: 10.1029/2001JD000343.
J. D. Hunter. Matplotlib: A 2D graphics environment. Computing In Science &
Engineering, 9(3):90–95, 2007.
H. E. Huppert and M. A. Hallworth. Bi-directional flows in constrained sys-
tems. Journal of Fluid Mechanics, 578(1):95–112, 2007. doi: 10.1017/
S0022112007004661.
K. Iacovino. An unexpected journey: Experimental insights into magma and volatile
transport beneath Erebus volcano, Antarctica. Ph.d. thesis, University of Cam-
bridge, 2014.
K. Iacovino, N. Peters, and C. Oppenheimer. Toward a unified method for the
quantification of volatiles in magmas via FTIR. Florence, Italy, 2013. doi: 10.
1180/minmag.2013.077.5.9. URL www.minersoc.org.
139
K. Iacovino, C. Oppenheimer, B. Scaillet, and P. Kyle. Erebus volcano as an
archetype for CO2-dominated rift volcanism: Experimental constraints on the
magma plumbing of Ross Island, Antarctic. Journal of Petrology (in review),
2014.
T. Ilanko, C. Oppenheimer, A. Burgisser, and P. Kyle. Cyclic degassing of Erebus
volcano, Antarctica. Bulletin of Volcanology (in review), 2014.
M. R. James, H. Pinkerton, and S. Robson. Image-based measurement of flux varia-
tion in distal regions of active lava flows. Geochemistry, Geophysics, Geosystems,
8(3):Q03006, Mar. 2007. ISSN 1525-2027. doi: 10.1029/2006GC001448.
E. Jones, T. Oliphant, and P. Peterson. SciPy: open source scientific tools for
Python, 2001. URL http://www.scipy.org/.
K. R. Jones, J. B. Johnson, R. Aster, P. R. Kyle, and W. McIntosh. Infrasonic track-
ing of large bubble bursts and ash venting at Erebus volcano, Antarctica. Journal
of Volcanology and Geothermal Research, 177(3):661–672, Nov. 2008. ISSN 0377-
0273. doi: 10.1016/j.jvolgeores.2008.02.001.
L. Jones. Terrestrial Laser Scanning (TLS) Observations of Erebus Volcano, Antarc-
tica: Insights into the Near-Surface Magmatic System. MSc thesis, New Mex-
ico Institute of Mining and Technology, 2013. URL http://www.ees.nmt.edu/
outside/alumni/papers/2013t_jones_lk.pdf.
L. K. Jones, P. R. KYLE, J. Frechette, M. H. Okal, and C. Oppenheimer. Terrestrial
Laser Scanning observations of geomorphic changes and lava lake levels at Erebus
volcano, Antarctica. Journal of Volcanology and Geothermal Research (in review),
2014.
H. Kantz and T. Schreiber. Nonlinear Time Series Analysis. Cambridge University
Press, 2nd edition, Nov. 2003. ISBN 0521529026.
E. P. Kantzas, A. McGonigle, G. Tamburello, A. Aiuppa, and R. G. Bryant.
Protocols for UV camera volcanic SO2 measurements. Journal of Volcanology
and Geothermal Research, 194(1–3):55–60, July 2010. ISSN 0377-0273. doi:
10.1016/j.jvolgeores.2010.05.003.
P. J. Kelly, P. R. Kyle, N. W. Dunbar, and K. W. Sims. Geochemistry and miner-
alogy of the phonolite lava lake, Erebus volcano, Antarctica: 1972-2004 and com-
parison with older lavas. Journal of Volcanology and Geothermal Research, 177
(3):589–605, Nov. 2008. ISSN 0377-0273. doi: 10.1016/j.jvolgeores.2007.11.025.
140
M. B. Kennel, R. Brown, and H. D. I. Abarbanel. Determining embedding dimension
for phase-space reconstruction using a geometrical construction. Physical Review
A, 45(6):3403–3411, Mar. 1992. doi: 10.1103/PhysRevA.45.3403.
C. Kern, F. Kick, P. Lu¨bcke, L. Vogel, M. Wo¨hrbach, and U. Platt. Theoretical
description of functionality, applications, and limitations of SO2 cameras for the
remote sensing of volcanic plumes. Atmospheric Measurement Techniques Discus-
sions, 3(1):531–578, Feb. 2010. ISSN 1867-8610. doi: 10.5194/amtd-3-531-2010.
C. Kern, C. A. Werner, M. P. Doukas, T. Elias, P. J. Kelly, and A. J. Sutton.
Imaging volcanic SO2 plumes with UV cameras. In EGU General Assembly
Conference Abstracts, volume 14, page 12596, Vienna, Apr. 2012. EGU. URL
http://adsabs.harvard.edu/abs/2012EGUGA..1412596K.
C. Kern, C. Werner, T. Elias, A. J. Sutton, and P. Lu¨bcke. Applying UV cameras for
SO2 detection to distant or optically thick volcanic plumes. Journal of Volcanology
and Geothermal Research, 262:80–89, July 2013. ISSN 0377-0273. doi: 10.1016/
j.jvolgeores.2013.06.009.
N. Kingsbury. Complex wavelets for shift invariant analysis and filtering of signals.
Applied and Computational Harmonic Analysis, 10(3):234–253, 2001. ISSN 1063-
5203. doi: DOI:10.1006/acha.2000.0343.
H. Knox. Eruptive characteristics and glacial earthquake investigation on Erebus
volcano, Antarctica. Ph.d. thesis, New Mexico Institute of Mining and Technology,
2012. URL http://www.ees.nmt.edu/outside/alumni/thesis.php.
S. G. Kraus. DOASIS, A Framework Design for DOAS. Ph.d. thesis, University of
Mannheim, 2006.
P. R. Kyle, K. Meeker, and D. Finnegan. Emission rates of sulfur dioxide, trace
gases and metals from Mount Erebus, Antarctica. Geophysical Research Letters,
17(12):PP. 2125–2128, Nov. 1990. doi: 199010.1029/GL017i012p02125.
A. La Spina, G. Salerno, and M. Burton. Recent SO2 camera and OP-FTIR field
measurements in Mexico and Guatemala. In EGU General Assembly Conference
Abstracts, volume 15, page 10478, Apr. 2013. URL http://adsabs.harvard.
edu/abs/2013EGUGA..1510478L.
J. Lawrence, M. Ashley, and J. Storey. A remote, autonomous laboratory for Antarc-
tica with hybrid power generation. Aust. Journal Elec. Electronics Engineering,
2(1):1–12, 2004.
141
E. Lev, M. Spiegelman, R. J. Wysocki, and J. A. Karson. Investigating lava flow
rheology using video analysis and numerical flow models. Journal of Volcanology
and Geothermal Research, 247–248:62–73, Dec. 2012. ISSN 0377-0273. doi: 10.
1016/j.jvolgeores.2012.08.002.
P. Lu¨bcke, N. Bobrowski, S. Illing, C. Kern, J. M. Alvarez Nieves, L. Vogel, J. Ziel-
cke, H. Delgado Granados, and U. Platt. On the absolute calibration of SO2
cameras. Atmos. Meas. Tech., 6(3):677–696, Mar. 2013. ISSN 1867-8548. doi:
10.5194/amt-6-677-2013.
J. Magarey and N. Kingsbury. Motion estimation using a complex-valued wavelet
transform. Signal Processing, IEEE Transactions on, 46(4):1069–1084, Apr. 1998.
ISSN 1053-587X. doi: 10.1109/78.668557.
T. Mather, A. Allen, B. Davison, D. Pyle, C. Oppenheimer, and A. McGonigle.
Nitric acid from volcanoes. Earth and Planetary Science Letters, 218(1–2):17–30,
Jan. 2004a. ISSN 0012-821X. doi: 10.1016/S0012-821X(03)00640-X.
T. A. Mather, V. I. Tsanev, D. M. Pyle, A. J. S. McGonigle, C. Oppenheimer, and
A. G. Allen. Characterization and evolution of tropospheric plumes from Lascar
and Villarrica volcanoes, Chile. Journal of Geophysical Research: Atmospheres,
109(D21):n/a–n/a, 2004b. ISSN 2156-2202. doi: 10.1029/2004JD004934.
A. J. McGonigle. Measurement of volcanic SO2 fluxes with differential optical ab-
sorption spectroscopy. Journal of Volcanology and Geothermal Research, 162(3–4):
111–122, May 2007. ISSN 0377-0273. doi: 10.1016/j.jvolgeores.2007.02.001.
A. J. S. McGonigle. Walking traverse and scanning DOAS measurements of volcanic
gas emission rates. Geophysical Research Letters, 29(20), 2002. ISSN 0094-8276.
doi: 10.1029/2002GL015827.
A. J. S. McGonigle, S. Inguaggiato, A. Aiuppa, A. R. Hayes, and C. Oppen-
heimer. Accurate measurement of volcanic SO2 flux: Determination of plume
transport speed and integrated SO2 concentration with a single device. Geo-
chemistry, Geophysics, Geosystems, 6(2):n/a–n/a, 2005. ISSN 1525-2027. doi:
10.1029/2004GC000845.
A. J. S. McGonigle, A. Aiuppa, M. Ripepe, E. P. Kantzas, and G. Tamburello.
Spectroscopic capture of 1 Hz volcanic SO2 fluxes and integration with volcano
geophysical data. Geophysical Research Letters, 36(21):n/a–n/a, 2009. ISSN 1944-
8007. doi: 10.1029/2009GL040494.
142
I. Molina, A. Burgisser, and C. Oppenheimer. Numerical simulations of convection
in crystal-bearing magmas: A case study of the magmatic system at Erebus,
Antarctica. Journal of Geophysical Research: Solid Earth, 117(B7), 2012. ISSN
2156-2202. doi: 10.1029/2011JB008760.
T. Mori and M. Burton. The SO2 camera: A simple, fast and cheap method for
ground-based imaging of SO2 in volcanic plumes. Geophysical Research Letters,
33(24):L24804, Dec. 2006. ISSN 0094-8276. doi: 10.1029/2006GL027916.
T. Mori and M. Burton. Quantification of the gas mass emitted during single
explosions on Stromboli with the SO2 imaging camera. Journal of Volcanol-
ogy and Geothermal Research, 188(4):395–400, Dec. 2009. ISSN 0377-0273. doi:
10.1016/j.jvolgeores.2009.10.005.
Y. Moussallam, C. Oppenheimer, A. Aiuppa, G. Giudice, M. Moussallam, and
P. Kyle. Hydrogen emissions from Erebus volcano, Antarctica. Bulletin of
Volcanology, 74(9):2109–2120, Aug. 2012. ISSN 0258-8900, 1432-0819. doi:
10.1007/s00445-012-0649-2.
Y. Moussallam, C. Oppenheimer, B. Scaillet, F. Gaillard, P. Kyle, N. Peters,
M. Hartley, K. Berlo, and A. Donovan. Tracking the changing oxidation state of
Erebus magmas, from mantle to surface, driven by magma ascent and degassing.
Earth and Planetary Science Letters, 393:200–209, May 2014a. ISSN 0012-821X.
doi: 10.1016/j.epsl.2014.02.055.
Y. Moussallam, N. Peters, C. Ramı´rez, C. Oppenheimer, A. Aiuppa, and G. Giu-
dice. Characterisation of the magmatic signature in gas emissions from Turrialba
volcano, Costa Rica. Solid Earth Discussions, 6(2):2293–2320, Aug. 2014b. ISSN
1869-9537. doi: 10.5194/sed-6-2293-2014.
P. A. Nadeau, J. L. Palma, and G. P. Waite. Linking volcanic tremor, degassing,
and eruption dynamics via SO2 imaging. Geophysical Research Letters, 38(1):
L01304, Jan. 2011. ISSN 0094-8276. doi: 10.1029/2010GL045820.
J. M. Oberhuber, M. Herzog, H.-F. Graf, and K. Schwanke. Volcanic plume simu-
lation on large scales. Journal of Volcanology and Geothermal Research, 87(1–4):
29–53, Dec. 1998. ISSN 0377-0273. doi: 10.1016/S0377-0273(98)00099-7.
C. Oppenheimer and P. R. Kyle. Probing the magma plumbing of Erebus volcano,
Antarctica, by open-path FTIR spectroscopy of gas emissions. Journal of Vol-
canology and Geothermal Research, 177(3):743–754, Nov. 2008. ISSN 0377-0273.
doi: 10.1016/j.jvolgeores.2007.08.022.
143
C. Oppenheimer, A. J. S. McGonigle, P. Allard, M. J. Wooster, and V. Tsanev.
Sulfur, heat, and magma budget of Erta ’Ale lava lake, Ethiopia. Geology, 32(6):
509–512, Jan. 2004. ISSN 0091-7613, 1943-2682. doi: 10.1130/G20281.1.
C. Oppenheimer, A. S. Lomakina, P. R. Kyle, N. G. Kingsbury, and M. Boichu.
Pulsatory magma supply to a phonolite lava lake. Earth and Planetary Science
Letters, 284(3-4):392–398, 2009. ISSN 0012-821X. doi: 10.1016/j.epsl.2009.04.043.
C. Oppenheimer, R. Moretti, P. R. Kyle, A. Eschenbacher, J. B. Lowenstern, R. L.
Hervig, and N. W. Dunbar. Mantle to surface degassing of alkalic magmas at
Erebus volcano, Antarctica. Earth and Planetary Science Letters, 306(3–4):261–
271, June 2011. ISSN 0012-821X. doi: 10.1016/j.epsl.2011.04.005.
T. R. Orr and J. C. Rea. Time-lapse camera observations of gas piston activity at
Pu’u ’O¯’o¯, Kı¯lauea volcano, Hawai’i. Bulletin of Volcanology, 74(10):2353–2362,
Oct. 2012. ISSN 0258-8900, 1432-0819. doi: 10.1007/s00445-012-0667-0.
E. Pacaud. Aravis - GNOME Live!, Jan. 2011. URL http://live.gnome.org/
Aravis.
M. R. Patrick, T. Orr, D. Wilson, D. Dow, and R. Freeman. Cyclic spattering,
seismic tremor, and surface fluctuation within a perched lava channel, Kı¯lauea
Volcano. Bulletin of Volcanology, 73(6):639–653, Aug. 2011. ISSN 0258-8900,
1432-0819. doi: 10.1007/s00445-010-0431-2.
M. R. Patrick, T. Orr, L. Antolik, L. Lee, and K. Kamibayashi. Continuous monitor-
ing of Hawaiian volcanoes with thermal cameras. Journal of Applied Volcanology,
3(1):1, 2014. ISSN 2191-5040. doi: 10.1186/2191-5040-3-1.
T. D. Pering, G. Tamburello, A. J. S. McGonigle, A. Aiuppa, A. Cannata, G. Giu-
dice, and D. Patane`. High time resolution fluctuations in volcanic carbon dioxide
degassing from Mount Etna. Journal of Volcanology and Geothermal Research,
270:115–121, Jan. 2014. ISSN 0377-0273. doi: 10.1016/j.jvolgeores.2013.11.014.
N. Peters. AvoPlot: An extensible scientific plotting tool based on matplotlib.
Journal of Open Research Software, 2(1), Feb. 2014. ISSN 2049-9647. doi: 10.
5334/jors.ai.
N. Peters, A. Hoffmann, T. Barnie, M. Herzog, and C. Oppenheimer. Use of mo-
tion estimation algorithms for improved flux measurements using SO2 cameras.
Journal of Volcanology and Geothermal Research, 2014a. ISSN 0377-0273. doi:
10.1016/j.jvolgeores.2014.08.031.
144
N. Peters, C. Oppenheimer, D. R. Killingsworth, J. Frechette, and P. Kyle. Corre-
lation of cycles in Lava Lake motion and degassing at Erebus Volcano, Antarc-
tica. Geochemistry, Geophysics, Geosystems, 15(8):3244–3257, Aug. 2014b. ISSN
15252027. doi: 10.1002/2014GC005399.
N. Peters, C. Oppenheimer, and P. Kyle. Autonomous thermal camera system
for monitoring the active lava lake at Erebus volcano, Antarctica. Geoscientific
Instrumentation, Methods and Data Systems, 3(1):13–20, Feb. 2014c. ISSN 2193-
0864. doi: 10.5194/gi-3-13-2014.
N. Peters, C. Oppenheimer, P. Kyle, and N. Kingsbury. Decadal persistence of cycles
in lava lake motion at Erebus volcano, Antarctica. Earth and Planetary Science
Letters, 395:1–12, June 2014d. ISSN 0012-821X. doi: 10.1016/j.epsl.2014.03.032.
U. Platt and J. Stutz. Differential Absorption Spectroscopy. In Differential Opti-
cal Absorption Spectroscopy, pages 135–174. Springer Berlin Heidelberg, Berlin,
Heidelberg, 2008. ISBN 978-3-540-21193-8, 978-3-540-75776-4. URL http:
//link.springer.com/content/pdf/10.1007/978-3-540-75776-4_6.
R. Polikar. The Wavelet Tutorial, Oct. 2010. URL http://users.rowan.edu/
~polikar/WAVELETS/WTtutorial.html.
N. Rappin and R. Dunn. wxPython In Action. Manning Publication Co., 2006.
M. Richter and T. Schreiber. Phase space embedding of electrocardiograms. Physical
Review E, 58(5):6392–6398, Nov. 1998. doi: 10.1103/PhysRevE.58.6392.
M. Ripepe, A. J. L. Harris, and R. Carniel. Thermal, seismic and infrasonic evidences
of variable degassing rates at Stromboli volcano. Journal of Volcanology and
Geothermal Research, 118(3–4):285–297, Nov. 2002. ISSN 0377-0273. doi: 10.
1016/S0377-0273(02)00298-6.
G. Salerno, M. Burton, C. Oppenheimer, T. Caltabiano, D. Randazzo, N. Bruno, and
V. Longo. Three-years of SO2 flux measurements of Mt. Etna using an automated
UV scanner array: Comparison with conventional traverses and uncertainties in
flux retrieval. Journal of Volcanology and Geothermal Research, 183(1–2):76–83,
May 2009. ISSN 0377-0273. doi: 10.1016/j.jvolgeores.2009.02.013.
J. Sanders. Veusz - A Scientific Plotting Package, 2013. URL http://home.gna.
org/veusz/.
145
G. Sawyer, G. Salerno, J. Le Blond, R. Martin, L. Spampinato, T. Roberts,
T. Mather, M. Witt, V. Tsanev, and C. Oppenheimer. Gas and aerosol emissions
from Villarrica volcano, Chile. Journal of Volcanology and Geothermal Research,
203(1–2):62–75, June 2011. ISSN 0377-0273. doi: 10.1016/j.jvolgeores.2011.04.
003.
G. M. Sawyer and M. R. Burton. Effects of a volcanic plume on thermal imag-
ing data. Geophysical Research Letters, 33:4 PP., July 2006. doi: 200610.1029/
2005GL025320.
T. Schreiber and A. Schmitz. Surrogate time series. arXiv e-print chao-dyn/9909037,
Sept. 1999. URL http://arxiv.org/abs/chao-dyn/9909037. Physica D 142
(2000) 346-382.
H. Shinohara. A new technique to estimate volcanic gas composition: plume
measurements with a portable multi-sensor system. Journal of Volcanology
and Geothermal Research, 143(4):319–333, May 2005. ISSN 0377-0273. doi:
10.1016/j.jvolgeores.2004.12.004.
H. Shinohara and J. B. Witter. Volcanic gases emitted during mild Strombolian
activity of Villarrica volcano, Chile. Geophysical Research Letters, 32(20):n/a–
n/a, 2005. ISSN 1944-8007. doi: 10.1029/2005GL024131.
J. Smekens, M. R. Burton, A. B. Clarke, A. Harijoko, H. Wibowo, and G. Sawyer.
High frequency SO2 flux measurements at Semeru volcano, Indonesia, using the
SO2 camera. AGU Fall Meeting Abstracts, -1:2882, Dec. 2013.
L. Spampinato, S. Calvari, C. Oppenheimer, and E. Boschi. Volcano surveillance
using infrared cameras. Earth-Science Reviews, 106(1–2):63–91, May 2011. ISSN
0012-8252. doi: 10.1016/j.earscirev.2011.01.003.
D. Sweeney, P. R. Kyle, and C. Oppenheimer. Sulfur dioxide emissions and degassing
behavior of Erebus volcano, Antarctica. Journal of Volcanology and Geothermal
Research, 177(3):725–733, Nov. 2008. ISSN 0377-0273. doi: 10.1016/j.jvolgeores.
2008.01.024.
G. Tamburello, E. Kantzas, A. McGonigle, A. Aiuppa, and G. Giudice. UV cam-
era measurements of fumarole field degassing (La Fossa crater, Vulcano Island).
Journal of Volcanology and Geothermal Research, 199(1–2):47–52, Jan. 2011. ISSN
0377-0273. doi: 10.1016/j.jvolgeores.2010.10.004.
146
G. Tamburello, A. Aiuppa, E. Kantzas, A. McGonigle, and M. Ripepe. Passive vs.
active degassing modes at an open-vent volcano (Stromboli, Italy). Earth and
Planetary Science Letters, 359–360:106–116, Dec. 2012. ISSN 0012-821X. doi:
10.1016/j.epsl.2012.09.050.
G. Tamburello, A. Aiuppa, A. J. S. McGonigle, P. Allard, A. Cannata, G. Giudice,
E. P. Kantzas, and T. D. Pering. Periodic volcanic degassing behavior: The Mount
Etna example. Geophysical Research Letters, 40(18):4818–4822, 2013. ISSN 1944-
8007. doi: 10.1002/grl.50924.
H. Tazieff. Permanent lava lakes: observed facts and induced mechanisms. Journal of
Volcanology and Geothermal Research, 63(1–2):3–11, Oct. 1994. ISSN 0377-0273.
doi: 10.1016/0377-0273(94)90015-9.
C. Textor, H. F. Graf, M. Herzog, J. M. Oberhuber, W. I. Rose, and G. G. J. Ernst.
Volcanic particle aggregation in explosive eruption columns. Part II: Numerical
experiments. Journal of Volcanology and Geothermal Research, 150(4):378–394,
Feb. 2006a. ISSN 0377-0273. doi: 10.1016/j.jvolgeores.2005.09.008.
C. Textor, H. F. Graf, M. Herzog, J. M. Oberhuber, W. I. Rose, and G. G. J. Ernst.
Volcanic particle aggregation in explosive eruption columns. Part I: Parameter-
ization of the microphysics of hydrometeors and ash. Journal of Volcanology
and Geothermal Research, 150(4):359–377, Feb. 2006b. ISSN 0377-0273. doi:
10.1016/j.jvolgeores.2005.09.007.
R. I. Tilling. Fluctuations in surface height of active lava lakes during 1972-
1974 Mauna Ulu Eruption, Kilauea Volcano, Hawaii. Journal of Geophysi-
cal Research: Solid Earth, 92(B13):13721–13730, 1987. ISSN 2156-2202. doi:
10.1029/JB092iB13p13721.
A. R. Van Eaton, M. Herzog, C. J. N. Wilson, and J. McGregor. Ascent dynamics
of large phreatomagmatic eruption clouds: The role of microphysics. Journal of
Geophysical Research: Solid Earth, 117(B3):B03203, Mar. 2012. ISSN 2156-2202.
doi: 10.1029/2011JB008892.
G. van Rossum. Python, 2013. URL http://www.python.org.
A. C. Vandaele, P. C. Simon, J. M. Guilmot, M. Carleer, and R. Colin. SO2 absorp-
tion cross section measurement in the UV using a Fourier transform spectrometer.
Journal of Geophysical Research: Atmospheres, 99(D12):25599–25605, 1994. ISSN
2156-2202. doi: 10.1029/94JD02187.
147
G. Ward and A. Baxter. Distributing Python Modules — Python v2.7.5 documen-
tation, 2013. URL http://docs.python.org/2/distutils/.
F. Witham and E. W. Llewellin. Stability of lava lakes. Journal of Volcanology
and Geothermal Research, 158(3-4):321–332, Nov. 2006. ISSN 0377-0273. doi:
10.1016/j.jvolgeores.2006.07.004.
F. Witham, A. W. Woods, and C. Gladstone. An analogue experimental model of
depth fluctuations in lava lakes. Bulletin of Volcanology, 69(1):51–56, Apr. 2006.
ISSN 0258-8900. doi: 10.1007/s00445-006-0055-8.
J. B. Witter, V. C. Kress, P. Delmelle, and J. Stix. Volatile degassing, petrology,
and magma dynamics of the Villarrica Lava Lake, Southern Chile. Journal of
Volcanology and Geothermal Research, 134(4):303–337, July 2004. ISSN 0377-
0273. doi: 10.1016/j.jvolgeores.2004.03.002.
R. Wright and E. Pilger. Satellite observations reveal little inter-annual variability
in the radiant flux from the Mount Erebus lava lake. Journal of Volcanology and
Geothermal Research, 177(3):687–694, Nov. 2008. ISSN 0377-0273. doi: 10.1016/
j.jvolgeores.2008.03.005.
T. Wright. Vispect, 2008. URL http://vispect.sourceforge.net/.
D. Zandomeneghi, R. Aster, P. Kyle, A. Barclay, J. Chaput, and H. Knox. Internal
structure of Erebus volcano, Antarctica imaged by high-resolution active-source
seismic tomography and coda interferometry. Journal of Geophysical Research:
Solid Earth, 118(3):1067–1078, 2013. ISSN 2169-9356. doi: 10.1002/jgrb.50073.
148