<figure><img src="image_70.png" /><figcaption>**Fig. 4** The response of synthetic (Webb 1991; Bagdassarov and Dingwell 1993a) and natural (Webb 1992; Bagdassarov and Dingwell 1993b) silicate liquids to low strain, high frequency oscillatory deformation tests. This is shown as a normalized imaginary component of the</figcaption></figure>

shear modulus $G''/G_\infty$ as a function of the frequency of applied oscillation normalized by the relaxation timescale $\omega\lambda_r$, which is equivalent to a Deborah number De in the frequency domain. We show the viscoelastic prediction for a Maxwell liquid for comparison (solid curve)

frequency-domain version of the Deborah number. The metric that we choose to track is the imaginary component of the complex elastic modulus $G''$ normalized by $G_\infty$. When this value $G''/G_\infty$ is close to zero in the low Deborah number limit, the system is dominated by liquid behaviour and when it rises from zero at increasing Deborah number, there is an increasing component of the response to the forced oscillation that is elastic. We present collated data for a sodium disilicate synthetic composition (Webb 1991), rhyolitic compositions (Webb 1992; Bagdassarov and Dingwell 1993b) and a synthetic composition used as an analogue for type calc-alkaline rhyolite systems (Bagdassarov and Dingwell 1993a). This imaginary component can be described by the generalized Debye model for viscoelasticity for systems with a single characteristic relaxation time as $G'' = \omega\lambda_r / [(\omega\lambda_r)^2 + 1]$. The data reproduce the broad shape of the Debye model, albeit with an

overprediction of $G''$ around $\omega\lambda_r \sim 1$ for the rhyolitic liquids, which is poorly understood. Nevertheless, we point out that the first onset of a measurable elastic component to the response of the liquid to deformation occurs at $De = \omega\lambda_r \approx 10^{-2}$. Similarly, the point above which the majority of the response to deformation is elastic (the peak of $G''$) occurs at $De = \omega\lambda_r \approx 1$ (Fig. 3). We use this observation to validate the two thresholds found in Fig. 3. This implies that fractures can propagate in silicate liquids when there is even a small component of elastic behaviour ($De \ge 10^{-2}$) and that those fractures can propagate vigorously when the elastic behaviour dominates over the viscous behaviour ($De \ge 1$). An interpretation might also be that $10^{-2} \le De \le 1$ is the range in which fracture propagation is competing with viscous relaxation of stress and therefore the fractures are unlikely to be long, sharp-tipped or pervasive. And that at $De \ge 1$, stress dissipation by fracture
propagation can be localized onto longer, sharp
and pervasive fracture networks.

The fact that a viscous limit to the expansion
of vesicular magma also scales with our Deborah
number criterion is tantalizing. Microphysically,
in this case it is not fracture of a deforming
homogeneous liquid, but fractures propagating at
bubble walls as the local rate of strain in the
liquid induced by expansion of the bubble meets
the Deborah number criterion for fracture prop-
agation (c.f. Kameda and Kuribara 2008). For
complete rupture of the vesicular material—
which is a fragmentation event or *violent rupture*
in Fig. 3—the fractures propagating from the
bubbles must interact. Presumably in the range
$10^{-2} \le \text{De} < 1$, fractures propagate but the frac-
ture tips are blunted during competing viscous
relaxation such that brittle behaviour can be
observed but is not catastrophic to the system.
Then, as the local strain rate at the bubble wall
increases further and De approaches 1, the frac-
tures propagate in a dominantly elastic medium
and can span the inter-bubble distances, inter-
acting to produce violent rupture. Therefore, it is
clear that a scaling of the Deborah number con-
cept to bubble wall dynamics would provide a
fragmentation criterion for viscoelastic vesicu-
lating magma (c.f. Namiki and Manga 2005;
Koyaguchi et al. 2008; Namiki and Manga
2008).

# 4 Laboratory-Scale Unrest Signals

When magma breaks in a laboratory experiment,
acoustic emissions—packets of acoustic energy
—are released and can be recorded (Benson et al.
2008; Lavallée et al. 2008). These signals appear
to represent large total released amounts of
energy when $\lambda$ is short (high $\dot{y}$), compared to
when $\lambda$ is long (low $\dot{y}$) (Lavallée et al. 2008).
This supports our posit that fracture networks are
more likely to be pervasive and large when De is
large than when De is small. For single phase
silicate liquids deformed at high De, Tuffen et al.
(2008) showed that the experimental acoustic
event frequency range and sample fracture length

scale scaled with natural constraints of in-conduit
volcanic fracture systems and natural frequencies
of low-frequency volcano-seismicity, indicating
that these viscoelastic fracturing events are
indeed the likely source mechanism for
low-frequency earthquakes at volcanoes. This
confirmed the conclusion of Neuberg et al.
(2006), who showed that these low-frequency
events were most likely to be associated with
repetitive fracturing of magma at a given depth in
the conduit during ascent. Neuberg et al. (2006)
further modelled magma ascent in confined
geometry and showed that a De ~ 1 is met at a
depth of 830 m below the surface using param-
eters typical for recent eruptions at Soufriere
Hills volcano and the magma thereof. Thomas
and Neuberg (2012) predicted a deeper source of
1500 m using the same model approach but by
invoking a conduit restriction (see Chap. 9),
consistent with inversions for low-frequency
sources using seismic data. Therefore, vis-
coelastic magma fracturing in the high Deborah
number regime appears to be a consistent model
for the source mechanics of low-frequency vol-
cano seismicity, often used for eruption
forecasting.

Vasseur et al. (2015) showed that the fore-
castability of full sample rupture scales with the
heterogeneity of the system—cast most simply as
a porosity (Fig. 5). The implication is that the
less vesicular the magma undergoing deforma-
tion, the less likely that accurate forecasts based
on accelerated precursory signals can be made,
with up to ~120% error on the timing of the
rupture event for single-phase homogeneous
liquids in the high De limit. It may be that the
magma vesicularity, pore-network structure
(Vasseur et al. 2017), crystallinity and textural
anisotropy, play key roles in determining fore-
casting success based on low-frequency earth-
quakes. More experimental work is clearly
required in this area, along with more rigorous
scaling between acoustic and seismic events that
originate from magma failure.

An implication of the models explored here,
encapsulated by Eq. 1, is that the peak stress
supported by a liquid $\sigma_m$ can be predicted as a
<figure><img src="image_71.png" /><figcaption>**Fig. 5** In the high-De regime, the error on a prediction of failure times scales with the porosity of the material such that low porosity magmas are unpredictable and high porosity magmas are predictable (Vasseur et al. 2015).</figcaption></figure>

The failure times are recorded as a large stress drop during uniaxial loading, and the approach to failure is monitored using acoustic emissions generated by pre-failure micro-fracturing events

<figure><img src="image_72.png" /><figcaption>**Fig. 6** The peak stress that can be supported by a shearing liquid in the low De limit is given by $\sigma_m = G_\infty De$. And in the high De limit, this value appears to asymptotically approach $\sigma_m = 10^{-2} G_\infty$ which is equivalent to the strength of glassy materials in that same regime (Simmons et al. 1982)</figcaption></figure>

function of De by $\sigma_m = G_\infty De$. Once $De \ge 10^{-2}$, however, this linear relationship appears to be invalid as fractures can form in the liquid. It is

perhaps significant that the average strength of glass in the high De regime is $\sim 10^8$ Pa (Simmons et al. 1982; Vasseur et al. 2013), consistent
with a rupture threshold of $10^{-2}G_∞$. These relationships are explored in Fig. 6.

Here we summarize just one example of how targeted laboratory experiments and dimensional analysis have demonstrated the most likely source mechanism of unrest signals commonly used to monitor magma movement and predict impending eruptions. Clearly, with scaling arguments for the applicability of an experimental set up—such as we show here with the Deborah number analysis—laboratory-based work can provide new insights into unrest mechanisms that facilitate more accurate forward-modelling of geophysical signals.

## References

Alidibirov M, Dingwell DB (1996) Magma fragmentation by rapid decompression. Nature 380:146–148

Bagdassarov N, Dingwell DB, Webb SL (1993) Effect of boron, phosphorus and fluorine on shear stress relaxation in haplogranite melts. Eur J Miner 3:409–425

Bagdassarov NS, Dingwell DB (1993) Frequency dependent rheology of vesicular rhyolite. J Geophy Res: Solid Earth 98(B4):6477–6487

Benson P, Vinciguerra S, Meredith P, Young R (2008) Laboratory simulation of volcano seismicity. Science (80-) 322:249–252

Chouet BA (1996) Long-period volcano seismicity: its source and use in eruption forecasting. Nature 380:309–316. doi:10.1038/380309a0

Chouet BA, Page RA, Stephens CD et al (1994) Precursory swarms of long-period events at Redoubt Volcano (1989–1990), Alaska: their origin and use as a forecasting tool. J Volcanol Geotherm Res 62:95–135. doi:10.1016/0377-0273(94)90030-2

Cordonnier B, Caricchi L, Pistone M et al (2012a) The viscous-brittle transition of crystal-bearing silicic melt: direct observation of magma rupture and healing. Geology 40:611–614

Cordonnier B, Schmalholz SM, Hess K, Dingwell DB (2012b) Viscous heating in silicate melts: an experimental and numerical comparison. J Geophys Res 117:B02203

Costa A, Macedonio G (2005) Viscous heating effects in fluids with temperature-dependent viscosity: triggering of secondary flows. J Fluid Mech 540:21–38

Costa A, Melnik O, Vedeneeva E (2007) Thermal effects during magma ascent in conduits. J Geophys Res 112: B12205. doi:10.1029/2007JB004985

Deubelbeiss Y, Kaus BJ, Connolly JA, Caricchi L (2011) Potential causes for the non-Newtonian rheology of crystal-bearing magmas. Geochem Geophy Geosys 12(5)

Dingwell DB (1995) Relaxation in silicate melts; some applications. Rev Mineral Geochemistry 32:21–66

Dingwell DB (1996) Volcanic dilemma: flow or blow? Science (80-) 273:1054–1055. doi:10.1126/science.273.5278.1054

Dingwell DB, Webb SL (1989) Structural relaxation in silicate melts and non-Newtonian melt rheology in geologic processes. Phys Chem Miner 16:508–516

Giordano D, Russell JK, Dingwell DB (2008) Viscosity of magmatic liquids: a model. Earth Planet Sci Lett 271:123–134

Gonnermann HM, Manga M (2003) Explosive volcanism may not be an inevitable consequence of magma fragmentation. Nature 426:432–435

Hess KU, Dingwell DB (1996) Viscosities of hydrous leucogranitic melts: a non-Arrhenian model. Am Mineral 81:1297–1300

Hess K-U, Cordonnier B, Lavallée Y, Dingwell DB (2008) Viscous heating in rhyolite: an in situ experimental determination. Earth Planet Sci Lett 275:121–126. doi:10.1016/j.epsl.2008.08.014

Hornby AJ, Kendrick JE, Lamb OD, Hirose T, De Angelis S, Aulock FW, Umakoshi K, Miwa T, Henton De Angelis S, Wadsworth FB, Hess KU (2015) Spine growth and seismogenic faulting at Mt. Unzen, Japan. J Geophys Res: Solid Earth 120(6):4034–4054

Ichihara M, Rittel D, Sturtevant B (2002) Fragmentation of a porous viscoelastic material: implications to magma fragmentation. J Geophys Res: Solid Earth 107(B10)

Kameda M, Kuribara H, Ichihara M (2008) Dominant time scale for brittle fragmentation of vesicular magma by decompression. Geophy Res Lett 35(14)

Kameda M, Ichihara M, Shimanuki S, Okabe W, Shida T (2013) Delayed brittle-like fragmentation of vesicular magma analogue by decompression. J Volcanol Geoth Res 258:113–125

Kendrick JE, Lavallée, Y, Hirose T, Di Toro G, Hornby AJ, De Angelis S, Dingwell DB (2014) Volcanic drumbeat seismicity caused by stick-slip motion and magmatic frictional melting. Nat Geosci 7(6):438

Koyaguchi T, Scheu B, Mitani NK, Melnik O (2008) A fragmentation criterion for highly viscous bubbly magmas estimated from shock tube experiments. J Volcanol Geotherm Res 178:58–71

Lavallée Y, Meredith PG, Dingwell DB et al (2008) Seismogenic lavas and explosive eruption forecasting. Nature 453:507–510. doi:10.1038/nature06980

Lensky NG, Lyakhovsky V, Navon O (2001) Radial variations of melt viscosity around growing bubbles and gas overpressure in vesiculating magmas. Earth Planet Sci Lett 186:1–6

Mader HM, Llewellin EW, Mueller SP (2013) The rheology of two-phase magmas: a review and analysis. J Volcanol Geoth Res 257:135–158

Mastin LG (2005) The controlling effect of viscous dissipation on magma flow in silicic conduits. J Volcanol Geoth Res 143(1):17–28

Miller AD, Stewart RC, White R et al (1998) Seismicity associated with dome growth and collapse at the
Soufriere Hills Volcano, Montserrat. Geophys Res Lett 25:3401–3404

Mueller S, Llewellin EW, Mader HM (2010) The rheology of suspensions of solid particles. Proc R Soc A Math Phys Eng Sci 466:1201–1228

Namiki A, Manga M (2005) Response of a bubble bearing viscoelastic fluid to rapid decompression: implications for explosive volcanic eruptions. Earth Planet Sci Lett 236(1):269–284

Namiki A, Manga M (2008) Transition between fragmentation and permeable outgassing of low viscosity magmas. J Volcanol Geoth Res 169(1):48–60

Neuberg JW, Tuffen H, Collier L et al (2006) The trigger mechanism of low-frequency earthquakes on Montserrat. J Volcanol Geotherm Res 153:37–50. doi:10.1016/j.jvolgeores.2005.08.008

Prousevitch AA, Sahagian DL, Anderson AT (1993) Dynamics of diffusive bubble growth in magmas: isothermal case. J Geophys Res Solid Earth 98:22283–22307

Salvage R, Neuberg JW (2016) Using a cross correlation technique to refine the accuracy of the Failure Forecast Method: application to Soufrière Hills volcano, Montserrat. J Volcanol Geotherm Res 324:118–133. doi:10.1016/j.jvolgeores.2016.05.011

Sammis CG, Ashby MF (1986) The failure of brittle porous solids under compressive stress states. Acta Metall 34:511–526

Simmons JH, Mohr RK, Montrose CJ (1982) Non-Newtonian viscous flow in glass. J Appl Phys 53:4075–4080

Sparks RSJ (1978) The dynamics of bubble formation and growth in magmas: a review and analysis. J Volcanol Geotherm Res 3:1–37

Stevenson RJ, Dingwell DB, Webb SL, Bagdassarov NS (1995) The equivalence of enthalpy and shear stress relaxation in rhyolitic obsidians and quantification of the liquid-glass transition in volcanic processes. J Volcanol Geotherm Res 68:297–306

Thomas ME, Neuberg J (2012) What makes a volcano tick-a first explanation of deep multiple seismic

sources in ascending magma. Geology 40:351–354. doi:10.1130/g32868.1

Tuffen H, Dingwell DB, Pinkerton H (2003) Repeated fracture and healing of silicic magma generate flow banding and earthquakes? Geology 31:1089–1092. doi:10.1130/g19777.1

Tuffen H, Smith R, Sammonds PR (2008) Evidence for seismogenic fracture of silicic magma. Nature 453:511–514

Vasseur J, Wadsworth FB, Lavallée Y et al (2013) Volcanic sintering: timescales of viscous densification and strength recovery. Geophys Res Lett 40:5658–5664

Vasseur J, Wadsworth FB, Lavallée Y, Bell AF, Main IG, Dingwell DB (2015) Heterogeneity: the key to failure forecasting. Sci Rep 5

Vasseur J, Wadsworth FB, Heap MJ, Main IG, Lavallée Y, Dingwell DB (2017) Does an inter-flaw length control the accuracy of rupture forecasting in geological materials? Earth Planet Sci Lett

Voight B (1988) A method for prediction of volcanic eruptions. Nature 332:125–130

Webb SL (1991) Shear and volume relaxation in Na₂ Si₂ O₅. Am Mineral 76(9–10):1449–1454

Webb SL (1992) Low-frequency shear and structural relaxation in rhyolite melt. Phys Chem Mineral 19 (4):240–245

Webb SL, Dingwell DB (1990a) The onset of non-Newtonian rheology of silicate melts. Phys Chem Miner 17:125–132

Webb SL, Dingwell DB (1990b) Non-Newtonian rheology of igneous melts at high stresses and strain rates: experimental results for rhyolite, andesite, basalt, and nephelinite. J Geophys Res Solid Earth 95:15695–15701

Zhang Y, Stolper EM, Wasserburg GJ (1991) Diffusion of a multi-species component and its role in oxygen and water transport in silicates. Earth Planet Sci Lett 103:228–240
**Open Access** This chapter is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, duplication, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, a link is provided to the Creative Commons license and any changes made are indicated.

The images or other third party material in this chapter are included in the work's Creative Commons license, unless indicated otherwise in the credit line; if such material is not included in the work's Creative Commons license and the respective action is not permitted by statutory regulation, users will need to obtain permission from the license holder to duplicate, adapt or reproduce the material.
<figure><img src="image_73.png" /></figure>

# Volcano Seismology: Detecting Unrest in Wiggly Lines

R.O. Salvage, S. Karl and J.W. Neuberg

## Abstract

Seismology is a useful tool to gain a better understanding of volcanic unrest in real time as it unfolds. The generation of seismic signals in a volcanic environment has been linked to a number of different physical processes occurring at depth, including fracturing of the volcanic edifice (producing high frequency seismicity) and movement of magmatic fluids (producing low frequency seismicity). Further classification of seismic signals according to their waveform similarity, in addition to their frequency content, allows greater detail in temporal and spatial changes of seismicity to be detected. At Soufrière Hills volcano, Montserrat, one of the target volcanoes of the VUELCO project, families of similar waveforms provided valuable insight into evaluating the significance of ongoing unrest. In June 1997 over 6000 more events were able to be identified over a 5 day period of interest (22 to 25 June) by using families of seismic events, rather than a standard amplitude-based detection algorithm. In total, 11 families were identified, with the events clustering into a number of swarms, suggesting a repeating and non destructive cyclic source mechanism. Since each family is believed to represent a distinct source location and mechanism, identifying 11 coexisting families

<footnote>R.O. Salvage (✉)
Observatorio Vulcanológico y Sismológico
de Costa Rica, Universidad Nacional,
Apartado Postal: 2386-3000, Heredia, Costa Rica
e-mail: beckysalvage@gmail.com</footnote>

<footnote>S. Karl
Shell, NAM offices, Schepersmaat 2,
9405 TA Assen, The Netherlands
e-mail: sandra.karl84@gmail.com</footnote>

<footnote>J.W. Neuberg
School of Earth and Environment, Institute of
Geophysics and Tectonics, University of Leeds,
Leeds LS2 9JT, United Kingdom
e-mail: j.neuberg@leeds.ac.uk</footnote>
reflects the complex diversity of physical processes which act simulta-
neous at this volcano. In July 2003, conditions at the volcano had clearly
changed since only one family of seismicity was identified. The source
location of this family appeared to shift with time from 8 July (when no
events from the family were identified) to 12 July (where most events had
a cross correlation coefficient over 0.9). In addition, the use of families
appears to greatly aid hindsight forecasting attempts for the large scale
dome collapses of 1997 and 2003 using the Failure Forecast Method.
Knowledge of the temporal and spatial extent of seismicity during periods
of unrest, its source mechanism and its relationship to physical processes
at depth is essential for decision and policy makers for risk mitigation.
However, the source mechanisms of such volcanic seismicity is still much
debated and appears to often be misinterpreted because of compromising
assumptions used in the numerical modelling of inverting such sources.
Use of a spatially extended source such as a ring fault structure, rather than
a single point for determining the origin of low frequency seismicity, is
now thought to be more realistic for the mechanism of such events since it
more accurately represents the movement of magma through a conduit.
However, use of this spatially extended source instead of a simple single
point results in a large underestimation of slip from P-wave amplitudes,
which may lead to an underestimation in magma ascent rates, with large
consequences for eruption forecasting. Additionally, the P-wave radiation
patterns exhibited by these two mechanisms are remarkably similar, and
can only be distinguished if the small radial radiation lobes can be
determined. In a volcanic environment this is extremely difficult due to
large uncertainties in earthquake source depth locations, and the
implementation of small aperture seismic networks.

Español

La sismología es una herramienta geofísica valiosa que brinda infor-
mación en tiempo real, permitiendo una mejor comprensión del compor-
tamiento de sistemas volcánicos que inician un proceso de reactivación o
de intensificación de la actividad. La generación de señales sísmicas en
ambientes volcánicos se ha relacionado con un número diverso de
procesos geofísicos que ocurren en el interior de los volcanes, incluyendo
fracturamiento del edificio volcánico (produciéndose sismicidad de alta
frecuencia) y movimiento de fluidos magmáticos (produciéndose sismi-
cidad de baja frecuencia). La clasificación de señales sísmicas basada en la
similitud de las formas de onda, además del contenido de frecuencias, ha
permitido detectar cambios temporales y espaciales de la sismicidad con
mayor detalle.

En el volcán Soufriere Hills en Monserrat, uno de los volcanes
investigados como parte del Proyecto VUELCO, el reconocimiento de
familias de señales sísmicas con formas de onda similares proveyó un
entendimiento valioso en la evaluación de la significancia de la
reactivación de su actividad volcánica. En junio de 1997, se pudieron identificar más de 6000 eventos ssmicos dentro de un periodo particular de 4 días (entre el 22 y el 25 de junio) mediante la determinación de familias de eventos sísmicos que fueron determinados por un algoritmo estándar de detección basado en la amplitud de la señal sísmica. En total, 11 familias de eventos sísmicos fueron identificadas, con grupos de eventos conformando enjambres sísmicos, lo que sugiere un mecanismo cíclico repetitivo de una fuente sísmica no destructiva. Ya que se considera que cada familia representa una fuente con ubicación espacial y mecanismo distinto, la identificación de 11 familias refleja la compleja diversidad de los procesos geofísicos que operan simultáneamente en este volcán en particular. En julio del 2003, el régimen del volcán cambió definitivamente ya que solamente 1 familia de eventos sísmicos fue identificada. La ubicación de la fuente sísmica de esta familia parece haber migrado con el tiempo entre el 8 de julio (cuando ningún evento sísmico de esta familia fue identificado) y el 12 de julio (cuando la mayoría de eventos tuvo un coeficiente de cros-correlación superior a 0,9). Además, el establecimiento de familias de sismos parece haber sido de gran ayuda en los intentos de pronosticar los colapsos del domo en gran escala de 1997 y el 2003 utilizando el método determinístico de pronóstico de rompimiento por fatiga (Failure Forecast Method, FFM). El conocimiento de la extensión temporal y espacial de la sismicidad durante periodos de reactivación volcánica, el mecanismo de la fuente, y su relación con los procesos geofísicos a niveles profundos son aspectos esenciales para los tomadores de decisiones y ejecutores de políticas relacionadas con la mitigación de riesgos.

En general, los mecanismos de la fuente generadora de la sismicidad volcánica continúan aún bajo gran debate y parecen ser con frecuencia mal interpretados debido a simplificaciones comprometedoras usadas en los modelos numéricos de inversión de dichas fuentes. El uso de una fuente espacial extendida tal como una falla estructural tipo anular, en lugar de un solo punto para determinar el origen de sismicidad de baja frecuencia, se piensa es más realista para visualizar el mecanismo de este tipo de eventos ya que representa con más exactitud el movimiento del magma a través de un conducto. Sin embargo, el uso de este modelo resulta en la subestimación significativa del deslizamiento en la falla determinada a través de las amplitudes de las ondas P, lo que podría llevar a una subestimación de las tasas de ascenso del magma, con consecuencias cruces para el pronóstico de erupción inminente. Por otra parte, los patrones de radiación de la onda P mostrados por estos dos mecanismos son marcadamente similares, y pueden ser solamente distinguidos si los pequeños lóbulos radiales de radiación pueden ser determinados. En un contexto volcánico, esta es una tarea extremadamente difícil debido a la gran incertidumbre inherente a la localización de la profundidad de la fuente sísmica, y a la implementación de redes sísmicas de pequeña extensión.
**Keywords**

Volcano seismology · Families of similar seismic events · Magmatic source mechanism · Eruption forecasting

**Index Terms**

Volcano seismology · Soufrière Hills · Chiles-Cerro Negro · Low frequency seismicity · Families · Cross correlation · Forecasting · Failure forecast method · Source mechanisms · Point source · CLVD · Ring fault · Spatially extended source · P-wave radiation patterns

## Volcanic Unrest

The monitoring of active volcanoes for the protection of society has excelled in recent decades fuelled by rapid technological advances, which have allowed the development and deployment of more cost-effective monitoring solutions. It is now possible to detect geophysical and geochemical signals from a volcano which previously would have been below the detection threshold. The routine monitoring of volcanoes during periods of quiescence is crucial, although not always feasible, in order to assess background levels of activity at volcanoes, and thus more rapidly detect the onset of future unrest. Most simply, volcanic unrest is defined as a deviation from background levels of activity towards a level which is a cause for concern over short time scales of hours to days. Volcanic unrest does not necessarily lead to eruption, although this is the most likely outcome, with 64% of 228 volcanoes that have experienced unrest since the year 2000 culminating in eruptive events (Phillipson et al. 2013).

Current monitoring efforts at volcanoes can be grouped into three categories: measurements of surface degassing; deformation; and seismic activity, which are all thought to result from the movement of magmatic fluids at depth, and therefore may provide key insights into an impending eruption. Unrest must be detected on short timescales which is appropriate for decision-making, and therefore seismicity has remained a primary monitoring tool, as it can be

remotely analysed in real-time, often by an automated system, if a number of sensors are placed around the volcano. Deviation from the background level is also often easier to determine for seismicity than for other signals.

## Seismic Event Characterisation

A wide variety of seismic signals exist within volcanic settings, associated with magmatic and hydrothermal fluid movement at depth, pressurization of the volcanic edifice, and/or the surface manifestation of the interaction of these processes. The variety of signals is a reflection of the number of different processes and the great structural heterogeneities found within this context. The characterisation of seismicity can be based upon waveform similarities, but is traditionally based upon the signals' time and frequency characteristics: different bands of frequency relate to different active source processes at depth, which can then be distinguished from one another, although the frequency bands associated with each process may overlap (Lahr et al. 1994).

The identification of seismic events in volcanic settings can assist with the detection of the onset and cessation of unrest, while the spatial and temporal patterns of occurrence may be informative of magma movement and changes in stress at depth, as well as the spatial extent of concern. An understanding of the physical processes occurring at depth is essential if accurate
<figure><img src="image_74.png" /><figcaption>**Fig. 1** Examples of waveforms and their frequency content seen in volcanic environments taken from Soufrière Hills Volcano, Montserrat in 1997. Soufrière Hills Volcano was a target volcano identified by the VUELCO project for investigation. **a**, **b** High frequency</figcaption></figure>

waveform with clear phase arrivals. **c**, **d** Low frequency waveforms with an emergent onset. Waveform filtered between 0.5 and 5 Hz. **e**, **f** Rockfall event with classic “cigar” shape

and timely forecasts of volcanic eruptions, and developing unrest scenarios, are to be made.

## Classification by Frequency Content

Seismic signals originating from processes occurring at depth within the volcanic system are usually split into high- and low-frequency end-members, although in reality a continuum across the spectrum exists between the two (Chouet and Matoza 2013). High frequency seismic signals (Fig. 1a, b), also known as Volcano-Tectonic (VT) events, have energy concentrated in the frequency range of 1 to 20 Hz, generally peaking between 6 and 8 Hz (Lahr et al. 1994). They are characterized by clear, impulsive P- and S-wave arrivals, followed by a short coda. High frequency seismicity is usually attributed to brittle failure within the volcanic edifice, where magmatic processes create enough elastic strain to force the surrounding rocks into failure (Arciniega-Ceballos et al.

2003), similar to the generation of tectonic earthquakes.

Low frequency (LF) seismic signals (Fig. 1c, d) usually occupy the spectral range of 0.2 to 5 Hz (Chouet and Matoza 2013), and are frequently characterized by emergent P-wave onsets and lack of S-wave arrivals. It has been suggested that the occurrence of low frequency waveforms is linked to resonance of seismic energy trapped at a solid-fluid interface either within a crack (e.g. Chouet 1988), or a volcanic conduit (e.g. Neuberg et al. 2000). The trigger mechanism of such seismic energy is further disputed, with suggestions that it may be generated by: (1) a stick-slip motion along the conduit walls as magma ascends (e.g. Iverson et al. 2006); (2) the brittle failure of magma itself either through an increase in viscosity and strain rates (Lavallée et al. 2008), which may be due to an increase in the ascent rate of magma through the conduit (Neuberg et al. 2006), changes in the crystal and/or bubble concentration in the magma (Goto 1999), or through a change in the
geometry of the conduit (Thomas and Neuberg 2012); (3) the interaction between the magmatic and hydrothermal system at depth (e.g. Nakano and Kumagai 2005); or (4) through slow rupture and failure of unconsolidated material on volcanic slopes (Bean et al. 2014).

Many volcanic seismic events fall between these two end-member categories and are termed “hybrid” events. Hybrid events are characterised by a high frequency onset with a long resonating low frequency coda, therefore distributing energy across a wider frequency spectrum (Chouet and Matoza 2013). Hybrid and LF events are often classified in the same group of volcanic seismicity, since source and path effects can result in a LF event recorded at one station being recorded as a hybrid event at another.

With the deployment of broadband sensors in many volcanic environments, it is now possible to detect seismicity within a much wider frequency band, up to 120 s periods (Chouet and Matoza 2013), known as Very Long Period (VLP) earthquakes. VLP events occupy the spectral range below 0.01 Hz. The generation of these waveforms is not yet fully understood, in particular how such long wavelengths can be generated in apparently small source volumes, although it has been linked to perturbations in the flow of fluid or gas within pressurized volcanic conduits or cracks (e.g. Dawson et al. 2011). Their large wavelength, sometimes over one hundred kilometers, means few path effects on the waveform and as such if identified, these waveforms provide an excellent choice for performing waveform inversion techniques to identify source characteristics (Chouet and Matoza 2013).

Furthermore, seismicity can be generated by surface processes, such as landslides, rockfall events, pyroclastic flows and lahars (Fig. 1e, f). These are particularly dominant during dome building eruptions and at volcanoes with glaciers during the spring and summer months due to partial melting of the ice (McNutt 2005). These signals can be exploited to determine the size and magnitude of such events, their location and their direction of travel (e.g. De Angelis et al. 2007). Typically rockfall events (small free falling rock

events) form a “cigar shaped” waveform with an emergent onset (Fig. 1e, f), whereby there is an initial increasing amplitude of the waveform as the amount of material falling down slope increases. Pyroclastic flow signals are distinguishable from rockfalls since their waveforms are at least an order of magnitude larger and they often occur over a longer duration since larger amounts of material are involved moving down slope (De Angelis et al. 2007), however the two are likely to exist on a continuum.

## Classification by Waveform Similarity

Seismic waveforms can also be classified according to their similarity with other detected seismic events. The frequency content of seismic waveforms is indicative of the active processes that may be occurring within the volcanic environment and the source mechanism involved in the generation of such seismicity. The further classification of seismic events into families which all have a similar waveform shape, as well as the same frequency content, allows the depiction of temporal and spatial changes in the source mechanism and the source location on a much smaller scale (e.g. Thelen et al. 2011; Salvage and Neuberg 2016). For example, the relative relocation of families of similar seismic events at Soufrière Hills volcano, Montserrat has produced very precise source locations (e.g. De Angelis and Henton 2011). By definition, families of seismic events must be generated by the same mechanism and at the same location in order for the detected waveforms to have the same shape at the seismometer, and therefore changes in either of these parameters affect the similarity of waveforms. In many instances it is assumed that families of seismic events are generated by the same mechanism and within a similar source location, estimated at between one quarter and one tenth of the wavelength (Geller and Mueller 1980; Neuberg et al. 2006).

Waveform similarity in terms of shape and duration can be evaluated by cross correlation. Identical signals will result in a cross correlation coefficient of 1 or -1, dependent upon their
<figure><img src="image_75.png" /><figcaption>**Fig. 2** Cross Correlation Matrix of events identified at Soufrière Hills volcano, Montserrat at a single station in June 1997. A total of 486 events were identified on 24 June 1997 and are shown in temporal order along the x and y axis. Events with a cross correlation coefficient of</figcaption></figure>

greater than 0.7 are shown on a colour scale, with those close to one being more similar. The autocorrelation of each event with itself is shown in *dark red* along the diagonal and is equal to a cross correlation coefficient of 1

relative polarity. Signals with no correlation result in a cross correlation coefficient of 0. A threshold must be chosen above which waveforms can be considered similar. Waveforms which are deemed similar can be grouped together into a family of events. The choice of similarity threshold is important: if it is too low there is a risk of placing events which are not similar into the same family; if it is too high similar events can be missed. Green and Neuberg (2006), Thelen et al. (2011) and Salvage and Neuberg (2016) suggest a cross correlation coefficient threshold of 0.7, since this is significantly above the correlation coefficient that can be produced from random correlations between noise and a waveform. Higher cross correlation coefficient thresholds can be used to identify families of almost identical waveforms, however Petersen (2007) suggests that this is probably not appropriate in volcanic settings due to additional noise in this environment.

The similarity between identified seismic events can be determined by cross correlating

each individual seismic event with every other seismic event. The result are typically presented as a similarity matrix, as seen in Fig. 2, where events which are deemed to be similar are shown on the colour spectrum. However, such a matrix may include a number of families of similar events since it only determines whether each event shows similarity to any of the other earthquakes analysed. In order to identify families of similar events, events with a high cross correlation coefficient as decided by the user are grouped together and removed from the matrix. This procedure is repeated across the entire investigated time period until all events have been classified into a family, or have been removed from the matrix. A master event is then determined from each family of events as the average of the stack of similar waveforms. This is representative of the family in terms of waveform shape.

Families of similar seismicity have been identified at a number of active volcanoes around the world, including Redoubt volcano, Alaska
<figure><img src="image_76.png" /><figcaption>**Fig. 3** Comparison of the timing and duration of swarms of families of events identified at a single station at Soufrière Hills volcano, Montserrat in June 1997. The timing of the dome collapse is represented by the vertical line on the 25 June 1997. The y axis is only an indication of each of the families present separated in space for the</figcaption></figure>

purpose of clarity on the plot and does not represent time or dominance; each master event is simply drawn below the last so that all can be compared. Each coloured rectangular box represents the times when the families were active during the 22–25 June analysis period

(e.g. Buurman et al. 2013); Mt. St. Helens, USA (Thelen et al. 2011); Colima, Mexico (Arámbula-Mendoza et al. 2011); Merapi, Indonesia (Budi-Santoso and Lesage 2016); Katla, Iceland (Svattoni et al. 2016); and Soufrière Hills volcano, Montserrat (Green and Neuberg 2006; Ottemöller 2008; Salvage and Neuberg 2016). The identification of families rather than simply detecting seismic events and classifying them according to their frequency content is advantageous as subtle temporal and spatial patterns can be identified, allowing detailed source information to be uncovered. In addition, this technique dramatically increases the number of identified events from the continuous seismic record, since low amplitude events and closely spaced events can still be identified. For example, using a standard amplitude-based detection algorithm, 1435 seismic events were identified at Soufrière Hills volcano, Montserrat between 22 and 25 June 1997, a period of interest due to increased seismicity

before a lava dome collapse. The cross correlation technique identified 7653 similar seismic events during the same time period, offering a five-fold increase in the number of detected earthquakes (Salvage and Neuberg 2016).

Low frequency families of seismicity identified during this unrest period at Soufrière Hills volcano were followed by a dome collapse on 25 June 1997. Soufrière Hills Volcano was chosen by the VUELCO project as a target volcano due to the longevity of its dome building and collapse cycles which have been ongoing since 1995, providing a wealth of associated geophysical data (Sparks and Young 2002; Wadge et al. 2014). In total, 11 distinct seismic sources (i.e. 11 families of seismicity) were identified during this period of unrest [Fig. 3; Green and Neuberg (2006); Salvage and Neuberg (2016)], which all broadly follow the same temporal pattern in the number of identified swarms present over this time period. However, the timing and duration of these
<figure><img src="image_77.png" /><figcaption>**Fig. 4** The evolution of the dominant cross correlation coefficient with time at Soufrière Hills volcano, Montserrat in July 2003. A dome collapse event occurred at the time of the *vertical line* (13:30 on 12 July 2003). The temporal gaps in the data represent drops in the seismometer recordings rather than a change in the cross correlation coefficient, indicated by the white space on the x-axis</figcaption></figure>

swarms can be seen to be different for each family of similar seismic events. Low frequency seismicity is associated with the movement of magnetic fluid at depth, and therefore in this case would suggest cyclic flow dynamics to generate such swarm-like behaviour. The source process for the generation of this seismicity must be stable and non-destructive in order to be repeatable (Green and Neuberg 2006; Petersen 2007), and must be able to occur at a number of different locations and/or by a number of different sources at the same time in order to generate a number of active families of events, reflecting the complex diversity of seismic sources and physical processes which act simultaneously at this volcano.

The identification of families can also be used to understand evolving seismicity with time. An evolving cross correlation coefficient with time, if not an artefact of data processing, may be indicative of a migrating source location or source mechanism. This was observed at Soufrière Hills volcano, Montserrat in July 2003 (Fig. 4; Salvage and Neuberg (2016)). The largest dome collapse to date observed at this volcano occurred on 12 July 2003, with removal of $2\times 10^6 \text{ m}^3$ of material (Herd et al. 2005), following a 4 day period from 8 to 12 July 2003

of heightened seismicity at Soufrière Hills. A migrating source mechanism can be identified from changing amplitudes in seismic events within the same family, however this characteristic cannot be identified from analysing cross correlation coefficients alone. The amplitudes of seismic events within the single family identified in July 2003 were relatively constant (Ottemöller 2008), suggesting the changing cross correlation coefficient is a consequence of a migrating source location at depth, rather than an evolving source mechanism. The generation of families ceased immediately prior to the dome collapse event, and no similar earthquakes were detected after the collapse (Fig. 4). This suggests that the physical conditions required for the generation of families were not met in the hours before, and after, the collapse event.

The analysis of families of seismicity in the time domain may also allow for the identification of spatial patterns in seismicity. Families detected at Chiles-Cerro Negro, a volcano within the Northern Andes on the border between Ecuador and Colombia in October 2014, suggests distinct temporal and spatial patterns of seismicity. The last eruption of the volcanic complex of Chiles-Cerro Negro is believed to have been
<figure><img src="image_78.png" /><figcaption>**Fig. 5** Cross Correlation Matrix of events identified using a simple amplitude based detection algorithm at Chiles-Cerro Negro on 20 October 2014, at a single station. A total of 597 events were identified and are shown in temporal order along the x and y axis. Events with a cross correlation coefficient of greater than 0.7 are shown on a colour scale, with those close to one (red)</figcaption></figure>

being more similar. The autocorrelation of each event with itself is shown in dark red along the diagonal and is equal to a cross correlation coefficient of 1. Distinct clusters of similar events can be identified, thought to suggest a temporal evolution in the dominant similar seismicity

3400 years ago, although seismicity has since been detected in the area, thought to be related to an active hydrothermal system (Ruiz et al. 2013). The volcanic complex is dissected by a large fault system, believed to have been active as recently as 1868, when two large seismic events occurred (Mw 6.6 and 7.2) (Beauval et al. 2010). Seismic activity increased in 1991 and then again in July 2013. However, a dramatic increase in seismicity from less than 50 events a day to over 150 events occurred in October 2014, concentrated beneath the summit of Chiles volcano at depths of less than 10 km (Ruiz et al. 2013). Although originally not a target volcano for the VUELCO project, the volcanic complex of Chiles-Cerro Negro is an excellent example of a re-awakening volcano, having shown no signs of magmatic unrest in recent history until the events of 2014. Analysis of seismicity identified in

October 2014 suggested not only the occurrence of families, but also their occurrence in distinct temporal patterns. The clustering of similar seismic events around the diagonal in a number of box-like formations within a similarity matrix suggests that a number of sources were active for discrete periods of time generating families of seismicity (Fig. 5). The distinct clusters of similar seismic waveforms may relate to a changing source location or mechanism at depth (Salvage 2015). Over this time period, no significant changes in the amplitude of events (indicative of a changing source mechanism) were identified. Since similar seismicity is thought to be generated through a similar source mechanism and a similar source location, distinct cluster of similar seismicity is most likely related to its own distinct spatial region, which generated seismicity during distinct periods of time.
## The Source Mechanisms of Low Frequency Earthquakes

Since low frequency events, and in particular families, appear to be important in detecting changes in unrest at volcanoes, it is important to understand their mechanism of generation. Synthetic modelling and moment tensor inversions of low frequency seismic wavefields are powerful tools for gaining information on the source mechanisms underlying volcanic earthquakes. Once instrument response and path effects have been accounted for, real data can be compared to synthetic models, and on the basis of a best-fit approach the obtained model parameters allow insights into the nature and geometry of the source (Chouet 1996; Shuler et al. 2013). However, as the fundamental assumptions behind the commonly used moment tensor inversions are based on plane surface geometries which are believed to be too simple to explain the generation of low frequency events in a volcanic environment, the application to more complex seismic sources has so far been inconclusive. In the framework of the VUELCO project, slip along bent surfaces (a complex source) was thought of as the underlying physical motion responsible for generating seismic energy. This novel way of investigating low frequency earthquakes can explain several features of the earthquakes under investigation without introducing compromising

assumptions such as slip along a single, unbent surface, which is believed to be unrealistic.

An example of a spatially extended source generating seismicity within a volcanic environment is a volcanic conduit through which magmatic fluids move. In these instances, the generation of low frequency seismicity may be related to the brittle failure of magma itself (Neuberg et al. 2006; Lavallée et al. 2008; Thomas and Neuberg 2012) or through a stick-slip motion at the conduit edge (Iverson et al. 2006). In either case, shallow source depths (1–2 km) and short epicentral distances to seismic receivers (a few kilometers) suggest that a spatially extended source is more realistic than a single point source.

The occurrence of slip (i.e. the generation of the seismic energy itself) of spatially extended sources may either be instantaneous along two or more slip surfaces, or may occur on different slip surfaces at different times, offset by a given time increment, delta *t*. A ring fault structure is a numerical description of seismogenic slip of magma along all of the conduit walls within a volcanic edifice, and can be numerically modelled by considering a cylinder representing the volcanic conduit with instantaneously slipping double couple (single point) sources bounding the circumference (Fig. 6). Upward movement inside the cylinder and downward movement outside represents the movement of magma

<figure><img src="image_79.png" /><figcaption>**Fig. 6** Schematic representation of a ring fault structure with movement directed upwards within the cylinder to represent the flow of magma through a conduit. Each</figcaption></figure>

planar surface is host to a single double couple source i.e. a point source (labelled 1–8)
through this channel. Other spatially extended sources which may evoke the generation of low frequency seismicity include: the upward movement of magma through a narrow dyke, numerically modelled by two oppositely directed double couple sources; slip along distinct segments of the volcanic conduit i.e. due to geometry changes (Thomas and Neuberg 2012), numerically modelled with double couple sources on distinct segments; or the generation of a number of seismic swarms of families of earthquakes occurring simultaneously (Salvage and Neuberg 2016), which can be numerically modelled by movement on two or more simultaneously acting ring fault structures.

Green and Neuberg (2006) suggested that accelerated magma movement at Soufrière Hills volcano can be linked to observed deformation cycles through low frequency seismic swarms, and that seismicity is only generated if significant magma movement takes place. Considering slip through extended sources brings us one step closer to estimating magma ascent rates. Once calibrated, the link between observed waveform amplitudes and the amount of seismogenic slip occurring during a seismic swarm will yield magma ascent rates and will ultimately contribute to forecasting volcanic eruptions more accurately. Point and extended source models yield great differences in observed P-wave amplitudes and waveforms, leading to remarkable differences when interpreting the amount of

slip occurring and the slip rates. In the case of slip along a ring fault, P-wave amplitudes are greatly reduced due to destructive interference, in comparison to simple double couple (point) sources. As a result, observed amplitudes yield an underestimation of actual slip by more than a factor of 3 if interpreted as a point source, when in reality a spatially extended source acts at depth. This underestimation in seismic moment consequently may lead to an underestimation of magma flow rate at depth, which in turn has severe implications for eruption forecasting (Karl 2014).

Furthermore, the P-wave radiation pattern for a spatially extended ring fault structure (i.e. for modelling the movement of magma within the entire conduit made up of a number of double couple sources) shows remarkable similarity to the P-wave radiation pattern for a compensated linear vector dipole (CLVD) source, which instead is a conservation of energy (Fig. 7). The derived ring fault radiation pattern shows rotational symmetry around the depth axis, and consists of a large compressional lobe directly above the source and an inversely polarised, dilatational lobe with the same amplitude below it. Only if the small radial P-wave radiation lobes can be determined for the ring fault structure (Fig. 8) is it possible to distinguish between these two source mechanisms. Seismic networks with small apertures typical in volcanic settings and

<figure><img src="image_80.png" /><figcaption>**Fig. 8** Zoom of smaller radial lobes of P-wave radiation pattern for a spatially extended ring fault structure. Compressional lobes are *red*, dilatational lobes are *blue*</figcaption></figure>

<figure><img src="image_81.png" /><figcaption>**Fig. 7** P-wave radiation pattern generated for a CLVD source, and for an extended ring fault structure source. The *red* lobes are compressional, the *blue* lobes are dilatational. In both radiation patterns a large compressional lobe is found above the source, with dilatational lobes extending radially which can lead to confusion for interpretation of first motion polarity patterns</figcaption></figure>
uncertainties in earthquake source depth locations will likely lead to difficulties in distinguishing between the two radiation patterns, since both can explain observed first motion polarity patterns of low frequency seismicity on volcanoes (Karl 2014).

## Forecasting Eruptive Activity

The ability to forecast the timing, intensity and type of volcanic activity is one of the key issues facing volcanologists today. The most notable instances of successful volcanic forecasting use precursory activity at andesitic-dacitic volcanoes. The cataclysmic eruption of Mt. Pinatubo, Philippines on 15 June 1991 was preceded by at least two months of heightened seismicity (Harrow et al. 1996). With increases in seismicity and an alarming sudden drop in SO₂, scientists were able to successfully evacuate over 45,000 local people and 14,500 military personnel to safety by 14 June, such that less than 300 people were killed in the ensuing volcanic activity on 15 June. More recently, the 2010 eruption of Merapi, Indonesia, on 26 September was preceded by approximately 6 weeks of precursory activity (Budi-Santoso et al. 2013): rates of seismicity and SO₂ during this time were comparable to, or higher than, the highest rates observed during previous (smaller) Merapi eruptions (1992–2007), and rapid deformation was observed. Consequently, one day prior to the explosive eruption, several tens of thousands of people were evacuated from a radius extending 10 km from the volcano, resulting in a greatly lowered death toll of 35.

Volcanic eruptions are often preceded by accelerating geophysical signals, associated with the movement of magma or other fluid towards the surface. Of these precursors, seismicity is at the forefront of forecasting volcanic activity since it is frequently observed and the change from background level can be observed in real time. Since forecasting of volcanic eruptions relies on the ability to determine the timing of magma reaching the surface, low frequency seismicity may act as a forecasting tool due to its

potential correlation with the movement of magmatic fluid at depth.

The Failure Forecast Method (FFM) is based on an empirical power-law relationship, which relates the acceleration of a precursor ($d^2\Omega/dt^2$) to the rate of that precursor ($d\Omega/dt$) (Voight 1988) method by:

$$ \frac{d^2\Omega}{dt^2} = K \left(\frac{d\Omega}{dt}\right)^\alpha \quad (1) $$

where K and α are empirical constants. Ω can represent a number of different geophysical precursors, for example low frequency seismic event rate (Salvage and Neuberg 2016), event rate of all recorded seismicity (Kilburn and Voight 1998), or the amplitude of seismic events (Ortiz et al. 2003). The parameter α is thought to range between 1 and 2 in volcanic environments (Voight 1988), or may even evolve from 1 towards 2 as seismicity proceeds (Kilburn 2003). α has also been calculated in hindsight as high as 3.3 for accelerating seismicity in 1991 at Mt. Pinatubo, although this extreme value appears rare and was calculated with only a small amount of seismic data (Smith and Kilburn 2010). An infinite $d\Omega/dt$ suggests an uncontrolled rate of change (a singularity) and in this environment is associated with an impending eruption. The inverse form of $d\Omega/dt$ is linear if α = 2, and therefore in this case the solution for the timing of failure is a linear regression of inverse rate against time, with the timing of failure relating to the point where the linear regression intersects the x-axis (Voight 1988).

Although assuming that α = 2 is the simplest method to estimate the timing of an eruption through a linear regression and therefore the most common application of the FFM in hindsight analysis, some authors have suggested that it may not be an appropriate assumption for use with the FFM (e.g. Bell et al. 2011). Additionally, some authors have argued that α may evolve with time as precursory sequences develop, which is not detailed in the FFM (Kilburn 2003). As the FFM follows a least squares regression analysis when α is equal to 2, the residual error between the observed event rate and the mean
event rate of seismicity should follow a typical Gaussian distribution (Bell et al. 2011). Greenhough and Main (2008) have suggested that since earthquake occurrence is a point process, the rate uncertainties are best described by a Poisson distribution. In this instance, a generalised linear model (GLM) where $\alpha = 1$, rather than a least squares regression model ($\alpha = 2$) may be more appropriate, since it can allow for a distribution of data that is non-Gaussian (Bell et al. 2011). At Soufrière Hills volcano, however, the use a GLM to forecast the timing of eruptive events in 1997 and 2003 failed to generate an appropriate forecast (Salvage and Neuberg 2016). Hammer and Ohrnberger (2012) suggested that this may be related to the fact that a Poisson process, and therefore the GLM, is a memoryless system, meaning that past events do not influence future patterns. A memoryless system is not consistent with the fundamental assumptions of the FFM, since previous geophysical observables form the basis of such a forecast.

One of the first instances of real time forecasting using the FFM was at Redoubt volcano, when the inverse average amplitude of seismic

events followed a linear regression trend for 4 days prior to a dome collapse event on 2 January 1990. Due to this trend, and the fact that the seismic intensity was far above background levels, the Alaskan Volcano Observatory issued a “formal warning” of an impending eruptive event on the morning of the 2 January, a few hours before the eruption began, although the FFM calculations suggested an eruption was likely within 0.5–2 days. A similar, if not clearer trend, that supported the forecast was found using the same precursory sequence but only using seismic events within the spectral range of 1.3–1.9 Hz (Cornelius and Voight 1994), suggesting an increased accuracy in forecasts when focusing on a single source process at depth.

Swarms of seismic events, i.e. a number of similar events within a short period of time, with typical swarm durations of hours to day, are not observed at all volcanoes, but have been commonly observed at Soufrière Hills volcano (e.g. Green and Neuberg 2006) and Redoubt volcano (e.g. Buurman et al. 2013). Using precursory seismicity and the FFM, Salvage and Neuberg (2016) forecast in hindsight the timing of a dome collapse event on 25 June 1997 at Soufrière

<figure><img src="image_82.png" /></figure>

**Fig. 9** Application of the FFM: the inverse average event rate per 10 min within swarms from 22 to 25 June 1997 at a single station at Soufrière Hills volcano. Each data point represents the inverse average event rate for each individual identified swarm of seismicity. The *vertical line* represents the known timing of dome collapse on the 25 June 1997 at 16:55 UTC. The

<figure><img src="image_83.png" /></figure>

graphical representation of the FFM is depicted by the linear regression (it is assumed that $\alpha = 2$ for simplicity) and the forecasted timing of failure can be read off the x-axis at the point where the linear regression crosses it.
**a** All triggered low frequency seismicity.
**b** Single family of similar seismicity.
Hills, based upon accelerating rates of low frequency earthquakes which occurred in swarms, rather than simply the number of low frequency events over the precursory time period (Fig. 9). Using the average event rate per swarm showed a clearer accelerating pattern over the entire seismic sequence, rather than using the traditional method of binning data in units of time. More accurate forecasts were determined when using only one single family of similar events to forecast the dome collapse, rather than all low frequency seismicity mixed together. A dome collapse on 12 July 2003 at Soufrière Hills volcano was also more accurately forecast when using a single family of similar events, rather than all low frequency seismicity, which occurred during the period of unrest (Salvage and Neuberg 2016). Consequently, the use of families of seismicity, and therefore concentration upon a single active system at depth, may allow a more accurate forecast of the timing of an eruptive event in these instances.

## Summary

Seismology is a powerful tool which can be used to understand processes occurring at depth, and their relationship to the surface at volcanoes, especially since seismic events are easily detected and the deviation from the background level is often notably pronounced. An increase in seismicity may be an indication of volcanic unrest, since it relates to a number of physical processes at depth including brittle failure and fracturing of the conduit or of the surrounding edifice (high frequency seismicity) and the movement of magmatic fluids at depth (low frequency seismicity). Using families of seismicity as an indicator of a single active system at depth, temporal and spatial patterns in the seismic events can be used to assess the potential migration of seismicity towards the surface, as well as to forecast the timing of volcanic eruptive events related to the acceleration of seismicity. However, a full understanding of the source mechanism of the generated seismic events is

essential to ensure that the magma flow rate is estimated accurately and therefore an accurate forecast can be generated for the timing of eruption. Using seismicity in combination with other monitoring tools, we are now closer to gaining a better understanding of evolving magmatic systems at depth.

**Acknowledgements** The past and present staff at the Montserrat Volcano observatory are fully acknowledged for their ongoing support in the upkeep and maintenance of the seismic network, and the sharing of data. All staff at the Instituto-Geofisico in Ecuador are also fully acknowledged for providing data from Chiles-Cerro Negro volcano, for useful discussions regarding the volcanic complex, and for their continued monitoring efforts of all volcanoes in Ecuador. We thank two anonymous reviewers for their detailed comments and suggestions which greatly enhanced the quality of this manuscript. Muchas gracias también a María Martínez-Cruz y Javier Pacheco por su ayuda con la traducción al español.

## References

Arámbula-Mendoza R, Lesage P, Valdés-González C, Varley N, Reyes-Dávila G, Navarro C (2011) Seismic activity that accompanied the effusive and explosive eruptions during the 2004–2005 period at Volcán de Colima, Mexico. *J Volcanol Geoth Res* 205(1):30–46

Arciniega-Ceballos A, Chouet B, Dawson P (2003) Long-period events and tremor at Popocatepetl volcano (1994–2000) and their broadband characteristics. *Bull Volc* 65(2):124–135

Bean CJ, De Barros L, Lokmer I, Métaxian J-P, O’Brien G, Murphy S (2014) Long-period seismicity in the shallow volcanic edifice formed from slow-rupture earthquakes. *Nat Geosci* 7(1):71–75

Beauval C, Yepes H, Bakun WH, Egred J, Alvarado A, Singaucho J-C (2010) Locations and magnitudes of historical earthquakes in the Sierra of Ecuador (1587–1996). *Geophys J Int* 181(3):1613–1633

Bell A, Naylor M, Heap M, Main I (2011) Forecasting volcanic eruptions and other material failure phenomena: an evaluation of the failure forecast method. *Geophys Res Lett* 38:L15304

Budi-Santoso A, Lesage P, Dwiyono S, Sumarti S, Jousset P, Metaxian J-P et al (2013) Analysis of the seismic activity associated with the 2010 eruption of Merapi Volcano, Java. *J Volcanol Geoth Res* 261:153–170

Budi-Santoso A, Lesage P (2016) Velocity variations associated with the large 2010 eruption of Merapi volcano, Java, retrieved from seismic multiplets and ambient noise cross-correlation. *Geophys J Int* 206 (1):221–240
Buurman H, West ME, Thompson G (2013) The seismicity of the 2009 Redoubt eruption. *J Volcanol Geoth Res* 259:16–30

Chouet B (1988) Resonance of a fluid-driven crack: radiation properties and implications for the source of long-period events and harmonic tremor. *J Geophys Res Solid Earth* (1978–2012) 93(B5), 4375–4400

Chouet BA (1996) New methods and future trends in seismological volcano monitoring. In: Monitoring and Mitigation of Volcano Hazards. Springer, pp. 23–97

Chouet BA, Matoza RS (2013) A multi-decadal view of seismic methods for detecting precursors of magma movement and eruption. *J Volcanol Geoth Res* 252:108–175

Cornelius R, Voight B (1994) Seismological aspects of the 1989–1990 eruption at Redoubt Volcano, Alaska: the materials Failure Forecast Method (FFM) with RSAM and SSAM seismic data. *J Volcanol Geoth Res* 62(1):469–498

Dawson PB, Chouet BA, Power J (2011) Determining the seismic source mechanism and location for an explosive eruption with limited observational data: Augustine Volcano, Alaska. *Geophys Res Lett* 38(3):L03302

De Angelis S, Henton S (2011) On the feasibility of magma fracture within volcanic conduits: constraints from earthquake data and empirical modelling of magma viscosity. *Geophys Res Lett* 38(19):L19310

De Angelis S, Bass V, Hards V, Ryan G (2007) Seismic characterization of pyroclastic flow activity at Soufrière Hills Volcano, Montserrat, 8 January 2007. *Nat Hazards Earth Syst Sci* 7:467–472

Geller R, Mueller C (1980) Four similar earthquakes in central California. *Geophys Res Lett* 7(10):821–824

Goto A (1999) A new model for volcanic earthquake at unzen volcano: melt rupture model. *Geophys Res Lett* 26(16):2541–2544

Green D, Neuberg J (2006) Waveform classification of volcanic low-frequency earthquake swarms and its implication at Soufrière Hills Volcano, Montserrat. *J Volcanol Geoth Res* 153(1):51–63

Greenhough J, Main I (2008) A poisson model for earthquake frequency uncertainties in seismic hazard analysis. *Geophys Res Lett* 35(19):L19313

Hammer C, Ohrnberger M (2012) Forecasting seismo-volcanic activity by using the dynamical behavior of volcanic earthquake rates. *J Volcanol Geoth Res* 229:34–43

Harlow DH, Power JA, Laguerta EP, Ambubuyog G, White RA, Hoblitt RP (1996) Precursory seismicity and forecasting of the June 15, 1991, eruption of Mount Pinatubo. Eruptions and lahars of Mount Pinatubo, Philippines, Fire and mud, pp 223–247

Herd RA, Edmonds M, Bass VA (2005) Catastrophic lava dome failure at Soufrière Hills volcano, Montserrat, 12–13 July 2003. *J Volcanol Geoth Res* 148(3):234–252

Iverson R, Dzurisin D, Gardner C, Gerlach T, LaHusen R, Lisowski M, Major J, Malone S, Messerich J, Moran S et al (2006) Dynamics of seismogenic volcanic extrusion at Mount St Helens in 2004–05. *Nature* 444(7118):439–443

Karl S (2014) The source mechanisms of low frequency seismic events on volcanoes. Ph.D. thesis, University of Leeds, UK, Online at: http://etheses.whiterose.ac.uk/id/eprint/8406

Kilburn C (2003) Multiscale fracturing as a key to forecasting volcanic eruptions. *J Volcanol Geoth Res* 125(3–4):271–289

Kilburn CR, Voight B (1998) Slow rock fracture as eruption precursor at Soufrière Hills volcano Montserrat. *Geophys Res Lett* 25(19):3665–3668

Lahr J, Chouet B, Stephens C, Power J, Page R (1994) Earthquake classification, location, and error analysis in a volcanic environment: implications for the magmatic system of the 1989–1990 eruptions at Redoubt Volcano, Alaska. *J Volcanol Geoth Res* 62(1):137–151

Lavallée Y, Meredith P, Dingwell D, Hess K-U, Wassermann J, Cordonnier B, Gerik A, Kruhl J (2008) Seismogenic lavas and explosive eruption forecasting. *Nature* 453(7194):507–510

McNutt SR (2005) Volcanic seismology. Annu Rev Earth Planet Sci 32:461–491

Nakano M, Kumagai H (2005) Response of a hydrothermal system to magmatic heat inferred from temporal variations in the complex frequencies of long-period events at Kusatsu-Shirane Volcano, Japan. *J Volcanol Geoth Res* 147(3):233–244

Neuberg J, Luckett R, Baptie B, Olsen K (2000) Models of tremor and low-frequency earthquake swarms on Montserrat. *J Volcanol Geoth Res* 101(1–2):83–104

Neuberg J, Tuffen H, Collier L, Green D, Powell T, Dingwell D (2006) The trigger mechanism of low-frequency earthquakes on Montserrat. *J Volcanol Geoth Res* 153(1):37–50

Ortiz R, Moreno H, García A, Fuentealba G, Astiz M, Peña P, Sánchez N, Tárraga M (2003) Villarrica volcano (Chile): characteristics of the volcanic tremor and forecasting of small explosions by means of a material failure method. *J Volcanol Geoth Res* 128(1):247–259

Ottemöller L (2008) Seismic hybrid swarm precursory to a major lava dome collapse: 9–12 July 2003, Soufrière Hills Volcano, Montserrat. *J Volcanol Geoth Res* 177(4):903–910

Petersen T (2007) Swarms of repeating long-period earthquakes at Shishaldin Volcano, Alaska, 2001–2004. *J Volcanol Geoth Res* 166(3):177–192

Phillipson G, Sobradelo R, Gottsmann J (2013) Global volcanic unrest in the 21st century: an analysis of the first decade. *J Volcanol Geoth Res* 264:183–196

Ruiz G, Cordova A, Ruiz M, Alvarado A (2013) Informe Tecnico de los volcanes Cerro Negro y Chiles. Tech. Rep., IG-EPN, in Spanish

Salvage RO (2015) Using seismic signals to forecast volcanic processes. Ph.D. thesis, University of Leeds, UK, Online at: http://etheses.whiterose.ac.uk/12268/. Restricted until April 2018

Salvage R, Neuberg J (2016) Using a cross correlation technique to refine the accuracy of the failure forecast method: application to Soufrière Hills volcano, Montserrat. *J Volcanol Geoth Res* 324:118–133
Sgattoni G, Jeddi Z, Gudmundsson Ó, Einarsson P, Tryggvason A, Lund B, Lucchi F (2016) Long-period seismic events with strikingly regular temporal patterns on Katla volcano's south flank (Iceland). J Volcanol Geoth Res 324:28-40

Shuler A, Ekström G, Nettles M (2013) Physical mechanisms for vertical-CLVD earthquakes at active volcanoes. J Geophy Res Solid Earth 118(4):1569-1586

Smith R, Kilburn C (2010) Forecasting eruptions after long repose intervals from accelerating rates of rock fracture: the June 1991 eruption of Mount Pinatubo, Philippines. J Volcanol Geoth Res 191(1):129-136

Sparks R, Young S (2002) The eruption of Soufrière Hills Volcano, Montserrat (1995-1999): overview of scientific results. Geol Soc Lond Mem 21(1):45-69

Thelen W, Malone S, West M (2011) Multiplets: their behavior and utility at dacitic and andesitic volcanic centers. J Geophys Res Solid Earth (1978-2012) 116 (B8):B08210

Thomas ME, Neuberg J (2012) What makes a volcano tick—a first explanation of deep multiple seismic sources in ascending magma. Geology 40(4):351-354

Voight B (1988) A method for prediction of volcanic eruptions. Nature 332:125-130

Wadge G, Voight B, Sparks R, Cole P, Loughlin S, Robertson R (2014) An overview of the eruption of Soufrière Hills Volcano, Montserrat from 2000 to 2010. Geological Society, London, Memoirs 39(1):1-40

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.
# The Ups and Downs of Volcanic Unrest: Insights from Integrated Geodesy and Numerical Modelling

J. Hickey, J. Gottsmann, P. Mothes, H. Odbert, I. Prutkin and P. Vajda

## Abstract

Volcanic eruptions are often preceded by small changes in the shape of the volcano. Such volcanic deformation may be measured using precise surveying techniques and analysed to better understand volcanic processes. Complicating the matter is the fact that deformation events (e.g., inflation or deflation) may result from magmatic, non-magmatic or mixed/hybrid sources. Using spatial and temporal patterns in volcanic deformation data and mathematical models it is possible to infer the location and strength of the subsurface driving mechanism. This can provide essential information to inform hazard assessment, risk mitigation and eruption forecasting. However, most generic models over-simplify their representation of the crustal conditions in which the deformation source resides. We present work from a selection of studies that employ advanced numerical models to interpret deformation and gravity data. These incorporate crustal heterogeneity, topography, viscoelastic rheology and the influence of temperature, to constrain unrest source parameters at Uturuncu (Bolivia), Cotopaxi (Ecuador), Soufrière Hills (Montserrat), and Teide (Tenerife) volcanoes. Such model complexities are justified by

J. Hickey and H. Odbert, previously at School of Earth Sciences, University of Bristol, UK.

J. Hickey (✉)
Camborne School of Mines, University of Exeter,
Exeter, UK
e-mail: J.Hickey@exeter.ac.uk

J. Gottsmann
School of Earth Sciences, University of Bristol,
Bristol, UK

J. Gottsmann
Cabot Institute, University of Bristol, Bristol, UK

P. Mothes
Instituto Geofisico, Escuela Politecnica Nacional,
Quito, Ecuador

H. Odbert
Met Office, Exeter, UK

I. Prutkin
Institute of Geosciences, Jena University, Jena,
Germany

P. Vajda
Earth Science Institute, Slovak Academy of
Sciences, Bratislava, Slovakia
geophysical, geological, and petrological constraints. Results highlight how more realistic crustal mechanical conditions alter the way stress and strain are partitioned in the subsurface. This impacts inferred source locations and magmatic pressures, and demonstrates how generic models may produce misleading interpretations due to their simplified assumptions. Further model results are used to infer quantitative and qualitative estimates of magma supply rate and mechanism, respectively. The simultaneous inclusion of gravity data alongside deformation measurements may additionally allow the magmatic or non-magmatic nature of the source to be characterised. Together, these results highlight how models with more realistic, and geophysically consistent, components can improve our understanding of the mechanical processes affecting volcanic unrest and geodetic eruption precursors, to aid eruption forecasting, hazard assessment and risk mitigation.

Extended Spanish Abstract

La deformación volcánica, caracterizada por pequeños cambios medibles en la morfología del volcán, a menudo, pero no siempre, precede a una erupción volcánica. Esta cuestión, sin embargo, se complica por el hecho de que los eventos de deformación (por ejemplo, inflación o deflación) pueden ser el resultado de una fuente magmática, no magmática o de fuentes mixtas/híbridas. Utilizando tanto la amplitud y patrones espacio-temporales de datos de deformación volcánica registrados, así como la utilización de modelos matemáticos, es posible inferir la ubicación y la fuerza de la fuente impulsora subyacente. Estos métodos pueden proporcionar información esencial para la evaluación y mitigación de riesgos, así como para el pronóstico de erupciones volcánicas. Sin embargo, la mayoría de los modelos genéricos son insatisfactorios en su representación de las condiciones de la corteza en las que reside la fuente. En este trabajo presentamos una selección de estudios que emplean modelos numéricos avanzados para la interpretación de datos de deformación y gravedad. Dichos datos incorporan la heterogeneidad de la corteza, la topografía, la reología inelástica y los efectos termo-mecánicos para constreñir los parámetros asociados a la fuente de perturbación en cuatro sistemas volcánicos. Las complejidades de estos modelos están justificadas por limitaciones geofísicas, geológicas y petrológicas. El estudio realizado en el volcán Uturuncu, localizado en Bolivia, destaca la importancia de la estructura sub-superficial y de los procesos dependientes de tiempo en la fuente para explicar los patrones de deformación espacial-temporal. La combinación de dichos resultados indica un ascenso del magma de tipo diapírico. En el volcán Cotopaxi, localizado en Ecuador, los nuevos modelos de inversión que emplean el Análisis por Elementos Finitos esclarecen la ubicación y el volumen de una intrusión magmática durante un episodio de actividad asísmico y no eruptivo con una baja tasa de suministro de magma. Estos modelos también proporcionan señales observables que podrían estar asociadas con futura actividad volcánica intrusiva o eruptiva. El análisis de la deformación intra-eruptiva en el volcán Soufrière Hills, en Montserrat,
mostró cómo la inflación registrada podría deberse a una serie de reservorios magmáticos apilados o un depósito alargado verticalmente, debido a respuestas termo-mecánicas de la corteza similares. Las condiciones de falla derivadas para las tasas de suministro de magma y reservorio son consistentes con las restricciones térmicas y mecánicas independientes. Utilizando datos gravimétricos y la curiosa falta de deformación asociada, el estudio en el volcán Teide, en Tenerife, se identifican, de forma separada, las contribuciones de fuentes poco profundas como de aquellas con más profundidad. La interpretación de ambos elementos sugiere que una intrusión magmática profunda activó un sistema hidrotermal superficial localizado encima de ésta, y por lo tanto demuestra un efecto causal de un origen mixto. Estos casos de estudio representan una contribución significativa para la comprensión de los procesos volcánicos durante los periodos de actividad intra-eruptiva y no eruptiva. La combinación de estos resultados enfatiza cómo una mejor parametrización de condiciones mecánicas de la corteza puede alterar fundamentalmente la forma en que el estrés y la tensión se reparten en la sub-superficie. Del mismo modo, una topografía compleja, como es el caso en los estratovolcanes con laderas empinadas, puede afectar la partición de la deformación superficial. Estos dos efectos impactan la inferencia de la localización de las fuentes y presiones magmáticas previo al fallo y la erupción, e indican cómo los modelos genéricos pueden conducir a interpretaciones engañosas debido a la simplificación de inferencias acerca de la corteza y a espacios planos y a medias. Resultados adicionales de la modelización son utilizados para inferir estimaciones cuantitativas y cualitativas de la tasa de suministro de magma y el mecanismo, respectivamente, contribuyendo así al entendimiento de la dinámica del transporte del magma. Además, la inclusión simultánea de datos de gravedad junto a mediciones de deformación, permiten la caracterización de la naturaleza magmática o no magmática de la fuente. Juntos, estos estudios destacan cómo los modelos que cuentan con componentes más plausibles y geofísicamente consistentes, pueden mejorar nuestro entendimiento de los procesos mecánicos que afectan la reactivación volcánica y de precursores geodésicos de erupciones. Estos también proporcionan un marco para ayudarnos a avanzar en el pronóstico de una erupción, así como en la evaluación y mitigación de los riesgos, proporcionando datos cuantitativos derivados de la modelización de mecanismos físicos adecuados y robustos.

**Keywords**

Volcano deformation • Gravity • Modelling • Crustal mechanics • Geodesy

**Palabras clave**

deformación volcánica • gravedad • modelización • mecánica de la corteza • geodesia
**Introduction**

Volcano deformation is a key observable during periods of volcanic unrest, and one of the main tools used to monitor developing crises (Sparks et al. 2012). Non-magmatic causes of volcanic deformation during an unrest period do not involve movement of new magma and are most commonly related to active hydrothermal systems (e.g., Fournier and Chardot 2012; Rouwet et al. 2014, and references therein). Magmatic causes reflect the active migration and accumulation of new magma, with usually deeper origins and wider deformation footprints. Hybrid mechanisms have also been proposed, where contributions from magmatic and hydrothermal systems combine to produce a single complex deformation pattern (Gottsmann et al. 2006a).

A variety of both ground and satellite based geodetic monitoring techniques are used to assess the spatial and temporal evolution of volcanic related surface deformation. The two most common techniques employed today are the Global Positioning System (GPS) and Interferometric Synthetic Aperture Radar (InSAR), which track changes in position due to the ground deforming and complement each other with advantages in their temporal and spatial coverage, respectively. Volcano gravimetry studies, the monitoring and interpretation of continuous or time-lapse spatio-temporal gravity variations, can also provide additional information on any density changes associated with the driving mechanism (Battaglia et al. 2008), and thus enable an estimate of its nature (e.g., magmatic or hydrothermal), which is particularly important in cases where there is no significant (observable) surface deformation.

The magnitude, spatial pattern, and temporal evolution of volcanic deformation can be used to infer the location and ‘strength’ of a causative subsurface source. These source parameters have important implications for volcanic hazards and risk mitigation. Geodetic data are also a key component in eruption forecasting efforts (e.g.,

Sparks et al. 2012), and constraints on source parameters (e.g., pressure or volume changes) from previous deformation episodes can help to quantitatively improve these forecasts. However, making the transition from surface deformation observations to subsurface source processes requires the use of geodetic models, which demand assumptions about crustal mechanics.

Generic analytical models of volcanic deformation (e.g., Mogi 1958) are often oversimplified and do not capture the intrinsic complexities of subsurface systems, which are essential for reliable assessment of causative processes, and thus eruption forecasting and hazard assessment. They usually represent the Earth’s crust as a homogeneous, isotropic, elastic, half-space with a flat, free surface. These assumptions can cause misleading interpretations when compared to numerical models that allow for more realistic crustal conditions (Masterlark 2007). Numerical models of volcanic deformation are thus becoming increasingly popular due to this ability to estimate source parameters in crustal conditions that are beyond the analytical realm. These improvements are synonymous with the recent advances in geodetic monitoring, which deserve a more in-depth analysis than generic analytical models can offer.

The extra complexities that numerical models can account for relate to both the source itself, as well as the crustal rocks in which it resides. The deformation source does not have to be represented as a point or cavity, but can contain its own material properties (e.g., Hickey et al. 2013; Gottsmann and Odbert 2014). Moreover, numerical models can incorporate inferences from other geophysical, geological and petrological observables and thus maintain a higher level of consistency by relaxing some of the assumptions that restrict generic analytical models. For example, this enables the inclusion of subsurface heterogeneity, topography, viscoelastic rheology, and temperature-dependent mechanics (e.g., Hickey et al. 2016, and references therein). Consequently, more robust source parameters and magma transport processes can
<table><caption>**Table 1** Summary of geodetic modelling</caption><thead><tr><th>Volcano</th><th>Unrest period</th><th>FWM</th><th>IVM</th><th>Data</th><th>Topo</th><th>CMX</th><th>TMX</th><th>Model</th><th>Ref</th></tr></thead><tbody><tr><td>Uturuncu</td><td>1992–2006</td><td>✓</td><td>–</td><td>InSAR</td><td>–</td><td>✓</td><td>–</td><td>2D</td><td>Hickey et al. (2013)</td></tr><tr><td>Cotopaxi</td><td>2001–2002</td><td>✓</td><td>✓</td><td>EDM</td><td>✓</td><td>✓</td><td>✓</td><td>3D</td><td>Hickey et al. (2015)</td></tr><tr><td>Soufrière hills</td><td>2003–2005</td><td>✓</td><td>–</td><td>cGPS</td><td>✓</td><td>✓</td><td>✓</td><td>2D</td><td>Gottsmann and Odbert (2014)</td></tr><tr><td>Central volcanic complex</td><td>2004–2005</td><td>–</td><td>✓</td><td>GPS & gravity</td><td>✓</td><td>–</td><td>–</td><td>3D</td><td>Prutkin et al. (2014)</td></tr></tbody></table>

FWM forward modelling, IVM inverse modelling, Topo topography, CMX crustal mechanics, TMX thermomechanics, 2D two-dimensional axisymmetric, 3D three-dimensional, Ref reference

be inferred, thereby improving the understanding of the links between deformation and eruption.

In this chapter we summarise the key findings from investigations of a variety of unrest episodes at a selection of the VUELCO target volcanoes; Cotopaxi (Ecuador), Soufrière Hills (Montserrat, British West indies) and the Central Volcanic Complex on Tenerife (Spain), as well as at Uturuncu volcano in Bolivia (Table 1). These examples highlight not only different timescales of unrest and spatial patterns of volcano deformation, but also the influence of crustal mechanical heterogeneity, thermal effects, and topography, on observed signals. The unrest episodes are investigated using both forward and inverse numerical modelling procedures and demonstrate the process of interpreting different geodetic data sets to constrain realistic source parameters. Thus, the purpose of this chapter is not to provide a detailed treatment of the mathematical and numerical approaches. Instead, it is intended to provide the reader with a general overview of the use of advanced numerical models to infer volcanic unrest driving mechanisms with realistic source characteristics.

## Implementing Complex Crustal Mechanics

One of the most important aspects of a deformation model is the representation of crustal mechanics. This has a fundamental control on the way in which stress and strain is distributed and

transferred through the Earth's crust, from the source to the surface. In this regard, as briefly mentioned above, generic analytical models are limited by their necessary assumptions of homogeneous and elastic conditions throughout the entire model domain. In reality, the Earth's crust is known to be layered, and volcanic regions in particular can have wide-ranging regions of stiff (high Young's Modulus) and soft (low Young's Modulus) rocks relating to the type of volcanic deposit that formed them, e.g., lava flows (stiff) compared to tuffs (soft) (Gudmundsson 2011). This is what we call subsurface heterogeneity. Where stiff and soft regions are adjacent in the crust, complex subsurface partitioning alters the way stress and strain are transferred to the surface. The outcome is a different surface deformation pattern compared to one that would be seen if the crust was in fact homogeneous. Consequently, generic analytical models restricted to homogeneous crustal mechanics will not adequately represent likely subsurface conditions, and the inferred source parameters from these models can be misleading.

Seismic studies are capable of delineating areas of relatively high and low seismic velocities. They can be used to estimate the dynamic Young's Modulus, $E_D$, Poisson's Ratio, $\nu$, and density, $\rho$, of the crust (Brocher 2005):

$$ \nu = 0.5 \times \left[ \left( \frac{V_p}{V_s} \right)^2 - 2 \right] / \left[ \left( \frac{V_p}{V_s} \right)^2 - 1 \right] \quad (1) $$
$$ \rho = 1.6612V_p - 0.4721V_p^2 + 0.067V_p^3 - 0.0043V_p^4 + 0.000106V_p^5 \quad (2) $$

$$ E_D = \frac{V_p^2 \rho (1 + v) (1 - 2v)}{(1 - v)} \quad (3) $$

where $V_p$ and $V_s$ are the primary and shear seismic velocities, respectively. The static Young's Modulus, $E$, is usually a factor of 2–9 smaller than $E_D$, and is more appropriate for deformation studies as the dominant processes taking place are significantly slower than the propagation of seismic waves (e.g., Gudmundsson 2011, and references therein). Numerical approaches can incorporate these spatially-variable calculations of Young's Modulus and Poisson's Ratio by allowing the mechanical parameterisation of the crust to vary within the model (Hickey and Gottsmann 2014). This can be achieved in one-dimension, where the values are only changed with depth, or in three-dimensions, where the mechanical representation additionally changes in the two horizontal axes. The Finite Element Method is the most common technique used when incorporating subsurface heterogeneity. With this approach, a heterogeneous model domain can be built using multiple different sized homogeneous ‘blocks’ of varying material properties (Hickey et al. 2013), with continuous interpolations as a function of depth or distance (Gottsmann and Odbert 2014), or with a three-dimensional coordinate interpolation (Hickey et al. 2016).

A further advantage of a numerical approach using the Finite Element Method is the ability to include a spatially variable inelastic rheology that is dependent on an estimated temperature distribution. Rocks do not always behave in an elastic manner (Ranalli 1995). Instead, where crustal conditions are hotter than the brittle-ductile transition, or subject to long-term loading, rocks can deform like very high viscosity fluids, and so viscoelastic effects are more likely (Ranalli 1995). This is especially important when a hot magmatic component is assumed to be present, or where elevated geothermal gradients are observed. Also, when models are restricted to

elastic mechanical behaviour it can be difficult to constrain realistic source processes, particularly if there is a complex temporal deformation pattern that can not be reproduced with the instantaneous elastic stress-strain relationship, and may be better represented by a non-instantaneous viscous stress response. A model that incorporates temperature-dependent mechanics can overcome these difficulties. This is achieved using the Finite Element Method with a two-step procedure (Hickey and Gottsmann 2014). First, a spatially-variable, stationary, temperature distribution is derived, primarily dependent on the local geothermal gradient(s) and an assumed (magmatic) source temperature. This temperature distribution, $T$, is then used to define a crustal viscosity, $\eta$, for example through an Arrhenius relationship, where:

$$ \eta = A_d \exp\left(\frac{H}{RT}\right) \quad (4) $$

and $A_d$ is the Dorn Parameter, $H$ is the activation energy, and $R$ is the universal gas constant. The viscosity is used in a viscoelastic material representation (Hickey and Gottsmann 2014), so where the viscosity is high the material response is dominated by elastic behaviour, and when the viscosity is lower the rocks show an increasingly higher tendency for viscous behaviour. To demonstrate the effect of temperature-dependent mechanics, and subsurface heterogeneity, the following case studies include examples of both.

## Case Studies

### Uturuncu

Uturuncu volcano is located in southern Bolivia, within the Altiplano-Puna Volcanic Complex. It has been steadily inflating since at least 1992 (Pritchard and Simons 2002), and combined with shallow seismicity and near-summit active fumaroles represents a volcanic system showing significant signs of unrest (Sparks et al. 2008).
Hickey et al. (2013) focused on the mechanism driving the 70 km wide region of ground uplift between 1992 and 2006. The aim was to constrain first-order source parameters that explain both the observed uplift rate of 1–2 cm/year and the large spatial deformation footprint (Pritchard and Simons 2002). Stress and strain from pressurised finite sources were solved numerically using Finite Element Analysis, accounting for both homogeneous and heterogeneous subsurface structure in elastic and viscoelastic rheologies. Crustal heterogeneity was constrained from seismic velocity data, which indicates a pervasive large low-velocity zone ~17 km below the surface. This is deduced to represent one of the world’s largest known regions of partial-melt: the Altiplano-Puna Magma Body (APMB) (Sparks et al. 2008).

The comparison between crustal heterogeneity and homogeneity highlights the significant effect of a mechanically weak source-depth layer (Fig. 1). The weak layer, with a lower Young’s Modulus, alters surface deformation patterns by accommodating more of the subsurface strain than its surrounding layers, thereby acting as a mechanical buffer. Continuous and regular time-dependent deformation, the long-lived nature of the source, and an anomalously high regional crustal heat-flux break the assumption of elastic conditions (e.g., Ranalli 1995), so a viscoelastic crustal rheology was tested, using the standard linear solid representation (e.g. Hickey and Gottsmann 2014). The elastic models could also only account for the spatial component of the observed uplift so their results were used solely to guide the parameters tested in the viscoelastic models. A range of possible source geometries were assessed, but spherical and oblate shapes were rejected on the grounds of their depth below the APMB and likely unsustainable pressurisation given the expected crustal mechanics. This left a prolate shaped source, whose minimum size was determined using maximum laboratory values for host-rock tensile strength. The final preferred model suggests that temporally-continuous pressurisation of a magma source protruding from the top of the APMB is

causing the observed spatial and temporal surface uplift, whereas the previous models could only infer that a simple-geometry source was being pressurised somewhere within the APMB (Pritchard and Simons 2002). Hence, this also demonstrates how a pressure-time function plays a first-order role in explaining time-dependent deformation.

Simultaneous work on an extended InSAR data set shows how the central uplift region at Uturuncu is surrounded by a ‘moat’ of subsidence (Fialko and Pearse 2012). To explain this observation they also constrained a model where magma rises out and up from the APMB with a diapiric-type ascent mechanism. Further evidence for this magmatic process is available through a complementary gravimetric study (del Potro et al. 2013). Therefore, this highlights how combinations of geodetic data and numerical models can not only constrain more plausible deformation source parameters, but can also infer magma transport dynamics.

## Cotopaxi

Cotopaxi is a large, glacier-clad stratovolcano situated in the Eastern Cordillera of the Ecuadorian Andes. The 1 km³ glacier presents a substantial lahar risk to people in the surrounding areas, and particularly to the 100,000 inhabitants that reside in the path of the 1877 lahar which descended the Inter-Andean Valley (Pistolesi et al. 2013).

Unrest was detected at Cotopaxi in 2001 and 2002. There was a significant increase in the amount of volcanic seismicity in the NE quadrant of the volcano, and this was originally interpreted to represent a dyke intrusion with subsequent gas-release and resonance of the crack (Molina et al. 2008). Hickey et al. (2015) revisited this unrest period, but approached it from a volcano deformation viewpoint. Analysis of an electronic distance meter (EDM) network over the 2001–2002 period also indicated an asymmetric inflation of the edifice that accompanied the recorded seismicity (Fig. 2). However, the irregular
<figure><img src="image_84.png" /><figcaption>**Fig. 1** A comparison of the effect crustal heterogeneity for the same prolate source geometry and depth. **a**, **b** The effect of a source depth soft layer on the subsurface deformation of a source. Both panels show the same source, embedded in a homogeneous (**a**) or heterogeneous domain (**b**). Colours relate to the radial displacement, and the white shape shows the exaggerated outline of the deformed source after the pressure is applied. In the heterogeneous model the source preferentially deforms</figcaption></figure>

into the softer layer (D), compared to the homogeneous medium, which exhibits a concentric deformation pattern. **c** Modelled surface displacement profiles from the homogeneous (**a**) and heterogeneous (**b**) models. The subsurface layering in **b** alters the displacement pattern produced at the surface as the soft layer modifies the subsurface strain partitioning. The *blue shaded area* represents the observed InSAR data and its estimated error bounds

acquisition in time of the EDM data prevented any systematic comparison between the two data sets. To solve for the optimum deformation source parameters, Hickey et al. (2015) implemented a novel numerical inversion procedure using Finite Element models. This is the first

volcanic deformation inversion study to explicitly account for both subsurface heterogeneity and surface topography while searching for a best-fit solution with a range of source shapes. The method works by solving for the predicted EDM deformation with an initial model and
source configuration. It then continually changes the source location and/or overpressure, within some predefined parameter limits, to minimise the misfit to the recorded EDM data (Fig. 2). After each inversion, the parameter limits (e.g., X and Y coordinates) are reduced around the previous solution to produce a set of decreasing size ‘Russian-doll-like’ parameter constraint grids that ensure a robust solution. Within this workflow, the Finite Element model geometry and mesh are automatically rebuilt, removing the need for repeated manual editing.

The inversion models converge on a shallow source beneath the SW flank. The individual best-fit model is inferred to represent a small oblate-shaped magmatic reservoir, approximately 4–5 km beneath the summit, with a volume increase of roughly $20 \times 10^6 \text{ m}^3$. A deformation source location in the SW is substantially different to the NE location proposed by Molina et al. (2008) when explaining the recorded seismicity. Despite this, when the deformation source was restricted to the NE quadrant the predicted

<figure><img src="image_85.png" /><figcaption>**Fig. 2** Deformation and modelling results from Cotopaxi volcano. a Map of the EDM network operational between 2001 and 2002 around the summit of the volcano, with the best-fit source from the Finite Element inversions indicated by a red star. The yellow squares represent the EDM base-stations and the orange triangles are the reflecting prisms. The inset map and arrow shows the location of Cotopaxi within Ecuador, South America.</figcaption></figure>

<figure><img src="image_86.png" /><figcaption>b The variation, and eventual reduction, of the misfit objective function with each iteration from the final, best-fit Finite Element inversion model. c The modelled values from the final inversion for the EDM baseline changes converge to their near-absolute values after approximately 25 iterations. d The best-fitting model (red circles) fits six out of seven of the EDM observations (black circles with error bounds)</figcaption></figure>
EDM measurements had a very poor fit to the observed data.

To clarify the difference in source location between the seismic and geodetic studies, Hickey et al. (2015) applied the best-fit source parameters from the elastic inversion models in a suite of temperature-dependent viscoelastic forward models to assess the rheology of the host-rock for a range of thermal parameters. The results indicated the most likely subsurface conditions would have promoted a large component of viscous deformation and that the deformation source in the SW would have consequently been pressurised aseismically. Fluid migration from the SW along existing NNE-SSW trending faults could have then caused the observed seismicity in the NE due to mass transport and excess pore pressures. The lack of eruption following this 2001–2002 unrest period and the aseismic nature of the event suggests that the magma supply rate for this period was low. A higher magma supply rate during a future unrest period would be more likely to produce seismicity around the reservoir, and could possibly indicate a level of unrest that signifies an increased likelihood of a forthcoming eruption.

## Soufrière Hills

Starting its latest eruptive episode in 1995 and displaying a remarkable diversity of eruptive activity, Soufrière Hills volcano (SHV) on Montserrat (British West Indies) is one volcano where a correlation between observed deformation and subsequent eruption can be directly established. Ground deformation data at SHV indicate cyclic behaviour of the andesitic magmatic system from periods of several hours to a few years (Odbert et al. 2014). This section focuses on the analysis of intra-eruptive unrest associated with island-wide ground uplift observed via a network of continuous GPS receivers between 28/07/2003 and 01/08/2005 (Odbert et al. 2014). After a major lava dome collapse in 2003, which marked the end of the second phase of dome extrusion at SHV, this

period preceded a restart of eruptive activity and renewed dome extrusion in August 2005.

Gottsmann and Odbert (2014) developed numerical models to test for the influence of temperature- and time-dependent stress evolution in a mechanically heterogeneous crust to explain the deformation data. Full details on the model setup and parameter derivation are given in Gottsmann and Odbert (2014) and not repeated here. They implemented two types of magma reservoir model. The first explored a series of pressurising, vertically-stacked reservoirs as proposed by Hautmann et al. (2010), while the second explored the time-dependent pressurisation of a single, vertically-elongated reservoir. Due to a similar temperature distribution from both the stacked and single reservoirs, and similar resultant rheological crustal properties, both suites of models provided equally good fits to the observed ground deformation. The study could hence not discriminate between pressurisation in a magmatic plumbing system consisting of either a single vertically elongated reservoir or a series of stacked reservoirs. Reservoir pressure changes between 4 and 7 MPa, for volumes between 60 and $100 \text{ km}^3$ and magma compressibility between $4 \times 10^{-11}$ and $1 \times 10^{-9} \text{ Pa}^{-1}$, provided plausible thermomechanical model parameters to explain the deformation data. The associated magma volume fluxes are between 0.015 and 0.021 $\text{km}^3$/year and match those derived from thermal modelling of active sub-volcanic systems (Annen 2009). Introducing a deep-crustal hot zone in the model, which modulates the partitioning of strain into the hotter underlying crust beneath the reservoir(s), promotes a further reduction in reservoir overpressures to values of around 1–2 MPa upon reservoir failure (Fig. 3). These pressure changes are significantly lower than those derived from models assuming a mechanically homogeneous and elastic crust. The deduced overpressures match those for sudden and rapid transcrustal reservoir activation prior to explosions at SHV from the analysis of volumetric strain data (Hautmann et al. 2014). The emerging eruption model at SHV hence involves the periodic failure of a compressible magma mush column beneath the volcano.
<figure><img src="image_87.png" /><figcaption>**Fig. 3** a The temperature distribution (in K) of the 2-D axisymmetric model domains, caused by a combination of a basal heat flux, the heat from the plumbing system (magma reservoir plus feeder pipe), and a deep-crusted hot zone. Background thermal distribution in the far field is caused exclusively by a basal heat flux of 0.09 W/m². Arrows indicate deformation of the medium by reservoir pressurisation. Substantial deformation is accommodated by hot and ductile parts of the mid and lower crust. b The model parameter space for reservoir volume, volume change, and magma compressibility (shown in log units on the *colour bar*) from the best fit model solutions to the 2003–2005 intra-eruptive deformation time series at Soufrière Hills Volcano. The overall fit quality decreases substantially toward higher compressibility (order 10⁻⁹ Pa⁻¹) and toward smaller reservoir volumes (<50 km³). The preferred parameter space is encompassed between the mapped upper limit of compressibility (ULC) and the lower limit of the magma reservoir volume (LLRV).</figcaption></figure>

## Las Cañadas

The Las Cañadas caldera (LCC) hosts the central volcanic complex (CVC) on Tenerife (Canary Islands), which includes the twin-volcanoes of Pico Viejo (PV) and Pico Teide (PT). Volcanic unrest on Tenerife started in April 2004 and ended a period of quiescence after the most recent magmatic eruption on the island in 1909. The unrest was geophysically characterised by anomalous seismic activity including a number of felt earthquakes and significant changes in the acceleration of gravity, yet, an absence of significant ground deformation (Gottsmann et al. 2006b). The initial interpretation of the gravity changes by Gottsmann et al. (2006b) pointed towards shallow (~2 km depth) fluid migration as the main source of the unrest, possibly related to a magma intrusion at greater depth. As part of the VUELCO project, Prutkin et al. (2014) revisited the microgravity data recorded between May 2004 and July 2005 and inverted the spatio-temporal residual gravity changes for the gravitational attraction of three-dimensional line segments.

The line segments are defined by 7 parameters: two end point coordinates (X, Y, Z) and a line strength (see Table 2). The latter is a proxy for the amplitude of sub-surface mass change, whereby a “weak” line segment represents a small amount of added mass compared to a “strong” line segment of greater added mass. An initial non-linear inversion yielded three line segments of different strengths at depths between 1 km a.s.l. and 2 km b.s.l. (i.e. between 1 and 4 km beneath the surface). Two of the segments were located to the NW of the PV-PT complex, while a third segment was located below the SW rim of the LCC (sources marked Sh 1–3 in Fig. 4). The locations of the segments correspond broadly to zones of heightened seismic activity and the unrest source location identified earlier by Gottsmann et al. (2006b).
<table><caption>**Table 2** UTM coordinates X (Easting in m), Y (Northing in m), and Z (elevation in m with respect to sea level) of the end points of the modelled shallow (Sh) and deep (Dp) sources represented by line segments (see Fig. 4)</caption><thead><tr><th>Line segment</th><th>X1</th><th>Y1</th><th>Z1</th><th>X2</th><th>Y2</th><th>Z2</th><th>ΔM</th></tr></thead><tbody><tr><td>Sh1</td><td>327,500</td><td>3,137,930</td><td>1280</td><td>330,330</td><td>3,137,260</td><td>600</td><td>0.31</td></tr><tr><td>Sh2</td><td>336,560</td><td>3,132,050</td><td>590</td><td>334,210</td><td>3,134,040</td><td>420</td><td>0.83</td></tr><tr><td>Sh3</td><td>335,460</td><td>3,120,730</td><td>1380</td><td>332,100</td><td>3,121,210</td><td>380</td><td>0.53</td></tr><tr><td>Dp1</td><td>334,890</td><td>3,133,910</td><td>-5440</td><td>334,630</td><td>3,134,200</td><td>-5680</td><td>7.25</td></tr><tr><td>Dp2</td><td>334,310</td><td>3,134,590</td><td>-5840</td><td>334,630</td><td>3,134,200</td><td>-5680</td><td>8.24</td></tr></tbody></table>

Mass additions (ΔM) are given in $10^{10}$ kg. Data Prutkin et al. (2014)

<figure><img src="image_88.png" /></figure>

**Fig. 4** Location of line segments from gravity data inversion. The figure shows the surface projections of three shallow sources (line segments Sh1–3 in *turquoise*) and the two deep sources (line segments Dp1 and 2 in *red*) superimposed over a Google Earth image of the Central Volcanic complex of Tenerife island (Spain). The source locations are derived by a nonlinear inversion of spatiotemporal residual gravity changes shown by coloured contours (in µGal) observed between 2004 and 2005 (see Prutkin et al. (2014) for details on inversion routine). The thickness of a line segment is indicative of the “strength” of the deduced mass change; i.e., the

amount of added mass. The mass added to the deep sources is more than 10 times higher than to the shallow sources. Line segments Sh1–3 are interpreted to represent near-surface sources in the NW and SW part of the Pico Teide (PT; 3718 m a.s.l.) and Pico Viejo (PV; 3135 m a.s.l.) volcanic complex, associated with fluid migration as a result of an intrusion of magma at around 5.8 km b.s.l. (Dp sources). See Table 2 for details on segments and associated mass changes. The southern caldera rim (CR) is shown for reference and the *inset* shows a digital elevation model of Tenerife, with the study area identified by a *white rectangle*.

However, closer inspection of the initial inversion results revealed that the line segments represent the superposition of deep and shallow seated sources. The decomposition of the gravity change data into shallow and deep fields (see Prutkin et al. 2014 for details) provided the basis

for the separate inversions of the two fields using the three-dimensional line segment approximation. The deep field inversion constrained two connected and strong line segments at a depth of about 5.8 km b.s.l. (marked Dp1 and Dp2 in Fig. 4), while the inversion of the shallow field
identified three similarly weak line segments situated at near-surface depths (<2 km from ground surface; see Fig. 4 and Table 2).

The most plausible interpretation of the inversion results is that the weak line segments represent sources dominated by hydrothermal fluids. In contrast, the deeper-seated sources can be interpreted as parts of an intrusion of new magma. This intrusion may have released fluids which consequently migrated towards shallower depths where they excited a shallow hydrothermal system. The geophysical signals resulting from these coupled magmatic-hydrothermal processes point towards a hybrid source nature for the unrest on Tenerife in 2004–2005. The identified link between deep and shallow unrest sources suggests the presence of permeable pathways for shallow fluid migration at the CVC.

Discussion

The Effect of Crustal Mechanics
on Stress, Strain and Pressure

The presented investigations at Uturuncu, Cotopaxi and SHV all incorporate subsurface heterogeneity, however the effects are somewhat different. At Cotopaxi the final inferred deformation source is shallow, at a level where the model does not incorporate substantial heterogeneity due to the limited amount of available seismic data. Hence, with this model configuration, the heterogeneity has not played a significant role in altering the location of the source. It is likely, however, that in reality the volcanic edifice and shallow subsurface does have a certain level of heterogeneity. It is therefore possible that had this been taken into account the inferred deformation source might be located more centrally beneath the edifice. Conversely, at Uturuncu and SHV subsurface heterogeneity played a crucial role in determining the deformation source locations. In both cases, the inclusion of vertical layering in the Young's Modulus distribution altered the inferred depth of the source; the heterogeneous models predict deeper sources than the generic homogeneous analytical models.

From this, it also follows that a horizontal vari-
ation in mechanical properties would influence
the horizontal location of a deformation source
(e.g., Hickey et al. 2016).

A second effect of subsurface heterogeneity relates to source pressure requirements. It is common with homogeneous crustal mechanics to require unrealistically high source pressurisation (>100's MPa) when attempting to fit an observed deformation signal, yet there is no petrological or mechanical evidence for such conditions. Theoretical work from a mechanical viewpoint suggests that the maximum overpressure a source can sustain without failing is equal to the tensile strength of the host rock (e.g., Gudmundsson 2011). Above this limit, a magmatic reservoir would trigger a dyke intrusion. Work at Uturuncu and SHV shows that elastic heterogeneous deformation models bring source over-pressure requirements more in line with both in situ and laboratory values of tensile strength. This is due to a relative reduction in the average Young's Modulus above the source, compared to a higher homogeneous Young's Modulus through the entirety of the crust, and thus a greater amount of deformation for a given pressure increment. Furthermore, with a maximum value for the source over-pressure, constraints can be placed on the minimum size of a deformation source, given the two are directly linked. This allows, for example, an estimate to be placed on the volume of a magma reservoir. When multiple magmatic sources might be present numerical models can additionally account for stress interactions between the two, something not considered by generic analytical models (Pascal et al. 2013).

As a next step to subsurface heterogeneity, the
case studies at Cotopaxi and SHV incorporated
temperature-dependent mechanics to evaluate the
effect of a viscoelastic rheology. At Uturuncu,
investigations were also carried out using a vis-
coelastic rheology without temperature-
dependence. In all three cases, a viscoelastic
medium reduced the over-pressure requirements
further when compared to the purely elastic
models. This is due to the viscous expansion that
follows an initial elastic inflation. The effect of
thermomechanics is greater than just modifying
source over-pressure requirements, however. Prolonged magma emplacement over thousands to millions of years in active volcanic areas builds up a significant thermal legacy within the crust. This results in elevated geothermal gradients, and in the case of continued active magmatism, deep crustal hot-zones (Annen 2009). These thermal perturbations are significant for the transfer of stress and strain. For example, at SHV, the combined thermomechanical effects of a deep-crustal hot zone and hot encasing rocks around a mid-crustal andesitic reservoir fundamentally alter the time-dependent subsurface stress and strain partitioning upon priming of the magma reservoir. These effects substantially influence surface strains recorded by volcano geodetic monitoring.

## Hybrid Unrest and Source Characterisation

An intrusion of magma commonly leads to the exsolution of fluids upon decompression, and fluid migration (and accumulation) in itself can produce measurable geodetic surface signals (e.g., Fournier and Chardot 2012; Rouwet et al. 2014, and references therein). A deep intrusion of magma, such as the 2004–2005 unrest on Tenerife, may not necessarily lead to observable surface deformation, meaning deformation data on its own can not provide any meaningful insights on the source process(es). In this case, the combination of deformation and gravity surveys allowed the characterisation of the unrest sources in much greater detail owing to their density contrasts, relative depths and mass additions. This highlights that, especially for the case of hybrid mechanisms where both magmatic and hydrothermal components are present, multi-disciplinary geodetic surveys can provide valuable information on source characterisation to help distinguish the processes driving unrest. This identification and discrimination of sources driving volcanic unrest via mathematical modelling of surface data is of vital importance for hazard and risk characterisation, given the different connotations associated with magmatic

and non-magmatic unrest, and plays a major role in probabilistic eruption forecasting (Rouwet et al. 2014).

## Application to Eruption Forecasting

The primary objective of the case studies was to develop advanced geodetic models to interpret spatial and temporal deformation monitoring signals, and provide better constraints on the subsurface processes causing volcanic unrest. This has been achieved by incorporating more plausible model components, such as subsurface heterogeneity, topography and temperature-dependent mechanics, to relax the assumptions that hinder analytical models and maintain consistency with inferences from geophysics, geology and petrology. Crucially, this highlighted the importance of pressure-time functions and inelastic rheology in deciphering temporal deformation patterns. In turn, a more thorough understanding of how a volcanic system behaves through time will benefit eruption forecasting, as quantitative estimates of key parameters such as magma supply rate and mechanism can be deduced.

When only considering the spatial deformation pattern, this work has further demonstrated some of the pitfalls associated with models assuming homogeneous, elastic, half-spaces (e.g., Mogi 1958). Crustal heterogeneity significantly effects the horizontal and vertical location of a deformation source by altering the subsurface strain distribution. Surface strain partitioning by complex topography is equally important, such as at steep-sided stratovolcanoes. The Cotopaxi case highlighted complex partitioning behaviour along deep ravines and adjacent lava flow ridges. If ground-based geodetic monitoring sites are positioned at localities with steep topographic gradients, surface strain partitioning plays a first-order role and needs to be accounted for during data modelling. Consequently, estimates of source parameters from generic, analytical models which cannot account for heterogeneities or complex topography may skew eruption forecasting and risk mitigation
efforts. On the other hand, more reliable deformation source locations from numerical models can be used to better estimate where and when an eruption might occur, and the associated hazard. It should also be pointed out that data limitations at some volcanoes may prevent complex numerical models being used, in which case simpler models may be the only option and their results should be carefully scrutinised and only applied cautiously to further analyses.

## Conclusions

The case studies provide new interpretations of volcanic processes during intra-eruptive and non-eruptive unrest. They also provide observations for distinct magmatic settings, thus contributing to the ongoing global comparisons of deformation and unrest, and the process of pattern recognition for identification of eruption precursors. The incorporation of multi-disciplinary data into integrated geodetic models has closed the gap between observations and interpretations of volcanic deformation and gravimetric changes. This will help to improve eruption forecasting as it moves away from a qualitative approach towards incorporation of more quantitative data derived from well-constrained physical mechanisms.

**Acknowledgements** Work presented herein has received funding by the European Commission (FP7; Theme: ENV.2011.1.3.3-1; Grant 282759: VUELCO). We are very grateful to Alvaro Guevara and Irving Munguía Gonzalez for their help translating Spanish, and to two anonymous reviewers for their constructive comments.

## Glossary Terms

**Deformation** the action of changing the shape of a physical object (e.g., a volcano) under the influence of some stress.

**Elastic** rheological behaviour in which an applied stress causes an immediate strain that is 100% recoverable when the stress is removed. See also rheology and viscoelastic.

**Finite Element Method/Analysis** a numerical modelling technique that subdivides an entire problem into a set of smaller 'elements'. Mathematical problems are then solved in each of the elements and combined together to calculate the response of the whole object.

**Geodesy** applied mathematics with the aim of measuring the geometric spatial representation of the Earth and its gravitational field in time-varying three-dimensions, as well as their orientation in space. Volcano geodesy specifically deals with recording the spatial and temporal patterns of crustal deformation and gravimetric changes in volcanic settings.

**Gravimetry** field of geophysics devoted to observing, processing and interpreting minute changes in the Earth's gravitational field.

**Hydrothermal system** an area in which heat and fluids from a partially molten magmatic body interact with a multi-phase groundwater system, causing chemical and thermal (heating) perturbations to the water.

**Mechanics (crustal, thermal)** the branch of physics that studies how stress can effect a physical object. Crustal mechanics relates to the values of the material parameters that describe the mechanical behaviour of the Earth's crust. Thermal mechanics relates to the variation in mechanical material properties due to changes in temperature.

**Model (inverse, forward)** a simulation of a process. Inverse models solve 'backwards' to determine optimal model parameters that fit a set of data. Forward models use predefined constant input parameters to calculate the expected model response.

**Rheology (crustal)** the behavioural response of the Earth's crust to forces that act upon or within it. See also elastic and viscoelastic.

**Strain** the change in dimension of an object (e.g., $\Delta X$) relative to the original dimension of the object (e.g., $X$). It has no units, as the units cancel: $\Delta X/X$.

**Stress** a measure of force per unit area, with the unit of pascals, Pa. $1 \text{ Pa} = 1 \text{ N/m}^2$.
**Viscoelastic** rheological behaviour in which an applied stress causes an elastic response, followed by a delayed (viscous) response. See also rheology and elastic.

**Index Terms**

Magma reservoir, deformation, gravity, GPS, InSAR, Cotopaxi, Soufrière Hills, Las Cañadas, Uturuncu, unrest, thermomechanics, modelling, finite element, inversion/inverse model, crustal mechanics, hybrid unrest, Young's Modulus, Poisson's Ratio, viscosity.

**References**

Annen C (2009) From plutons to magma chambers: thermal constraints on the accumulation of eruptible silicic magma in the upper crust. Earth Planet Sci Lett 284(3-4):409-416

Battaglia M, Gottsmann J, Carbone D, Fernández J (2008) 4D volcano gravimetry. Geophysics 73(6):WA3-WA18

Brocher TM (2005) Empirical relations between elastic wavespeeds and density in the Earth's crust. Bull Seismol Soc Am 95(6):2081-2092

del Potro R, Dez M, Blundy J, Camacho AG, Gottsmann J (2013) Diapiric ascent of silicic magma beneath the Bolivian Altiplano. Geophys Res Lett 40(10):2044-2048

Fialko Y, Pearse J (2012) Sombrero uplift above the Altiplano-Puna magma body: evidence of a ballooning mid-crustal diapir. Science 338(6104):250-252

Fournier N, Chardot L (2012) Understanding volcano hydrothermal unrest from geodetic observations: insights from numerical modeling and application to White Island volcano, New Zealand. J Geophys Res 117(B11):B11208

Gottsmann J, Odbert HM (2014) The effects of thermo-mechanical heterogeneities in island-arc crust on time-dependent pre-eruptive stresses and the failure of an andesitic reservoir. J Geophys Res. doi:10.1002/2014JB011079

Gottsmann J, Rymer H, Berrino G (2006a) Unrest at the Campi Flegrei caldera (Italy): a critical evaluation of source parameters from geodetic data inversion. J Volcanol Geoth Res 150(1-3):132-145

Gottsmann J, Wooller L, Marti J, Fernandez J, Camacho AG, Gonzalez PJ, Garcia A, Rymer H (2006b) New evidence for the reawakening of Teide volcano. Geophys Res Lett 33(20)

Gudmundsson A (2011) Rock fractures in geological processes. Cambridge University Press, New York

Hautmann S, Gottsmann J, Sparks RSJ, Mattioli GS, Sacks IS, Strutt MH (2010) Effect of mechanical heterogeneity in arc crust on volcano deformation with application to Soufrière Hills volcano, Montserrat, West Indies. J Geophys Res 115(B9):B09203

Hautmann S, Witham F, Christopher T, Cole P, Linde AT, Sacks IS, Sparks RSJ (2014) Strain field analysis on Montserrat (W.I.) as tool for assessing permeable flow paths in the magmatic system of Soufriere Hills Volcano. Geochem Geophys Geosyst 15:2013G. doi:10.1002/C005087

Hickey J, Gottsmann J (2014) Benchmarking and developing numerical finite element models of volcanic deformation. J Volcanol Geoth Res 280:126-130

Hickey J, Gottsmann J, del Potro R (2013) The large-scale surface uplift in the Altiplano-Puna region of Bolivia: a parametric study of source characteristics and crustal rheology using finite element analysis. Geochem Geophys Geosyst 14(3):540-555

Hickey J, Gottsmann J, Mothes P (2015) Estimating volcanic deformation source parameters with a finite element inversion: the 2001-2002 unrest at Cotopaxi volcano, Ecuador. J Geophys Res Solid Earth 120(3):1473-1486

Hickey J, Gottsmann J, Nakamichi H, Iguchi M (2016) Thermomechanical controls on magma supply and volcanic deformation: application to Aira caldera, Japan. Sci Rep 6(August):32691

Masterlark T (2007) Magma intrusion and deformation predictions: sensitivities to the Mogi assumptions. J Geophys Res 112(B06419)

Mogi K (1958) Relations between the eruptions of various volcanoes and the deformations of the ground surfaces around them. Bull Earthq Res Inst 36:99-134

Molina I, Kumagai H, Garca-Aristizábal A, Nakano M, Mothes P (2008) Source process of very-long-period events accompanying long-period signals at Cotopaxi Volcano, Ecuador. J Volcanol Geothermal Res 176(1):119-133

Odbert HM, Stewart RC, Wadge G (2014) Cyclic Phenomena at the Soufrière Hills Volcano, Montserrat, The Eruption of Soufrière Hills Volcano, Montserrat from 2000 to 2010, vol 39. The Geological Society, London, pp 41-60. doi:10.1144/M39.2

Pascal K, Neuberg J, Rivalta E (2013) On precisely modelling surface deformation due to interacting magma chambers and dykes. Geophys J Int 196(1):253-278

Pistolesi M, Cioni R, Rosi M, Cashman KV, Rossotti A, Aguilera E (2013) Evidence for lahar-triggering mechanisms in complex stratigraphic sequences: the post-twelfth century eruptive activity of Cotopaxi Volcano, Ecuador. Bull Volcanol 75(3):698

Pritchard ME, Simons M (2002) A satellite geodetic survey of large-scale deformation of volcanic centres in the central Andes. Nature 418(6894):167-171

Prutkin I, Vajda P, Gottsmann J (2014) The gravimetric picture of magmatic and hydrothermal sources driving hybrid unrest on Tenerife in 2004/5. J Volcanol Geoth Res 282:9-18
Ranalli G (1995) Rheology of the Earth. Chapman and Hall, London

Rouwet D, Sandri L, Marzocchi W, Gottsmann J, Selva J, Tonini R, Papale P (2014) Recognizing and tracking volcanic hazards related to non-magmatic unrest: a review. J Appl Volcanol 3:1–17

Sparks RSJ, Folkes CB, Humphreys MC, Barfod DN, Clavero J, Sunagua MC, McNutt SR, Pritchard ME (2008) Uturuncu volcano, Bolivia: volcanic unrest due to mid-crustal magma intrusion. Am J Sci 308(6):727–769

Sparks RSJ, Biggs J, Neuberg JW (2012) Monitoring volcanoes. Science 335(6074):1310–1311

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

<figure><img src="image_89.png" /></figure>
<figure><img src="image_90.png" /></figure>

# Fluid Geochemistry and Volcanic Unrest: Dissolving the Haze in Time and Space

Dmitri Rouwet, Silvana Hidalgo, Erouscilla P. Joseph and Gino González-Ilama

## Abstract

The heat and gas released by a degassing magma affects the overlying predominantly meteoric aquifers to form magmatic-hydrothermal systems inside the solid body of a volcano. This chapter reviews how fluid geochemical signals help to track the evolution throughout the various stages of volcanic unrest. A direct view into a degassing magma is possible at open-conduit degassing volcanoes. Nevertheless, in most cases gas is trapped (i.e. scrubbed) by abundant water, leading to the loss of the pure signal the magma ideally provides. Deciphering how magmatic gas rises through, reacts, and re-equilibrates with the liquids in the magmatic-hydrothermal system in time and space is the only way to trace back to the pure signal. The most indicative magmatic gas species ($CO_2$, $SO_2-H_2S$, HCl and HF) are released as a function of their solubility in magma. The less soluble gas species are released early from a magma at higher pressure conditions ($CO_2$) (deeper), whereas the more soluble species are released later, at lower pressures ($SO_2$, HCl and HF) (shallower depth). When these gases hit the water during their rise towards the surface, they will be more or less scrubbed. Depending on the chemical equilibria inside the magmatic-hydrothermal system (e.g. $SO_2-H_2S$ conversion, acidity), the gas that eventually reaches the surface will carry the history of its rise from bottom to top. Tracking volcanic unrest implies a time frame; the kinetics of magma degassing throughout the liquid cocktail inside the volcano impose the maximum resolution the volcano provides and hence the monitoring time window to be adopted for

<footnote>D. Rouwet (✉)
Istituto Nazionale di Geofisica e Vulcanologia,
Sezione di Bologna, Via Donato Creti 12, 40128
Bologna, Italy
e-mail: dmitri.rouwet@ingv.it</footnote>

<footnote>S. Hidalgo
Instituto Geofisico, Escuela Politécnica Nacional,
Quito, Ecuador</footnote>

<footnote>E.P. Joseph
Seismic Research Centre, University of the West
Indies, St. Augustine, Trinidad & Tobago</footnote>

<footnote>G. González-Ilama
Centro de Investigaciones en Ciencias Geológicas,
Universidad de Costa Rica, San José, Costa Rica</footnote>
each volcano. Gas-dominated systems are “faster” and require a higher monitoring frequency, water-dominated systems are slower and require a lower monitoring frequency.

## Resumen

El calor y gas liberados por la desgasificación del magma afecta los acuíferos de origen predominantemente meteórico para formar sistemas magmático-hidrotermales dentro el cuerpo sólido del volcán. Este capítulo revisa como la geoquímica de fluidos puede ayudar a trazar la evolución a través de las varias etapas de “unrest” volcánico. Una visión directa dentro de un magma en desgasificación es solo posible para volcanes de conducto abierto. Sin embargo, en la mayoría de las situaciones el gas queda atrapado (i.e. “scrubbing”) en el agua, que conduce a la pérdida de la señal pura que el magma idealmente puede proporcionar. Descifrar como el magma sube a través de los líquidos, y reacciona y re-equilibra con ellos dentro el sistema magmático-hidrotermal, en un marco de tiempo y espacio, es la única manera para rastrear el origen de la señal del magma. La desgasificación de magma se da por cuatro procesos: (1) durante la subida de magma, (2) por la descompresión debido al eliminar una porción del edificio volcánico, (3) debido a la convección interna dentro la cámara magmática, o (4) después de “ebullición secundaria” siguiendo el enfriamiento y consecuente cristalización. Las especies gaseosas magmáticas más indicativas ($CO_2$, $SO_2-H_2S$, HCl y HF) se liberan en función de su solubilidad en el magma. Las especies menos solubles se liberan antes del magma, bajo regímenes de presiones más altas ($CO_2$), mientras que las especies más solubles se liberan después, bajo regímenes de presiones más bajas ($SO_2$, HCl y HF). En términos espaciales, $CO_2$ se libera a lo largo de una área espacial más amplia (desde lo más profundo). La presencia de $SO_2$ es una indicación clara de un magma que sube hacia un ambiente más somero. La llegada de HCl en la superficie generalmente indica la presencia de una remesa de magma somera (cientos de metros hasta pocos kilómetros). Especialmente un aumento en la proporción $CO_2/SO_2$ es indicativo para elucidar un estado de “inquietud” (“unrest”). Una disminución consecutiva en el $CO_2/SO_2$, después de un aumento, es una indicación de que el magma está cerca de la superficie y es propenso a una erupción. Cuando estos gases alcanzan el agua durante su ascenso hacia la superficie, serán más, o menos, absorbidos. Dependiendo de los equilibrios químicos dentro el sistema magmático-hidrotermal (e.g. conversión de $SO_2-H_2S$, acidez), el gas que al final llega a la superficie lleva consigo la historia de su ascenso desde el fondo hasta la superficie. La desgasificación magmática es un proceso más rápido, mientras que la dinámica hidrotermal en el sistema rocoso FeO-FeO$_{1.5}$ es más lenta. Por eso, el $H_2S$ se suele llamar un “gas hidrotermal”, y el $SO_2$ un “gas magmático”. Trazar “unrest volcánico” implica un encuadramiento de tiempo más específico. Si la ventana de tiempo de monitoreo es más largo
que el tiempo definido por la cinética de la migración de gas, los detalles de la dinámica de desgasificación se perderán inevitablemente. Contrariamente, si la ventana de tiempo de monitoreo es más corto que la ventana definida por la cinética de la migración de gas, en tal caso, resultará en una visión demasiado detallada de lo que el sistema magmático-hidrotermal puede proveer. Estudios recientes han demostrado que acuíferos extremadamente ácidos pueden "desacelerar" la señal dejada por el sistema magmático-hidrotermal dominado por el gas (e.g. fumarolas), pero pueden "acelerar" la señal dejada por el sistema magmático-hidrotermal dominado por agua (e.g. lagos cratéricos ácidos). Estos hallazgos tienen implicaciones significativas para el encuadramiento de tiempo en reconocer la desgasificación magmática, y por tanto, para la frecuencia del monitoreo.

**Keywords**

Fluid geochemistry • Magmatic-hydrothermal systems • Volcanic unrest • Volcano monitoring

**Palabras clave**

Geoquímica de fluidos • Sistemas magmático-hidrotermales • Inquietud volcánica • Monitoreo volcánico

**Introduction**

Fluid geochemical monitoring tracks variations in gaseous species and fluid phases released in various manners from a volcano. Fluids infiltrate, move, migrate, rise, react and re-equilibrate in the water- and vapor-filled solid body of the volcano. These processes are invisible as they occur in the subsurface, and can only be deduced from measurements at the surface. Nevertheless, gas and water are more mobile than rock and, when a volcano shifts into the gear of unrest, a change in degassing is often the first sign to be detected. As such, fluid geochemistry offers a crucial means to recognize unrest in a timely matter.

Most volcanic edifices store a large volume of water, as cold or thermal aquifers, with various hydrogeological architectures. When gas hits water in its rise towards the surface, the original signature of the gas is mostly lost. This disadvantage can however be overcome.

Understanding how gas absorbs and solutes react in the liquid phase, and how it is eventually released to the atmosphere, is key to linking surface manifestations to magma dynamics.

Volcano monitoring largely focuses on how, where and when magma migrates towards the surface, as, intrinsically, every magmatic eruption is anticipated by magma rise. Nevertheless, volcanoes can become hazardous even without an eruption. A volcano is in a state of magmatic unrest if we recognize signals of a magma-on-the-move; in any other situation, given volcanic unrest, the volcano is in a state of non-magmatic unrest (hydrothermal or tectonic) (Rouwet et al. 2014a). To recognize volcanic unrest, the background behavior of a volcanic system should be tracked for a sufficiently long period, in order to know when a deviation from this background becomes a cause for concern (i.e. unrest, Phillipson et al. 2013). This background behavior and deviations from it are volcano-specific and can be monitored in several
ways, besides the geochemical approach reviewed here.

The aim of this chapter is to scan through magmatic-hydrothermal systems during the process of magmatic degassing, from bottom to top, and describe how fluids behave in time and space. What are the lessons learned from fluid geochemistry throughout the evolution of volcanic quiescence, re-awakening, volcanic unrest, magmatic unrest and non-magmatic/hydrothermal unrest?

## Magmatic-Hydrothermal Manifestations

One of the first signs of re-awakening after prolonged volcanic quiescence to a state that eventually causes concern is often the appearance of fumarolic exhalations from a crater. This happened, for instance, in 1994 during the reawakening of Popocatépetl (Mexico), and recently, at Cotopaxi (Ecuador), both VUELCO target volcanoes (De la Cruz-Reyna and Tilling 2008; Hall and Mothes 2008) (Fig. 1a). For open-conduit volcanoes the presence of a plume (i.e. a visible gas-vapor cloud originating from an open volcano crater) can become the prominent manifestation of degassing.

Some volcanoes are characterized by decade to century long high-temperature fumarolic degassing in a closed-conduit setting (Fig. 1b), suggesting the presence of a stable, but shallow magma chamber. This constantly high-temperature background degassing is generally no sign of unrest (e.g. Momotombo's >700 °C, Satsuma-Iwojima's >900 °C and Kudryavy's >700 °C; Menyailov et al. 1986; Shinohara et al. 1993; Taran et al. 1995), whereas the increase of fumarolic temperatures from low (boiling point of water at a given altitude, hence atmospheric pressure) to high (above boiling to magmatic temperature) can be a sign of resumed unrest (e.g. the 1980–1990s crisis at Vulcano, Italy, Capasso et al. 1997).

Magmatic-hydrothermal systems are aquifers inside a volcano or beneath a volcanic area, heated by a magma, at an unspecified depth.

The origin of the water is generally meteoric (i.e. rain, snow and its melt water). How hot, and how gas-rich such magmatic-hydrothermal systems are depend on the proportion of the water volume with respect to the heat and gas provided by the magma. The latter depends on the residing depth of the magma. Boiling-temperature fumaroles are a common manifestation at magmatic-hydrothermal systems (Fig. 1c). When the water table of such systems intersects the surface, in craters or volcano flanks, boiling pools appear (Fig. 1d). Such pools can manifest bubbling degassing, and are nothing less than a water-rich fumarole (Fig. 1e). Depending on the dominant gas they exhale, paired with water vapor, the manifestations are called solfataras (S-rich gases) or mofettes (CO₂-rich gases); depending on the temperature and vapor/water proportion they emit they are called fumaroles (boiling or above boiling steam vents), thermal springs (liquid water emission) or geysers (water + vapor jets with a cyclic behavior). Thermal springs can discharge inside active craters or on volcano flanks in a degassed state, without bubbling or boiling (Fig. 1f). Heated water can fill (parts of) craters and form volcanic lakes (Fig. 1g). Depending on the degassing state and depth of the underlying magma, degassing features (bubbling or diffuse degassing) and evaporation can occur at the lake surface (Fig. 1h).

The pictures of the degassing manifestations in Fig. 1 show a trend from gas-dominated, active plume degassing in an open-conduit setting towards more water-dominated, hydrothermal, passive degassing in a closed-conduit setting. These visual observations only give a first glance of magmatic-hydrothermal activity, and do not reflect the state of unrest of a volcano. Throughout the life-time of magmatic-hydrothermal systems (centuries to millenia) volcanoes can evolve from gas-dominated to water-dominated, and vice versa. The next sections present what we know on the theoretical level, following the laws of chemistry. Despite these classic rules, it will become clear that the range of manifestations and variations in fluid signatures is wide, and volcano-dependent.
<figure><img src="image_91.png" /></figure>

**Fig. 1** Degassing manifestations at magmatic-hydrothermal systems. **a** Open-conduit degassing at Popocatépetl. **b** High-temperature fumarolic degassing at Vulcano. **c** Boiling-temperature fumarolic degassing at the Arbol Quemado fracture in the Turrialba crater area. **d** Boiling pools inside the El Chichón crater. **e** Bubbling

thermal spring at the SE flank of El Chichón. **f** Thermal flank spring at El Chichón. **g** The El Chichón crater lake. **h** Evaporative degassing from Laguna Caliente, Poás. The cyan arrow points towards water-dominance; the yellow arrow towards more gas-dominance

## Magma Degassing from Bottom to Top

### Magma Degassing

The degassing of a magma increases when the confining pressure in the magma decreases. Magma decompression can occur in four ways (Fig. 2): (1) magma rise, (2) decompression by uncovering a portion of the volcanic edifice, (3) internal convection in a magma chamber, or (4) secondary boiling upon cooling and consequent crystallization. The first process is often induced by the input of a deeper magma, into a shallower magma chamber. Magma rise towards the lower pressure regime results from the buoyancy difference between the stagnant

magma, the rising new melt and the surrounding rocks. The second degassing process can be triggered by the mass removal from part of the volcanic edifice. This superficial process can even trigger eruptions (e.g. 1980 Mt St. Helens). Internal magma convection causes degassing of a less dense, gas-rich magma batch. Once degassed, the now denser magma batch sinks (e.g. Stromboli, Aiuppa et al. 2009). The fourth process increases the gas/melt ratio in the magma, due to the loss of crystals from the cooling magma, hence favoring degassing.

Based on gas geochemistry only it is hard to rule out which process actually occurs at depth, and thus distinguish between magmatic and non-magmatic unrest. Geophysical signals are needed. Following the definition of magmatic unrest by Rouwet et al. (2014a) (i.e. the
recognition of a magma-on-the-move), only the first process is initially consistent with the requisite of magmatic unrest.

Variations in magma degassing can be detected both qualitatively and quantitatively. Detailed insights in the degassing state of a magma can only be obtained if both are measured contemporaneously. Which gas species are released, how much and when? When a magma starts to degas, by any of the above processes, the less soluble species is released first (i.e. at higher confining pressure in the magma chamber). The order in solubility of indicative magmatic gas species is CO₂ < SO₂ < HCl < HF; the order of release when a magma progressively degasses is “CO₂-first till HF-last” (Giggenbach 1987). Hence, tracking variations in ratios between

these species gives qualitative insights into the degassing state of a magma. A consecutive increase with time in first CO₂/SO₂, then SO₂/HCl, then HCl/HF ratios reflects the evolution in degassing state from a magma moving from depth towards the surface. Especially an increase in the CO₂/SO₂ is indicative of the state of unrest, pointing to an input of poorly degassed magma at great depths. A consecutive decrease in CO₂/SO₂ ratio, after the increase, is an indication of magma moving towards the surface. The latter two ratios come into play when eruption of magma is imminent, or even ongoing: the highly soluble species HCl and HF are released from a highly degassed magma, a situation that reflects near-surface degassing (Aiuppa et al. 2002). The arrival of HCl at the surface (e.g. in

<figure><img src="image_92.png" /><figcaption>**Fig. 2** Sketch of the four mechanisms that instigate magma degassing (not to scale)</figcaption></figure>
fumaroles) generally indicates the presence of a shallow magma batch (hundreds of meters, or less).

The “purest” proxy of the magmatic gas is provided by direct sampling of a fumarole with a near-magmatic temperature, followed by the analysis of its chemical composition. During the past five-six decades fumarolic gases have been extensively sampled and analyzed representing a wide range of temperatures and states of volcanic activity. Table 1 presents a compilation of chemical compositions of fumarolic gases for the high- and low-temperature ranges.

Before dealing with the abundant water vapor in fumaroles, we tackle the other gas species in the “dry-gas phase”. Regardless of the fumarole temperature, CO₂ is the most abundant gas species (Table 1). CO₂ is often released at the surface across a wider spatial extent, as it degasses from greater depth (Fig. 3). Old and deep magma bodies can continue to release CO₂ for tens to hundreds of millenia.

The second most abundant “dry-gas” species are the sulphur species SO₂ and H₂S. SO₂ is more soluble in magma than CO₂, and will thus be released at lower pressure. As SO₂ degasses at lower depth, the degassing tends to be more centralized along the central conduit (open or closed) of the volcano (Fig. 3). As such, SO₂ is often measured in volcanic plumes. An increase of SO₂ has often been interpreted as an indicator of magma rise into the shallower environment.

At high-temperature magmatic conditions the following reaction applies (Giggenbach 1987; Delmelle and Bernard 2015):

$$3\text{SO}_2 + 7\text{H}_2 = \text{H}_2\text{S} + 2\text{S}^\circ + 6\text{H}_2\text{O} \quad (1)$$

This explains the presence of both SO₂ and H₂S in high-temperature fumaroles. The SO₂/H₂S ratio is sensitive to temperature and the oxidation state (i.e. the role of H₂ in Eq. 1 reflects the oxidation vs reduction state, Giggenbach 1987).

On the other hand, it is noted that the low-temperature fumaroles lack SO₂, instead, the dominant sulphur species is H₂S (Table 1). Low-temperature hydrothermal conditions (T < 300 °C) favor the reduced S-species H₂S (Table 1), equilibrated by the rock phase, following the reaction:

$$\text{SO}_2 + 6(\text{FeO})_{\text{rock}} + \text{H}_2\text{O} = \text{H}_2\text{S} + 6(\text{FeO}_{1.5})_{\text{rock}} \quad (2)$$

or simply through reduction of SO₂ in the gaseous environment:

$$\text{SO}_2 + \text{H}_2 = \text{H}_2\text{S} + \text{O}_2 \quad (3)$$

To convert SO₂ into H₂S in a rock matrix (FeO–FeO₁.₅-system), following reaction (2), a major constraint is “sufficient time”. Magmatic degassing (Eq. 1) is a faster process, while hydrothermal dynamics (Eq. 2) are slower. This

<table><caption>**Table 1** Chemical composition of fumarolic gases (concentrations in micromol/mol), for high (magmatic, 280°–1130 °C) and low (hydrothermal, 83°–160 °C) temperature conditions, expressed as the minimum and maximum measured concentrations for 12 and 59 samples, respectively</caption><thead><tr><th></th><th>T(°C)</th><th>H<sub>2</sub>O</th><th>CO<sub>2</sub></th><th>SO<sub>2</sub></th><th>H<sub>2</sub>S</th><th>HCl</th><th>HF</th><th>H<sub>2</sub></th><th>N<sub>2</sub></th><th>CH<sub>4</sub></th><th>CO</th></tr></thead><tbody><tr><td>High-T minimum (#12)</td><td>280</td><td>311000</td><td>1200</td><td>320</td><td>4</td><td>275</td><td>21</td><td>3</td><td>40</td><td>0.1</td><td>0.2</td></tr><tr><td>High-T maximum (#12)</td><td>1130</td><td>993000</td><td>672000</td><td>67800</td><td>21460</td><td>14200</td><td>2500</td><td>14900</td><td>1800</td><td>7.1</td><td>4600</td></tr><tr><td>Low-T minimum (#59)</td><td>83</td><td>638200</td><td>2655</td><td>0</td><td>0</td><td>0</td><td>0</td><td>16</td><td>32</td><td>0.4</td><td>0.011</td></tr><tr><td>Low-T maximum (#59)</td><td>160</td><td>997200</td><td>355000</td><td>0</td><td>3700</td><td>0</td><td>0</td><td>1220</td><td>6800</td><td>5330</td><td>1.6</td></tr></tbody></table>
<figure><img src="image_93.png" /><figcaption>**Fig. 3** Sketch of the passing “gas train” during magma degassing, resulting from the difference in solubility of the various gas species (not to scale)</figcaption></figure>

is why H₂S is often called a “hydrothermal gas” and SO₂ a “magmatic gas” (Table 1).

As noted in Table 1, the major gas species in fumaroles is often water, regardless of their temperature. Even high-temperature fumaroles are water-dominated.

Knowing that water boils at 100 °C (at atmospheric pressure, at sea level), in theory, the temperature of a fumarole is buffered at 100 °C and cannot rise until water is exhausted in the underlying plumbing system. Moreover, the critical temperature of water is 374 °C (i.e. the temperature at which vapor and liquid water cannot coexist anymore), imposing a second temperature buffer. This implies that (1) water directly originates from a high-temperature magma under super-critical conditions (“andesitic water” with T > 374 °C, Taran et al. 1989), and/or (2) water is excessively present in the fumarole plumbing system with respect to the gas, and will hardly ever exhaust.

## When the Gas Hits the Water

The shallow subsurface environment of the different sections of the Earth's crust hosts numerous aquifers at various depths originating from the infiltration and storage of meteoric water, or seawater in the case of low-lying islands; volcanic edifices are no different, being “small dots”

at the Earth's surface. As described before, a magma degasses and heats the space between the magma and the surface, and will inevitably heat and modify the volcanic aquifers. As the magma heats the overlying aquifers, between the magma and the surface, from bottom to top, a gas-only, vapor + gas zone, a vapor + liquid zone and a liquid-only zone can be found (Fig. 4). When a magma rises, or the heat input from the magma increases, the vapor + liquid zone or vapor + gas zone will be pushed upwards until intersecting the surface (e.g. in a volcano crater), manifested at the surface as boiling or bubbling pools and fumaroles (Fig. 4). If the distance between the magma and surface is larger, the thermal aquifer will intersect the surface and create thermal springs (Fig. 4). Facing unrest, a liquid to vapor transition at a surface manifestation reflects heating of the hydrothermal system.

When the above “gas train” (Fig. 3) consecutively reaches the liquid-only zone, gas species will be absorbed and react depending on their specific chemical properties in water. The capacity to absorb magmatic gases in the liquid phase is called “scrubbing” (Symonds et al. 2001). The CO₂ that reaches the water from greater depths during magma degassing (Fig. 3), will create CO₂-dominated bubbling thermal springs, and HCO₃-rich slightly acidic springs (pH 5–7) (Fig. 4). The second least soluble gas species that hits the water is SO₂ that will be
<figure><img src="image_94.png" /><figcaption>**Fig. 4** Sketch of the gas that hits the water in an hypothetical “wet volcano” (not to scale). L = liquid only zone (turquoise area), V + L = vapor + liquid zone, (yellow area) V + G = vapor + gas zone (grey area), G = gas only zone (white area around the degassing</figcaption></figure>

magma, dark grey). Red dotted lines are isotherms. Cl-geothermal water are deep remnant waters, not of interest for geochemical monitoring of volcanoes for being “old and slow”

hydrolyzed as sulphuric acid (H₂SO₄) and dissociate into HSO₄⁻ or SO₄²⁻, depending on the pH of the water. The dissolved H⁺ creates the high acidity or low pH, following pH = -log [H⁺]. Volcanic environments are renowned for being acidic; the acid is generated by the scrubbing of acidic gases into the water.

The following SO₂ disproportionation reactions occur (Kusakabe et al. 2000):

$$3\text{SO}_2 + 3\text{H}_2\text{O} = 2\text{HSO}_4^- + \text{S}^\circ + 2\text{H}^+ \quad (4)$$

or

$$4\text{SO}_2 + 4\text{H}_2\text{O} = 3\text{HSO}_4^- + \text{H}_2\text{S} + 3\text{H}^+ \quad (5)$$

Reaction (4) occurs under relatively oxidizing conditions, low temperature and high total sulphur concentrations, whereas reaction (5) occurs under relatively reducing (i.e. oxygen poor) conditions, high temperature and low total sulphur concentrations. “To disproportionate SO₂” means that for each three or four moles of SO₂ that hit the water, one mole of S°, or one mole of

H₂S is given in return (reaction 4 and 5, respectively).

As SO₂ degasses at shallower depth in a magma system than CO₂, the resulting acid SO₄-rich springs are found near the central degassing conduit (e.g. inside active craters). For SO₂-dominated magmatic systems a pH near 2 or less is common; for H₂S-dominated hydrothermal systems a pH of 2–2.5 is the most acidic water can get (i.e. hydrothermal or “steam heated” waters).

The next gas species to be added to the “liquid cocktail” are HCl and HF (Fig. 3). HCl is highly volatile, but also highly hydrophile. This *contradictio in terminis* means that when HClgas reaches water it will be trapped in the liquid phase as Cl⁻ and H⁺. Generally, Cl⁻ is considered “conservative” in the liquid phase, and is therefore often used as a tracer in the hydrothermal environment (see Section “Rock Leaching Upon Weathering”). Conservative means that Cl⁻ can be hardly lost from the solution as Cl-salts are highly soluble in acidic
and hot environments, and HCl should not degas from “high pH water” (>2). Sulphate minerals in their turn can be lost from solutions by precipitation, demonstrating their non-conservative character (Rouwet and Ohba 2015). For extremely acidic environments (pH 0 or <0) the reaction

$$HCl = H^+ + Cl^- \quad (6)$$

moves to the left due to an H⁺ excess with respect to Cl⁻ (i.e. HCl degassing). The same counts for HF. Bearing in mind the abundance of the acid SO₂ in the magmatic-hydrothermal environment, providing a large part of the acidity, this “secondary HCl degassing from the liquid phase” is less unexpected than previously thought. This implies that HCl can degas from a low-temperature aquifer, as long as the aquifer is extremely acid (Capaccioni et al. 2016). Moreover, as hot water releases vapor, this loss enriches the remnant liquid in solutes, including H⁺ (i.e. salinity and acidity increase), leading to the fact that even SO₂ tends to degas from the liquid, instead of remaining in the water phase as SO₄²⁻ or HSO₄⁻.

Acid water (pH < 3.8) is completely transparent for the omni-present CO₂. This implies that CO₂, one of the “deepest signals” available for a degassing volcano, will completely outgas from bottom to top. In active craters, often underlain by acidic thermal aquifers, the release of CO₂ thus behaves as though there is no water present. This is a great advantage to monitoring volcanoes and tracking unrest, especially to detect the onset of unrest.

In conclusion, the above insights demythologize two generally accepted facts: (1) high-temperature fumaroles cannot contain water vapor (Table 1), and (2) low-temperature fumaroles cannot release acidic gases.

Within the scope of this book, tracking unrest using fluid geochemistry requires the introduction of a time frame, or a monitoring time window and frequency. Does the fumarole reflect the exact moment of degassing, or is it rather an average degassing for the longer previous period stored and steadily released from the

magmatic-hydrothermal plumbing system? What is the time delay between the moment the gas hits the water and the eventual release at the surface? The kinetics (i.e. “speed”) of the gas migration from a magma towards the surface are still poorly constrained. If we can estimate the residence time of gas and water in the magmatic-hydrothermal system, we are able to define a monitoring time window, and hence adopt an adequate monitoring frequency. As explained, the acidity of the feeding aquifer plays a role. For extreme acidic conditions, less gas scrubbing occurs and the fumarolic system will react faster, and hence, shorter monitoring time windows can be adopted. If the monitoring time window is longer than the window defined by the kinetics of gas migration, details in degassing dynamics will be lost. On the contrary, the monitoring time window should not be a lot shorter than the time window defined by the kinetics of gas migration, if so, it will provide a too detailed view of what the magmatic-hydrothermal system can maximally provide.

## The Other Liquid: Elemental Sulphur

Whereas water melts at 0 °C and boils at 100 °C, sulphur melts at 119 °C and boils at 444 °C (Oppenheimer and Stevenson 1989, and references therein). This physical fact on phase transitions implies that in the hydrothermal environment (boiling water) elemental sulphur is solid, and that during the initial phase of transition towards a more magmatic, high-temperature regime sulphur will become liquid. Molten sulphur in the hydrothermal plumbing system can be remobilized, clear vugs and vents and eventually be expelled as a liquid sulphur flow from fumarole mouths, or “fill-and-freeze” pores in the shallower hydrothermal system. The first process opens up degassing pathways; the second process can decrease rock porosity and permeability near the surface, thus sealing a hydrothermal system. During the evolution from low-temperature (>119 °C, unrest) towards high-temperature (occasionally magmatic unrest), the viscosity of the liquid sulphur increases
2000-fold (>160 °C) to become an extremely efficient sealer of a magmatic-hydrothermal system. Pressure build-up beneath seals in hydrothermal systems can lead to phreatic eruptions (Rouwet and Morrissey 2015; Rouwet et al. 2016). Monitoring fumarolic temperatures is thus essential, and probably one the simplest methods to apply.

# Tracking Hydrothermal Unrest and Related Hazards: Methods from Case-Studies

## From Quiescence to Unrest, to Phreatic Eruptions, to Magmatic Eruptions

### Turrialba, Costa Rica (2001–2016)

A transition from the stable, passively degassing hydrothermal system of Turrialba volcano (Costa Rica), to hydrothermal unrest, to phreatic eruptions, to magmatic eruptions, is a recent example of an evolution from volcanic quiescence heading towards eruption.

In 2001, increased fumarolic activity (appearance of SO₂ in late 2001) was paired to seismic swarms and ground deformation (Martini et al. 2010; Vaselli et al. 2009) (i.e. unrest). In 2007, the increased SO₂/H₂S molar ratio in fumaroles (>100), combined with an increase in exhalation temperature up to 282 °C (in early 2008, Martini et al. 2010; Vaselli et al. 2009), point to more oxidized and magmatic conditions (i.e. magmatic unrest). Clear plume degassing resumed in early 2007 (Fig. 5a), and SO₂ fluxes reached 740 t/d in January 2008 (Martini et al. 2010). In late 2009, fumarolic degassing was vigorous and extended into the Arbol Quemado fracture, newly formed in 2002. The first phreatic eruption occurred on 5 January 2010, in the inner crater wall of the actively degassing SW crater. Strong jet-like degassing occurred afterwards from this new vent (Fig. 5b), while diffuse fumarolic degassing diminished in the SW crater. A second phreatic eruption occurred on 12 January 2012, from a vent inside the Arbol Quemado fracture. The day before, this eruption was preceded by liquid sulphur flowing out of the

Arbol Quemado fumaroles (González et al. 2015). A third phreatic eruption episode involved the 2010 and 2012 vents simultaneously (21 May 2013). During this 3.5-year long phreatic cycle the SO₂ flux from Turrialba's plume was high: from 2500 to 4300 t/d (Campion et al. 2012; Moussallam et al. 2014). CO₂/SO₂ molar ratios in March 2013 were relatively low (2.6), hinting at a CO₂-depleted and SO₂-rich magma (Moussallam et al. 2014). From 2001 to 2013, Moussallam et al. (2014) suggest the progressive “drying-out” of the underlying hydrothermal system.

The first magmatic eruption at Turrialba since the 1864–1866 phreatomagmatic activity occurred during the night of 29 and 30 October 2014 (Mora-Amador et al. 2015). At the time of writing, the last magmatic eruptions took place in September 2016.

Despite the well-monitored and tracked evolution from volcanic quiescence to magmatic eruption it remains unclear why some volcanoes quickly evolve from quiescence to eruption, while at Turrialba it took 14 years from quiescence to magmatic eruption, passing the complete range of unrest manifestations during this relatively long time span.

### Cotopaxi, Ecuador (2015–2016)

Another example of volcanic quiescence to magmatic unrest is the one of Cotopaxi volcano in 2015. As Cotopaxi is a dangerous volcano whose activity would potentially affect densely populated areas its monitoring network has been continuously improved since the late 70s. After 73 years of quiescence, the first sign of unrest at Cotopaxi was a progressive increase in the amplitude of transient seismic events in April 2015. SO₂ is measured at Cotopaxi by DOAS stations installed on the flanks of the volcano since 2008. The permanently measured SO₂ emissions showed an increase on May 20 from almost non detectable up to ~3000 t/d. The fumaroles showed increased activity and a gas plume from the crater was usually observed on clear days. By early June SO₂ emissions yielded up to 5000 t/d. On July 20 a green lake was observed filling the crater of the volcano,
<figure><img src="image_95.png" /><figcaption>**Fig. 5** Hydrothermal unrest at Turrialba, Poás and Irazú volcanoes. a Vegetation die-back due to resumed plume degassing. b The 2010 phreatic eruption vent (picture by S. Calabrese). c Lack of vegetation at Cerro Pelón (Poás)</figcaption></figure>

nevertheless no significant changes in the SO₂
emission was observed, suggesting that the lake
was of high acidity and/or too small to be an
efficient scrubber. The first phreato-magmatic
explosions occurred on August 14 and produced
ash and gas columns reaching up to 9 km above
the crater. The satellite-borne instruments such as
OMI and OMPS reported 16,400 and 12,500 t/d
of SO₂ released to the atmosphere on that day
(http://so2.gsfc.nasa.gov/pix/daily/0815/ecuador_
0815z.html). Continuous ash emissions followed
the initial explosions producing a permanent gas
and ash plume drifting westward. SO₂ measured
in this permanent plume, by mobile-DOAS tra-
verses or by the permanent stations, reached

due to plume degassing and acid rain fall downwind.
d Efficient rock mass removal through Río Sucio, 30 km
downstream Irazú volcano

24,000 t/d and decreased progressively until the
end of the activity in late November 2015.

Since June and more consistently since 14
August 2015, BrO was also detected in the plume
(Dinger et al. 2016). Airborne Multi-GAS mea-
surements showed that the plume had a low CO2/
SO2 ratio, and that SO2 was >99% of total sulfur
(SO2 + H2S), indicating a shallow magmatic
origin for the gas. At the time of writing
(September 2016), SO2 emissions decreased to
background levels. SO2 permanent monitoring
proved to be a useful tool at Cotopaxi providing
real time data contributing, together with other
geophysical methods, to better evaluate volcanic
unrest scenarios.
## Gas Impact and Acid Rain

Persistent high-temperature fumarolic and plume degassing impact volcano flanks down-wind (Fig. 5c). This can happen when it rains through a plume, generating acid rain. Many volcanoes are subject to such long-term, non-eruptive hazard, making the surrounding ground harsh living environments. Other volcanoes suffer vegetation die-back, when activity resumes after prolonged quiescence (Fig. 5a). This was clearly visible during the increased activity at Turrialba volcano (Section “Turrialba, Costa Rica (2001–2016)”, Fig. 5a, b; González et al. 2015) and downwind Poás’ western flank (Cerro Pelón, Fig. 5c). To assess such hazards, a meteorological station (wind direction, speed, air humidity and rainfall), and a DOAS device to measure SO₂ fluxes are valuable tools.

## Rock Leaching upon Weathering

Absorption of magmatic gases into aquifers creates acidic magmatic-hydrothermal systems that, sooner or later, will exit the volcano. Acidic water attacks the wall rock and becomes loaded with solutes (Delmelle et al. 2015). If (1) meteoric recharge is high, (2) acid input is high, (3) wall rock is fresh, and (4) (thermal) spring water discharge is high, rock leaching capacity can reach thresholds to even mechanically destabilize volcano flanks. Enhanced chemical leaching for long periods favors physical rock removal, causing rock fall, landslides, or even flank and sector collapses of volcanic edifices. Even if magmatic degassing seems absent and a volcano may be long dormant, hazards loom due to the scrubbing capacity of deeper aquifers (Delmelle et al. 2015, and references therein).

The best suited method to quantify and track rock mass removal from a volcanic edifice is by monitoring the discharge rates from thermal springs, and their dissolved solutes and solids. The “Cl-inventory” uses Cl as the conservative tracer. As mentioned earlier, this is true for less

acidic magmatic-hydrothermal systems (Ingebritsen et al. 2001; Taran and Peiffer 2009; Chiodini et al. 2014; Collard et al. 2014). Measuring the Cl-release from rivers draining thermal springs, and knowing the Cl-content and Cl/solute ratios in thermal spring waters, the rock mass removal rate can be estimated by:

$$Q_r * C_r = Q_s * C_s \quad (7)$$

where Qr and Qs are the discharges of rivers and thermal springs, respectively, in L/s; Cr and Cs are the concentrations of Cl of rivers and thermal springs, respectively, in mg/L. Measuring the river discharge (Qr) and analyzing the river and spring waters for its Cl content (Cr and Cs, respectively) thus enables to estimate the spring discharge (Qs), which would otherwise be impossible to directly measure in the field (e.g. numerous spring discharges, too irregular spring mouths). This method combines gas-water-rock interaction and hydrology of magmatic-hydrothermal systems in order to assess indirect hazard. Volcanoes with high rock mass removal rates are e.g. Irazú (Costa Rica, Fig. 5d); extremely acidic magmatic-hydrothermal rock removers are Kawah Ijen (Java), Poás (Costa Rica), and Copahue (Argentina-Chile). In the most extreme cases, the acidity and toxic metal load affects agricultural activities and human health.

Moreover, through the same Cl-inventory approach, the geothermal potential (i.e. heat output) from springs can be estimated, by multiplying the enthalpy of discharged spring waters, often based on geothermometric temperatures of the deep system, with the spring discharge rates. Such estimates were obtained for the active magmatic-hydrothermal systems of El Chichón and Tacaná (both in Chiapas, Mexico; Taran and Peiffer 2009; Collard et al. 2014), and Domuyo (Argentina; Chiodini et al. 2014), and originally of the Cascades Volcanic Range by Ingebritsen et al. (2001). Understanding the state of unrest on the long term of a specific volcano is needed to rule out if the volcano would be a feasible target for geothermal exploitation, or not.
# Volcanic Lakes

## Acid Peak-Activity Lakes in a State of Unrest

Volcanic lakes are the intersection of the crater surface and the underlying aquifer (Fig. 1g, h). Hence, they become “windows” into the depths of magmatic-hydrothermal systems. While fumaroles directly lose their signal from depth to the atmosphere as a “snapshot” (but see Section “When the Gas Hits the Water”), volcanic lakes preserve a past gas marker for a certain period. The time we can track back by studying the water chemistry depends on the duration the water resides in the lake (i.e. residence time). The residence time (RT) is estimated by dividing the lake volume (V in m³) by the input or output rate (Q in m³/s) of fluids, assuming steady-state conditions (Varekamp 2003; Rouwet et al. 2014b):

$$RT = V/Q \qquad (8)$$

Small lakes with high rates of fluid flushing offer a better time-resolution than large lakes with slow rates of fluid in- and output. The monitoring frequency (e.g. lake water sampling) should be tuned to this residence time; a higher monitoring frequency will oversample the lake chemistry, while a lower monitoring frequency will lead to the loss of information the lake potentially provides (Rouwet et al. 2014b).

Volcanic lakes are excellent gas scrubbers. Nevertheless, recent studies quantify the gas release from the lake surface of the most acidic lakes (e.g. Aso, Copahue, Poás, Kawah Ijen; Shinohara et al. 2015; Tamburello et al. 2015; Capaccioni et al. 2016; de Moor et al. 2016; Gunawan et al. 2016). Under the most extreme pH conditions (<0) HF, HCl and even SO₂ can degas freely from the lake. This means that acidic lakes are more sensitive than thought before, as acid gas flashes through the water body with only minor scrubbing. Monitoring frequency can thus increase and, hence, spectroscopic or electrochemical sensor tools become extremely useful (DOAS to measure SO₂ fluxes from volcanic

plumes, Multi-GAS to measure ratios between gas species, e.g. CO₂/SO₂).

Considering acidic lakes as “open-air” fumaroles (and fumaroles as “buried” acidic lakes) has changed the monitoring time frame, which might lead to better chances to forecast phreatic eruptions (Rouwet et al. 2014a, b; Rouwet and Morrissey 2015).

## Volcanic Lake Response to External Triggers in the Absence of Magmatic Unrest

The Boiling Lake, in Dominica of the Lesser Antilles (a VUELCO target volcano), is a high-temperature volcanic crater lake, that is believed to be formed as a result of phreatic or phreatomagmatic activity (Lindsay et al. 2005). It is approximately 50 m × 60 m in size and ca. 12–15 m deep, with an estimated volume of ~1.2 × 10⁴ m³ when filled (Fournier et al. 2009). Over the last 150 years, temperatures taken at the edge of the lake have generally ranged between 80 and 90 °C, and the pH from 4 to 6 (Joseph et al. 2011), however, the lake experiences periods of instability where water level, temperature and state of hydrothermal activity fluctuate suddenly. The long term stability of the Boiling Lake is attributed to the hypothesis that the lake is suspended above the water table by the buoyancy of rising steam bubbles from the underlying hydrothermal system (Fournier et al. 2009). Perturbations of this interaction as a result of landslides into the lake or regional seismicity is attributed to the cause of the sudden periods of instability (i.e. lake drainage and refill cycles), rather than changes in underlying magmatic activity. This is mainly because no coincident anomalous hydrothermal activity was observed in the Valley of Desolation geothermal area, which is located nearby the Boiling Lake.

The most recent episode of instability occurred during December 2004 to April 2005 (Fig. 6), when the lake water level dropped by ~8–10 m, and the temperature at the water’s edge decreased to <30 °C. Water acidity went from the usual acidic pH of 4–6 to neutral, while
<figure><img src="image_96.png" /><figcaption>**Fig. 6** The refill of an empty Boiling Lake during a period of unusual activity where water levels and geothermal activity were rapidly fluctuating. **a** An empty Boiling Lake on 7 April 2005. **b** A full Boiling Lake on 13 April 2005. Pictures by Arlington James, Forestry Officer, Forestry and Wildlife Division, Dominica (used with permission)</figcaption></figure>

Cl concentration dropped from the typical 2000–6000 to 29–50 mg/L, and SO₄ concentration dropped from 1500–4000 to 100–270 mg/L, indicating a drastic decrease of hydrothermal fluid input into the lake (Joseph et al. 2011). Additionally, total dissolved solid content decreased from 13,400 to 4500 mg/L, suggesting strong dilution by fresh water. Measurements of temperature, pH and chemical composition taken in August 2006 indicate that the lake had

returned to its normal steady-state of activity (Joseph et al. 2011). This episode is reported to have been triggered as a result of the extensional strain induced by a regional Mw 6.3 earthquake that occurred on November 21, 2004 offshore Les Saintes (Guadeloupe) (Feuillet et al. 2011), which may have contributed to diminished water inflow.

It should be noted, however, that a phreatic explosion and gas release occurred at an “empty”
<table><caption>**Table 2** Geochemical signals and what they indicate with respect to volcanic unrest</caption><thead><tr><th>Type of unrest</th><th>Geochemical signal</th><th>Indication</th></tr></thead><tbody><tr><td rowspan="4">Unrest</td><td>CO<sub>2</sub> flux above background</td><td>Changes in deep degassing dynamics</td></tr><tr><td>Increase in T of hot springs and/or fumaroles</td><td>Increased heat input</td></tr><tr><td>Changes in H<sub>2</sub>O/CO<sub>2</sub> ratios in fumaroles</td><td>Changes in water/gas ratio</td></tr><tr><td>Appearance of new fumaroles and/or hot springs</td><td>Aerial extension of activity</td></tr><tr><td rowspan="5">Magmatic unrest</td><td>Appearance of acidic gases (SO<sub>2</sub>, HCl, HF)</td><td>Changes in mid- to shallow magma dynamics</td></tr><tr><td>T fumarole &gt;119 °C</td><td>Remobilisation of sulphur</td></tr><tr><td>SO<sub>2</sub> flux &gt; X t/d</td><td>SO<sub>2</sub> flux above background, volcano-dependent</td></tr><tr><td>Increase in CO<sub>2</sub>/SO<sub>2</sub> ratio</td><td>Arrival of an undegassed magma at depth</td></tr><tr><td>Extreme increase in T fumaroles (&gt;300 °C)</td><td>Towards magmatic T</td></tr><tr><td rowspan="4">Magmatic eruption</td><td>Decreasing CO<sub>2</sub>/SO<sub>2</sub> ratios after increase</td><td>More superficial magma degassing</td></tr><tr><td>Increase in Cl, Br, F concentrations in hot springs/pools</td><td>Input of highly soluble acidic gases</td></tr><tr><td>Decrease in H<sub>2</sub>O/CO<sub>2</sub> and/or H<sub>2</sub>S/SO<sub>2</sub> and/or SO<sub>2</sub>/HCl ratios</td><td>More gas with a more magmatic signature</td></tr><tr><td>New fumaroles</td><td>Aerial extension of activity</td></tr><tr><td rowspan="6">Hydrothermal unrest</td><td>Anomalous glacier defrosting</td><td>Sudden removal of water mass... lahars</td></tr><tr><td>Water to vapour transition</td><td>Pushing vapour front from below</td></tr><tr><td>Changes in hydrothermal features</td><td>Variations or aerial extension of activity</td></tr><tr><td>Increase in B and/or NH<sub>4</sub> in waters</td><td>Input of vapour</td></tr><tr><td>Increase in CH<sub>4</sub>/CO<sub>2</sub> in fumaroles</td><td>A more hydrothermal signature in fumaroles</td></tr><tr><td>Variations in phreatic level in aquifers</td><td>Pushing vapour front from below</td></tr><tr><td rowspan="4">Hydrothermal eruption</td><td>120 °C &lt; T fumarole &lt;200 °C</td><td>Self-sealing by a change in S viscosity</td></tr><tr><td>Extension of alteration areas or fumarolic fields</td><td>Aerial extension of activity</td></tr><tr><td>Appearance of muddy pools</td><td>Clearing bugs and vents, unplugging</td></tr><tr><td>Boiling/bubbling of pools that previously didn't</td><td>Rising vapour front and/or extra heating and degassing</td></tr></tbody></table>

Boiling Lake on 10 December 1901 that resulted in the deaths of two individuals (Elliot 1938; Bell 1946). This suggests that hazards related to volcanic lakes such as the Boiling Lake, may occur without magmatic input.

## Take-Home Ideas: Implications for Geochemical Monitoring

Over the past five-six decades, gas geochemistry at magmatic-hydrothermal systems has mainly focussed on chemical equilibria and kinetics in

the subsurface environment. Over the last 10–15 years more attention has been paid to remote sensing of volcanic gas plumes (DOAS, Multi-GAS) with the obvious advantage of increased safety and frequency of data gathering. Nevertheless, the relationship between the fumarole and plume has yet to be better constrained. The best proxy of a magmatic gas remains a direct sample of a high-temperature fumarole, although such target fumaroles are often inaccessible, especially during eruptive phases. Compromises between data fidelity and safety of the operators, and the frequency of data
gathering should be framed in terms of what we want and maximally can unravel. The advantage of fluid geochemistry in volcano monitoring arrises from the fact that volatiles are mobile and thus reach the surface often before physical changes manifest. Timely recognition of unrest and especially hydrothermal unrest is often possible. Table 2 summarizes geochemical signals and how they relate to the various states of volcanic unrest, useful for monitoring based on deterministic research and probabilistic modeling.

Future research should focus on better constraining degassing dynamics at the surface-atmosphere boundary. Recent studies have demonstrated that extremely acidic aquifers can “slow down” the signal released from gas-dominated magmatic-hydrothermal systems (fumaroles), but “speed up” the signal released from water-dominated systems (e.g. acidic crater lakes). These findings have strong implications for the time frame of magma degassing, and hence for the monitoring frequency.

## References

Aiuppa A, Federico C, Paonita G, Valenza M (2002) S, Cl and F degassing as an indicator of volcanic dynamics: the 2001 eruption of Mount Etna. Geophys Res Lett 29(11):1559. doi:10.1029/2002GL015032

Aiuppa A, Federico C, Giudice G, Giuffrida G, Guida R, Gurrieri S, Liuzzo M, Moretti R, Papale P (2009) The 2007 eruption of Stromboli volcano: insights from real-time measurements of the volcanic gas plume CO₂/SO₂ ratio. J Volcanol Geoth Res 182(221):230. doi:10.1016/j.jvolgeores.2008.09.013

Baldoni E, Rouwet D, Mora-Amador R, Ramírez C, González-Ilama G, Lucchi F, Capaccioni B, Tranne CA, Pecoraino G (submitted) Hydrogeochemical model of the Irazú-Turrialba volcanic complex (Costa Rica) and implications for hazard assessment and volcanic surveillance. In: Caudron C, Capaccioni B, Ohba T (eds) GSL special volume geochemistry and geophysics of volcanic lakes

Bell H (1946) Glimpses of a governor's life—Dominica. Samson Low, London

Campion R, Martínez-Cruz M, Lecocq T, Caudron C, Pacheco J, Pianrdi G, Hermans C, Carn S, Bernard A (2012) Space- and ground-based measurements of sulphur dioxide emissions from Turrialba Volcano (Costa Rica). Bull Volcanol. doi:10.1007/s00445-012-0631-z

Capaccioni B, Rouwet D, Tassi F (2016) HCl degassing from extremely acidic crater lakes: empirical results from experimental determinations and implications for geochemical monitoring. In: Caudron C, Capaccioni B, Ohba T (eds) GSL special publications 437 geochemistry and geophysics of volcanic lakes. doi:10.1144/SP437.12

Capasso G, Favara R, Inguaggiato S (1997) Chemical features and isotopic composition of gaseous manifestations on Vulcano Island, Aeolian Islands, Italy: an interpretative model of fluid circulation. Geochim Cosmochim Acta 61(16):3425-3440

Chiodini G, Liccioli C, Vaselli O, Calabrese S, Tassi F, Caliro S, Caselli A, Agusto M, D'Alessandro W (2014) The Domuyo volcanic system: an enormous geothermal resource in Argentine Patagonia. J Volcanol Geoth Res 272:71-77. doi:10.1016/j.jvolgeores.2014.02.006

Dinger F, Arellano S, Battaglia J, Bobrowski N, Galle B, Hernández S, Hidalgo S, Hörmann C, Lübcke P, Platt U, Ruíz M, Warnach S, Wagner T (2016) Variations of the BrO/SO₂ molar ratios during the 2015 Cotopaxi eruption. EGU 2016-1001

Collard N, Taran Y, Peiffer L, Campion R, Jacome Paz MP (2014) Solute fluxes and geothermal potential of Tacana volcano-hydrothermal system, Mexico-Guatemala. J Volcanol Geoth Res 288:123-131. doi:10.1016/j.jvolgeores.2014.10.012

De la Cruz-Reyna S, Tilling RI (2008) Scientific and public responses to the ongoing volcanic crisis at Popocatépetl Volcano, Mexico: importance of an effective hazards-warning system. J Volcanol Geoth Res 170:121-134

Delmelle P, Bernard A (2015) The remarkable chemistry of sulfur in hyper-acid crater lakes: a scientific tribute to Bokuichiro Takano and Minoru Kusakabe. In: Rouwet D, Christenson B, Tassi F, Vandemeulebrouck J (eds) Book chapter in volcanic lakes, Springer, Heidelberg, pp 239-260. doi:10.1007/978-3-642-36833-2_10

Delmelle P, Henley RW, Opfergelt S, Detienne M (2015) Summit acid crater lakes and flank instability in composite volcanoes. In: Rouwet D, Christenson B, Tassi F, Vandemeulebrouck J (eds) Book chapter in volcanic lakes. Springer, Heidelberg, pp 289-306. doi:10.1007/978-3-642-36833-2_12

de Moor JM, Aiuppa A, Pacheco J, Avard G, Kern C, Liuzzo M, Martinez M, Giudice G, Fischer TP (2016) Short-period volcanic gas precursors to phreatic eruptions: insights from Poás Volcano, Costa Rica. Earth Planet Sci Lett 442:218-227. doi:10.1016/j.epsl.2016.02.056

Elliot EC (1938) Boiling lake—the 1900 story. In: Broken Atoms. Unpublished report presented to the Government of Dominica

Feuillet N, Beauducel F, Jacques E, Tapponnier P, Delouis B, Bazin S, Vallée M, King GCP (2011) The Mw = 6.3, November 21, 2004, Les Saintes earthquake (Guadeloupe): Tectonic setting, slip model
and static stress changes. J Geophys Res 116(B10):1-25. doi:10.1029/2011JB008310

Fournier N, Withal F, Moreau-Fournier M, Bardou L (2009) The Boiling Lake of Dominica, West Indies: high temperature volcanic crater lake dynamics. J Geophys Res 114:B02203. doi:10.1029/2008JB005773

Giggenbach WF (1987) Redox processes governing the chemistry of fumarolic gas discharges from White Island, New Zealand. Appl Geochem 2:143-161

González G, Mora-Amador R, Ramirez C, Rouwet D, Alpizar Y, Picado C, Mora R (2015) Actividad historica y analysis de la amenaza del volcano Turrialba, Costa Rica. Rev Geol Am Centr 52:129-149. doi:10.15517/rgac.v0i52.19033

Gunawan H et al (2016) New insights into Kawah Ijen's volcanic system from the wet volcano workshop experiment. In: Caudron C, Capaccioni B, Ohba T (eds) GSL special volume geochemistry and geo-physics of volcanic lakes

Hall M, Mothes P (2008) The rhyolitic-andesitic eruptive history of Cotopaxi volcano, Ecuador. Bull Volcanol 70:675-702. doi:10.1007/s00445-007-0161-2

Ingebritsen SE, Galloway DL, Colvard EM, Sorey ML, Mariner RH (2001) Time variation of hydrothermal discharge at selected sites in the western United States: implications for monitoring. J Volcanol Geoth Res 111:1-23

Joseph EP, Fournier N, Lindsay J, Fischer T (2011) Gas and water geochemistry of geothermal systems in Dominica, Lesser Antilles island arc. J Volcanol Geoth Res 206(1-2):1-14. doi:10.1016/j.volgeores.2011.06.007

Kusakabe M, Komoda Y, Takano B, Abiko T (2000) Sulfur isotopic effects in the disproportionation reaction of sulfur dioxide in hydrothermal fluids: implications for the d34S variations of dissolved bisulphate and elemental sulfur from active crater lakes. J Volcanol Geoth Res 97:287-307

Lindsay J, Robertson R, Shepherd J, Ali S (2005) Volcanic hazard atlas of the Lesser Antilles. In: 279. St. Augustine: Seismic Research Unit, University of the West Indies

Martini F, Tassi F, Vaselli O, Del Potro R, Martínez M, Van der Laat R, Fernandez E (2010) Geophysical, geochemical and geodetical signals of reawakening at Turrialba volcano (Costa Rica) after 150 years of quiescence. J Volcanol Geoth Res 198:416-432. doi:10.1016/j.jvolgeores.2010.09.021

Menyailov IA, Nikitina LP, Shapar VN, Pilipenko VP (1986) Temperature increase and chemical changes of fumarolic gases at Momotombo volcano, Nicaragua, in 1982-1985: are these indicators of a possible eruption? J Geophys Res 187:12199-12214

Mora-Amador R, Ramírez-Umaña C, González G, Rouwet D, Lucchi F, Forni F, Sulpizio R, Baldoni E, Alpízar-Segura Y, Tranne CA (2015) First documentation of the ongoing phreatic-strombolian eruptions of Turrialba volcano (Costa Rica). IUGG General Assembly, Prague 3746

Moussallam Y, Peters N, Ramirez C, Oppenheimer C, Aiuppa A, Giudice G (2014) Characterisation of the magmatic signature in gas emissions from Turrialba Volcano, Costa Rica. Solid Earth 5:1341-1350. doi:10.5194/se-5-1341-2014

Oppenheimer C, Stevenson D (1989) Liquid sulphur lakes at Poás volcano. Nature 342:790-793

Phillipson G, Sobradelo R, Gottsmann J (2013) Global volcanic unrest in the 21st century: an analysis of the first decade. J Volcanol Geoth Res 264:183-196

Rouwet D, Morrissey M (2015) Mechanisms of crater lake eruptions: physical and numerical modeling. In: Rouwet D, Christenson B, Tassi F, Vandemeulebrouck J (eds) Book chapter volcanic lakes. Springer, Heidelberg, pp 73-92. doi:10.1007/978-3-642-36833-2_3

Rouwet D, Ohba T (2015) Iotope fractionation and HCl partitioning during evaporative degassing from active crater lakes. In: Rouwet D, Christenson B, Tassi F, Vandemeulebrouck J (eds) Book chapter volcanic lakes. Springer, Heidelberg, pp 179-200. doi:10.1007/978-3-642-36833-2_7

Rouwet D, Sandri L, Marzocchi W, Gottsmann J, Selva J, Tonini R, Papale P (2014a) Recognizing and tracking hazards related to non-magmatic unrest: a review. J Appl Volcanol 3:17. doi:10.1186/s13617-014-0017-3

Rouwet D, Tassi F, Mora-Amador R, Sandri L, Chiarini V (2014b) Past, present and future of volcanic lake monitoring. J Volcanol Geoth Res 272:78-97. doi:10.1016/j.jvolgeores.2013.12.009

Rouwet D, Mora-Amador R, Ramírez-Umaña C, González G, Inguaggiato S (2016) Dynamic fluid recycling at Laguna Caliente (Poás, Costa Rica) before and during the 2006-ongoing phreatic eruption cycle (2005-2010). In: Caudron C, Capaccioni B, Ohba T (eds) Geochemistry and geophysics of volcanic lakes. Geological Society London Special Publications. doi:10.1144/SP437.11

Shinohara H, Giggenbach WF, Kazahaya K, Hedenquist JW (1993) Geochemistry of volcanic gases and hot springs of Satsuma-Iwojima, Japan: following Matsuo. Geochem J 27:271-285

Shinohara H, Yoshikawa S, Miyabuchi Y (2015) Degasbing activity of a volcanic crater lake: volcanic plume measurements at the Yudamari crater lake, Aso volcano, Japan. In: Rouwet D, Christenson B, Tassi F, Vandemeulebrouck J (eds) Book chapter volcanic lakes. Springer, Heidelberg, pp 73-92. doi:10.1007/978-3-642-36833-2_8

Symonds RB, Gerlach TM, Reed MH (2001) Magmatic gas scrubbing: implications for volcano monitoring. J Volcanol Geoth Res 108:303-341

Tamburello G, Agusto M, Caselli A, Tassi F, Vaselli O, Calabrese S, Rouwet D, Capaccioni B, Di Napoli R, Cardellino C, Chiodini G, Bitetto M, Brusca L, Bellomo S, Aiuppa A (2015) Intense magmatic degassing through the lake of Copahue volcano, 2013-2014. J Geophys Res doi:10.1002/2015JB012160
Taran YA, Hedenquist JW, Korzhinsky M, Tkachenko SI, Shmulovich KI (1995) Geochemistry of magmatic gases from Kudryavy volcano, Iturup, Kuril Islands. *Geochim Cosmochim Acta* 59:1749–1761

Taran YA, Peiffer L (2009) Hydrology, hydrochemistry and geothermal potential of El Chichón volcano-hydrothermal system, Mexico. *Geothermics* 38:370–378

Taran YA, Pokrovsky BG, Dubik YM (1989) Isotopic composition and origin of water from andesitic magmas. *Dokl Acad Sci* 304:440–443

Varekamp JC (2003) Lake contamination models for evolution towards steady state. *J Limnol* 62(1):67–72

Vaselli O, Tassi F, Duarte E, Fernández E, Poreda RJ, Delgado Huertas A (2009) Evolution of fluid geochemistry at the Turrialba volcano (Costa Rica) from 1998 to 2008. *Bull Volcanol*. doi:10.1007/s00445-009-0332-4

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

<figure><img src="image_97.png" /></figure>
<figure><img src="image_98.png" /></figure>

# Geophysical Footprints of Cotopaxi's Unrest and Minor Eruptions in 2015: An Opportunity to Test Scientific and Community Preparedness

Patricia A. Mothes, Mario C. Ruiz, Edwin G. Viracucha, Patricio A. Ramón, Stephen Hernández, Silvana Hidalgo, Benjamin Bernard, Elizabeth H. Gaunt, Paul Jarrín, Marco A. Yépez and Pedro A. Espín

## Abstract

Cotopaxi volcano, Ecuador, experienced notable restlessness in 2015 that was a major deviation from its normal background activity. Starting in April and continuing through November 2015 strong seismic activity, infrasound registry, hikes in SO₂ degassing and flank deformation with small displacements were some of the geophysical anomalies that were registered. Obvious superficial changes, such as small hydromagmatic eruptions, emission of vapor and ash columns, thermal hotspots around the crater and in nearby orifices and exacerbated glacier melting were also observed. Our contribution provides an overview of the 2015 Cotopaxi unrest by presenting the patterns of geophysical data and the sequence of events produced by the volcano. Cotopaxi's last important VEI 4 eruption was in 1877. Then it had devastating effects because of the transit of huge lahars down 3 major drainages. Comparatively, the 2015 activity never surpassed a magnitude VEI 2 and principally produced limited hydro-magmatic explosions and semi-continuous low energy emissions and light ashfalls. Given the potential of major destruction from a large Cotopaxi eruption it is important to understand the geophysical fingerprints that characterized the 2015 episode with an eye to identifying onset of future restless periods. Overall, the monitoring activities, the data interpretation, formulation of reasonable eruptive scenarios, and finally, the preparation of a stream of constant information being relayed to concerned authorities and the public, was a real test of the IGEPN's capacity to deal with a complicated eruption situation whose outcome was not apparent at the beginning, but which concluded in a very small eruptive episode.

P.A. Mothes (✉) · M.C. Ruiz · E.G. Viracucha ·

P.A. Ramón · S. Hernández · S. Hidalgo · B. Bernard ·

E.H. Gaunt · P. Jarrín · M.A. Yépez · P.A. Espín

Instituto Geofísico, Escuela Politécnica Nacional,
Quito, Ecuador

e-mail: pmothes@igepn.edu.ec
## Resumen

En 2015 el volcán Cotopaxi, Ecuador experimentó un notable cambio, que fue una desviación importante de su actividad normal de base. A partir de abril y hasta noviembre de 2015 fuerte actividad sísmica, registros de infrasonido, incremento en la desgasificación de SO₂ y pequeños cambios en la deformación de los flancos fueron algunas de las anomalías geofísicas registradas. Evidentes cambios superficiales también fueron observados como pequeñas erupciones hidromagmaticas, emisión de vapor, columnas de cenizas, puntos calientes alrededor del cráter y el deshielo de los glaciares. Nuestra contribución proporciona una visión general de las anomalías del Cotopaxi en el 2015, mediante la presentación de patrones de los datos geofísicos y la secuencia de eventos producidos por el volcán. La última erupción importante del Cotopaxi fue un VEI 4 en 1877. Esta tuvo efectos devastadores debido al descenso de enormes lahares por sus tres drenajes mayores. Comparativamente, la actividad del año 2015 nunca superó una magnitud VEI 2, principalmente produciendo explosiones hidromagmaticas, escasas emisiones y leves caídas de ceniza. Debido a la potencial destrucción por una eventual erupción grande del Cotopaxi es importante entender los registros geofísicos que caracterizó el episodio de 2015 para poder identificar el inicio de futuros períodos eruptivos. En general, las actividades de vigilancia, la interpretación de datos, formulación de escenarios eruptivos razonables y por último, la preparación de un flujo de información constante que llegue a las autoridades interesadas y el público, fue una verdadera prueba de la capacidad del IGEPN para hacer frente a una situación de erupción cuyo resultado no era evidente al principio, pero que finalizó como una erupción pequeña.

## Keywords

Volcanic unrest • Precursory geophysical patterns • Precursory LPs and VLPs • Volcano monitoring • Cotopaxi volcano-Ecuador • State of preparedness

## Introduction

Long dormant volcanoes that begin to awaken may have start and stop activity that has to be evaluated with respect to the volcano's known past. A volcano's past activity is known by study of its stratigraphy, other physical evidence and possibly historical chronicles (Tilling 1989). Many uncertainties preclude knowing the final outcome of a restless volcano (Newhall 2000; Sparks and Aspinall 2004), since a volcano may

awaken for short term low-level activity, then resume repose until additional magma inputs herald an episode of further unrest (Phillipson et al. 2013). Eruptions that barely bring magma to the surface may be classified as “failed”, since so little juvenile magma erupts (Moran et al. 2011).

At Cotopaxi ample geological and historical information exists with regard to its past activity (Hall and Mothes 2008; Garrison et al. 2011; Pistolesi et al. 2012). Formulation of eruptive
scenarios with respect to the 2015 unrest period were based on our collective knowledge of the volcano's geology and eruptive history and published information as well as interpretation of the abundant geophysical data streams available through instrumental and observational networks operated by the Instituto Geofísico of the Escuela Politécnica Nacional (IGEPN)-Quito, Ecuador, the entity in charge of volcano and tectonic monitoring in Ecuador. The combination of these inputs allowed scientists at the IGEPN to transmit a coherent image of the evolving unrest presented during 2015 and to indicate the most likely eruption/activity scenarios. Two earlier unrest periods are known: 1975-1976, when the IGEPN had limited seismic equipment operating on the volcano and then in 2001-2002. Both periods were comprised of increased fumarolic activity both inside and outside of the crater and a hike in seismicity for the 2001 period (Molina et al. 2008). The 2015 unrest displayed important changes in seismic and deformation patterns, gas output and superficial activity, compared to background, whose level was established since around 1990. In sum, Cotopaxi's 2015 unrest displayed a progressive crescendo of geophysical signals, then minor hydromagmatic explosions, followed by overall seismic energy decrease at the end of 2015, which was accompanied by fewer superficial manifestations. Like many other volcanoes that have displayed unrest, this recent episode did not culminate in a full-fledged eruption with large volumes of juvenile pyroclastics (Phillipson et al. 2013). Moran et al. (2011) maintain that an eruption is “failed” when magma reaches but does not pass the “shallow intrusion” stage, i.e., the magma gets close to, but does not reach the surface. In the actual case, the amount of erupted material was minor, and had a dense rock equivalent volume of ~0.5 Mm³ (Bernard et al. 2016).

## Cotopaxi Volcano

Cotopaxi volcano, located in central Ecuador atop the Eastern Cordillera, is a large, symmetrical stratocone with a basal diameter of 18 km

and an altitude of 5897 m (Hall and Mothes 2008). Its actual glacier cap of 10.49 km² is rapidly diminishing due to climatic change Cáceres et al. (2016) (Fig. 1). The volcano's last important VEI = 4 eruption was on 26 June 1877. Then it generated highly erosive pyroclastic flows that melted glaciers and triggered voluminous lahars (~100 Mm³ per drainage). Each lahar traveled hundreds of kilometers down several drainages enroute to the Amazon basin, Pacific Ocean and to the Atlantic, respectively (Mothes et al. 2004; Mothes and Vallance 2015). These past lahar routes now host sprawling suburbia, important economic activities and vital infrastructure. Ecuador's second most visited national park (Parque Nacional Cotopaxi-PNC) is centered on the volcano and draws some 200,000 tourists a year.

The volcano's five most important eruptive episodes during the historical period (since 1532) have been of andesitic composition and ranged between VEI = 3 and 4 (Pistolesi et al. 2012). Nonetheless, the volcano is bi-modal and produces VEI = 5 rhyolitic eruptions about every 2000 years (Hall and Mothes 2008). The last important rhyolitic eruption, the Peñas Blancas event, occurred about 2800 years BP (Mothes et al. 2015a). The youngest andesitic eruptive products contain intergrowths of plagioclase and pyroxene and four different populations of plagioclase crystals which indicate pervasive magma mixing (Garrison et al. 2011).

Given the high probability for the generation of long-distance lahars, wide dispersal of pyroclastic fall, and the consequential negative impact on many economic activities and the compromise of critical infrastructure should an eruptive period last for months to years, Cotopaxi is considered a “National” volcano, located in the center of Ecuador, near to Quito, the capital and other populated areas. Even a short-lived VEI 4 eruption (VEI = Volcano Explosivity Index) (Newhall and Self 1982) has the potential to gravely affect Ecuador's overall productivity and major transport lines. Lastly, the volcano is considered one of the most dangerous in Ecuador, given the possible exposure of >300,000 residents to primary lahars and ashfalls during
<figure><img src="image_99.png" /><figcaption>**Fig. 1** Location of Cotopaxi volcano in Ecuador’s Inter-Andean Valley and major cities</figcaption></figure>