Structure of Masaya and Momotombo volcano, Nicaragua, investigated with a temporary seismic network

Journal of Volcanology and Geothermal Research 379 (2019) 1–11

Contents lists available at ScienceDirect

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

Structure of Masaya and Momotombo volcano, Nicaragua, investigated with a temporary seismic network Anne Obermann a,⁎, Irene Molinari b,c, Jean-Philippe Métaxian d,e, Francesco Grigoli a, Wilfried Strauch f, Stefan Wiemer a a

Swiss Seismological Service, ETH, Zurich, Switzerland Institute of Geophysics, ETH, Zurich, Switzerland Istituto Nazionale di Geofisica e Vulcanologia, Bologna, Italy d Institut de Physique du Globe de Paris, University Paris Diderot, Paris, France e ISTerre, IRD R219, CNRS, Université de Savoie Mont Blanc, Le Bourget-du-Lac, France f Instituto Nicaragüense de Estudios Territoriales, Managua, Nicaragua b c

a r t i c l e

i n f o

Article history: Received 5 February 2019 Received in revised form 21 April 2019 Accepted 23 April 2019 Available online 4 May 2019 Keywords: Volcano Tomography Ambient seismic noise Lava lake, Masaya, Momotombo, volcanic unrest

a b s t r a c t Since the end of 2013, the region around the two volcanoes Masaya and Momotombo, which includes the Nicaraguan capital Managua, has shown an unusually high seismic and volcanic activity. In December 2015, the Momotombo volcano erupted after 110 years of quiescence. Since mid-December 2015, the Masaya volcano has also shown gradually increasing activity, including the formation of a lava lake in its main crater. By adding 30 broadband stations, we had temporarily (December 2016–March 2017) densified the permanent Nicaraguan seismic network around these volcanoes to study the local seismicity and image the subsurface structure. During the observation period, we observed an overall low level of seismicity. Recorded events around Momotombo likely consist of aftershocks of the M5.5 earthquake that struck this area on September, 26th, 2016. At Masaya, we did not observe volcano-tectonic events. Using the continuous waveform recordings, we perform a 3D ambient seismic noise tomography that reveals a first image of the subsurface velocity structure below the Masaya and Momotombo volcanoes. While Momotombo shows a typical elongated low shear-wave velocity anomaly that reaches depths of about 8 km, Masaya does not show indications of a deep plumbing system. At Masaya, we have indications of a shallow (0–3 km) magmatic chamber, offset to the west and not directly below the active Santiago vent, where the crater lake is located At greater depth (3–8 km) a low velocity anomaly towards the northeast coincides in location with a modelled positive gravity anomaly and could indicate the presence of a former intrusive body. With this study we want to trigger further interest in the diverse tectonic and volcanic features of Nicaragua. Future, long-term seismic imaging and monitoring projects are of critical interest for the estimation of seismic and volcanic risks in Managua and the surroundings. © 2019 Elsevier B.V. All rights reserved.

1. Introduction Nicaragua lays at the Central American Convergence margin, where the Cocos plate is subducting below the Caribbean plate at a rate of about 8 cm/yr (DeMets, 2001). The convergence azimuth along the Nicaraguan margin is up to 25° oblique to the trench (LaFemina et al., 2009). This convergence obliquity leads to slip partitioning along the Central American margin and trench-parallel motion of the fore arc at rates as high as 14–15 cm/yr (DeMets, 2001; Turner et al., 2007). The Central American volcanic arc that formed as a result of the subduction, counts eighteen Holocene volcanoes in Nicaragua (Siebert et al., 2011). ⁎ Corresponding author. E-mail address: [email protected] (A. Obermann).

https://doi.org/10.1016/j.jvolgeores.2019.04.013 0377-0273/© 2019 Elsevier B.V. All rights reserved.

The Momotombo and Masaya volcanoes are part of this arc. Momotombo is positioned at the northern tip and Masaya at 40 km south of Lake Managua (Xolotlan). Lake Managua lies along the Nicaraguan Depression, a large NW-SE trending graben that crosses the Western part of the country. The origin of the depression as an intra-arc basin is explained by the oblique subduction process combined with the clockwise rotation of the Caribbean plate and counterclockwise rotation of the fore-arc sliver (Schliz, 2011). Here, the mafic igneous basement of oceanic lithosphere is overlaid by Cretaceous to Miocene marine sediments (Walther et al., 2000). Nicaragua's capital Managua (N1 Million inhabitants), is situated along this tectonic structure at the Southern shore of Lake Managua. The city was destroyed twice in the last century (1931, 1972) by disastrous earthquakes (e.g. Brown Jr et al., 1974; Sultan, 1931).

2

A. Obermann et al. / Journal of Volcanology and Geothermal Research 379 (2019) 1–11

Since the end of 2013, the Masaya and Momotombo volcanoes (Figs. 1, 2) have shown an unusually high seismic and volcanic activity. On April 10th, 2014, a M6.3 earthquake occurred close to the Momotombo volcano, followed by intense aftershock activity and a migration of seismicity towards Managua (Strauch, 2014). Over the following two years, the seismic activity has remained considerably higher than in the previous observation period from 1975 to 2013 (Fig. 1), with three additional large earthquakes ranging from M5.6–6.0 in the region. The last M6.0 earthquake (15/09/2016) near the Momotombo volcano and its aftershocks caused destructions in rural areas near the epicenter, in the village of Puerto Momotombo and the town La Paz Centro. Additionally, at Momotombo (Fig. 2C,D), an eruption occurred between December 2015 and January 2016, after 110 years of quiescence (Tenorio et al., 2015). At the Masaya volcano (Fig. 2A, B), a lava lake appeared on December 11th, 2015, along with a modest associated increase in seismicity in the surroundings of the volcano. As of April 2019 the lava lake is still present. Nicaragua has a permanent seismic network (Strauch et al., 2018) consisting of around 90 short period, broadband and strong motion stations with real time communication. The network is run by Instituto Nicaraguense de Estudios Territoriales (INETER), the Nicaraguan national geosciences institute (http://www.ineter.gob.ni). However, only a handful of stations operate in the direct vicinity of both volcanoes and most with only short-period sensors. With 30 broadband stations

from the instrumental pool of the Geo-Forschungs-Zentrum Potsdam (GFZ), we temporarily densified the seismic network for a period of three months (December 2016–March 2017). While this operation time is too short to gain a better understanding of the seismological/tectonic parameters of the recent activity, it is long enough to image the depth and spatial extent of the magma plumbing systems using ambient seismic noise tomography. After a brief review of the geological and tectonical background in the area of interest (Section 2), we discuss the temporary seismic network around the Momotombo and Masaya volcanoes (Section 3). We then present the results from a 3D ambient seismic-noise Rayleighwave tomography (Section 4) below the volcanoes; and discuss the regional seismicity during the observation period (Section 5). The imaged extent of Masaya's magma plumbing system is then compared with results from an earlier gravimetric survey (Métaxian, 1994) and discussed in the context of the tectonics of the region (Section 6). 2. Tectonic and volcanological background 2.1. Tectonics in the area Along much of the Middle America Trench, oblique subduction at a high rate of convergence results in northwest-directed trench-parallel block motion. However, in Nicaragua such faults are not well developed.

Fig. 1. Seismic events larger than ML 4.3 in the volcanic chain of Nicaragua A) from 1975 to 2013. B) from 2014 to 2016. The magnitudes of the largest events are denoted in black. A schematic representation of bookshelf faulting along Central American volcanic arc in Nicaragua is added following LaFemina et al. (2002). C) Zoom into the region of interest showing the broadband seismic network used in this study. Stations from the temporary MOMANIC network are marked in red, stations from the permanent network run by INETER in yellow. With black dots, we plot the seismicity observed from December 2016 to March 2017 using automated seismic event location and waveform stacking. The uncertainty of the location is about 0.5–1 km. Green and blue lines mark clusters discussed in the main text in more detail. D) Zoom into the Masaya caldera. Santiago and Masaya craters are indicated.

A. Obermann et al. / Journal of Volcanology and Geothermal Research 379 (2019) 1–11

3

Fig. 2. A) Lava lake in the Santiago crater at the Masaya volcano. B) Station deployment on a former lava stream at the Masaya volcano. C, D) View on the Momotombo volcano.

LaFemina et al. (2002) suggested instead that this motion is accommodated by bookshelf faulting that includes northeast-striking left-lateral faults (schematic representation in Fig. 1A). Previous works presented earthquake epicenter and focal mechanism data. These mapped fracture and fault data consistent with this model. Trench-ward migration of the volcanic arc since the Miocene and reactivation of northeast-striking Miocene structures may have led to the development of this arc- and trench-normal fault system. Changes in regional tectonic stress may initiate small-volume eruptions. Cailleau et al. (2007) investigated tectonic stress in Nicaragua and stated that earthquakes are encouraged on faults orientated normal to the volcanic arc and located between active volcanic centers that cause significant crustal weakening. Others, such as Alvarez et al. (2018), Funk et al. (2009), and Suárez et al. (2015) question this model and interpret the structural patterns near Managua as a time-transgressive rift opening. They suggest that the oldest extension (late Oligocene/early Miocene) began in the southeast and migrated to the northwest earlier phase of intra-arc normal rifting presently being superimposed by arc-parallel, right-lateral shear related to the northwestward transport of the Central America fore-arc sliver. 2.2. Momotombo volcano The Momotombo volcano is a quaternary stratovolcano with a height of 1297 m at the north shore of Lake Managua. In 1905 Momotombo erupted basaltic material and has remained in quiescence until December 2015, while maintaining a persistent state of hightemperature fumarolic activity (Menyailov et al., 1986). The high fumarole temperatures suggest high proportions of magmatic compounds and limited interaction with meteoric water or the hydrothermal system (Frische et al., 2006). In the Momotombo region, intensive faulting of the down-going slab facilitates deep penetration of seawater into the subduction zone Ranero et al. (2000). As a consequence, the magmatic plumbing system is associated with thick deep- and shallow-water sediment deposits (Ranero et al., 2000; Snyder and Fehn, 2002; Walther et al., 2000). At the southern flank of the volcano is a geothermal field, which was discovered in 1970. The field drilling program was initiated in 1974. In

2016, forty-seven wells have been drilled (depths of 310–2839 m) with a combined wellhead energy capacity of 35–77 MW (Kaspereit et al., 2016). Since 2000, seven wells in the eastern part of the wellfield have been used for injection (previously Lake Managua was used to dispose of the brine). Kaspereit et al. (2016) constructed a conceptual model that proposes that Momotombo possesses many of the characteristics of a typical volcanic-hosted geothermal system: hot geothermal fluids ascend around the volcanic intrusion and due to the hydraulic gradient, then migrate downward and laterally through permeable flanks of the volcano. 2.3. Masaya volcano The Masaya volcano is an active volcanic complex, about 25 km southeast of Managua, on the western edge of the Nicaraguan depression. The inner caldera (560 m above sea level) measures 6 × 11.5 km2 and was formed by Plinian eruptions of about 14 km3 of material, b6000 years go (Williams, 1983). It is occupied by several pit craters, San Pedro, Nindiri and the dominant craters Santiago and Masaya (Fig. 1D). The inner caldera is enclosed by a 25 km wide older structure, the Las Sierras caldera. The Masaya volcano has been the site of frequent periods of unrest over the past 500 years; voluminous degassing, lava lakes, minor explosions from the Santiago cone and summit pit craters (Rymer et al., 1998; Stoiber et al., 1986; Viramonte and Incer-Barquero, 2008). Since the 1990s, strong degassing of up to 600–1800 metric tons day−1 in 1999 has occurred with high SO2, CO2 concentrations causing intense environmental degradation and health hazard (Burton et al., 2000; Delmelle et al., 2002). Although situated in a subduction regime, the Masaya caldera shows basaltic and tholeitic activity (Bice, 1985; Carr, 1984; Williams, 1983) likely because of the low-pressure intra-crustal fractionation due to a thin crust (Carr, 1984). This thin crust allows the formation of superficial magmatic intrusions and lava lakes. The visible part of the intermittent lava lakes is situated at the bottom of the Santiago crater. The actual extent of the lava lake might be much larger. The most recent appearances of lava are documented from 1989 and 1993 (Smithsonian Institution, 1993). In the 1990s,

A. Obermann et al. / Journal of Volcanology and Geothermal Research 379 (2019) 1–11

Group velocity (km/s)

4

there was probably a shallow lava reservoir at the bottom of Santiago crater. This theory is supported by vents appearing for a few weeks each time, which permitted incandescences to be seen. However, a lava lake was not directly visible. Viramonte and Incer-Barquero, 2008 collected historic records that mention the appearance and disappearance of lava lakes at Masaya in the past centuries. The lava lake that appeared on December 11th, 2015 (also observed by satellite-based MODIS thermal observations (Aiuppa et al., 2018)) was still present in April 2019 and is clearly visible from the crater rim.

4

A

3 2 1 0 0

# measurements

In a benchmark study to obtain more information about the sources of unrest at the MOmotombo and MAsaya Volcanoes in NICaragua (MOMANIC Project), we set up a temporary seismic network of 30 broadband stations (Trillium-Compact equipped with CubeDataLoggers). These additional stations around the Momotombo and Masaya volcanoes temporarily densified the permanent network (Fig. 1C). We acquired data for three months from the end of December 2016 to the end of March 2017 at a sampling rate of 200 Hz. The stations were either placed in the free field or, if possible, in a quiet corner of a local farm. All stations were buried at about 30 cm depth, placed on a cobblestone and covered with a bucket. However, particularly at the Masaya volcano, the ground conditions varied largely. The superficial volcanic structure is composed of a very heterogeneous, poorly consolidated sequence of recent tephra and lava flow layers, rendering proper station placement often challenging or impossible. Stray animals and failing batteries unfortunately reduced the overall data completeness to 80%.

10

15

20

15

20

Period (s) 300

3. The seismic network of the temporary MOMANIC Project

5

B

200 used for analysis

100 0 0

5

10

Period (s) Fig. 3. a) Group-velocity dispersion curves of vertical-vertical component CCs obtained using a FTAN analysis. b) Number of group-velocity measurements as a function of period. We restrict our analysis to periods with at least 150 measurements around each volcano, which limits us to a frequency band from 1 to 15 s.

velocities in the period range from 1 to 5 s. Fig. 3b shows the number of measurements as a function of period. We limit our analysis to period ranges with at least 50 measurements around each volcano (N150 measurements in total), which restricts us to periods between 1 and 15 s. 4.3. Inversion for 2D group velocity maps

4. Ambient seismic noise tomography For the ambient noise tomography, we first compute the crosscorrelations CC between the continuous vertical component recordings of all station pairs (Section 4.1). Then, we calculate the frequency dependent group travel times (Section 4.2) and invert them to construct 2D group velocity maps at different frequencies (Section 4.3). To obtain the depth structure, we invert the regionalized dispersion curves for local 1D shear velocity models in every cell of the grid using a neighborhood algorithm (Section 4.4). 4.1. Pre-processing and computation of noise cross-correlations Prior to the calculation of the CCs, we apply the following processing steps: instrumental correction; resampling of the data to a sampling frequency of 10 Hz; band-pass filtering between 0.2 and 30 s; elimination of 2 h signal segments that show amplitudes N3 times the standard deviation of the daily trace (e.g. local earthquakes), spectral whitening of the amplitude from 0.2 to 30 s. We then calculate the CCs between all station pairs for the remaining two-hour segments and stack them over the three months. We average positive and negative lag-times to enhance the part of the signal that is symmetric. This procedure slightly increases the signal to noise ratio and helps to homogenize the noise sources distribution. 4.2. Rayleigh-wave group velocity measurements We estimate the group velocity dispersion curves with a frequencytime analysis from 0.5 to 20 s (Levshin et al., 1989; Ritzwoller and Levshin, 1998). We use a graphical user interface that involves analyst validation of the dispersion curves (Fallahi et al., 2017; Mordret et al., 2014). Group velocities from inter-station distances smaller than 1.5 wavelength and SNRb5 are not considered. Fig. 3a shows the final set of group velocity dispersion curves calculated for this study. The dispersion curves are relatively flat, with velocities between 0.7 and 3 km/s. In particular, dispersion curves from the Masaya volcano show low

We perform tomographic inversions of the group-velocity measurements for periods between 1 and 15 s using the algorithms described by Barmin et al. (2001) and implemented by Mordret et al. (2013) in a Cartesian version. The algorithms are based on ray theory involving a regularization function composed of a spatial Gaussian smoothing function and a constraint on the amplitude of the perturbation depending on local path density. For the 2-D models we use 56 × 37 = 2072 cells with a grid size of 0.03°×0.03° (4 km EW, 4 km NS). Please note that the small grid-size was chosen according to the main area of interest (region around the Masaya and Momotombo volcanoes) which has a dense seismic network. We are aware that the tomography is not well constrained outside of this area, but still think that the data from the far away stations can add information about the structure at depth. The initial model for the inversion has a constant velocity that is taken as the mean group velocity for each period. We perform the inversion in two steps (Moschetti et al., 2007; Obermann et al., 2016). First, we invert a smooth map that is used to identify and reject measurements with travel time residuals greater than two standard deviations. Then, we use the remaining measurements to produce the final group velocity maps. The topography is not taken into account during the inversion procedure. Brenguier et al. (2007) estimated the error of this approximation on the velocity as b5% for Piton de la Fournaise with an elevation of 2600 m. Momotombo has an elevation of 1297 m and Masaya 635 m. The topography effect should hence be negligible compared to the group velocity variations of about 5–20% at each period (Fig. 3). In Fig. 4 the ray paths with respective group velocities (A) and the inverted 2D group-velocity maps (B) are shown at periods of 1, 4, and 10 s, corresponding to an increasing depth penetration. The mean velocity obtained from the inversion for the group-velocity maps, slightly increases with the period from 1.3 to 2.5 km/s. At short periods (1–4 s), we observe a low velocity zone at the location of the volcanic cones Masaya and Momotombo. With increasing period this low velocity zone appears elongated between the two volcanoes. Poorly resolved areas with b4 ray paths per grid cell are not shown. A variance reduction of the travel-time residuals (norms in L2) exhibits values ranging from

A. Obermann et al. / Journal of Volcanology and Geothermal Research 379 (2019) 1–11

5

Fig. 4. A) Rayleigh-wave group velocities associated with each path. The black square marks the inset for which the 2D-Rayleigh-wave group velocity maps are shown (B). The results are presented at periods of 1, 4 and 10 s. Areas that are not well-resolved are masked in the group velocity maps.

68 to 83%, indicating that the retrieved velocity maps fit the data well. In Fig. 5 we show the corresponding spatial resolution. A resolution shift b4 km indicates that the centre of the retrieved anomaly is in the same cell as the input spike, which is the case for the central areas around Momotombo and Masaya, at all frequencies. Fig. 4A shows that with increasing period the path density is dominated by SE-WN striking rays and cross-firing is practically absent for periods below 5 s. For this reason, we also perform a detailed tomography of Masaya volcano alone, using periods from 1 to 5 s. Here, we use a smaller grid size of 0.013° × 0.013°. The results are shown in Fig. 6. We observe a low velocity zone in the south-west part of the Santiago crater that is most pronounced at 2–3 s. At larger periods the velocity anomaly disappears, indicating that there is a low velocity body at shallow depth of about 1.5–3 km. This low velocity zone extends with less strength towards the NW. The star shape of the anomaly is caused by the distribution of stations. 4.4. Depth inversion To obtain the depth structure beneath the volcano, we average the group velocities in each model cell, as obtained for the 2D inversions

at periods from 1 to 15 s (with steps of 0.1 s between 1 and 10 s and steps of 1 s afterwards). From the averages, we construct local dispersion curves for each cell that we subsequently invert to assess a 1-D local layered velocity model. For this inversion we use a Monte-Carlo global direct-search technique (Molinari et al., 2015; Sambridge, 1999; Wathelet, 2008). The solution space is sequentially and non-uniformly sampled, taking into account the data fit achieved by old samples. To reduce the number of parameters to be inverted, we discretize the continuous functions into constant thickness layers with constant velocity. In particular, the vertical parametrization consists, from top to bottom, of: 4 layers with 0.5 km thickness, 4 layers with 1 km thickness, 2 layers with 1.5 km, 3 of 2 km, and 5 layers of 4 km thickness. For each layer, the allowed parameter ranges for the shear-wave velocity (Vs) are defined based on the plausible crustal models in volcanic region, allowing a large variability to resolve both low and high velocity body. In total 34,000 Vs models are sampled for each location and the best-fitting 1000 models are kept to calculate a mean model for each cell location. The similarity between best and mean models is very high, meaning that the solution is generally robust (similar to the mean). In Fig. 7, we show horizontal slices that we extracted from the best-fitting 3D Vs model at various depths ranging from 0.75 to 14 km. The shear-wave

6

A. Obermann et al. / Journal of Volcanology and Geothermal Research 379 (2019) 1–11

Fig. 5. Spatial resolution at 1, 4 and 10s. A) Resolution shift. A resolution shift b4 km indicates that the centre of the retrieved anomaly is in the same cell as the input spike. B) Spatial resolution maps.

velocity variations are shown with respect to the mean velocity in each layer. Areas that are not well constrained are masked. Around the larger Masaya area, we observe a low velocity anomaly of up to −15 to −20%. At shallow depth (0.75 km, Fig. 7), the low velocity anomaly is particularly strong towards the SW and likely related to surface geology consisting of former lava flows with large voids. At greater depth, the low velocity zone below Masaya concentrates NE/NW of the crater area (3.5 km, Fig. 7). The low velocity anomaly at Momotombo is about −5%. Both anomalies are no longer visible in the 6.5 km depth slice (Fig. 7). At greater depth, the low velocity anomaly concentrates along a SE-WN line. This could be due to the subduction trench, but can equally well be a bias, as we are missing cross-firing at periods below 5 s (Fig. 4A) and most of the ray paths fall on this SE-WN line. Vertical cross-sections along the Nicaraguan volcanic chain (Fig. 8A) and across Masaya (Fig. 8B) and Momotombo (Fig. 8C) volcano allow a better overview of the depth extent of the anomalies, as here the shearwave velocity variations are plotted with respect to the mean velocity in each layer. The location of the cross-sections is shown in Fig. 7. We notice the very shallow (1.5 km depth) extended low velocity anomaly below Masaya (Fig. 8A) which is offset with increasing depth to the NW/NE (Fig. 8A,B) and disappears at a depth of 6 km. At Momotombo (Fig. 8C), we observe an extension of the low velocity zone down to about 8 km depth. We notice a small low velocity imprint of Momotombito Island. 5. Regional seismicity Using all available stations, from December 2016 to March 2017, we detected about 300 earthquakes within the region of interest using the

Seiscomp3 module Scanloc developed for local-regional seismicity monitoring operations (Grigoli et al., 2018). Earthquakes locations are then refined using a waveform stacking based method and are shown in Fig. 1C. Since a local magnitude scale is not available for the region, magnitude estimation was performed by using the original relation as defined by Richter (Richter, 1935). The estimated magnitude of the events range from 1.5 to 2.5. During our investigation period, we did not detect a seismic cluster at the Masaya volcano; most events were located around Momotombo. We therefore focus on the Momotombo area. 5.1. Detection and location with waveform stacking methods Seismicity datasets in volcanic environments are often characterized by a massive number of small seismic events which are strongly noise contaminated. For these reasons, standard techniques for seismic data analysis based on automatic picking are affected by various problems reducing their performances (Grigoli et al., 2018). To overcome these problems alternative waveform-based methods for seismicity characterization have been proposed recently. Waveform based detection and location methods exploit the coherence of the waveforms recorded at different stations and do not require any automated picking procedure. The main advantage of these methods relies on their robustness even when the recorded waveforms are very noisy. On the other hand, like any other location method, the location performance strongly depends on the accuracy of the available velocity model. This is particularly true in volcanic environments, where strong topography effects are combined with a poorly known and highly heterogeneous 3D velocity structure. In this work, we used a location method which combines

A. Obermann et al. / Journal of Volcanology and Geothermal Research 379 (2019) 1–11

7

Fig. 6. Focus on the Masaya volcano. 2D velocity maps at different periods of 1, 2, 3 and 4 s. The star shape is an artefact from the station geometry.

Fig. 7. Horizontal cross-sections of the obtained 3D shear-wave velocity models at various depths. Lateral velocity perturbations are shown relative to the mean velocity. The significance level of velocity variations is approximately ±3%. Additionally, the slices show the relocated seismic events around Momotombo using a master-event waveform stacking locator.

8

A. Obermann et al. / Journal of Volcanology and Geothermal Research 379 (2019) 1–11

Fig. 8. Shear-wave velocity variations with respect to the mean velocity in each layer: A) along the Nicaraguan volcanic chain; B) SW/NE profile across Masaya volcano together with the approximate depth location of the modelled density anomalies from Métaxian (1994). C) WN/ES profile across the Momotombo volcano together with the local seismicity. Areas that are not well constrained are masked. The locations of the vertical cross sections are indicated in Fig. 7. The dotted lines indicate the modelled gravimetric anomalies from Métaxian (1994).

features of relative location techniques, such as the source specific station correction term (Richards-Dinger and Shearer, 2000) with a waveform based location method. This location approach, which inherits typical features of relative location methods (Grigoli et al., 2016), is less dependent on the velocity model and topography effects introduced by the volcanic structure and presents several benefits, which improve the location accuracy: 1) it accounts for phase delays due to local site effects, e.g. surface topography or variable sediment thickness 2) theoretical velocity models are only used to estimate travel time within the source volume, and not along the entire source-sensor path. In addition, this approach inherits the main advantage of waveform based location techniques, the noise robustness. The idea behind this approach is quite simple. If the distance between the target events and the master event is much smaller than the average source-station distance, then their ray-paths can be considered approximatively identical, thus the travel time correction for the master event is also valid for the target events. This location approach reduces the dependence on the velocity model and is useful if, for a particular region, a detailed velocity model is not available. Since location uncertainties related to the master event are propagated to the target ones it is necessary to choose the best located event (those with the best signal to noise ratio at all stations of the network) to be the master event. For the Momotombo area, earthquake locations were initially performed using waveform stacking, with a velocity model extracted from the CRUST2.0 database (Bassin et al., 2000). In a second step, we used the best located events (i.e. those with large magnitude) and calculated the source specific correction terms at each station. It is important to note that, even after the correction, the coherence level does not change significantly. This is mainly due to the fact that many stations of the network deployed at the Momotombo volcano are characterized by a very low signal-to-noise ratio, therefore they do not contribute to the coherence enhancement during the waveforms stacking process. With this method, we could identify two distinct clusters of seismicity in the area of Momotombo (Fig. 1C).

6. Discussion 6.1. Masaya volcano 6.1.1. Spatial extent and depth of the magma plumbing system At the Masaya volcano, the tomography results give no indication of a deep-rooted magma plumbing system. The low velocities in the caldera region are constrained to the upper 3 km of the subsurface. In the detailed tomography (Fig. 6), we observed that the largest negative shear-wave anomaly is offset to the west and not directly below the active Santiago vent (location of the lava lake). The absolute Vs velocities range from 1.1 to 1.4 km/s. This result corresponds to the observation of (partially) molten crustal magma from other studies (Mordret et al., 2015; Obermann et al., 2016; Stankiewicz et al., 2010). This observation is also in agreement with observations from satellite geodesy by Stephens and Wauthier (2018), who observed up to 8 cm ground inflation in the Masaya caldera from November 2015 to September 2016. Stephens and Wauthier (2018) explain their data with a spherical magma reservoir offset to the NW from the active Santiago vent, extending to a depth of 3 km. The observed offset is suggested to be a result of preexisting caldera structures. The presence of a very shallow magma chamber at Masaya is also in agreement with gravimetric studies from the nineties, as well as past and recent gas flux observations (Aiuppa et al., 2018; Connor and Williams, 1989; Métaxian, 1994; Stoiber et al., 1986). Métaxian (1994) showed that the central part of the caldera including the Santiago and Masaya crater has a lower density (2.05–2.15 gcm−3) than the surrounding areas (Figs. 8B, 9B, anomaly c). This region of lower density values corresponds well with the largest negative shear-wave anomaly (10%) that we obtained from the detailed tomography around Masaya at 3 s. We therefore conclude that Masaya has a very shallow magma chamber that is offset to the west from the active Santiago vent. At lower periods, in the depth range of 3–6 km, we observe a shift of the low velocity anomaly towards the NE of the caldera (Fig. 8B). Gravimetric studies (Connor and Williams, 1989; Métaxian, 1994) showed

A. Obermann et al. / Journal of Volcanology and Geothermal Research 379 (2019) 1–11

9

Fig. 9. A) Bouguer anomaly calculated for a density of 2.15 gcm−3 from Métaxian (1994). B) Lateral projection of the inverted density anomalies (a,b,c) and the 2D group velocities at periods of 3 s for Masaya (from Fig. 6). The area with lowest density (c) corresponds to the area with the largest low velocity anomaly (−10%).

the existence of a positive gravity anomaly of 0.5 gcm−3, with densities of 2.55–2.65 gcm−3 towards the calderas northeastern corner (Fig. 9B, anomaly a). The main component of this anomaly has been modelled by Métaxian (1994) as a dense body centered at 6 km depth and 6 km thickness with a density contrast N0.5. g·cm−3 compared to the surrounding areas (Fig. 8B, anomaly a). Métaxian (1994) interpreted this anomaly as being related to a former intrusive body that appeared prior to the caldera formation, as the anomaly shape does not follow the caldera morphology. During their modelling, Métaxian (1994) introduced a second intrusive body at depths of 3–8 km below the Western flank of the caldera that is, however, poorly constrained (Fig. 8A, anomaly d). The anomalies approximate location corresponds to the positive shear-wave anomaly location that we observe beneath the South-Western part of the caldera. In his model, Métaxian (1994) additionally introduced thin (200 m) positive gravity anomalies close to the surface (Fig. 9B, anomaly b) that are thought to be due to the contribution of dense basaltic lava flows. They are easily identified in the morphology and their thickness is estimated at 300–350 m. Our results agree with the gravimetric study which suggested that the Masaya caldera was probably not formed by a collapse triggered by the emptying of a large magma chamber, but more likely by explosive dynamics, repeated over time. This eruptive process has been mentioned by several authors (Bice, 1980; Kieffer and Creusot-Eon, 1992; Williams, 1983). 6.1.2. Seismicity at Masaya During our seismic campaign from December 2016 to March 2017, we observed a low level of seismicity spread around the Masaya region. This could be coincidence, as increased seismic unrest had been observed prior to the appearance of the lava lake and throughout 2016 in the Masaya region (Tenorio, 2015; Tenorio, 2016). Contrary to Métaxian et al. (1997) who studied the seismicity of the Masaya volcano in the nineties and could show a low-level volcano-tectonic activity (1 VT per day on average) located at shallow depth in the eastern part of Santiago crater, we did not observe any volcano-tectonic events that locate in the vicinity of the crater. Power-spectral density noise plots show significant ‘high frequency’ noise between 2 and 2.5 Hz close to the Masaya crater that diminishes with increasing distance from the crater and is barely visible from a 3 km distance. This signal could be related to the continuous degassing, as it has been observed in geysering systems (Karyono et al., 2017). From the crater edge, we can observe that the lava is moving very fast on the surface of the lava lake. These large moving masses should generate an oscillation that could be related to the high frequency noise.

Stations at the crater rim showed two prominent frequency peaks at around 1.5 and 2 Hz that might be related to tremor signals as postulated by Métaxian et al. (1997). A much more in-depth analysis would be needed to confirm this hypothesis, however, which is beyond the scope of this paper. Longer-term seismic studies would be of great interest to conclude on the seismic activity related to the lava lake. 6.2. Momotombo volcano 6.2.1. Spatial extent and depth of the magma plumbing system Below the cone of the Momotombo volcano, we observe a vertically elongated, 4 km wide, low shear-wave velocity anomaly that reaches depths of about 8 km (Fig. 8A,C). This is likely the main magma chamber of the volcano. At shallow depths of up to 3 km, the low velocity anomaly extends laterally below the Southern flank of Momotombo (Fig. 8A). This low velocity anomaly is likely linked to the volcanic hostedgeothermal system that is exploited in this area since the seventies. This observation agrees with the conceptual model proposed by Kaspereit et al. (2016): hot geothermal fluids ascend around the volcanic intrusion and due to the hydraulic gradient migrate downward and laterally through permeable flanks of the volcano. Further towards the south, we see a shallow circular low velocity anomaly that is constrained to the top 2–3 km. This anomaly occurs below Momotombito, a small stratovolcano 15 km off the coast from Momotombo (Fig. 8A). Likely the region between Momotombo and Momotombito is relatively weak, giving rise to the occurrence of active fluid flow, that is absent below the other flanks of Momotombo. 6.2.2. Seismicity around Momotombo Most of the seismic activity at Momotombo is concentrated approximately 10 km NW from the volcanic complex (Figs. 1C, 7, 8A). Here, the seismicity likely consists of aftershocks of the M5.5 earthquake that struck this area on September 28th, 2016, about 15 km WNW from the top of the volcanic complex (Fig. 1B). Epicentral locations in this area show a trend along an ENE-WSW direction reflecting the strike of the hosting faults system (Fig. 1C). A second cluster of seismicity can be found between the volcano and the Momotombito Island (Figs. 1C, 8A), where the geothermal plant is located and the M6.3 took place in 2014. At this location, seismicity is more clustered and shows a cylindrical shaped spatial distribution, which may be linked to non-tectonic processes. In both areas the hypo-central depth of seismicity is quite shallow, mostly between 3 and 10 km (Fig. 8A,C). The seismicity tends to avoid

10

A. Obermann et al. / Journal of Volcanology and Geothermal Research 379 (2019) 1–11

the areas of pronounced low shear-velocity and is more abundant in the better consolidated rock volumes. 7. Conclusion From this short-term seismic installation we gained a first glimpse of the magma plumbing systems at the Masaya and Momotombo volcanoes. While Momotombo shows a typical elongated low shear-wave velocity anomaly that reaches depths of about 8 km, Masaya does not show indications of a deep plumbing system. At Masaya, we observe the largest negative shear-wave anomaly at shallow depth (0–3 km) offset to the west and not directly below the active Santiago vent, where the crater lake is contained. We interpret this anomaly as a shallow magma chamber. At greater depth (3–8 km) a low velocity anomaly towards the northeast coincides in location with a modelled positive gravity anomaly and could indicate the presence of a former intrusive body. Seismicity around Momotombo shows a trend along an ENE-WSW direction reflecting the strike of the hosting faults system. A cluster below located between the volcano and Momotombito Island could be linked to the geothermal exploitation, but cannot be definitely confirmed given the short duration of our seismic installation. With this study, we hope to trigger further interest in the diverse tectonic and volcanic features of Nicaragua. Future, long-term seismic imaging and monitoring projects are of critical interest for the estimation of seismic and volcanic risks in Managua and the surrounding areas. Acknowledgments This work is supported by the government of Nicaragua on behalf of the Instituto Nicaraguense de Estudios Territoriales (INETER). The GeoForschungs-Zentrum Potsdam, Germany (GFZ-Potsdam) provided the 30 mobile seismic broad band stations from its geophysical instruments pool. The cooperation between SED/ETHZ and INETER is promoted and supported by the Swiss Agency for Development and Cooperation DEZA. The research leading to these results has received funding from the European Community's Seventh Framework Programme under grant agreement No. 608553 (Project IMAGE) and the Swiss Federal office of Energy with the project GEOBEST. F.G. is funded by the European Union's Horizon 2020 research and innovation programme under the Marie Skodowska-Curie grant agreement No. 790900. We would like to thank Ivan Koulakov and an anonymous reviewer for their comments that greatly helped to improve this manuscript. References Aiuppa, A., de Moor, J.M., Arellano, S., Coppola, D., Francofonte, V., Galle, B., Giudice, G., Liuzzo, M., Mendoza, E., Saballos, A., et al., 2018. Tracking formation of a lava lake from ground and space: Masaya volcano (Nicaragua), 2014–2017. Geochem. Geophys. Geosyst. 19 (2), 496–515. Alvarez, L., Rodriguez, A., Gonzalez, O., Moreno, B., Cabrera, A., 2018. Seismotectonics of the Nicaraguan depression from recent seismicity. J. Geol. Geophys. 7 (446), 2. Barmin, M., Ritzwoller, M., Levshin, A., 2001. A fast and reliable method for surface wave tomography. Monitoring the Comprehensive Nuclear-Test-Ban Treaty: Surface Waves. Springer, pp. 1351–1375. Bassin, C., Laske, G., Masters, G., 2000. The current limits of resolution for surface wave tomography in North America. Eos 81. Bice, D.C., 1980. Origin of the Masaya Caldera, Nicaragua. Geological Society of America Abstract program 12 (3), 98. Bice, D.C., 1985. Quaternary volcanic stratigraphy of Managua, Nicaragua: correlation and source assignment for multiple overlapping plinian deposits. Geol. Soc. Am. Bull. 96 (4), 553–566. Brenguier, F., Shapiro, N.M., Campillo, M., Nercessian, A., Ferrazzini, V., 2007. 3-D surface wave tomography of the Piton de la Fournaise volcano using seismic noise correlations. Geophys. Res. Lett. L02305. Brown Jr., R., Ward, P.L., Plafker, G., 1974. Geologic and seismologic aspects of the Managua, Nicaragua, earthquakes of December 23, 1972: US geological survey professional paper 838, 1973. Bull. Seismol. Soc. Am. 4 (1031–1031). Burton, M.R., Oppenheimer, C., Horrocks, L.A., Francis, P.W., 2000. Remote sensing of CO2 and H2O emission rates from Masaya volcano, Nicaragua. Geology 28 (10), 915–918.

Cailleau, B., LaFemina, P., Dixon, T.H., 2007. Stress accumulation between volcanoes: an explanation for intra-arc earthquakes in Nicaragua? Geophys. J. Int. 3, 1132–1138. Carr, M.J., 1984. Symmetrical and segmented variation of physical and geochemical characteristics of the Central American volcanic front. J. Volcanol. Geotherm. Res. 20 (3–4), 231–252. Connor, C.B., Williams, S., 1989. Interpretation of gravity anomalies, Masaya Caldera Complex, Nicaragua. Preprint: Transaction of the 12th Caribbean Conference, St Croix, Virgin Islands. Delmelle, P., Stix, J., Baxter, P., Garcia-Alvarez, J., Barquero, J., 2002. Atmospheric dispersion, environmental effects and potential health hazard associated with the lowaltitude gas plume of Masaya Volcano, Nicaragua. Bull. Volcanol. 64 (6), 423–434. DeMets, C., 2001. A new estimate for present-day Cocos-Caribbean plate motion: Implications for slip along the Central American volcanic arc. Geophys. Res. Lett. 28 (21), 4043–4046. Fallahi, M.J., Obermann, A., Lupi, M., Karyono, K., Mazzini, A., 2017. The plumbing system feeding the Lusi eruption revealed by ambient noise tomography. Journal of Geophysical Research: Solid Earth (10), 8200–8213. Frische, M., Garofalo, K., Hansteen, T.H., Borchers, R., 2006. Fluxes and origin of halogenated organic trace gases from Momotombo Volcano (Nicaragua). Geochem. Geophys. Geosyst. (5). Funk, J., Mann, P., McIntosh, K., Stephens, J., 2009. Cenozoic tectonics of the Nicaraguan depression, Nicaragua, and Median Trough, El Salvador, based on seismic-reflection profiling and remote-sensing data. Geol. Soc. Am. Bull. (11−12), 1491–1521. Grigoli, F., Cesca, S., Krieger, L., Kriegerowski, M., Gammaldi, S., Horalek, J., Priolo, E., Dahm, T., 2016. Automated microseismic event location using master-event waveform stacking. Sci. Rep. 6, 25744. Grigoli, F., Scarabello, L., Böse, M., Weber, B., Wiemer, S., Clinton, J.F., 2018. Pick-and waveform-based techniques for real-time detection of induced seismicity. Geophys. J. Int. 213 (2), 868–884. Karyono, K., Obermann, A., Lupi, M., Masturyono, M., Hadi, S., Syafri, I., Abdurrokhim, A., Mazzini, A., 2017. Lusi, a clastic-dominated geysering system in Indonesia recently explored by surface and subsurface observations. Terra Nova (1), 13–19. Kaspereit, D., Rickard, W., Osborn, W.L., Mann, M., Perez, M., 2016. Field management and expansion potential of the momotombo geothermal field using numerical simulation and conceptual modeling. Proceedings, 41st Workshop on Geothermal Reservoir Engineering. Kieffer, G., Creusot-Eon, A., 1992. La Caldeira de Masaya (Nicaragua): une dépression polyphasée de type maar. Comptes rendus de l'Académie des sciences. Série 2, Mécanique, Physique, Chimie, Sciences de l'univers, Sciences de la Terre. vol. 315 (11), pp. 1403–1409. LaFemina, P., Dixon, T.H., Govers, R., Norabuena, E., Turner, H., Saballos, A., Mattioli, G., Protti, M., Strauch, W., 2009. Fore-arc motion and Cocos ridge collision in Central America. Geochem. Geophys. Geosyst. 10 (5). LaFemina, P.C., Dixon, T.H., Strauch, W., 2002. Bookshelf faulting in Nicaragua. Geology 30 (8), 751–754. Levshin, A., Yanocskaya, T., Lander, A., Bukchin, B., Barmin, M., Ratnikova, L., Its, E., 1989. Seismic Surface Waves in a Laterally Inhomogeneous Earth. Springer-Verlag. Menyailov, I., Nikitina, L., Shapar, V., Pilipenko, V., 1986. Temperature increase and chemical change of fumarolic gases at Momotombo Volcano, Nicaragua, in 1982–1985: are these indicators of a possible eruption? Journal of Geophysical Research: Solid Earth 91 (B12), 12199–12214. Métaxian, J.-P., 1994. Etude sismologique et gravimétrique d'un volcan actif: dynamisme interne et structure de la Caldeira Masaya, Nicaragua. Ph.D. thesis. Université de Savoie. Métaxian, J.-P., Lesage, P., Dorel, J., 1997. Permanent tremor of Masaya Volcano, Nicaragua: wave field analysis and source location. Journal of Geophysical Research: Solid Earth 102 (B10), 22529–22545. Molinari, I., Verbeke, J., Boschi, L., Kissling, E., Morelli, A., 2015. Italian and Alpine threedimensional crustal structure imaged by ambient-noise surface-wave dispersion. Geochem. Geophys. Geosyst. 16 (12), 4405–4421. Mordret, A., Landès, M., Shapiro, N.M., Singh, S., Roux, P., Barkved, O., 2013. Near-surface study at the Valhall oil field from ambient noise surface wave tomography. Geophys. J. Int. 193 (3), 1627–1643. Mordret, A., Landès, M., Shapiro, N.M., Singh, S., Roux, P., 2014. Ambient noise surface wave tomography to determine the shallow shear velocity structure at Valhall: depth inversion with a Neighbourhood Algorithm. Geophys. J. Int. 198 (3), 1514–1525. Mordret, A., Rivet, D., Landès, M., Shapiro, N.M., 2015. 3-D shear-velocity anisotropic model of Piton de la Fournaise volcano (La Réunion Island) from ambient seismic noise. Journal of Geophysical Research: Solid Earth 120 (1), 406–427. Moschetti, M.P., Ritzwoller, M.H., Shapiro, N.M., 2007. Surface wave tomography of the western United States from ambient seismic noise: Rayleigh wave group velocity maps. Geochem. Geophys. Geosys. 8 (Q08010). Obermann, A., Lupi, M., Mordret, A., Jakobsdóttir, S.S., Miller, S.A., 2016. 3D-ambient noise Rayleigh wave tomography of Snæfellsjökull volcano, Iceland. J. Volcanol. Geotherm. Res. 317, 42–52. Ranero, C.R., Huene, R., Flueh, E., Duarte, M., Baca, D., McIntosh, K., 2000. A cross section of the convergent Pacific margin of Nicaragua. Tectonics 19 (2), 335–357. Richards-Dinger, K., Shearer, P., 2000. Earthquake locations in southern California obtained using source-specific station terms. Journal of Geophysical Research: Solid Earth 105 (B5), 10939–10960. Richter, C.F., 1935. An instrumental earthquake magnitude scale. Bull. Seismol. Soc. Am. 25 (1), 1–32. Ritzwoller, M., Levshin, A., 1998. Eurasian surface wave tomography: group velocities. J. Geophys. Res. 103, 4839–4878.

A. Obermann et al. / Journal of Volcanology and Geothermal Research 379 (2019) 1–11 Rymer, H., de Vries, B.v.W., Stix, J., Williams-Jones, G., 1998. Pit crater structure and processes governing persistent activity at Masaya Volcano, Nicaragua. Bull. Volcanol. 59 (5), 345–355. Sambridge, M., 1999. Geophysical inversion with a neighborhood algorithm. I. Searching a parameter space. Geophys. J. Int. 138 (2), 479–494. Schliz, M., 2011. Analogue Models of the Nicaraguan Depression: Intra-Arc Rifting and Volcanic Activities. Ph.D. thesis. Shimane Univ., Interdisciplinary Faculty of Science and Engineering, Japan. Siebert, L., Simkin, T., Kimberly, P., 2011. Volcanoes of the World. Univ of California Press. Smithsonian Institution, 1993. Bull. Global Volcan. Network 18 (6), 4. Snyder, G.T., Fehn, U., 2002. Origin of iodine in volcanic fluids: 129i results from the Central American volcanic arc. Geochim. Cosmochim. Acta 66 (21), 3827–3838. Stankiewicz, J., Ryberg, T., Haberland, C., Natawidjaja, D., et al., 2010. Lake Toba volcano magma chamber imaged by ambient seismic noise tomography. Geophys. Res. Lett. 37 (17). Stephens, K., Wauthier, C., 2018. Satellite geodesy captures offset magma supply associated with lava lake appearance at Masaya Volcano, Nicaragua. Geophys. Res. Lett. 45 (6), 2669–2678. Stoiber, R.E., Williams, S.N., Huebert, B.J., 1986. Sulfur and halogen gases at Masaya Caldera Complex, Nicaragua: total flux and variations with time. Journal of Geophysical Research: Solid Earth 91 (B12), 12215–12231. Strauch, W., 2014. Resumen del Terremoto del 10 de Abril de 2014 y sus réplicas? Impacto en los alrededores del Lago Xolotlán y en Managua. Boletín mensual Sismos y Volcanes de Nicaragua? Abril de 2014 (Technical Report). Strauch, W., Talavera, E., Tenorio, V., Ramírez, J., Argüello, G., Herrera, M., Acosta, A., Morales, A., 2018. Toward an earthquake and tsunami monitoring and early warning system for Nicaragua and Central America. Seismol. Res. Lett. 89 (2A), 399–406.

11

Suárez, G., Muñoz, A., Farraz, I.A., Talavera, E., Tenorio, V., Novelo-Casanova, D.A., Sánchez, A., 2015. The 10 April 2014 Nicaraguan crustal earthquake: evidence of complex deformation of the Central American volcanic arc. Geodynamics of the Latin American Pacific Margin. Springer, pp. 3305–3315. Sultan, D.I., 1931. The Managua earthquake. The Military Engineer 23 (130), 354–361. Tenorio, V., 2015. Sismicidad? Volcán Masaya, Boletín mensual Sismos y Volcanes de Nicaragua Diciembre 2015 (Technical Report). Tenorio, V., 2016. Sismicidad? Volcán Masaya, Boletín mensual Sismos y Volcanes de Nicaragua Enero 2016 (Technical Report). Tenorio, V., Navarro, M., Strauch, W., 2015. Actividad Sísmica del 2013 al 2015 y Actividad eruptiva del volcán Momotombo en diciembre de 2015. Boletín mensual Sismos y Volcanes de Nicaragua Septiembre 2015 (Technical Report). Turner, I.I.I., Henry, L., LaFemina, P., Saballos, A., Mattioli, G.S., Jansma, P.E., Dixon, T., 2007. Kinematics of the Nicaraguan forearc from GPS geodesy. Geophys. Res. Lett. 34 (2). Viramonte, J.G., Incer-Barquero, J., 2008. Masaya, the ‘Mouth of Hell’, Nicaragua: volcanological interpretation of the myths, legends and anecdotes. J. Volcanol. Geotherm. Res. 176 (3), 419–426. Walther, C., Flueh, E.R., Ranero, C.R., Von Huene, R., Strauch, W., 2000. Crustal structure across the Pacific margin of Nicaragua: evidence for ophiolitic basement and a shallow mantle sliver. Geophys. J. Int. 141 (3), 759–777. Wathelet, M., 2008. An improved neighborhood algorithm: parameter conditions and dynamic scaling. Geophys. Res. Lett. 35 (9). Williams, S.N., 1983. Plinian airfall deposits of basaltic composition. Geology 11 (4), 211–214.