Signs of magma ascent in LP and VLP seismic events and link to degassing: An example from the 2010 explosive eruption at Merapi volcano, Indonesia

Journal of Volcanology and Geothermal Research 261 (2013) 171–192

Contents lists available at SciVerse ScienceDirect

Journal of Volcanology and Geothermal Research journal homepage: www.elsevier.com/locate/jvolgeores

Signs of magma ascent in LP and VLP seismic events and link to degassing: An example from the 2010 explosive eruption at Merapi volcano, Indonesia Philippe Jousset a, b,⁎, 1, Agus Budi-Santoso c, d, 2, Arthur D. Jolly e, 3, Marie Boichu f, g, 4, 5, Surono h, 6, S. Dwiyono c, Sri Sumarti c, Sri Hidayati h, Pierre Thierry a a

BRGM, 3 Avenue Claude Guillemin, BP36009, 45060 Orléans Cedex 2, France Helmholtz Center GFZ, Telegrafenberg, 14473 Potsdam, Germany BPPTK (Balai Penyelidikan dan Pengembangan Teknologi Kegunungapian), Jalan Cendana 15, Yogyakarta 55166, Indonesia d ISTerre, CNRS, Université de Savoie, 73376 Le Bourget du Lac cedex, France e GNS Sciences, Wairakei Research Centre, Private Bag 2000, Taupo, 3352, New Zealand f The University of Cambridge, Department of Geography, Downing Place, Cambridge CB23EN, United Kingdom g Institut Pierre Simon Laplace, Laboratoire de Météorologie Dynamique, Ecole Polytechnique, 91128 Palaiseau cedex, France h Center of Volcanology and Geological Hazard Mitigation, Jalan Diponegoro 57, 40122 Bandung, Indonesia b c

a r t i c l e

i n f o

Article history: Received 15 March 2012 Accepted 15 March 2013 Available online 27 March 2013 Keywords: Volcano-seismology Long-period seismicity Complex frequency evolution Fluids and gas Eruption triggering VLP waveform inversion Merapi

a b s t r a c t The link between seismicity and degassing is investigated during the VEI 4 eruptions of Merapi volcano (Indonesia) in October and in early November 2010. Seismicity comprised a large number and variety of earthquakes including Volcano-Tectonic events, a sustained period of Long Period Seismicity (LPS), i.e., Long-Period events (LP), Very Long Period events (VLP) and tremor. LPS seismicity is ascribed to the excitation of fluid-filled cavity resonance and inertial displacement of fluids and magma. We investigate here LPS that occurred between 17 October and 4 November 2010 to obtain insights into the volcano eruption processes which preceded the paroxysmal phase of the eruption on 4–5 November. We proceed to the moment tensor inversion of a well-recorded large VLP event during the intrusion phase on 17 October 2010, i.e., before the first explosion on 26 October. By using two simplified models (crack and pipe), we find a shallow source for this VLP event at about 1 km to the south of the summit and less than 1 km below the surface. We analyze more than 90 LP events that occurred during the multi-phase eruption (29 October–4 November). We show that most of them have a dominant frequency in the range 0.2–4 Hz. We could locate 48 of those LP events; at least 3 clusters of LP events occurred successively. We interpret these observations as generated by different fluid-filled containers in the summit area that were excited while magma rose. We also observe significant variations of the complex frequency during the course of the eruption. We discuss these changes in terms of a variable ratio of fluid to solid densities and/or by possible conduit geometry change and/or permeability of the conduit or the edifice and/or by resonance of different fluid-containers during the release of more than 0.4 Tg of SO2 and large but unknown masses of other volcanic gases. Finally, we also discuss how the major explosions of the eruption were possibly triggered by passing waves resulting from regional tectonic earthquakes on 3 and 4 November. © 2013 Elsevier B.V. All rights reserved.

1. Introduction

⁎ Corresponding author. Tel.: +49 30 288 1299. E-mail addresses: [email protected] (P. Jousset), [email protected] (A. Budi-Santoso), [email protected] (A.D. Jolly), [email protected] (M. Boichu), [email protected] (Surono), [email protected] (P. Thierry). 1 Tel.: +49 30 288 1299. 2 Tel.: +62 274 514192. 3 Tel.: +64 7 374 8211. 4 Tel.: +44 1223 333386. 5 Tel.: +33 169335171. 6 Tel.: +62 22 727 26 06. 0377-0273/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.jvolgeores.2013.03.014

The understanding of eruption mechanisms is one of the main goals of volcano research. Volcano-seismology attempts to constrain key parameters of volcanic processes, such as the geometry of the magmatic conduit, the physical properties of fluids within the edifice and the dynamics of volcanic processes, including degassing and eruption. Seismic signals observed on volcanoes and hydrothermal systems are classified according to their frequency content (e.g., McNutt, 2000) and interpreted with the help of complementary geophysical and geochemical ground and satellite observations (Surono et al., 2012) and magma modeling (e.g., Neuberg and O'Gorman, 2002). Volcano-tectonic earthquakes (VT) contain frequencies ranging 2–50 Hz and are widely accepted to result from brittle fracture of

172

P. Jousset et al. / Journal of Volcanology and Geothermal Research 261 (2013) 171–192

rocks under the stress of volcanic activity, including fluid movements, stress changes, and temperature gradients evolution. They indicate zones where stress concentrates at depth and how those stress concentration evolves during eruptions (e.g., Roman et al., 2006; Hidayati et al., 2008). They do not indicate the dynamic state and properties of magmatic fluids, which are the driving force for eruptions. Long-Period Seismicity (LPS) include tremor, Long Period (LP) also called Low-Frequency (LF) or Multiphase (MP) and Very-Long Period (VLP) events. Their analysis and modeling have been one of the main focuses of volcano seismology for more than 35 years. LP events contain frequencies ranging from 0.2 to ~5 Hz. LP seismic models associate LP events with mechanical interactions between fluids and rocks in volcanic media. A widely accepted mechanism for LP events is the resonance of fluid-filled cavities (Aki et al., 1977; Chouet, 1985, 1986; Ferrazzini and Aki, 1987; Chouet, 1996a, 1996b; Neuberg, 2000; Jousset et al., 2003). Other source models have also been proposed for certain types of LP events, such as oscillations due to fluid-driven flow (Julian, 1994; Rust et al., 2008) or abrupt mass shifts of solidified domes, conduit magma or magma pads (Johnson et al., 2008) or stick–slip models (Harrington and Brodsky, 2007; Pallister et al., 2013a). Very-Long Period (VLP) events contain frequencies lower than 0.2 Hz; they are linked to inertial displacement of material such as magma or gas (Ohminato et al., 1998; Legrand et al., 2000; Chouet et al., 2005; Waite et al., 2008; Jolly et al., 2012), or pit-crater or cavern collapse (Jousset and Rohmer, 2012). In this work, we assume that LPS seismicity is associated with resonance: the LPS signals are caused by oscillations of fluid-filled conduits or cracks driven by magma or fluid movement. However, we do not rule out other mechanisms that may also be involved. We focus on LPS recorded before and during the M ~ 4 explosive 2010 eruption (Surono et al., 2012) at Merapi volcano. It is the largest eruption in more than a century (Hartmann, 1934, 1935). The digital seismic records during the 2010 eruption revealed all types of

a

earthquakes previously identified on this volcano (Shimozuru et al., 1969; Hidayat et al., 2000; Ratdomopurbo and Poupinet, 2000; Hidayat et al., 2002; Surono et al., 2012; Budi-Santoso et al., 2013; Ratdomopurbo et al., 2013; see also Fig. 1). Merapi stratovolcano is located 25–30 km north of the metropolitan area of Yogyakarta, Indonesia with a population of 1.6 million. It overlies the Java subduction zone (Luehr et al., 2013) and is composed mainly of basaltic–andesite tephra, pyroclastic flow, lavas, and lahar deposits (Camus et al., 2000; Newhall et al., 2000; Voight et al., 2000; Charbonnier et al., 2013; Costa et al., 2013; Cronin et al., 2013; Komorowski et al., 2013; Preece et al., 2013). Surono et al. (2012) distinguish four main phases during the 2010 explosive eruption; accurate study of deposits (Komorowski et al., 2013) indicate additional details of this chronology: (1) From 31 October 2009 to 26 October 2010, a new magmatic intrusion was revealed by inflation of the volcano, increased number of earthquakes (VT seismic swarms and VLP events), and changes in fumaroles gas composition. (2) A phreato-magmatic explosive phase started with the violent eruption of the 26 October, which generated 8-km-long pyroclastic density currents. The explosion was followed by the emission of a sustained white/gray eruptive plume and the occurrence of moderate-size explosive eruptions on 29 and 31 October; this phase ended with an increased number of LP events. (3) The magmatic phase started on 1 November with the rapid (up to 35 m3s−1) growth of a new summit lava dome observed from satellite imaging (Pallister et al., 2013b). Komorowski et al. (2013) proposed that the magmatic phase started as soon as 30 October. A swarm of LP events or explosions began suddenly on 3 November at about 4:05 UTC and ended with a large explosion on 3 November at 08:05 UTC, shortly after a regional M ~ 4.2 tectonic earthquake. The explosion generated a 12-km high

b

d

c

e

Fig. 1. Several types of seismic events recorded during 2010 eruption (vertical component at DEL station – see Fig. 2 for location) and their corresponding spectrograms. (a) Volcano-Tectonic (VT) earthquake, resulting from brittle rupture of rocks; frequency content 1–25 Hz. (b) Long-Period (LP) or Low-Frequency (LF) events, related to resonance of fluid-filled containers; frequency content 0.2–5 Hz. (c) Pyroclastic flow, corresponding to dome collapse: signal with cigar shape, large amplitude, long duration (up to more than 2 h at Merapi volcano) and saturating short-period sensors with high gain. (d) “Multiphase” (MP) event, also called “hybrid event”, corresponding to dome growth; frequency content 0.2–20 Hz; (e) “Guguran” corresponding to rock falls with shape similar to Pyroclastic flow signals, but lower amplitude and shorter duration. Other seismic events, e.g., tremor and VLP events are not shown here (See Fig. 4 and S2).

P. Jousset et al. / Journal of Volcanology and Geothermal Research 261 (2013) 171–192

The monitoring network at Merapi in 2010 comprised 4 shortperiod (Mark Products L-4 seismometers) permanent stations (PUS, KLA, DEL, and PLA, Fig. 2). Seismic monitoring and analysis were carried out in real time during the eruption and used qualitatively to infer magmatic and eruptive processes (Surono et al., 2012). As with previous eruptions at Merapi, most VT earthquakes hypocenters were localized in the upper part of the edifice below the 2006 dome in two separate clusters between 2.5 and 5 km, and 0 and 1.5 km below the summit (Budi-Santoso et al., 2013). In this work, we use data from the short-period monitoring network and from a temporary broadband (BB) network of 5 stations. All of these stations were telemetered to BPPTK using wireless transmission. The BB stations include: 1 STS2 Strekeisen seismometer (station LBH) and 4 Güralp CMG40T seismometers installed in July 2009 as part of the MIAVITA (MItigate and Assess risk from Volcanic Impact on Terrain and human Activities) European research project (Thierry et al., 2008). The BB network comprised stations GMR, GRW, PAS, WOR and operated from July 2009 to September 2010. The L56 station was operated beginning in September 2010. We established seismometers close to the summit (Fig. 2) to capture near-field effects of the volcano-seismic sources in the dome and conduit and to attempt to constrain accurately the source

600

900 1200 1500 1800 2100 2400 2700 3000

15 00

9170

Longitude UTM (km)

GRW

KLA

PAS 0

200

L56 WOR

2500

PUS

00 15

2. Monitoring of the 2010 seismicity at Merapi volcano

300

00

Prior to and simultaneously with the eruptions, an unprecedented number of earthquakes and an increased earthquake rate were recorded, in comparison with the previously instrumented eruptions of Merapi (Budi-Santoso et al., 2013; Ratdomopurbo et al., 2013). We briefly describe the seismic monitoring system that was used during the course of the eruption in Section 2. We draw inferences about the physical status, behavior and geometry of fluid-filled conduits prior to the eruption on the basis of moment tensor mechanism inversion of one large amplitude VLP event clearly recorded at 4 broadband stations and we discuss its significance in Section 3. In Section 4, we analyze the complex frequencies of more than 90 LP events preceding the main eruptive phases. Our study reveals variations in the complex frequencies which are interpreted as changes of (1) the conduit geometry and/or (2) fluid compositions and/or (3) permeability of the conduit or the edifice during the eruption and/or (4) changes in the resonator as the eruption proceeds. On the basis of these LPS observations, we propose mechanisms to explain the eruption dynamics, and we develop and discuss an integrated model of the eruption, which is based on influx of new magma and gas (Section 5) resulting in an unstable magmatic mixture (Hill et al., 2002) that was excited by small perturbations such as seismic waves from regional earthquakes (e.g., Moran et al., 2002; Davis et al., 2007; Cannata et al., 2010; Takekawa et al., 2013).

m 0

20

convective column and pyroclastic density currents that extended to 9 km from the summit. Multiple explosions with increasing intensity and gradually shorter delays between explosions occurred on 4 November (see Fig. 16), following another regional tectonic ML ~ 4 earthquake. The paroxysmal event produced a 17 km high eruptive column which disturbed seriously air traffic (Picquout et al., 2013) and a series of dome collapse and pyroclastic density currents that extended 16 km from the summit vent (Komorowski et al. 2013; Jenkins et al., 2013). This VEI ~4 (or M ≥4.0) explosion sequence was recorded on seismograms up to >45 km from the vent and contributed to the evacuation warning for the main explosion during the late night and early morning of 4–5 November, local time (Surono et al., 2012). (4) Lava dome growth stopped by 8 November and was followed by decreasing volcanic tremor interspersed with moderate to small explosions and emissions of ash and gas from vents that penetrated the newly formed lava dome. Lahars became the main post-eruption hazard at this stage (de Bélizal et al., 2013), as rainfall mobilized the new pyroclastic deposits.

173

GMR 9165 1500

1000

LBH

DEL

PLA

9160 435

440

Latitude UTM (km) Fig. 2. Simplified map (digital elevation model) of Merapi volcano used for the Finite Difference (FD) modeling. Location of seismic network stations during the 2010 Merapi eruption: empty triangles = short-period sensors; filled triangles = broadband stations. The black box indicates the inspected volume for wave form inversion (see Section 3). The range of the FD computation grid is latitude 430 km to 450 km and longitude 9155 km to 9175 km (UTM grid Zone 49) and from the surface down to 8 km depth and grid step of 150 m (table 2). Locations for best-fit inverted cylinder and crack are shown by the gray disk and the small gray rectangle, respectively.

processes of volcanic earthquakes (Hidayat et al., 2000; Lokmer et al., 2007); we also established a remote station (CRM at 46 km from the summit) to study possible link between seismic tectonic and volcanic activities. Summit BB seismometers (L56, PAS and GMR) were destroyed during the eruption. The remaining BB stations were included in the monitoring routine during the course of the eruption, after short-period sensors KLA, DEL and PUS were also destroyed by pyroclastic density currents on 4 November at 11:46:08, 17:07:03 and 21:25:07, respectively (Fig. 16; all times are cited in UTC unless otherwise noted). Removal of the instrumental response and calibration were performed before signal analysis (Jousset et al., 2011a). Technical problems including poor time synchronization at some stations (due to loss of GPS time signal) prevented a full analysis in real-time during the eruption (Surono et al., 2012). In this work, we resynchronized those records by adapting a cross-correlation technique (Sens-Schönfelder, 2008) to our network. The technique consists in determining clock differences between stations by pair-wise correlation of a random seismic wave field and then in calculating clock errors by a generalized inversion of the matrix of the cross-correlation and the clock difference (Vasseur, 2009; Jousset et al., submitted for publication; see BudiSantoso et al., 2013; Fig. S1; and Appendix A for details). Fig. 3 illustrates steps for clock corrections to resynchronize records on 17 October 2010, when a large VLP event occurred. We checked that the accuracy of

174

P. Jousset et al. / Journal of Volcanology and Geothermal Research 261 (2013) 171–192

a

14:16

14:17

14:18

14:19

GMR LBH

L56 PAS

PLA

DEL KLA PUSZ

b

14:16

14:17

14:18

14:19

GMR ~6s

LBH

~-6 s

~-170 s ~

L56 PAS

~ 110 s

~-61 61 s

PLA DEL KLA PUSZ

c

14:18

14:19

GMR LBH L56 PAS PLA DEL KLA PUSZ

Fig. 3. Signals recorded at several stations, including a Very Long Period event. (a) Before time synchronization and unfiltered. (b) before time synchronization and band-pass filtered between 0.02 and 0.15 Hz. Note that we observe time shift between records of about 6 s, 170 s and 61 s for stations GMR, LBH and L56, respectively. These rough values can be compared to the shifts retrieved from ambient noise correlation values (Table 6).

resynchronized signal is better than several tenths of ms for events on 17 October, by analyzing the standard deviation of cross-correlation

lags between stations having good GPS clock signals (see Appendix A and Fig. S1). This level of uncertainty is small enough to be used to

P. Jousset et al. / Journal of Volcanology and Geothermal Research 261 (2013) 171–192

pick phases and locate VT earthquakes (Budi-Santoso et al., 2013) and/ or to perform waveform inversion to constrain locations and source mechanism of volcano-seismic events, including VLP events. 3. Moment-tensor inversion of a precursory Very-Long Period event on 17 October 2010 3.1. Analysis of VLP seismicity and methodology More than 200 VLP signals (0.01–0.2 Hz) were recognized at the summit BB stations (L56 or WOR and PAS) during the intrusion phase before the 26 October eruption (Fig. S2). VLP events have been associated with MP events (as previously observed by Hidayat et al., 2000 at Merapi), and at other times associated with VT earthquakes. We noted swarms of VLP events that merged into VLP tremor on several occasions (e.g., Fig. S2c). Qualitative comparison between isolated VLP events reveals that two main waveform types can be recognized (Fig. 4) at the summit stations (L56 or WOR), both hidden

175

within higher frequency signals: the vertical component is either upward or downward. From 5 until 13 September 2010, we counted about 100 VLP events at summit station WOR, 56% of them had an upward first motion and 44% of them had a downward first motion. From 14 until 21 October 2010, we counted at about 90 VLP events at summit station L56, 63% of them displaying a downward first motion and 37% of them an upward first motion. This preliminary observation suggests that a limited number of trigger mechanisms from repeating process at stationary source locations produced the VLP signals and that the mechanism may have evolved with time during the intrusion. Most stations were too far from the summit area and did not record VLP waveforms. On few occasions, however, more powerful VLP events occurred and were recorded at several distal broadband stations. A strong VLP event occurred on 17 October 2010 at 14:18 and was recorded at all broadband stations (Fig. 3). It has a downward first motion at the summit station L56. Note that VLP signals are also recorded at short-period sensors. Because the resolution of VLP signals in records

15/10/2010 – 07:29 EW

a Frequency (Hz)

b

NS

Z 40 20 0

EW

NS

c

GAS RELEASE

Z 60 s

Filter = 0.02-0.2 Hz

21/10/2010 – 00:50 EW

a

NS

b

Frequency (Hz)

Z 40 20 0

EW

MAGMA INTRUSION

NS

c Z 60 s Filter = 0.02-0.2 Hz Fig. 4. Two types of Very Long period events recorded at the summit station L56. (a) and (d) Unfiltered signals for the 3 components of each VLP type. (b) and (e) Corresponding spectrograms (c) and (b) Same signals filtered in 0.01–0.5 Hz band. See Fig. 3 and S2 for additional VLP examples.

176

P. Jousset et al. / Journal of Volcanology and Geothermal Research 261 (2013) 171–192

of short-period instrument is poor, we did not use short-period data for VLP quantitative analysis. In order to locate VLP seismic events, several methods have been proposed, including the semblance method (Kawakatsu et al., 2000; Almendros and Chouet, 2003) or full or constrained wave-form inversion (e.g., Chouet et al., 2005). For the VLP recorded on our network, particle motions are not very clear (possibly affected by the strong topography at Merapi; Neuberg and Pointer, 2000), reducing the effectiveness of the semblance method. Instead, in order to obtain an understanding of the seismic source processes, we followed the procedure outlined in Jolly et al., 2012. We used constrained waveform inversion to solve the discrete form of the ground displacement at seismic stations (Eq. (1)) and obtain source locations and mechanisms. We calculated the Green's functions (see Section 3.2) and convolved them with appropriate moment components to obtain synthetic seismograms (Section 3.3), which can be compared to observed data (Section 3.4). 3.2. Green's functions We use the notation and theoretical development of Aki and Richards (2002) and we exclude the single force term. The ground displacement at a seismometer site can be seen as a convolution (noted *) between the perturbation at the source and elastic Green's functions, as given by un ðt Þ ¼ Mpq ðt ÞGnp;q ðt Þ;

ð1Þ

with p,q = x,y,z and where the displacement response un(t) can result from the elastic response to a set of unidirectional (single forces) or bidirectional (dipoles or moments) impulses Gnp,q(t) of the Green tensor, and Mpq(t), is the time history of the pq elements of the moment tensor, n is the component of displacement at the receiver, and q implies spatial differentiation with respect to the q coordinate. To generate the Green's functions, we use a visco-elastic finite difference approach (Jousset et al., 2004).

(Tessmer et al., 1992). There is no accurate velocity model at Merapi volcano, despite tomography experiments performed to date at regional scale (Wagner et al., 2007; Bohm et al., 2013; Luehr et al., 2013). Ratdomopurbo and Poupinet (2000) used Vp = 3 km s−1 and Vp/Vs = 1.86; Wassermann and Ohrnberger (2001) used Vp = 2.8 km s −1. We applied a 1D velocity model (Table 1) for our initial modeling. The 1D layer model is distorted by the topography transformation method, and produces a suitable approximation of the wave field in the volcanic environment (e.g. Jousset et al., 2004; Jolly et al., 2010, 2012). Since we have little information about the attenuation structure of the volcano, we set the attenuation parameters to a low value. Hence, the visco-elastic system reduces to a fully elastic case for our present modeling purposes. To obtain stable waveforms within the frequency of interest (b 1 Hz), we used a grid spacing of 150 m laterally and 75 m in depth. This choice provided stable waveforms throughout the model area and allowed inversion of signals observed for the example events. Given the velocity model used in this study and the frequency of interest, we obtain the wavelength as >2 km (for 1D model velocities), giving more than 12 grid points per wavelength. Our method models topography with stretched grids and is less sensitive to instabilities caused by staircase topography, as used in other studies (e.g., Ohminato and Chouet, 1997; Ripperger et al., 2003) which require at least 25 points per minimum wavelength to ensure computational stability. In addition, our modeling propagates waves to short distances, so numerical dispersion is negligible. The Green's functions produced by the finite difference implementation would ideally be generated from a delta pulse consisting of all frequencies. However, the finite grid size dictates an upper frequency limit. To obtain stable Green's functions, we used a Ricker pseudo-impulse having duration of 12.5 s within a 25-s window with a band limited frequency range. This choice matches closely the lowest frequencies observed in the data and is within the stability requirements dictated by the finite difference method at the chosen grid spacing. A summary of the computational setup is shown in Table 2.

3.3. Modeling seismic wave propagation in complex media with topography

3.4. Moment tensor inversion

We use a 3-D finite-difference velocity–stress time domain formulation of the viscoelastic wave equations (Robertsson et al., 1994; Hestholm and Ruud, 1998), in a heterogeneous medium with non-flat free-surface topography (Hestholm and Ruud, 2000). This formulation (Hestholm, 1999) comprises (1) the full wave equation of motion (including P-SV waves, SH waves, and surface waves); (2) constitutive laws; and (3) equations for memory functions. The incorporation of surface topography is achieved by using a linear mapping transformation (Fig. 5) of the rectangular computational grid into a curved grid

Our inversion approach follows closely that of Nakano and Kumagai (2005a, 2005b) who found that stable moment tensor inversion results could be obtained with as few as 2–3 seismic stations if the point source models were limited to the simplest yet realistic geometries that are particularly relevant in the volcanic environment. We used the visco-elastic method outlined above (Section 3.3) to calculate Green's functions for moment tensor components within the Cartesian coordinate system described in Fig. 6. The 9 moment tensor components (Mxx, Myy, Mzz, Mxy = Myx, Mxz = Mzx, Myz = Mzy) can be found

Fig. 5. Curved grid in the (x,y,z) system and rectangular grid in the (ζ,ξ,η) system (Jousset et al., 2004).

P. Jousset et al. / Journal of Volcanology and Geothermal Research 261 (2013) 171–192 Table 1 Velocity model used for initial waveform inversion, Vp and Vs are primary and secondary wave velocities, rho is the density and Qp and Qs are the quality factors describing the attenuation of P and S waves respectively. Layer

z(m)

Vp

Vs

rho

Qp

Qs

1 2

0 20,000

3000 5200

1750 2920

2600 2650

950 900

920 900

for simple geometric systems (e.g., Chouet, 2003; Nakano and Kumagai, 2005b) which include: Crack geometry
 M xx ¼ Mo λ=μ þ 2sin2 θ cos2 Φ
 M xy ¼ Myx ¼ M o 2sin2 θ sin Φ cos Φ M xz ¼ Mzx
¼ M o ð2sin θ cos θ cosΦÞ M yy ¼ Mo λ=μ þ 2sin2 θ cos2 Φ

ð2Þ

T

Pipe geometry
 M xx ¼ Mo λ=μ þ cos2 θ cos2 Φ þ sin2 Φ
 2 M xy ¼ Myx ¼ −Mo sin θ sin Φ cos Φ M xz ¼ Mzx
¼ −M o ðsin θ cos θ cos ΦÞ 2

2

2

M yy ¼ Mo λ=μ þ cos θ sin θ þ cos Φ



ð3Þ

M yz ¼ Mzy ¼ −Mo ðsin θ cos θ sin ΦÞ
 M zz ¼ Mo λ=μ þ sin2 θ where λ and μ are elastic Lame parameters, and the angles θ and Φ are shown in Fig. 6. We assume λ = μ (Poisson ratio of 0.25, a value commonly used in crustal modeling studies, e.g. Nakano and Kumagai, 2005b) and tested a range of values consistent with the velocity model (Table 1). We solve for the source time function by inverting waveform data and the analytical solutions in Eq. (2) or (3) using the matrix inversion scheme given by: d ¼ Gm

ð4Þ

where d is waveform data, G contains a matrix including Green's functions, and m contains the moment tensor components. We can invert for the unknown scalar moment Mo (Lokmer et al., 2007), by expressing Eq. (4) as d ¼ GMo f ðλ=μ; θ; ΦÞ

ð5Þ

Table 2 Parameters used for all Green's function calculations. Nodes for damping Latitude (x direction, m) Longitude (y direction, m) Height (z direction, m) Space sampling (x, m) Space sampling (y, m) Space sampling (z, m) x position of source (m) y position of source (m) z position of source (m) Duration of source (s) Computation duration (s) Numerical time step (s) Sampling rate of seismogram (Hz)

where f is a column vector containing particular values of λ, μ, θ and Φ obtained from the analytical solutions in Eqs. (2)–(4). We use a standard inverse solver to obtain the moment estimates for each time step, and solve for each position, geometry and orientation separately as a grid search. The grid search component of our analysis is conducted over 200 m intervals. The search space ranges (Fig. 2) from UTM grid longitudes 437 km to 441.5 km in the x component (east-west) and latitudes 9163.5 km to 9168.5 km in the y component (south-north), and from 200 m depth to 3700 m depth below surface. The LBH gain factor is unknown and hence could not be used for the inversion portion of the VLP analysis. Instead, the inversion for the moment is computed using three well calibrated seismic stations (GMR, L56 and PAS) and station LBH waveforms are introduced only for the subsequent residual calculation (see Section 3.5). The inversion is computed in the time domain. To obtain an estimate of fit between observed and computed waveforms, we calculated residual estimates using the method outlined in Lokmer et al. (2007), where T

R ¼ ðd−GmÞ ðd−GmÞ=d d

M yz ¼ Mzy ¼ M o ð2sin θ cos θ sin ΦÞ
 M zz ¼ Mo λ=μ þ 2cos2 θ

15 10,000.0 (435 km to 445 km UTM) 10,000.0 (9160 km to 9170 km UTM) 8000.0 150.0 150.0 75.0 437 km to 441.2 km (UTM) 9163.8 km to 9168.8 km (UTM) 200 to 3700 25.0 25.0 0.01 100

177

ð6Þ

In our case, the residual is determined in a grid search for 10° increments in Φ and θ ranging 1° to 360° and 90° to 180°, respectively, for both pipe and crack models and also for a range of gain values for station LBH. To compare unambiguously the results for crack and pipe geometries and orientations, we would ideally need to account for different number of free parameters for the different geometries (e.g., Akaike, 1974; Nakano and Kumagai, 2005a, 2005b). However, as we perform a grid search, this is not required (O'Brien et al., 2010). With a network configuration like ours, we are aware that it is very difficult to distinguish between the source types and location (Lokmer and Bean, 2010) with the limited data and their geometric relationship to the source region. This produces an unconstrained tradeoff between the source parameters, as expressed in the Table 3. 3.5. Results of the inversion The location of the modeled crack is found at longitude = 438.2 km, latitude = 9164.8 km and depth = 200 m below surface (Fig. 1 and Table 3). The crack is oriented θ = 150°, Φ = 70°, and has a residual of 0.5094. The best pipe model is located at longitude = 439.7 km, latitude = 9165.3 km, and depth = 200 m below the surface. The pipe is oriented with θ = 130° and Φ = 340°, and has a residual of 0.5075 (Fig. 7). Fig. 8 shows the inverted time-independent moment components and Fig. 9 shows results of the waveform inversion for both model types. It is important to stress that due to the many sources of errors involved in data processing (resynchronization of stations and unknown amplitude of the signals recorded at LBH) our results may include substantial uncertainties. This is reflected in the wide range of possible locations with similar residual values, so the source position is approximate. The minimum solutions for crack and pipe models in our inversion are similar as expected and therefore we do not prefer one model over the other. However, the inversion method allows us to examine the range and location of plausible source models and to estimate the source volume that may be excited. We find that whatever the source mechanism, the source is found at shallow depth (b1 km) and in the southern unstable flank, as opposed to the northern stable flank of Merapi; and it is certainly located below the 2006 dome (Walter et al., 2013). We find moment values of the order of 1.5 to 4.5 × 10 12 Nm for the crack and the pipe, respectively. The principal axes of the moment tensor for a tensile crack have amplitudes λΔV, λΔV, and (λ + 2μ)ΔV, where ΔV represents the volume change associated with the crack opening or closure and λ and μ are the Lame coefficients of the elastic rock matrix (Chouet, 2003). The principal axes of the moment tensor for a pipe have amplitudes of (λ + μ)ΔV, (λ + μ)ΔV, and μΔV, where ΔV represents the volume

178

P. Jousset et al. / Journal of Volcanology and Geothermal Research 261 (2013) 171–192

Fig. 6. Coordinate system used for testing model geometries for Merapi volcano VLP events. Model geometries include volumetric cracks (left) and pipes (right).

change associated with the pipe inflation or deflation. Assuming λ = μ and ν = 0.25 with μ = 1 × 10 9 Pa to 1 × 10 10 Pa, we obtain ΔV = 150–1500 m 3 and ΔV = 200–2250 m 3 for the crack and for the pipe, respectively. Owing to the high uncertainties of our modeling, we conservatively conclude that these values are similar and of the order of 1000 m3. This value is comparable to volumes (3500 m 3) found at Kilauea, Hawaii (Ohminato et al., 1998), but one to two orders of magnitude larger than found at Merapi previously; Hidayat et al. (2002) estimated volumes of 10–60 m 3 for moment values of one order of magnitude less, during a slow dome-growth period. 3.6. VLP seismicity: Evidence of magma intrusion and overpressured gas release from a closed system In principle, the volume change represents the volume difference between the repose state of the source and the increase or decrease of the source volume, due to mass displacement within the source volume; it can be flow (in or out) of magma or gas, due to a pressure increase or decrease. We tentatively discriminate between magma and gas as a source of the VLP event. Considering a gas volume change in 200–2500 m 3 at about 200 m depth (in agreement with the inversion of the VLP signals), the corresponding release of SO2 at surface would be 0.1–1.2 tons, assuming that the gas is at lithostatic pressure, edifice density is 2600 kg/m 3, proportion of sulfur in the gas phase is ~ 1% and fumarolic gas temperature is 575 °C (Surono et al., 2012). Therefore, even for 200 VLP events, this represents a very low amount of erupted SO2 compared to the amounts of SO2 expelled during the full eruption (0.44 × 109 tons from Surono et al., 2012) and consistent with the low amount of gas expelled by satellites before the eruption (Fig. 14). We interpret this observation as due to an accumulation of gas in the upper edifice before the first explosive eruption. The discrimination between magma or gas release can be done on the basis of duration of the signal and its relation to the fluid viscosity. Hidayat et al. (2002) interpreted a VLP source at Merapi as being due to hot gas release and not magma, because their seismic VLP signal

Table 3 Summary of inversion minimum parameters. VLP

Crack Pipe

Location (UTM km)

Orientation

X (E-W)

Y (N-S)

Z (meters)

Θ

Φ

438.2 439.7

9164.2 9165.3

200 200

40 20

70 340

was too short (6 s) to allow movement of a viscous magma. The volume we find is larger (1000 m 3) and the duration is somewhat longer, i.e. 10–15 s for the downward first motion VLP phase and 20–30 s for the upward first motion VLP phase, see Figs. 3 and 4. The average magma effusion rate during the eruption was 25 to 35 m 3 s −1 (Surono et al., 2012; Pallister et al., 2013b). A displacement of 1000 m 3 of material at 35 m 3 s −1 requires ~ 30 s. This analysis suggests that the longer VLP signals would be related to magma movement, whereas the shorter duration VLP would be compatible with gas. The downward first movement signal (with shorter duration) would be associated with a loss of mass (gas release) and would decrease the volume below the summit. The source location of this VLP is found at very shallow subsurface and the first eruption was mainly phreato-magmatic (Surono et al., 2012). This is an indication that this shallow VLP was probably not due to magma flow in the crack (which would have produced a new dome at the start of the eruption, as observed for the recent eruptions in the past), but rather to highly-pressurized fluid escaping from below the dome. From the violence of the eruption on 26 October that destroyed the 2006 dome (Pallister et al., 2013b), pressure must have built progressively until the strength of the dome root (a volcanic plug) was overcome. However, Fig. 3 shows that the duration of the VLP with upward first motion is longer, on the order of a minute. Although we cannot not locate all VLP events, this result suggests that VLPs with upward first movement and long duration may carry signatures of magma rise into the edifice. Considering these arguments, we suggest that exsolved gas from the newly intruded magma was not sufficiently released though porosity and cracks and the 2006 dome and its solidified root acted as a sealed plug, building pressure in the underlying magmatic conduit and causing local deformation (Surono et al., 2012; Saepuloh et al., 2013). The numerous MP events observed before the 26 October explosion may have been expressions of the 2006 dome responding to pressure of the new intruded degassing magma. We conclude that the shallow part of the subsurface conduit was effectively “closed” before the first explosion on 26 October. This result suggests that magma started to rise before the first VLP occurred, i.e., before 9 September 2010, i.e., less than 2 months before the first eruption.

Residual

4. Long-Period eruption-related seismicity

0.5094 0.5075

We analyzed LP seismicity associated with the Merapi 2010 eruption to help constrain mechanisms for the 2010 eruption. We recognized more than 100 Long-Period (LP) events, including some

P. Jousset et al. / Journal of Volcanology and Geothermal Research 261 (2013) 171–192

179

Fig. 7. Crack (left) and pipe (right) model residuals with respect to position for 200 m, 700 m, and 1200 m depths. Diamonds are the three stations PAS, WOR and LBH (see Fig. 1).

180

P. Jousset et al. / Journal of Volcanology and Geothermal Research 261 (2013) 171–192 12

4

x 10

PIPE

3 CRACK

Moment (N.m)

2 1 0 −1 −2 −3 −4 −5 0

10

20

30

40

50

60

Time (s) Fig. 8. Time-independent moment time function for the two model types (crack and pipe). Dark lines are the best fit moment for the crack and dashed line for the pipe. The moment calculations shown assume λ = μ = 5.0 × 1010 N/m2.

monochromatic LP events. Fifty LP events were reported by the monitoring, including 27 LP events on 23–24 October; 1 LP event on 27 October and 22 LP events on 31 October (Budi-Santoso et al., 2013).

4.1. LP events characteristics and classification Among the numerous LP events that occurred before and during the eruption (Budi-Santoso et al., 2013), we focus here on 90 LP events that occurred during the most intensive activity (from 23 October to 5 November) (Fig. 10a). We utilize signals filtered in the 0.5–4 Hz band. This frequency band is subject to lower scattering effects, which are known to be prominent at Merapi (Wegler and Luehr, 2001).

a

The correlation matrix between waveforms has been used as a basis for event classification (e.g., Green and Neuberg, 2006; Saccorotti et al., 2007). We compute the correlation matrix of the vertical component of motion from station DEL for all pairs of detected filtered signals. Each element Sjk in the correlation matrix (Fig. S3) represents the maximum value of the cross-correlation between a pair (j,k) of signals, thereby yielding a quantitative measure of the degree of similarity between waveform pairs. Events were classified based on an equivalent class algorithm derived from Eardley's algorithm (Press et al., 1986). This algorithm identifies groups of elements satisfying a commutative “sameness” condition and sorts them into similar trees. In our case, the sameness condition is selected as Sjk > Smax, where Smax ranges from 0 to 1. When Smax is close to 1, the equivalence class algorithm leads to the identification of a small number of clusters (or none) consisting of seismograms with extremely similar signals (note that a cluster consists of at least two events). When Smax is reduced, new event pairs become eligible for consideration, and either augment existing clusters or form new clusters. Reducing Smax further may cause two or more clusters to coalesce into a single cluster. As Smax approaches the cross-correlation noise value, coalescence between event traces eventually dominates, resulting in a small number of very large clusters containing almost all the events. Because of the competing processes of new cluster formation and cluster coalescence, some value of Smax will yield a maximum number of clusters. Our dataset yields a maximum of 4 clusters corresponding to Smax ~ 0.8 (Fig. 10b). One of those clusters has only 2 events and is discarded here. Interestingly, LP events also happen to occur clustered in time; i.e., all events belonging to one cluster occur before the next LP cluster starts (Fig. 10c). The first cluster of LP events occurred before the fast growth of the new dome on 1–2 November (Surono et al., 2012; Pallister et al., 2013b). The second cluster comprises events with less clear correlation and occurs after the new dome growth had begun on 1 November. We also note that, after the dome appeared, the correlation coefficient decreased with time (Fig. 10c), while tremor amplitude increased (Surono et al., 2012), suggesting an evolving medium toward less organized structure as the eruption progressed until the time of the paroxysmal explosion.

b L56 EW

L56 EW

L56 NS

L56 NS

L56 Z

PAS EW

CRACK −4

L56 Z

PAS EW

PIPE

1 x 10 m/s PAS NS

−4

1 x 10 m/s

PAS NS PAS Z

PAS Z LBH EW

LBH EW LBH NS

LBH NS

LBH Z

GMR EW

GMR EW

GMR NS

GMR NS

GMR Z

GMR Z

Fig. 9. Example of model inversion results for (a) the best fit crack model and (b) the best fit pipe model. Black lines show data and gray lines indicate the best fit model. Parameters for crack and pipe models are indicated in Table 3.

P. Jousset et al. / Journal of Volcanology and Geothermal Research 261 (2013) 171–192

181

a 30−Oct−2010, 19:29:58

PUS

DEL

KLA

PLA

10 s

PUS

DEL

KLA

PLA 1

2

3

5

4

Frequency [Hz]

Number of clusters

b

4

3

2

1 0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

Smax Fig. 10. Families of Long-Period events. (a) Examples of LP events typical for each cluster. (Fig. S3 shows the correlation matrix of the correlation coefficient between events); (b) number of clusters found by the Earley classification algorithm as a function of the maximum correlation coefficient Smax taken to define families (see Table 4 for further details). (c) Correlation coefficients for a representative event from each cluster with time. This figure demonstrates that the resonators are excited one after another.

182

P. Jousset et al. / Journal of Volcanology and Geothermal Research 261 (2013) 171–192

c

1

Coefficient Correlation max

0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 30/10

31/10

01/11

02/11

03/11

Fig. 10 (continued).

4.2. Long-Period seismic event locations Table 4 and Fig. 11 give the list and location of LP events for each cluster. LP events within each cluster are expected to share a common location; however, different locations may be expected for different clusters. In order to locate the LP events, we did not use stations with unsynchronized signals; we used 4 to 5 stations (PUS, KLA, DEL, PLA and GRW) with accurate timing from GPS clock. LP events cannot be picked accurately, due to the emergent character of their signal amplitude (no clear P- and no S-wave). We did not stack LP events within each cluster (Green and Neuberg, 2006), but we located each single LP individually and compared their positions. To that aim, we combined an improved picking technique and a non-linear localization algorithm (Jousset et al., 2011b). We implemented a picking algorithm that computes and displays picking criteria (Akaike, 1978) directly on the amplitude trace of the seismogram (Maeda, 1985), without using autoregressive coefficients as in methods such as that of Leonard and Kennett (1999). For a digital seismogram x[1,N] of length N samples, the Akaike Information Criteria (AIC) value is defined as: AICðkÞ ¼ klog½varðx½1; kÞ þ ðN−k−1Þlog½varðx½k þ 1; NÞ

ð7Þ

in which k ranges through all the seismogram samples, var means the variance, and log is the natural logarithm. When selecting windows where a wavelet is located, this AIC picker defines the onset point as the global minimum (Fig. 11a). This picking method minimizes

operator errors. We also implemented an estimate of picking uncertainty, which is determined by the use of a pick interval, rather than the traditional single pick time. The pick interval is defined by two picked times, during which the wave's true arrival time (unknown) might be observed. This pick uncertainty is then taken into account in the non-linear inversion localization algorithm. The non-linear inversion localization algorithm performs a grid search within the volcanic edifice space of the minimum for differences between synthetic travel times computed at the nodes of the model parameterized grid and observed travel times in the least square sense. The result is a probability density function (pdf) of the difference in travel times being at its minimum, which yields the most probable hypocenter. When the pick interval is small (low pick uncertainty), then the pdf is rather narrow and location is more accurate, whereas when the pick interval is large (high pick uncertainty), the pdf is large and the location is less accurate. With this technique, the uncertainty of trigger onset of LP events due to emergent character of waves is incorporated into the final location error and yields reasonably consistent locations. We applied several iterations for the search grid by reducing the grid search volume and refining the scanned grid (see Fig. S4). The first grid is coarse (1 km spacing) and samples a large volume including the whole edifice (typically a cube of 10 km edge length). The second search volume is reduced is size (we use a volume defined by the 68% isovalue of the pdf from the previous iteration) and samples the medium with denser points (twice as many, e.g., 500 m spacing for the second iteration). In addition, we perform the grid search only within the edifice, i.e., we take into account the topography of the volcano, unlike

Table 4 Classification of LP events using the Eardley algorithm and setting the maximum correlation value to 0.8. Each cluster contains signal having a crosscorrelation coefficient larger than 0.8. They are clustered in space and time. This analysis demonstrates that the LP events reflect the resonance of several fluid containers excited successively as magma rises.

Cluster 1 (18 LP events)

Cluster 2 (6 LP events) Cluster 3 (3 LP events)

Times

Main frequencies (Hz)

Average Location (UTM km)

30.10: 13:39:40; 16:47:30; 19:07:39; 19:29:50; 20:19:20; 31.10: 00:19:45; 00:22:05; 00:23:55; 00:31:55; 01:16:09; 03:20:35; 03:58:19; 05:29:20; 07:01:10; 07:22:10; 08:13:35; 01.11: 02:21:59;02:48:52 02.11: 01:25:55; 02:54:45; 04:56:20; 10:18:50; 12:16:10; 15:46:14 03.11: 06:12:50; 06:23:15; 07:13:35

f = 1–2 Hz Q = 30–100

Lon. = 439.5 +/− 0.7 Lat. = 9167.5 +/− 0.5 Elev. = 2.35 +/− 0.3

f = 3–4 Hz Q = 170–250

Lon. = 439.1 +/− 0.3 Lat. = 9167.0 +/− 0.35 Elev. = 2.45 +/− 0.2

P. Jousset et al. / Journal of Volcanology and Geothermal Research 261 (2013) 171–192

183

a

b Cluster 1

Cluster 2

1 −> 30−Oct−2010 13:39:51 2 −> 30−Oct−2010 16:47:37 3 −> 30−Oct−2010 19:07:53 4 −> 30−Oct−2010 19:30:13 5 −> 30−Oct−2010 20:19:25 6 −> 31−Oct−2010 00:19:59 7 −> 31−Oct−2010 00:24:03 8 −> 31−Oct−2010 00:32:05 9 −> 31−Oct−2010 03:20:45 10 −> 31−Oct−2010 05:29:29 11 −> 31−Oct−2010 08:13:44 12 −> 01−Nov−2010 02:22:13 13 −> 01−Nov−2010 02:47:47

Cluster 3

1 −> 02−Nov−2010 01:26:02 2 −> 02−Nov−2010 02:54:50 3 −> 02−Nov−2010 04:56:26 4 −> 02−Nov−2010 10:18:59 5 −> 02−Nov−2010 12:16:14

1 −> 03−Nov−2010 06:13:04 2 −> 03−Nov−2010 06:23:23 3 −> 03−Nov−2010 07:13:42

9169 GRW

10

9167

Latitude (UTM km)

8 6 3 1 11 9 2 13 1 3 12 2 5 2 1 1 PUS 43 4

9168

KLA

9166

9165

9164

DEL

9163

9162

PLA

9161

434

435

436

437

438

439

440

441

442

443

Longitude (UTM km) Fig. 11. LP event location using double picking approach and non-linear localization algorithm that includes topography. a. Example of the picking technique using AIC criteria. The time window must be chosen to include the segment of seismogram of interest only. For a very clear P-wave onset in the seismogram, AIC values have a very clear global minimum at the sample where the P-wave is thought to begin. For a low signal to noise ratio, the global minimum does not give satisfactory results. b. Map of Merapi topography indicating locations of each cluster of LP events.

184

P. Jousset et al. / Journal of Volcanology and Geothermal Research 261 (2013) 171–192

Amplitude (ms−1)

Amplitude (ms−1)

many other hypocenter location linear algorithms (e.g., Lahr, 1999). Consequently, our procedure does not produce erroneous locations above the ground surface. Fig. 11b shows epicentral locations of 21 LP event hypocenters for the 3 clusters found in Section 4.1 (see also Table 4). LP events from the first cluster (14 events) are found at longitude 439.5 +/− 0.7 km, latitude 9167.5 +/− 0.5 km and elevation = 2350 +/− 300 m, which is close to the surface on the northern flank of the volcano. LP events from the second cluster (6 events, 5 locatable) are found at longitude 439.1+/− 0.25 km, latitude 9167.0 +/− 0.35 km, elevation 2450 +/− 200 m, which is higher in the edifice and closer to the summit. LP events from the third cluster (3 events) are more scattered and also in the north part of the edifice. LP events which were not classified within clusters locate in the same general area of the clusters. Globally, LP event locations are scattered within the north part of the edifice. They always locate at shallow depth (e.g., a few hundred meters deep), except for two deeper events that took place on 2 November at 15:47 (elevation 1750 m = ~2 km below the summit), and at 16:47 (elevation −2875 m = ~5 km below the summit). The successive locations of the LP main clusters suggest a temporal migration of hypocenters in the edifice toward the summit. The shallow locations and upward trend are similar to results for VT hypocenter locations and migration (Budi-Santoso et al., 2013). These results suggest that these migrations are associated with ascent of the magma that fueled the eruption. The migration of VT events has also been observed at basaltic volcanoes, but more rarely at silicitic volcanoes (e.g., Thelen et al., 2008; Budi-Santoso et al., 2013).

1

4.3. LP event complex frequency analysis It is now widely accepted that the frequency and time properties of LP events are directly related to the characteristics of a fluid-filled container, i.e., the container geometry, fluid/rock physical properties and the nature of the fluid (Chouet, 1996b). The characteristic waveform of an LP event is described by the decay of the oscillating acoustic response of a fluid-filled resonator (Kumagai and Chouet, 2000, 2001) and may be triggered by a time-localized excitation (e.g., Neuberg et al., 2006). Properties of LP events can be inferred through the Sompi analysis method, i.e. the analysis of the complex frequencies of the recorded signals, especially in their coda (Kumagai and Chouet, 2000). The complex frequency is defined as f - ig, where f is the frequency and g is the growth rate. We define the quality factor Q in terms of the complex frequency as Q = −f/(2g). The acoustic properties of the fluid and solid and the geometry of the resonator determine the complex frequencies. The acoustic properties of a crack containing magmatic or hydrothermal fluids can be described using the Sompi analysis method on synthetic LP events with different fluid composition and fluid container geometry. Kumagai and Chouet (2000) described comprehensively the acoustic properties of a crack with various types of fluids similar to volcanic fluids. They found that the complex frequency content, i.e., frequency and Q of LP events generated with such a model is directly dependent on the crack geometry and the type of fluid within the container. Their analysis is based on the transverse mode with wavelength 2/5 of the width of the crack (Chouet et al., 1994). Following

DEL−30−Oct−2010, 08:05:50.00

x 10−6

a

0

−1 0

1

10

20

2

4

30

40

50

60

6

8

10

12

x 10−6

b

0

−1 0

Gradient (sec−1)

Amplitude (ms−1Hz−1)

Time (s) 2

x 10−4

c

1

0 0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

d 0

1000 250

−0.02

100 50

−0.04 0.5

1

1.5

2

2.5

3

3.5

4

4.5

Frequency (Hz) Fig. 12. Complex frequency analysis performed on the coda of a Long-Period event recorded on 31 October 2010 at 00:20 at PUS station (1 km from the summit, see Fig. 2). a. Vertical component of the LP event. b. Signal used for the analysis; the box indicates the portion of the signal used for analysis. c. Corresponding Fourier spectral amplitude lower than 5 Hz. d. Complex frequencies of the individual wave elements (frequencies lower than 5 Hz) for all trial of the autoregressive orders (5–40). The clusters of points within rectangles represent clear signals, and the scattered points represent incoherent noise. The dotted lines represent lines along which the factor Q is constant (values indicated along the dotted lines).

P. Jousset et al. / Journal of Volcanology and Geothermal Research 261 (2013) 171–192

4.4. Spatio-temporal variations of the complex frequencies during the eruption

a

5 S5f

Unclustered Cluster 1 Cluster 2 Cluster 3

4

Frequency (Hz)

Kumagai and Chouet (2000), we assume that the LP events at Merapi were produced by the resonance of fluid-containers and we assume that the resonance mode they use is the active mode during the eruptive period for each resonator associated with each LP cluster. Checking for all resonance modes is still to be performed accurately; however, Kumagai and Chouet (2000) suggest that their main conclusions remain valid for other resonance modes. In our case, we consider our results of Sompi analysis to be only indicative of possible mechanisms. Surono et al. (2012) found Q values in the range 20–30 by analyzing an LP event on 31 October at 00:20. We find that this LP event belongs to the first LP cluster. In this study, we analyzed systematically and in detail seismograms and spectrograms to identify all LP events, using the method illustrated in Fig. 1. We also analyzed the complex frequency of all detected events (Figs. 12 and S5). Through analysis of many LP events, we confirm that Q values are in the 20–30 range, as suggested by Surono et al. (2012). In neglecting intrinsic attenuation effects, these Q values may be interpreted as resulting from resonance of cracks containing either mixtures of CO2–H2O, bubbly water, or basaltic magma with bubbles (Kumagai and Chouet, 2000; Kumagai, pers. comm., 2011).

185

S5g

3 S5a

2

S5d

S5e

S5b

1

S5c

New dome E

0 30.10

b

31.10

E

E

E

1.11

2.11

3.11

4.11

300 250 200 S5g

Q

150

It has been shown at many volcanoes that the complex frequencies of LP events may show spatial as well as temporal variations of both their main frequency and quality factor, with the latter ranging from 1 to several hundred (Lesage and Surono, 1995; Nakano et al., 1998; Kumagai et al., 2002; Nakano and Kumagai, 2005c; Lokmer et al., 2008). By analyzing several LP events with time, Kumagai and Chouet (2000) provided a basis for the interpretation of temporal variations in the observed complex frequencies of LP events in terms of the nature of the fluids beneath volcanoes. At Merapi, temporal analysis of the complex frequency for several LP events helps us in discriminating between the fluid mixture possibilities suggested by the Sompi analysis described in Section 4.3. We look at the temporal evolution of the principal frequency and corresponding Q value of more than 90 LP events spanning 29 October 2010 to 02 November 2010, during the course of the eruption, when most LP events were recognized (Surono et al., 2012; Budi-Santoso et al., 2013). For this analysis, we are interested in frequencies less than ~4 Hz recorded at the PUS station, which is the closest station to the sources. The maximal distance between the LP hypocenters determined previously (Section 3) and PUS is b 1 km. Waves traveling a shorter distance are less affected by scattering within the Merapi media (Wegler and Luehr, 2001). Our results show that the main frequency lies between ~1 and ~4 Hz and Q = 10 to 250. We observe variations with time (Fig. 13). From 29 October to 30 October the average frequency was 4.5 +/− 2 Hz. A decrease occurred on late 30 October resulting in a minimum with mean of ~1.5 +/− 1 Hz. A quiet period then ensued and lasted several hours. After 1 November, the frequency increased again to more than 4 Hz. On 2 November, a swarm of intense LP events started suddenly preceding the 3 November large explosion (this sequence is discussed later in the paper). All of these LP events have their main frequencies around 3–4 Hz, but with a second peak at ~ 1.5 Hz in the run-up to the 3 November explosion. Although the time-series is noisier, a similar evolution holds for the Q value defined by the complex frequency analysis (Fig. 13b). The global main trend may be described by a positive correlation between decreasing main frequencies and decreasing Q value. Assuming that the LP events observed at Merapi are generated by the resonance of fluid filled containers, we analyze qualitatively which model best explains our observations. The parameters of the model are α/a, b/μ and the crack stiffness C = bL/μd where α is the compression wave velocity of the rock matrix, a is the sound speed of the fluid, b is bulk modulus of the fluid, μ is the rigidity of the solid and L and d

S5f

100

Unclustered Cluster 1 Cluster 2 Cluster 3

S5a S5d S5e

50 E S5b

0 30.10

E

E

E

New dome

S5c

31.10

1.11

2.11

3.11

4.11

Fig. 13. Temporal evolution of a. the principal frequency and b. Q value derived from Sompi analysis of the vertical component record from station PUS between 29 October and 3 November 2010 at Merapi.

are the crack length and aperture, respectively. In this model, Q is mainly dependent on α/a, and the dimensionless frequency υ = fL/α is dependent on both α/a and ρf/ρs, where f is the frequency, and ρf and ρs are fluid and solid densities, respectively. A qualitative analysis of Plate 1 of Kumagai and Chouet (2000) and our values suggests that a simultaneous decrease of both Q (from ~250 to ~20) and of the dimensionless frequency ν (from ~0.12 to 0.1) may be associated with an increase of the ratio ρf/ρs (from ~0.75 to ~1) and a decrease in the ratio α/a (from ~30 to ~5). The dimensionless frequency is related to the dimensionless phase velocity of the crack wave, the wavelength of the wave and geometrical parameters of the resonator, such as its length. Within the limitations of the hypotheses of this modeling, there are several possible mechanisms to explain the frequency and resonance variations: (1) a change in the magma properties (increase in density through degassing or new intrusion of denser magma or spatial variations of the bubble content in the magma) and/or (2) a change in the geometrical properties of the resonator (Jousset et al., 2003), and/or (3) a change of resonator location (different resonators excited successively), as suggested by the waveform similarity analysis (Section 4.1). 5. Discussion on the eruption model and monitoring implications 5.1. Summary of the eruptive model Seismic features outlined in this paper (large number of earthquakes, elevated rate of seismicity, numerous VLP of two types, and long-period events clustered in few successive families) support the conclusion that the eruption was fed by a large volume of volatile-rich magma, as proposed by Surono et al. (2012) and

186

P. Jousset et al. / Journal of Volcanology and Geothermal Research 261 (2013) 171–192

106

SO2 flux determined from: Ground DOAS

SO2 flux (ton/day)

AIRS (12h mean)

10

5

10

4

IASI (12h mean) OMI (24h mean) Over-or under-estimation of the true flux Overview of the mean SO2 flux combining all observations

103

? ? ? Maximum flux during 1992-2007 eruptions Flux between eruptions, 1992-2007 (COSPEC)

102

RSAM Amplitude

PHASE IV

7 6

2.0 106

1.0 10

6

5

L L

E E E

4 3 2

E E

0

E

1

20 21 22 23 24 25 26 27 28 29 30 31 01 02 03 04 05 06 07 08 09 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24

October 2010

LP principal frequency (Hz)

PHASE III E

PHASE II

PHASE I

OMI detection limit

0

November 2010

Fig. 14. Comparison between SO2 flux and RSAM data (from Surono et al., 2012) with superimposed principal frequencies of the LP events from this study. a) Overview of SO2 degassing during the 2010 Merapi volcano eruption (UTC time). SO2 fluxes were determined from ground-based scanning DOAS measurements (hourly means) and satellite images, from IR IASI and AIRS sensors (12-h means) and the UV OMI sensor (24-h means). Question marks indicate gaps in IASI or OMI data coverage, or interference from SO2 plumes emitted by other Indonesian volcanoes. b, RSAM computed for the Plawangan station (6 km from the summit). PHASE I to IV refer to eruptive phases as defined in Surono et al. (2012). E stands for explosion; L for Lahar.

confirmed by several detailed analysis (e.g., Bignami et al., 2013; Borisova et al., 2013; Budi-Santoso et al., 2013; Charbonnier et al., 2013; Costa et al., 2013; Komorowski et al., 2013; Pallister et al., 2013b; Troll et al., 2013). A large volume of new magma and gas intruded and breached the edifice (generating many VT earthquakes (Budi-Santoso et al., 2013). The magma and gas pressurized the 2006 dome producing many MP events and summit deformation (Surono et al., 2012). The magma, volatiles, and hydrothermal fluids also pressurized the volcano's conduit which produced VLP events with upward first motion (as signature of the magma rise) and VLP events with downward first motion, as signature of gas release through the permeable upper conduit and 2006 dome. Gas release through cracks and conduit permeability did not balance gas accumulation due to the large magma flux, which powered the highly explosive phases of the eruption. After the first eruption on 26 October, the system was open and degassing occurred with oscillations as shown by satellite monitoring (Fig. 14). As shown in Fig. 14, as degassing increases the principal LP frequency increases. Magma heated the hydrothermal system, producing numerous LP events that occurred at different locations in the edifice, possibly due to magma migration towards the summit. Several mechanisms may explain the frequency and attenuation variations of the LP events including (1) geometry change of the resonator after violent explosions; (2) magma degassing from the intruded body; (3) successive intrusions of denser magma pockets; (4) different resonating areas within the volcano or, (5) a combination of those. 5.2. Seismic interpretation with respect to volcanology, petrology and gas emissions Accurate moment tensors can be recovered when the number of observed waveform allows good sampling of the radiation pattern (e.g., De Barros et al., 2011). In our study, location of sources is hindered by the low number of stations, by errors due to our resynchronization

processing, by the poor knowledge of the velocity model, and possibly also by scattering of seismic waves at shallow depth. However, our seismic observations confirm the high pressurization of the system before the first eruption and that the shallow parts of the subsurface feeder system were effectively closed before the first explosion on 26 October. Independently, EDM observations prior to the eruption revealed large deformations only of the southern flank of the edifice (Surono et al., 2012; Budi-Santoso et al. 2013). A displacement rate of about 10 cm/day was recorded on 17 October, concurrently with the largest VLP. However, Saepuloh et al. (2013) showed that deformation at Merapi was limited at the regional scale. We suggest that the deformation may have been triggered by pressurized magmatic fluids opening new fractures and moving into pre-existing planes of structural weakness (e.g., unconformities between eruptive episodes). Such a process would both inflate this part of the flank and effectively “lubricate” these discontinuities, thereby enabling gravitational creep. The bulk composition of erupted magmas did not show any strong change between 2006 and 2010 eruption (Borrisova et al., 2013; Costa et al., 2013) and therefore, composition alone cannot explain the changes between explosive and effusive phases. Instead, a much larger volume of rapidly ascending magma and gas has been proposed (Surono et al., 2012; Pallister et al., 2013b). Permeability has been shown to be important parameter to explain explosivity of volcanoes (Collinson and Neuberg, 2012; Degruyter et al., 2012). A large volume of gas was released during the Merapi eruption (Fig. 14) which confirms the existence of a large proportion of gas in the magma (large number of bubbles) at the beginning of the eruption. After pressure decrease associated with dome collapse, some delay is necessary for new bubbles to nucleate and existing bubbles to expand (Melnik and Sparks, 2002), offering an explanation for the lack of vesicular juvenile clasts in products from the paroxysmal explosive eruption on 5 November (Costa et al., 2013) and suggesting nearly complete separation of magma and gas phases by that time. However, Komorowski et al. (2013) found vesicular fragments in pyroclastic density currents emplaced between 29 October and 4 November. During this

P. Jousset et al. / Journal of Volcanology and Geothermal Research 261 (2013) 171–192

CRM (46 km)

Lempungearthquake M5.5, 700 km

Regional earthquake M4.2, <200 km

Large explosion

Swarm of LP events/explosions

DEL (3 km)

Fig. 15. Swarm of LP events triggered by a local M4.2 tectonic earthquake.

DEL 2.1 km

stop at 17:07:03.1

PUS <1km

stop at 21:25:08.8

PLA 5 km

CMR >46 km

(stop due to low battery) Series of explosions, Increasing intensity, Decreasing interexplosive duration

VEI4 Tectonic earthquake explosion

Series of domes collapses Pyroclastic flows and sub-plinian column

DEL 2.1 km

PUS <1km

PLA 5 km

CMR >46 km Fig. 16. Short-period sensor records during the main eruption and corresponding deposits features from Komorowski et al. (2013).

187

188

P. Jousset et al. / Journal of Volcanology and Geothermal Research 261 (2013) 171–192

period of time, large episodic releases of SO2 could also be observed at surface (Fig. 14). Consequently, the correlation between peak LP activity, large SO2 release at the surface and the only presence of vesicular juvenile components (in several days leading up to 4 November eruption) is circumstantial evidence for the proposed seismic model, in which major degassing was ongoing during this time interval, activating resonance of one or more fluid-filled conduits. Therefore, magma degassing provides a plausible mechanism to explain at least qualitatively the frequency variations observed during the eruption sequence. This model would result in a simultaneous increase of the ratio of fluid density to the solid density of about 25% (e.g., by degassing from the magma mixture) and/or an increase in the velocity of the seismic wave in the fluid, assuming a constant velocity of the rock matrix. This analysis neglects the possibility of a velocity increase of Merapi edifice media prior to or during eruptions (Ratdomopurbo and Poupinet, 2000). Given such an increase in media velocity, the increase in the wave velocity in the fluid would be even larger, which would imply even more degassing. If the frequency changes are due to density changes only, one may invoke a crystallization of the new intruded magma – a process which produces denser magma and releases gas (Tait et al., 1989). However, in Fig. 14, we also observe rapid variations of the LP principal frequency, changes that are too rapid to be explained by crystallization (Tait et al., 1989). Variations in magma density at these frequencies could rather be linked to variations in the magma volatile content associated with ascent-driven degassing (Johnson et al., 2008) or to the rise of a less-degassed magma batch (Jousset et al., 2000). An additional reason to favor variation in abundance of gas bubbles in the magma is that it rose very quickly as suggested by the LP events sequence, the VT migration and the large EDM rates. Strong gas release from a large magma volume coupled with insufficient permeability of the conduit and overlying dome may have contributed to the accumulation of gas and resulted in the explosive eruption of 26 October. Permeability changes also occur during the magma ascent (Costa et al., 2013), as isolated exsolved bubbles in deep part of the conduit start to coalesce and increase permeability, leading as a consequence, to a change of the geometry or the location of the resonator. 5.3. LP models and MP events To explain LP events, the resonance model for LP events (Chouet, 1996a; Neuberg, 2000) has been proposed at volcanoes having similar activity to Merapi (e.g. Montserrat; Neuberg et al., 2000). This model has been challenged by several models including fluid flow (Rust et al., 2008) or stick–slip models (Harrington and Brodsky, 2007; Pallister et al., 2013a) or inertial displacement of the whole column of magma (Johnson et al., 2008). Hybrid earthquakes have been shown to belong to a continuum with LF events at Merapi (Hidayat et al., 2000) as well as at Montserrat (Neuberg, 2000). They may therefore be special cases of LP events (Neuberg et al., 2000). At Merapi volcano, MP (Multiphase) events are also called hybrid events; they are clearly associated to dome growth (Ratdomopurbo and Poupinet, 2000). Their number decreased substantially after the 26 October eruption after the 2006 dome exploded (Budi-Santoso et al., 2013). At Montserrat series of similar hybrid events (Green and Neuberg, 2006) sometimes merged into tremor before dome collapses (Jousset et al., 2003). At Mount St Helens, similar exceptionally regular “drumbeat” series of hybrid events occurred before tremor and faster dome growth (Moran et al., 2008). Similar drumbeat activity occurred at Merapi during dome growth in October 1996, before the eruption eventually became explosive in January 1997. Although systematic analysis has not been performed on the frequency content of MP events at Merapi, the similarity of such drumbeat seismic activity at Merapi and Mount St Helens suggest that either a stick–slip model (Iverson et al., 2006; Moran et al., 2008; Pallister et al., 2013a) or a different type of resonating fracture (Waite et al., 2008) could explain the Merapi MP events. Note that

Table 5 Delays computed by cross correlation of ambient noise during 5 days before 18 October (values in seconds). Stations

GMR

LBH

L56

PAS

DEL

PUS

KLA

GMR LBH L56 PAS DEL PUS KLA PLA

160.0 60.7 12.6 5.97 6.15 2.2 6.46

−114.0 −170.0 −172.35 −172.5 −173.0 −175.8

−60.69 −61.4 −61.42 −75.85 −60.2

−2.3 1.52 7.5 6.91

0.02 0 −0.01

−0.01 −0.02

0.01

PLA

Pallister et al. (2013a) tend to favor stick–slip, but do not rule out the resonating fracture model of Waite et al. (2008) as an alternative to stick–slip for the Mount St. Helens drumbeats. At Merapi, we suggest that MP events would be stick–slip with low fluid involvement and LP events would be resonance of cracks with high fluid involvement. 5.4. Dynamic triggering of main eruptions in 2010 The magma volatile content varied during the eruption, which could explain transitions between explosive and effusive phases of the eruption as revealed from satellite observations (Pallister et al., 2013b) and deposits analysis (Komorowski et al., 2013). In addition, increasing gas pressures may cause a state of instability such that little external triggering force is needed to trigger an eruption (Walter et al., 2007; Eichelberger, 2010). We suggest that this was the case on 3 November 2010 when a ML ~ 4.2 tectonic earthquake occurred b200 km south of Yogyakarta at 4:05:21.17 (latitude 8.94°S, longitude 111.48°E, depth 61 km) (Fig. 15) and triggered a LP explosion sequence at Merapi. We suggest that the volcanic system had reached a critical state in which the small perturbation from passing seismic waves induced rapid degassing, magmatic fragmentation and eruption (Davis et al., 2007; Hill et al., 2002). Statistical studies of potential triggered eruptions (Marsan and Nalbant, 2005) show that the evidence for tectonic-earthquaketriggered eruptions is limited. We concur, and stress that the main trigger for the eruption of Merapi was pressurization due to magma ascent. However, increasing numbers of examples in which distant tectonic earthquakes induce seismic swarms and changes in eruptive activity (including at Merapi volcano; Walter et al., 2007) suggest that under appropriate conditions, dynamic triggering of eruptions does take place. 6. Summary and concluding remarks We analyzed seismic features of the 2010 eruption at Merapi volcano using advanced seismological techniques. More than 200 VLP events occurred before the first eruption and more than 90 LP signals were found during the eruption. Those numbers, larger than during any previously instrumentally monitored eruption at Merapi, were indicative of the much larger-than-normal eruption (Surono et al., 2012). We found two types of VLP events. The first type has upward first motion and longer signal duration, indicative of mass addition and upward movement of magma within the conduit. The second type has downward first motion and shorter signal duration, indicative of Table 6 Clock errors on 17 October 2010 obtained from cross-correlation of ambient noise. Errors are estimated from the standard deviation of delays for stations with synchronized GPS (see Fig. S1c). Station

GMR

LBH

L56

PAS

DEL

PUS

KLA

PLA

Clock (s)

3.9

−170.0

−61.1

1.13

0.93

0.44

2.07

0

P. Jousset et al. / Journal of Volcanology and Geothermal Research 261 (2013) 171–192

mass loss and gas release. Using records of cross-correlation of ambient noise we were able to resynchronize the signals from a clear VLP event with downward first motion that was recorded at all stations. We then conducted waveform inversion and demonstrated that our resynchronisation technique is accurate enough for localization and moment tensor wave-form inversion. We found that short duration deformations produce VLP signals that are localized within the southern flank of the volcano and at shallow depth, in the area of maximum deformation as determined by EDM measurements. Two mechanisms for the VLP events were tested, but we could not discriminate between them. However, a volume change of (~ 1000 m 3) at shallow depth is favored and likely records repetitive build up and release of pressurized gas in the region of the north flank beneath the 2006 dome. From this result, we suggest that the permeability of the 2006 dome and its underlying plug were inadequate to release gas pressure and prevent the initial eruption from being explosive. We analyzed 90 LP events and found 3 clusters in time and space suggesting that several fluid-filled resonators were successively excited as magma rose and the eruption proceeded. We also analyzed the complex frequency content, evolution of the principal frequencies and Q factors with time. Large variations were observed and interpreted as resulting from a resonating crack system and a change in the resonance pattern with time. We offer several hypotheses for this change in resonance: a change in the density ratio of fluid to solid, a change in the geometry of a single resonator, or excitation of different resonators. The latter hypothesis seems more plausible as it is corroborated by the evolution of LP cluster hypocenters, and VT earthquake migration (Budi-Santoso et al., 2013). We present evidence that Merapi volcano is subject to dynamic triggering. Consequently, in the future, whenever an eruption starts at Merapi volcano, careful analysis of the regional seismicity should be performed as well. As Surono et al. (2012) suggested, we stress again the value of distal stations within volcano monitoring networks; as such stations do not become saturated during eruptions such that they can provide critical information during crises. In addition, distal stations provide the data needed to evaluate the role of eruption triggering. During the most intense period of the eruption, the short-period stations were saturated as the eruption was very intense (Surono et al., 2012). However, low pass filtering between 0.01 and 0.1 Hz of signals from short-period instruments reveals VLP signals that were modulated by amplitude variations of several minutes to hours (Fig. 16). At first inspection, one would think that filtering the recorded clipped signal from overgained digitizers would not result in an interpretable signal. But surprisingly, the filtered Merapi records are extremely consistent and correspond well with the detailed volcanologic chronology of the eruption. The largest amplitudes in the filtered seismic signal correspond to the longest pyroclastic density currents, as identified in the chrono-stratigraphic record of deposits (Komorowski et al., 2013). The largest explosions contributed more seismic energy and kinematic energy to the pyroclastic density currents. At those frequencies well outside the classical use of short-period instruments, the instrumental response is nonlinear and it is unrealistic to derive any more sophisticated quantitative analysis from those records. However, this interesting finding may open new perspectives in terms of monitoring of large eruptions with inexpensive short-period sensors. This result suggests that even when signals are clipped, filtered seismic observations may be used to gain information during a powerful eruption on the size and nature of embedded eruptive phases.

189

seismic stations. S., P.J. and P.T. acknowledge the discussions with Amélie Vagner, Gonéry Le Cozannet, Franck Lavigne, Jean-Christophe Komorowski, Jean-Philippe Metaxian, and Philippe Lesage. Reviews from Ivan Lokmer, an anonymous reviewer and of Guest Editor John Pallister greatly improved the manuscript. P.J.'s work as Guest Editor for the Special Issue on Merapi was aided by wise advice from co-Guest Editors John Pallister and Surono, and Editor-in-Chief Jurgen Neuberg. Thank you for this invaluable experience, which provided also great pleasure. P.J., S., M.B. and P. T. received financial support from the European MIAVITA project. The MIAVITA project was financed by the European Commission under the 7th Framework Programme for Research and Technological Development, Area “Environment”, Activity 6.1 “Climate Change, Pollution and Risks”. AJ was funded by NZ Core funded research. The revision of the paper was carried out at Helmholtz Center, GFZ Potsdam. Appendix A. Resynchronization of seismic signals with lost GPS by ambient noise correlation (adapted from Vasseur, 2009 and Jousset et al., submitted for publication) For earthquake localization and tomographic experiments, the synchronization of a seismic network has to be very accurate (b 10 ms). The seismometer records should therefore be locked to an external clock via continuous GPS or radio signal. When the connection to the external clock is lost, the less accurate internal clock of the instrument is used and data may no longer be synchronized between instruments. In order to check the precision of the network's timing and optimally to correct clock errors, we determine clock differences between stations by pair-wise correlation of the observed random seismic wavefield. Using this procedure, clock errors were determined at chosen times by generalized inversion of the matrix of the cross-correlation and the clock difference (Sens-Schönfelder, 2008). The first-order estimation of the time difference between the clocks of two receivers is the shift that can be directly measured (Budi-Santoso et al., 2013). Clock differences between two stations correspond to the time lag with maximum correlation. More precise estimations of clock errors can be measured following a multistep procedure: 1. Computation of all Green's functions on a 1-day-long basis; 2. Calculation of the reference Green's functions by stacking all-day-long Green's functions over the reference duration (which can be the whole studied period); 3. Measurement of the time difference for each set of 1-day-long Green's functions, with respect to the reference Green's functions; 4. Calculation of the final time differences for each day's Green's functions by summing the previous measurement to the clock difference measured for the reference Green's function; 5. Construction of the vector Δ that contains clock differences obtained for all pairs of stations at each time; 6. Generalized inversion of the matrix G from Eq. (A1) which leads to the clock error's δ vector with respect to a reference station that was never unlocked from its external clock.

0

1 B1 B B1 G:Δ ¼ δ⇔B B0 B @0 0

−1 0 0 1 1 0

0 1 1 0 0 δLB1−LB2 0 1 B δLB1−LB3 C 1 0 C B C C ΔLB1 C B C B 0 −1 C C:B ΔLB2 C ¼ B δLB1−LB4 C @ ΔLB3 A B δLB2−LB3 C −1 0 C B C C @ δLB2−LB4 A 0 −1 A ΔLB4 1 −1 δLB3−LB4

ðA1Þ

Acknowledgments We acknowledge the efforts of the BPPTK staff for their careful monitoring of Gunung Merapi for many years. We thank J. Subandryio for his help in providing access to BPPTK facilities and I. Nurdin, Sadana and F. Purwoto for their assistance during field work to establish the

Eq. (A1) is expressed here for 4 stations arbitrarily named LB1, LB2, LB3 and LB4. The maximum of the reference Green's function provides the mean of the clock differences over the whole studied period and for each sensor pair. The inverse problem on Eq. (A1) implies that the absolute network time cannot be determined. Thus,

190

P. Jousset et al. / Journal of Volcanology and Geothermal Research 261 (2013) 171–192

the network time has to be adjusted by setting the setting the clock offset of a station that remained locked to its external clock. We applied this technique to a VLP recorded on 17 October 2010 (see Section 3). This analysis is restricted to the vertical component. The frequency band between ~ 0.017 Hz (minimal frequency of the receivers) and b 0.5 Hz proved to be useful at Merapi network. Fig. 3 shows 5 minutes of raw data records at selected stations of the network. A large event is recorded at all stations. Note that the VLP event is visible when records are filtered in the range [0.01–0.1 Hz]. We see that GMR, LBH and L56 have erroneous GPS clocks. We computed the clock delay by cross-correlation of ambient noise for several days before this event for 8 stations (GMR, LBH, L56, PAS, DEL, PUS, KLA and PLA). Some delays between stations are shown in Fig. S1c. The vector δ is built from the delays computed from the maximum of the cross correlation for each station couple (see Table 5), as calculated by generalized inversion of G for 8 stations with respect to PLA. The clock errors are given for this period in Table 6. The uncertainties of the method are quantified by computing standard deviations of time differences measurements throughout the period of interest for stations with synchronized GPS clocks. The average uncertainty for each set corresponds to the root mean square of sliding standard deviations. Results are shown for GRW (synchronized clock) with respect to the short-period stations (Fig. S1c). For the period considered, standard deviation are 44 ms, 29 ms, 57 ms and 55 ms, for station couples GRW-DEL, GRW-KLA, GRW-PAS and GRW-PUS, respectively. For the 17 October VLP event (Fig. S1b), we have standard deviations of 93 ms, 59 ms, 36 ms, 1013 ms, and 29 ms for station couples GMR-L56, L56-DEL, GMR-PLA, GMR-DEL, and L56-KLA, respectively. The results are independent of station type (short-period or broad band). The values obtained show that the accuracy of the method can exceed the accuracy of internal clocks after only few days of uncorrected instrumental drift and they follow a Gaussian probability distribution. Standard theoretical drift of instruments provided by the manufacturer is 8.10−7, i.e. 69.12 ms/day. The error of arrival-wave picking to locate events is about 10 ms to 100 ms; the analysis above indicates that at Merapi with cross-correlation stacked for several days, the clock error reduces to ~50 ms, implying an error on travel path of 150 m for a velocity of 3000 ms−1. This confirms that the method is very powerful. Appendix B. Supplementary data Supplementary data to this article can be found online at http:// dx.doi.org/10.1016/j.jvolgeores.2013.03.014. References Akaike, H., 1974. A new look at the statistical model identification. IEEE Transactions on Automatic Control 9, 716–723. Akaike, H., 1978. A Bayesian analysis of the minimum AIC procedure. Annals of the Institute of Statistical Mathematics 30, 9–14. Aki, K., Richards, P., 2002. Quantitative Seismology, Second edition. University Science Books, Sausalito, California, USA (704 pp.). Aki, K., Fehler, M., Das, S., 1977. Source mechanism of volcanic tremor: fluid-driven crack models and their application to Kilauea eruption. Journal of Volcanology and Geothermal Research 2, 259–287. Almendros, J., Chouet, B.A., 2003. Performance of the Radial Semblance Method for the location of Very Long Period volcanic signals. Bulletin of the Seismological Society of America 93 (5), 1890–1903. Bignami, C., Ruch, J., Chini, M., Neri, M., Buongiorno, M.F., Hidayati, S., Sayudi, D.S., Surono, 2013. Pyroclastic density current volume estimation after the 2010 Merapi volcano eruption using X-band SAR. Journal of Volcanology and Geothermal Research 261, 236–243. Bohm, M., Haberland, C., Asch, G., 2013. Imaging fluid related subduction processes beneath Central Java (Indonesia) using seismic attenuation tomography. Tectonophysics 590, 175–188. Borisova, A.Y., Martel, C., Gouy, S., Pratomo, I., Sumarti, S., Toutain, J.-P., Bindeman, I.N., de Parsival, P., Metaxian, J.-P., Surono, 2013. Highly explosive 2010 Merapi eruption: Evidence for shallow-level crustal assimilation and hybrid fluid. Journal of Volcanology and Geothermal Research 261, 193–208.

Budi-Santoso, A., Lesage, P., Dwiyono, S., Sumarti, S., Subandriyo, Surono, Jousset, P., Metaxian, J.-P., 2013. Analysis of the seismic activity associated with the 2010 eruption of Merapi Volcano, Java. Journal of Volcanology and Geothermal Research 261, 153–170. Camus, G., Gourgaud, A., Mossand-Berthommier, P.-C., Vincent, P.-M., 2000. Merapi (Central Java, Indonesia): an outline of the structural and magmatological evolution, with a special emphasis to the major pyroclastic events. Journal of Volcanology and Geothermal Research 100, 139–163. Cannata, A., Di Grazia, G., Montalto, P., Aliotta, M., Patane, D., Boschi, E., 2010. Response of Mount Etna to dynamic stresses from distant earthquakes. Journal of Geophysical Research 115, B12304. http://dx.doi.org/10.1029/2010JB007487. Charbonnier, S.J., Germa, A., Connor, C.B., Gertisser, R., Preece, K., Komorowski, J.-C., Lavigne, F., Dixon, T., Connor, L., 2013. Evaluation of the impact of the 2010 pyroclastic density currents at Merapi volcano from high-resolution satellite imagery, field investigations and numerical simulations. Journal of Volcanology and Geothermal Research 261, 295–315. Chouet, B.A., 1985. Excitation of a buried magmatic pipe: a seismic source model for volcanic tremor. Journal of Geophysical Research 90, B2. http://dx.doi.org/ 10.1029/JB090iB02p01881. Chouet, B.A., 1986. Dynamics of a fluid-driven crack in three dimensions by the finite difference method. Journal of Geophysical Research 91, 13967–13992. http://dx.doi.org/ 10.1029/JB091iB14p13967. Chouet, B.A., 1996a. Long-period volcano seismicity: its source and use in eruption forecasting. Nature 380, 309–316. http://dx.doi.org/10.1038/380309a0. Chouet, B.A., 1996b. New methods and future trends in seismological volcano monitoring. In: Springer-Verlag (Ed.), Monitoring and Mitigation of Volcano Hazards Scarpa/ Tilling, pp. 23–97 (Berlin). Chouet, B.A., 2003. Volcano-seismology. Pure and Applied Geophysics 160, 739–788. Chouet, B.A., Page, R.A., Stephens, C.D., Lahr, J.C., Power, J.A., 1994. Precursory swarms of long-period events at Redoubt volcano (1989–1990), Alaska: their origin and use as forecasting tool. Journal of Volcanology and Geothermal Research 62, 95–135. Chouet, B.A., Dawson, P., Arciniega-Ceballos, A., 2005. Source mechanism of Vulcanian degassing at Popocatépetl Volcano, Mexico, determined from waveform inversions of very long period signals. Journal of Geophysical Research 110, B07301. http://dx.doi.org/10.1029/2004JB003524. Collinson, A.S.D., Neuberg, J.W., 2012. Gas storage, transport and pressure changes in an evolving permeable volcanic edifice. Journal of Volcanology and Geothermal Research 243–344, 1–13. Costa, F., Andreastuti, S., Bouvet de Maisonneuve, C., Pallister, J.S., 2013. Petrological insights into the storage conditions, and magmatic processes that yielded the centennial 2010 Merapi explosive eruption. Journal of Volcanology and Geothermal Research 261, 209–235. Cronin, S.J., Lube, G., Dayudi, D.S., Sumarti, S., Subandriyo, S., Surono, 2013. Insights into the October–November 2010 Gunung Merapi eruption (Central Java, Indonesia) from the stratigraphy, volume and characteristics of its pyroclastic deposits. Journal of Volcanology and Geothermal Research 261, 244–259. Davis, M., Koenders, M.A., Petford, N., 2007. Vibro-agitation of chambered magma. Journal of Volcanology and Geothermal Research 167, 24–36. De Barros, L., Lokmer, I., Bean, C.J., O'Brien, G.S., Saccorotti, G., Métaxian, J.-P., Zuccarello, L., Patanè, D., 2011. Source Mechanism of Long Period events recorded by a high density seismic network during the 2008 eruption on Mt Etna. Journal of Geophysical Research 116 (B01304), 1–17. de Bélizal, E., Lavigne, F., Hadmoto, D.S., Degeai, J.-P., Dipayana, G.A., Mutaqin, B.W., Marfai, M.A., Coquet, M., Le Mauf, B., Robin, A.-K., Vidal, C., Cholik, N., Aisyah, N., 2013. Rain-triggered lahars following the 2010 eruption of Merapi volcano, Indonesia: A major risk. Journal of Volcanology and Geothermal Research 261, 330–347. Degruyter, W., Bachman, O., Burgisser, A., Manga, M., 2012. The effects of outgassing on the transition between effusive and explosive silicic eruptions. Earth and Planetary Science Letters 349–350, 161–170. Eichelberger, J., 2010. Messy magma mixture. Nature Geoscience 3, 593–594. http://dx.doi.org/10.1038/ngeo951. Ferrazzini, V., Aki, K., 1987. Slow waves trapped in fluid-filled infinite crack: implications for volcanic tremor. Journal of Geophysical Research 92 (B9), 9215–9223. Green, D., Neuberg, J., 2006. Waveform classification of volcanic low-frequency earthquake swarms and its implication at Soufrière Hills Volcano, Montserrat. Journal of Volcanology and Geothermal Research 153 (1–2), 51–63. Harrington, R.M., Brodsky, E.E., 2007. Volcanic hybrid earthquakes that are brittlefailure events. Geophysical Research Letters 34, L06308. http://dx.doi.org/ 10.1029/2006GL028714. Hartmann, M.A., 1934. Der grosse Ausbruch des Vulkanes G. Merapi Mittel Java im Jahre 1872. Natuurkundig Tijdschrift voor nederlandsch-Indie 94, 189–209. Hartmann, M.A., 1935. Die Ausbruche des G. Merapi (Mittel Java) bis zum Jahre 1883. Neues Jahrbuch für Mineralogie, Geologie, und Paleontologie 75 (B), 127–162. Hestholm, S., 1999. Three-dimensional finite difference viscoelastic wave modeling including surface topography. Geophysical Journal International 139, 852–878. Hestholm, S.O., Ruud, B.O., 1998. 3D finite difference elastic wave modeling including surface topography. Geophysics 63, 613–622. Hestholm, S., Ruud, B.O., 2000. 2D finite-difference visco-elastic wave modeling including surface topography. Geophysical Prospecting 48, 341–373. Hidayat, D., Voight, B., Langston, C., Ratdomopurbo, A., Ebeling, C., 2000. Broadband seismic experiment at Merapi Volcano, Java, Indonesia: very-long-period pulses embedded in multiphase earthquakes. Journal of Volcanology and Geothermal Research 100, 215–231.

P. Jousset et al. / Journal of Volcanology and Geothermal Research 261 (2013) 171–192 Hidayat, D., Chouet, B.A., Voight, B., Dawson, P., Ratdomopurbo, A., 2002. Source mechanism of very-long-period signals accompanying dome growth activity at Merapi volcano, Indonesia. Geophysical Research Letters 29 (23), 2118. http://dx.doi.org/ 10.1029/2002GL015013. Hidayati, S., Ishihara, K., Iguchi, M., Ratdomopurbo, A., 2008. Focal mechanism of volcano-tectonic earthquakes at Merapi Volcano, Indonesia. Indonesian Journal of Physics 19 (3), 75–82. Hill, D.P., Pollitz, F., Newhall, C., 2002. Earthquake–volcano interactions. Physics Today 41–47 (S-0031-9228-0211-020-9). Iverson, R.M., Dzurizin, D., Gartner, C.A., Gerlach, T.M., LaHusen, R.G., Lisowski, M., Major, J.J., Malone, S.D., Messerich, J.A., Moran, S.C., Pallister, J.S., Qamar, A.I., Schilling, S.P., Vallance, J.W., 2006. Dynamics of seismogenic volcanic extrusion at Mount St. Helens in 2004–2005. Nature 444, 439–443. http://dx.doi.org/10.1038/nature05322. Jenkins, S., Komorowski, J.-C., Baxter, P.J., Spence, R., Picquout, A., Lavigne, F., Surono, 2013. The Merapi 2010 eruption: An interdisciplinary impact assessment methodology for studying pyroclastic density current dynamics. Journal of Volcanology and Geothermal Research 261, 316–329. Johnson, J.B., Lees, J.M., Gerst, A., Sahagian, D., Varley, N., 2008. Long-period earthquakes and co-eruptive dome inflation seen with particle image velocimetry. Nature 456, 20. http://dx.doi.org/10.1038/nature07429. Jolly, A.D., Sherburn, S., Jousset, P., Kilgour, G., 2010. Eruption source processes derived from seismic and acoustic observations of the 25 September 2007 Ruapehu eruption—North Island, New Zealand. Journal of Volcanology and Geothermal Research 191, 33–45. Jolly, A.D., Neuberg, J., Jousset, P., Sherburn, S., 2012. New source process for evolving repetitious earthquakes at Ngauraohoe volcano, New Zealand. Journal of Volcanology and Geothermal Research 215–216, 26–36. http://dx.doi.org/10.1016/ j.jvolgeores.2011.11.010. Jousset, P., Rohmer, J., 2012. Evidence for remotely triggered micro-earthquakes during salt cavern collapse. Geophysical Journal International 191 (1), 207–223. http://dx.doi.org/10.1111/j.1365-246X.2012.0o55598.x. Jousset, P., Mori, H., Okada, H., 2000. Possible magma intrusion revealed by temporal gravity, ground deformation and ground temperature observations at Mount Komagatake (Hokkaido) during the 1996–1998 crisis. Geophysical Journal International 143, 557–574. Jousset, P., Neuberg, J., Surton, S., 2003. Modelling the time-dependent frequency content of low-frequency volcanic earthquakes. Journal of Volcanology and Geothermal Research 128, 201–223. http://dx.doi.org/10.1016/S0377-0273(03)00255-5. Jousset, P., Neuberg, J., Jolly, A., 2004. Modelling low frequency volcanic earthquakes in a viscoelastic medium with topography. Geophysical Journal International 159, 776–802. http://dx.doi.org/10.1111/j.1365-246X.2004.02411.x, 2004. Jousset, P., Budi-Santoso, A., Faria, B., Kouokam, E., Sincioco, J., 2011a. Catalogue at each of the volcano target (Merapi, Indonesia; Fogo, Cape Verde; Kanlaon, Philippines and Mount Cameroon). Report FP7-MIAVITA Deliverable 3.h. (63 pp.). Jousset, P., Haberland, C., Bauer, K., Arnason, K., 2011b. Hengill geothermal volcanic complex (Iceland) characterized by integrated geophysical observations. Geothermics 40, 1–24. Jousset, P., Delatre, M., Bouchot, V., Vasseur, J., Bitri, A., Sanjuan, B., 2013. Ambient seismic noise techniques: a new tool for geothermal exploration and monitoring — example at Bouillante geothermal field, French West Indies. Geothermics (submitted for publication). Julian, B., 1994. Volcanic tremor: non-linear excitation by fluid flow. Journal of Geophysical Research 99 (B6), 11859–11877. Kawakatsu, H., Kaneshima, S., Matsubayashi, H., Ohminato, T., Sudo, Y., Tsutsui, T., Uhira, K., Yamasato, Y., Ito, H., Legrand, D., 2000. Aso94: Aso seismic observation with broadband instruments. Journal of Volcanology and Geothermal Research 101, 129–154. Komorowski, J.-C., Jenkins, S., Baxter, P.J., Picquout, A., Lavigne, F., Charbonnier, S., Gertisser, R., Cholik, N., Budi-Santoso, A., Surono, 2013. Paroxysmal dome explosion during the Merapi 2010 eruption: Processes and facies relationships of associated high-energy pyroclastic density currents. Journal of Volcanology and Geothermal Research 261, 260–294. Kumagai, H., Chouet, B.A., 2000. Acoustic properties of a crack containing magmatic or hydrothermal fluids. Journal of Geophysical Research 105 (B11), 25493–25512. Kumagai, H., Chouet, B.A., 2001. The dependence of acoustic properties of a crack on the mode and geometry. Geophysical Research Letters 28, 3325–3328. Kumagai, H., Chouet, B.A., Nakano, M., 2002. Temporal evolution of a hydrothermal system in Kuratsu–Shirane volcano, Japan, inferred from the complex frequencies of longperiod events. Journal of Geophysical Research 107 (B10), 2236. http://dx.doi.org/ 10.1029/2001JB000653. Lahr, J.C., 1999. HYPOELLIPSE: a computer program for determining local earthquake hypocentral parameters, magnitude and first-motion patterns. U.S. Geological Survey Denver Federal Center, Denver, USA, Open-file report 99-23, paper. ((112 pp.) and online (http://greenwood.cr.usgs.gov/pub/open-file-reports/ofr-99-0023)). Legrand, D., Kaneshima, S., Kawakatsu, H., 2000. Moment tensor analysis of near-field broadband waveforms observed at Aso Volcano, Japan. Journal of Volcanology and Geothermal Research 101 (1–2), 155–169. http://dx.doi.org/10.1016/S03770273(00)00167-0. Leonard, M., Kennett, B., 1999. Multi-component autoregressive techniques for the analysis of seismograms. Physics of the Earth and Planetary Interiors 113, 247–264. Lesage, P., Surono, 1995. Seismic precursors of the February 10, 1990 eruption at Kelut volcano, Java. Journal of Volcanology and Geothermal Research 65 (135–146), 1995. Lokmer, I., Bean, C.J., 2010. Properties of the near-field term and its effect on polarization analysis and source locations of long-period and very-long period seismic events at volcanoes. Journal of Volcanology and Geothermal Research 192, 35–47.

191

Lokmer, I., Bean, C., Saccorotti, G., Patane, D., 2007. Moment tensor inversion of LP events recorded on Etna in 2004 using constraints obtained from wave simulation tests. Geophysical Research Letters 34, L22316. http://dx.doi.org/10.1029/ 2007GL031902. Lokmer, I., Saccorotti, G., Di Lieto, B., Bean, C., 2008. Temporal evolution of long-period seismicity at Etna volcano, Italy, and its relationships with the 2004–2005 eruption. Earth and Planetary Science Letters 266, 205–220. Luehr, B.-G., Koulakov, I., Rabbel, W., Zschau, J., Ratdomopurbo, A., Brotopuspito, K.S., Fauzi, S., Sahara, D.P., 2013. Fluid ascent and magma storage beneath Gunung Merapi revealed by multi-scale seismic imaging. Journal of Volcanology and Geothermal Research 261, 7–19. Maeda, N., 1985. A method for reading and checking phase times in auto-processing system of seismic wave data. Zizin 38, 365–379. Marsan, D., Nalbant, S.S., 2005. Methods of measuring seismicity rate changes: a review and a study of how the Mw 7.3 Landers affected the aftershock sequence of the Mw 6.1 Joshua Tree earthquake. Pure and Applied Geophysics 162, 5–25. http://dx.doi.org/ 10.1007/s00024-004-2665-4. McNutt, S.R., 2000. Volcanic Seismicity. In: Sigurdsson, H., Houghton, B., McNutt, S.R., Rymer, H., Stix, J. (Eds.), Chapter 63 of Encyclopedia of Volcanoes. Academic Press, San Diego, CA, pp. 1015–1033. Melnik, O., Sparks, R.S.J., 2002. Dynamics of magma ascent and lava extrusion at Soufriere Hills Volcano, Montserrat. In: Druitt, T.H., Kokelaar, B.P. (Eds.), The eruption of Soufriere Hills Volcano, Montserrat from 1995 to 1999: Geological Society, London, Memoirs, 21, pp. 153–171. Moran, S.C., Stihler, S.D., Power, J.A., 2002. A tectonic earthquake sequence preceding the April–May 1999 eruption of Shishaldin Volcano, Alaska. Bulletin of Volcanology 64, 520–524. Moran, S.C., Malone, S.D., Qamar, A.I., Thelen, W.A., Wright, A.K., Capan-Auerbach, J., 2008. Seismicity associated with renewed dome building at Mount St. Helens, 2004–2005. Chap. 2 of In: Sherrod, D.R., Scott, W.E., Stauffer, P.H. (Eds.), A Volcano rekindled; the renewed eruption of Mount St Helens, 2004–2006: U.S. Geological Survey Professional Paper 1750 (856 pp.). Nakano, M., Kumagai, H., 2005a. Waveform inversion of volcano-seismic signals assuming possible source geometries. Geophysical Research Letters 32, L12302. http://dx.doi.org/10.1029/2005GL022666. Nakano, M., Kumagai, H., 2005b. Correction to “Waveform inversion of volcano-seismic signals assuming possible source geometries”. Geophysical Research Letters 32, L12302. http://dx.doi.org/10.1029/2005GL02484. Nakano, M., Kumagai, H., 2005c. 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. Journal of Volcanology and Geothermal Research 147, 233–244. Nakano, M., Kumagai, H., Kumazawa, M., Yamaoka, K., Chouet, B.A., 1998. The excitation of and characteristic frequency of the long-period volcanic event: an approach based on an inhomogeneous autoregressive model of a linear dynamic system. Journal of Geophysical Research 103, 10031–10046. Neuberg, J., 2000. Characteristics and causes of shallow seismicity in andesite volcanoes. Philosophical Transactions of the Royal Society of London, Series A: Mathematical, Physical and Engineering Sciences 358, 1533–1546. Neuberg, J., O'Gorman, C., 2002. The seismic wavefield in gas-charged magma. In: Druitt, T.H., Kokelaar, B.P. (Eds.), The Eruption of Soufrière Hills Volcano, Montserrat, from 1995 to 1999. : Geological Society of London, Memoir, 21. Geological Society, London, pp. 603–609. Neuberg, J., Pointer, T., 2000. Effects of volcano topography on seismic broad-band waveform. Geophysical Journal International 143, 239–248. Neuberg, J., Luckett, R., Baptie, B., Olsen, K., 2000. Models of tremor and low-frequency earthquake swarms on Montserrat. Journal of Volcanology and Geothermal Research 101, 83–104. Neuberg, J., Tuffen, H., Collier, L., Green, D., Powell, T.W., Dingwell, D., 2006. The trigger mechanism of low-frequency earthquakes on Montserrat. Journal of Volcanology and Geothermal Research 153, 37–50. http://dx.doi.org/10.1016/ j.jvolgeores.2005.08.008. Newhall, C., Bronto, S., Alloway, B., Banks, N.G., Bahar, I., del Marmol, M.A., Hadisantono, R.D., Holcomb, R.T., MCGeehin, J., Miksic, J.N., Rubin, M., Sayudi, S.D., Sukhyar, R., Andreastuti, S., Tilling, R.I., Torley, R., Trimble, D., Wirakusumah, A.D., 2000. 10,000 years of explosive eruptions of Merapi Volcano, Central Java: archaeological and modern implications. Journal of Volcanology and Geothermal Research 100, 9–50. O'Brien, G., Lokmer, I., Bean, C., 2010. Statistical selection of the ‚best' source mechanisms from inversions of synthetic volcanic long-period events. Journal of Volcanology and Geothermal Research 115, B09303. Ohminato, T., Chouet, B.A., 1997. A free-surface boundary condition for including 3D topography in the finite-difference method. Bulletin of the Seismological Society of America 87 (2), 494–515. Ohminato, T., Chouet, B.A., Dawson, P., Kedar, S., 1998. Waveform inversion of very long period impulsive signals associated with magmatic injection beneath Kilauea Volcano, Hawaii. Journal of Geophysical Research 103 (B10), 23839–23862. Pallister, J.S., Cashman, K.V., Hagstrum, J.T., Beeler, N.M., Moran, S.C., Denliger, R.P., 2013a. Faulting within the Mount St. Helens conduit and implications for volcanic earthquakes. Geological Society of America Bulletin 125, 359–376. http://dx.doi.org/ 10.1130/B30716.1. Pallister, J.S., Schneider, D.J., Griswold, J.P., Keeler, R.H., Burton, W.C., Noyles, C., Newhall, G.C., Ratdomopurbo, A., 2013b. Merapi 2010 eruption—Chronology and extrusion rates monitored with satellite radar and used in eruption forecasting. Journal of Volcanology and Geothermal Research 261, 144–152.

192

P. Jousset et al. / Journal of Volcanology and Geothermal Research 261 (2013) 171–192

Picquout, A., Lavigne, F., Mei, E.T.W., Grancher, D., Noer, C., Vidal, C.M., Hadmoko, D.S., 2013. Air traffic disturbance due to the 2010 eruption of Merapi volcano. Journal of Volcanology and Geothermal Research 261, 366–375. Preece, K., Barclay, J., Gertisser, R., Herd, R.A., 2013. Textural and micro-petrological variations in the eruptive products of the 2006 dome-forming eruption of Merapi volcano, Indonesia: implications for sub-surface processes. Journal of Volcanology and Geothermal Research 261, 98–120. Press, W., Flannery, B., Teukolsky, S., Vetterking, W., 1986. Numerical Recipes. Cambridge University Press, Cambridge. Ratdomopurbo, A., Poupinet, G., 2000. An overview of the seismicity of Merapi volcano, (Java, Indonesia), 1983–1995. Journal of Volcanology and Geothermal Research 100, 193–214. Ratdomopurbo, A., Beauducel, F., Subandriyo, J., Nandaka, I.G.M.A., Newhall, C.G., Suharna, Sayudi, D.S., Suparwaka, H., Sunarta, 2013. Overview of the 2006 eruption of Mt. Merapi. Journal of Volcanology and Geothermal Research 261, 87–97. Ripperger, J., Igel, H., Wasserman, J., 2003. Seismic wave simulation in the presence of real volcano topography. Journal of Volcanology and Geothermal Research 128, 31–44. Robertsson, J., Blanch, J., Symes, W., 1994. Visco-elastic finite-difference modeling. Geophysics 59 (9), 1444–1456. Roman, D.C., Neuberg, J., Luckett, R.R., 2006. Assessing the likelihood of volcanic eruption through analysis of volcano-tectonic earthquake fault-plane solutions. Earth and Planetary Science Letters 248, 229–237. http://dx.doi.org/10.1016/j.epsl.2006.05.029. Rust, A.C., Balmforth, N.J., Mandre, S., 2008. The feasibility of generating low-frequency volcano seismicity by flow through a deformable channel. In: Lane, S.J., Gilbert, J.S. (Eds.), Fluid Motions in Volcanic Conduits: A Source of Seismic and Acoustic Signals: Geological Society of London, Special Publication, 307, pp. 45–56. Saccorotti, G., Lokmer, I., Bean, C., Grazia, G.D., Patane, D., 2007. Analysis of sustained long-period activity at Etna volcano, Italy. Journal of Volcanology and Geothermal Research 160, 340–354. Saepuloh, A., Urai, M., Aisyah, N., Sunarta, Widiwijayanti, C., Subandryio, Jousset, P., 2013. Interpretation of ground surface changes prior to the 2010 large eruption of Merapi volcano using ALOS/PALSAR, ASTER TIR and gas emission data. Journal of Volcanology and Geothermal Research 261, 130–143. Sens-Schönfelder, C., 2008. Synchronizing seismic networks with ambient noise. Geophysical Journal International 174, 966–970. http://dx.doi.org/10.1111/j.1365246X.2008.03842.x. Shimozuru, D., Miyazaki, T., Gioda, N., Matahelemual, J., 1969. Seismic observation at Merapi volcano. Volcanological Survey of Indonesian volcanoes: Bulletin of Earthquake Research Institute of Japan, 14, 1, p. 36. Surono, P., Jousset, J., Pallister, M., Boichu, M.F., Buongiorno, A., Budi-Santoso, F., Costa, S., Andreastuti, F., Prata, D., Schneider, L., Clarisse, H., Humaida, S., Sumarti, C., Bignami, J., Griswold, S., Carn, C., Oppenheimer, 2012. The 2010 explosive eruption of Java's Merapi volcano — a ‘100-year’ event. Journal of Volcanology and Geothermal Research 241–242, 121–135. Tait, S., Jaupart, C., Vergniole, S., 1989. Pressure, gas content and eruption periodicity of shallow crystallizing magma chamber. Earth and Planetary Science Letters 92, 107–123.

Takekawa, J., Mikada, H., Goto, T., 2013. Effect of heterogeneities on evaluating earthquake triggering of volcanic eruptions. Natural Hazards and Earth System Sciences 13, 231–237. http://dx.doi.org/10.5194/nhess-13-231-2013. Tessmer, E., Kosloff, D., Behle, A., 1992. Elastic wave propagation simulation in the presence of surface topography. Geophysical Journal International 108, 621–632. Thelen, W.A., Crosson, R.S., Kreager, K.C., 2008. Absolute and relative locations of earthquakes at Mount St Helens, Washington, using continuous data: implication for magmatic processes. Chap. 4 of In: Sherrod, D.R., Scott, W.E., Stauffer, P.H. (Eds.), A Volcano rekindled; the renewed eruption of Mount St Helens, 2004–2006: U.S. Geological Survey Professional Paper 1750 (856 pp.). Thierry, P., Jousset, P., Le Cozannet, G., 2008. MIAVITA: MItigate and Assess risk from Volcanic Impact on Terrain and human Activities. General Assembly IAVCEI, 17–22 August 2008, Reykjavik. Troll, V.R., Deegan, F.M., Jolis, E.M., Harris, C., Chadwick, J.P., Gertisser, R., Schwartzkopf, L.M., Borisova, A.Y., Bindeman, I.N., Sumarti, S., Preece, K., 2013. Magma differentiation processes at Merapi Volcano: inclusion petrology and oxygen isotopes. Journal of Volcanology and Geothermal Research 261, 38–49. Vasseur, J., 2009. Eléments de compréhension du champ géothermique de Bouillante (Guadeloupe). Master Degree report.IPG Paris, University of Paris 7 (38 pp.). Voight, B., Constantine, E.K., Sismowidjoyo, S., Torley, R., 2000. Historical eruptions of Merapi Volcano, Central Java, Indonesia, 1768–1998. Journal of Volcanology and Geothermal Research 100, 69–138. Wagner, D., Koulakov, I., Rabbel, W., Luehr, B.-G., Wittwer, A., Kopp, H., Bohm, M., Asch, G., MERAMEX Scientists, 2007. Joint inversion of active and passive seismic data in Central Java. Geophysical Journal International 170 (2), 923–932. http://dx.doi.org/ 10.1111/j.1365-246X.2007.03435.x. Waite, G.P., Chouet, B.A., Dawson, P.B., 2008. Eruption dynamics at Mount St. Helens imaged from broadband seismic waveforms: interaction of the shallow magmatic and hydrothermal systems. Journal of Geophysical Research 113, B02305. http://dx.doi.org/10.1029/2007JB005259. Walter, T.R., Wang, R., Zimmer, M., Grosser, H., Lühr, B., Ratdomopurbo, A., 2007. Volcanic activity influenced by tectonic earthquakes: static and dynamic stress triggering at Mt. Merapi. Geophysical Research Letters 34, L05304. http://dx.doi.org/ 10.1029/2006GL028710. Walter, T.R., Ratdomopurbo, A., Subandriyo, Aisyah, N., Brotopusptio, K.S., Salzer, J., Lühr, B., 2013. Dome growth and coulée spreading controlled by surface morphology, as determined by pixel offsets in photographs of the 2006 Merapi eruption. Journal of Volcanology and Geothermal Research 261, 121–129. Wassermann, J., Ohrnberger, M., 2001. Automatic hypocenter determination of volcano seismic transients based on wavefield coherence — an application to the 1998 eruption of Mt Merapi, Indonesia. Journal of Volcanology and Geothermal Research 110 (1–2), 57–77. Wegler, U., Luehr, B., 2001. Scattering behaviour at Merapi volcano (Java) revealed from an active seismic experiment. Geophysical Journal International 145 (3), 579–592. http://dx.doi.org/10.1046/j.1365-246x.2001.01390.x.