Improving Seismic Surveillance at Mt. Etna Volcano by Probabilistic

Earthquake Location in a 3D Model

by Antonino Mostaccio, Tiziana Tuvè, Domenico Patanè,* Graziella Barberi,and Luciano Zuccarello

Abstract The increasing accuracy of 3D velocity models developed recently forMt. Etna has enabled their use today in routine earthquake locations. In this work, wetested the potential and performance of a global-search probabilistic earthquake lo-cation method (NonLinLoc) in a 3D velocity model, to improve earthquake locationsfor seismic surveillance. In addition, NonLinLoc hypocenter locations and those ob-tained by standard iterative-linear 3D locations, SimulPS-14, have also been com-pared. To this end, a dataset of 328 selected earthquakes, occurring during the2002–2003 Etna flank eruption, and the recent highest resolution 3D velocity model,have been used. The results revealed that the differences in hypocentral coordinatesbetween the two methods are typically of the same order or smaller than the spatiallocation uncertainty. To evaluate the consistency of results between the two 3Dlocation algorithms, synthetic datasets with real source–receiver configuration are alsoconsidered. Furthermore, by using NonLinLoc we estimated the influence of thesource–receiver geometry on the quality of hypocenter locations. If we vary the net-work geometry in a dense and well-distributed network like at Etna, reducing thenumber of stations (by 20% and 50%), it is significant that no large systematic hypo-central shifts of the relocated earthquakes are observed if they occur within the net-work. NonLinLoc is a fast and promising approach for automatic earthquake locationsand surveillance purposes at Mt. Etna, because (1) it works well with a reduced num-ber of seismic pickings, which are usually available in the automatic locations; (2) it isnot particularly sensitive to tolerable levels of random noise in arrival times; and (3) itproduces full location uncertainty and resolution information with respect to standarditerative-linear 3D locations.

Introduction

Seismology and geodesy are generally viewed as themost reliable diagnostic tools for monitoring highly activeor erupting volcanoes. Many of the seismic observationsperformed at active volcanoes in the 1960s and 1970s wereobtained using only a small number of instruments. Sincethen, nearly all well-monitored volcanoes have beenequipped with at least a few (4–6) instruments and some evenwith dozens of instruments. The success of the most recentseismological studies relies on the accuracy of earthquakelocations. The goodness of hypocenter locations is governedby several factors, including the network geometry, availablephases, arrival-time reading accuracy, and knowledge of thecrustal structure (Poupinet et al., 1984; Got et al., 1994;Richards-Dinger and Shearer, 2000).

A layered velocity model that does not represent a 3Dvelocity structure can introduce significant errors into thecalculated travel times for earthquake location. In the past,in order to optimize earthquake locations, researchersdeveloped various techniques to reduce the effects of a3D velocity structure.

More recently, new approaches have been developedthat are not restricted to 1D or 2D velocity models and whichseek to locate the sources in a nonlinear, probability-basedmanner (Lomax et al., 2000). Unfortunately, there are nogood velocity models for monitoring and surveillance atmost volcanoes, and the inversion for the source location isstill frequently done using the Hypo family of computercodes and 1D velocity models, such as the Hypo71 (Leeand Lahr, 1972), Hypoinverse (Klein, 1978), or Hypoellipse(Lahr, 1989) programs. Most of the computed source coor-dinates, especially when focusing on shallow events, must beseen solely as an approximation of the true hypocenter.

*Also at Instituto Andaluz de Geofísica, Universidad de Granada,Granada, Spain.

2447

Bulletin of the Seismological Society of America, Vol. 103, No. 4, pp. 2447–2459, August 2013, doi: 10.1785/0120110202

From the early 1980s, seismic activity has been moni-tored at Etna volcano by a permanent seismic network thathas been progressively improved (Patanè et al., 2004). Since2001, the Sezione di Catania of Istituto Nazionale diGeofisica e Vulcanologia (INGV) has managed this network,which was initially made up of analog single-componentshort-period (1 s) stations. Of about 40 operating stations(Fig. 1; Patanè et al., 2004, 2006), 5 were analog 3-Cshort-period (1 s) and 3 were digital 3-C broadband (40 s)stations. As a rule, a 3-C broadband temporary array com-prising from 5 to 10 stations, depending on seismicity loca-tion and the advent of eruptive crises, is deployed to integratethe permanent network, and during the study period 10 sta-tions were used. The permanent network has been consider-ably enhanced since 2005 by 24-bit digital stations equippedwith broadband (40 s) sensors, and today, thanks to a con-figuration of 33 BB and 12 short-period stations, there is verygood coverage of the volcanic area and the collected data isof high quality.

Over the last thirty years, a number of investigations intothemain structural characteristics beneathEtna have been per-formed by several authors (i.e., Sharp et al., 1980; Hirn et al.,1991; De Luca et al., 1997; Chiarabba et al., 2000; Laigleet al., 2000; Patanè et al., 2002), and crustal velocity modelsusing data obtained with seismic refraction experiments andtomographic inversions of local and teleseismic earthquakeshave been developed. In recent years, more robust 3D velocitymodels have been proposed for Etna (Patanè, De Gori, et al.,2003; Chiarabba et al., 2004; Patanè et al., 2006), yieldingbetter constraints on earthquake locations.

These studies and insights allow better constraining ofthe seismic activity at almost all depths and enable tacklingthe relationship between seismicity and eruptive activity.

Notwithstanding these improvements, at the Sezione diCatania of the INGV, the automatic and refined (offline) hy-pocentral locations of earthquakes reported in the catalog arestill obtained using the Hypoellipse algorithm (Lahr, 1989)and 1D velocity model, taking into account the difference inaltitude of the seismic stations, and allowing hypocenterswithin the volcanic cone.

In this paper, we present the hypocentral results obtainedin a 3D model, considering the probabilistic nonlinear earth-quake location method (NonLinLoc) proposed by Lomaxet al. (2000) applied to the seismicity recorded during the2002–2003 eruptive period. We also compare these hypo-central results with the routine locations performed forsurveillance purposes by using the Hypoellipse algorithm(Lahr, 1989) and with those obtained with the standarditerative-linear 3D location method (SimulPS-14, Thurber,1983; Eberhart-Phillips, 1990), to evaluate if NonLinLoc isbetter in terms of both location accuracy and computationalefficiency, and more importantly, if it is less sensitive to noise.

Volcano Tectonic Earthquake Features at Mt. Etna

Mt. Etna is one of the most active basaltic volcanoes inthe world, with a historical record of documented eruptionsgoing back over 2000 years. It is located in eastern Sicily in acomplex geodynamic framework, where major regionalstructural lineaments play a key role in the dynamic proc-esses of the volcano (e.g., Ferrucci and Patanè, 1993; Patanèet al., 2004).

Similarly to other volcanic areas, the large variety ofseismovolcanic signals observed at Etna volcano (e.g., DelPezzo et al., 1993; Patanè et al., 2004) can be grouped into:(1) volcano tectonic earthquakes (VT), which are tectonic-like events, generated by regional tectonic stresses and/or bylocal stresses deriving from magma migration within theEarth’s crust; and (2) volcanic events, which are seismicmanifestations of fluid dynamics, including volcanic signalsdirectly connected with an eruption itself, magma outflow tothe surface, and explosive activity at the craters.

The VT earthquakes usually recorded at Etna are:(1) earthquakes located within the deep (ca. 15–30 km) andintermediate crust (ca. 5–15 km), with high-frequencycontent (mostly above 4–5 Hz), and (2) shallow earthquakeslocated at depths of less than 5 km below the sea level (b.s.l.),producing arrivals with medium-to-low frequency content(mostly between 4 and 2 Hz) and/or complex signatures atstations just a few kilometers away from the epicentralarea (Patanè et al., 1997; Patanè and Giampiccolo, 2004).At Etna, VT earthquakes mainly occur in swarms, whereasforeshock–mainshock–aftershock sequences are rarely re-corded. In terms of magnitude, this VT seismicity is madeup of micro or small earthquakes, which seldom exceed mag-nitude 4.5.

Figure 1. Map of the Mt. Etna area showing the seismic stationsoperating during the study period. The concentric white curves areelevation contours at 1000 m intervals. In the upper left inset is aschematic tectonic map of Sicily with the main active faults.

2448 A. Mostaccio, T. Tuvè, D. Patanè, G. Barberi, and L. Zuccarello

All these earthquakes are believed to result fromregional tectonics and/or local stresses determined by excessmagmatic pressure on the surrounding rocks (e.g., Musumeciet al., 2004; Alparone et al., 2012), providing the necessaryenergy for rock failure, while no fluids are generally involvedin the source processes.

For example, the July–August 2001 flank eruption washeralded by one of the strongest seismic swarms recorded atEtna in the last 30 years. This intense VT seismicity, begin-ning on 12 July, was associated with a rapidly rising dike andopening of a complex fracture system that preceded the erup-tion on 17 July (Patanè et al., 2002; Patanè, De Gori, et al.,2003, Bonaccorso et al., 2004). In total, from 12 July to 9August, when the eruption ended, 2694 low-energy earth-quakes (1 ≤ MD ≤ 3:9) were recorded by the local perma-nent seismic network, but most (2645) occurred in the123 hours before the eruption onset (Patanè, Privitera, et al.,2003). Fifteen months after this eruptive episode ended, anew eruption began abruptly on 26 October 2002. Only afew hours of premonitory seismicity preceded and accompa-nied the opening of eruptive fissures along a biradial direc-tion, on both the northeastern and southern flanks of thevolcano. The local permanent network was integrated withfurther digital 3-C broadband (60 s) stations immediatelyafter the eruption onset (Fig. 1). 874 earthquakes were re-corded (1:0 ≤ MD ≤ 4:4) until the eruption end on 28 Janu-ary 2003. All these events included strike-slip and dip-slipearthquakes on normal and thrust faults (Gambino et al.,2004; Patanè et al., 2005).

Seismic Data and Hypocentral Locations withHypoellipse and SimulPS

In this study, we examined the seismicity during the2002–2003 Etna flank eruption.

Given the low magnitude and shallow depths, a goodsignal-to-noise ratio is mostly available only at the closeststations. First P-wave arrivals were also often masked by thecoda waves of a previous event and/or by the high tremoramplitude, which increases markedly at the eruption onsetin concomitance of the lava outflow (e.g., Patanè, Privitera,et al., 2003). To tackle the problem, the arrival times and S-phase were picked manually together with P-onset polaritydeterminations. Then, a dataset of 328 earthquake locations(MD ≥ 1:1), obtained by the Hypoellipse code (Lahr, 1989)with the 1D velocity model derived from Hirn et al. (1991),was selected for the study.

Several previously performed tomographic studies atEtna (i.e., Patanè et al., 2002; Patanè, De Gori, et al., 2003;Chiarabba et al., 2004), highlighted the pronounced lateralheterogeneity of the crust beneath the volcano that can stronglyinfluence earthquake locations. This prompted the need to re-locate Etna’s seismicity in a 3D model. To obtain the high-resolution 3D velocity model of the volcanic area in the uppercrust, at depths of less than 10 km b.s.l., Patanè et al. (2006)considered a large dataset (2001–2003), which naturally

already included the 2002–2003 seismicity. A total of 712well-constrained earthquakes, occurring from 10 August 2001to 18 January 2003, were used, and 8587 P-wave arrivals, and2293 S-wave arrivals were inverted by using the SimulPS-14software (hereafter named SimulPS) with a grid node spacingof 2 × 2 × 1 km within a volume of 36 × 36 × 18 km.

The inversion method developed originally by Thurber(1983) and modified by Eberhart-Phillips (1993) andEberhart-Phillips and Reyners (1997), solves the nonlinearcoupled hypocenter-velocity problem by a linearized itera-tive damped least-squares scheme. Each iteration consists ofan inversion for VP and optionally VP=VS variations, and forhypocenter locations.

The several tests performed at Etna demonstrated that agrid spacing of 2 × 2 × 1 km is the best compromise amongmodel parameterization, spatial resolution, and a reliable rep-resentation of the velocity field down to about 10 km depth.Conversely, finer grid spacing causes the introduction ofartifacts and false anomalies.

3D velocity model developed by Patanè et al. (2006)(Fig. 2) is an improvement on previous tomographic modelsfor Etna and provides better definition of the shape andgeometry of the upper portion of the high-velocity VP vol-ume, located below the south and southeastern part of thevolcano, and interpreted as a solidified intrusive body (seeChiarabba et al., 2004, and references therein). This high-VP body, already recognized by the first tomographic inves-tigation at Etna (Hirn et al., 1991), represents the main struc-tural feature beneath the volcano, which is embedded in thecontinental crust together with thick sedimentary rocks.

If we compare 3D SimulPS and 1D Hypoellipse loca-tions, the former show an increase in clustering of seismicityand some clear linear features are also better resolved. More-over, the root mean square (rms) residuals and hypocentralerrors with SimulPS are improved, resulting less than 0.15 sand 0.5 km, respectively. In Figure 3a,b, maps and west–eastsections of the selected earthquakes located using Hypoel-lipse and SimulPS are shown.

The epicenters of events were spread over a wide area,mostly in the eastern part of the volcano, limited by the North-east Rift and Pernicana fault system to the north and by theSanta Venerina–Santa Tecla fault system (Timpe fault) tothe south. The focal depths are constrained within 5 km b.s.l.

Although the 1D locations seem to bewell constrained, itis evident that the 1D model is unable to represent the volca-no’s complex structure and many events are forced to the topof the model (Fig. 3a). Conversely, using a 3D model allowsbetter constraining of the locations, which are now releasedand the depths of the events reach up to 3 km above thesea level (a.s.l.) (Fig. 3b).

Probabilistic Earthquake Hypocenter Relocations inthe 3D Velocity Model

The 328 earthquake dataset was further relocated usingthe NonLinLoc algorithm (NLLoc hereafter, Lomax et al.,

Improving Seismic Surveillance at Mt. Etna Volcano by Probabilistic Earthquake Location in a 3D Model 2449

2000). The NLLoc program is able to manage the aforemen-tioned 3D VP model in the relocation process. The posteriorprobability density function (PDF) as computed by NLLocrepresents a complete probabilistic solution to the earthquakelocation problem, including information on uncertainty andresolution (Lomax et al., 2000).

NLLoc follows the probabilistic formulation of inver-sion given in Tarantola and Valette (1982) and Tarantola(1987) and the equivalent methodology for earthquake loca-tion (i.e., Tarantola and Valette 1982; Moser et al., 1992;Wittlinger et al., 1993; Lomax et al., 2000). The software(Lomax et al., 2000) can be used with any available velocitymodel (1D, 2D, or 3D) and the results represented by a pos-teriori PDF that can be computed via: (1) a systematic gridsearch, (2) a stochastic Metropolis-Gibbs search, or (3) ahybrid “Oct-Tree” method.

First, we applied the Oct-Tree sampling algorithm(Lomax and Curtis, 2001), which uses recursive subdivisionand sampling of rectangular cells in a 3D space to generate acascade structure of sampled cells. This algorithm was usedto determine a location PDF by first taking a set of sampleson a coarse regular grid of cells throughout the searchvolume. This was followed by a recursive process that takesthe cell k having the highest probability Pk of containing theevent location, and subdividing this cell into eight child cells(hence the name Oct-Tree), from which eight new samples ofthe PDF were obtained. These samples were added to a list ofall previous samples and then the highest probability cell wasonce again sought. This process was continued until a pre-determined number of samples was obtained, or until anothertermination criterion was reached.

Figure 2. VP velocity model (Patanè et al., 2006) in the well-resolved layers. In the 2 km layer the nodes (black cross) and thecenter (blue circle) of the horizontal grid used are shown. The white-gray lines are elevation contours (every 1000 m), and the red tri-angle indicates the central craters of the volcano. The west–eastcross section of the VP model crossing the grid center (red horizon-tal line) is shown for the −7 km layer. Adapted with permissionfrom Patanè et al., 2006.

Figure 3. (a) Epicentral map and west–east cross section ofevents located using Hypoellipse code and the 1D model. (b) Epi-central map and west–east cross section of events located by usingSimulPS code and 3D velocity model.

2450 A. Mostaccio, T. Tuvè, D. Patanè, G. Barberi, and L. Zuccarello

The results give efficient, accurate, and complete map-ping of earthquake location PDFs in a 3D space, and thealgorithm is roughly 100 times faster than the systematic gridsearch. The number of sampled cells follows the values of thePDF, thus leading to a higher density of smaller cells in areasof higher PDF (lower misfit). All the events were relocatedusing the same parameters, such as number of samples drawnfrom PDF (i.e., 1000).

NonLinLoc uses a gridded representation of the 3Dmodel consisting of constant velocity, with cubic cells of1 km side length, while the SimulPS uses a smooth cubicB-spline parameterization on a 2 × 2 × 1 km grid in the x,y, and depth directions, respectively.

The locations were performed within a volume of36 × 36 × 18 km, centered at 37°73 in latitude and 15°00in longitude located about 1 km south of the crater area. Bothin the horizontal and vertical components, a 1 km step wasused. The grid nodes in depth direction were positioned from3.0 km a.s.l. to 15 km b.s.l.

We used all available P phases (4251 observations),including the arrival times of lower quality (e.g., low signal-to-noise ratio recorded during the eruptive periods).

Concerning the S phases, as 3D S-wave velocity has notbeen considered and the use of a global constant VP=VS ratiomay lead to biased absolute hypocenter locations (Maurerand Kradolfer, 1996; Husen et al., 2003), we selected onlythe best quality (good signal-to-noise ratio) S-phase picks(1024 observations).

In Figure 4a, we plot the map and sections of NLLocrelocated earthquakes that, similarly to those obtained bySimulPS-14, prove more clustered compared to the 1D loca-tions (Fig. 3a).

Because the Oct-Tree algorithm may not detect narrowlocal minima in the PDF, we used the grid-search algorithm toobtain a complete mapping of the multiple minima in thePDF of the earthquake locations. Finally, the results fromindividual locations were stacked to obtain a compositelocation probability for all the relocated earthquakes (Saccor-otti et al., 2001). In Figure 4b, we show the composite mar-ginal probability distributions projected over the horizontaland the west–east vertically oriented planes.

Comparison between NLLoc, SimulPS, andHypoellipse Earthquake Locations

In this section, we compare the relocation results by theglobal-search probabilistic method (NLLoc) in the 3D modelwith those obtained by the standard iterative-linear 3DSimulPS, as well as by those using a standard algorithm suchas Hypoellipse in a 1D model.

Figure 5a–e shows the differences in depth, latitude,longitude, rms, and GAP between 1D Hypoellipse versus3D NLLoc and 3D SimulPS versus 3D NLLoc.

Comparing 1D Hypoellipse and 3D NLLoc, we note thesystematic shift of NLLoc hypocenters toward deeper foci(Fig. 5a). The median value is 330 m (standard deviation

1820 m), even though larger individual shifts can also beobserved (30% of the events show a >1:0 km shift). Nosignificant shifts are observed in latitude (Fig. 5b), whilelongitudinal movements of foci westwards (the median valueis 110 m with a standard deviation of 897 m; Fig. 5c) areclearly recognizable. Overall, the observed shift in epicenterlocations is less than�1:0 km (both in latitude and longitudefor about 80% of the events). By using NLLoc, there is aclear improvement in rms and GAP calculation (Fig. 5d,e).GAP and rms carried out by Hypoellipse are greater thanthose using NLLoc for about 80% of the hypocenters.

A systematic shift to greater depths is also observedcomparing 3D SimulPS hypocenter locations and 3D NLLoc(Fig. 5a). Even if the median value is 155 m and the standarddeviation is 1906 m, 25% of hypocenters show movements>1:0 km in the focal depth. Regarding latitude (Fig. 5b),about 60% of the events move southwards (the median valueis 260 m with standard deviation 1650 m), while in longitude(Fig. 5c) about 60% of the events move westwards, with amedian value of 170 m and standard deviation 1300 m. Thismeans that, overall, the epicenters move toward the south-west inside the seismic network. This can clearly be observedin Figure 5e, where the GAP values of relocated earthquakesare regularly lower than those obtained using SimulPS-14code. Finally, the values of rms (Fig. 5d) are systematicallygreater than those obtained with SimulPS.

Overall, on comparing 3D SimulPS and 3D NLLochypocenters, there is a relatively small systematic shift of focitoward greater depth, even if the maximum depth of seismic-ity remains unchanged. The relocated shallow seismicitymoves to greater depth because such events are more difficult

Figure 4. (a) Epicentral map and west–east cross section ofevents located by using NLLoc code and 3D velocity model.(b) Composite marginal probability distribution projected overthe horizontal plane and on the west–east vertical plane. Accordingto the column scale shown in the figure, probability increases fromlighter to darker colors. In particular, the green area marks the 90%confidence region for source location (i.e., the region inside whichthe probability of finding a hypocenter is 90%).

Improving Seismic Surveillance at Mt. Etna Volcano by Probabilistic Earthquake Location in a 3D Model 2451

to locate with linear methods, due to the lack of stations closeto the focal depth.

With respect to the 1D Hypoellipse locations, the relo-cated hypocenters have improved accuracy using the 3Dmodel and the probabilistic method. Indeed, we note that a3D probabilistic location is more stable (as suggested bycomplete solution information), even near velocity gradients

and outside of but not too far from the network. At Mt. Etna,3D NLLoc and 1D Hypoellipse locations show significantdifferences mainly in the depth determination and for earth-quakes occurring in proximity of the wide high-velocity VP

volume, located in the central-southeastern part of thevolcano. This suggests that the locations obtained by usinga layered model can give acceptable location and uncertaintyresults for those earthquakes lying within a seismic network,if the medium is not strongly heterogeneous. If we considerlocations obtained by NLLoc with the same layered modelused by Hypoellipse and compare them with the 3D NLLoc,significant differences in depths and minor variations in thelatitude and longitude are still observed.

3D Velocity-Model Performance

The difference between nonlinear and linear earthquakelocation is that linearized algorithms are more likely to pro-duce unsatisfactory results when the location problem ispoorly conditioned. Linearized locations are known to gen-erate accurate hypocenter estimates only when network cov-erage around an event is good, reading errors are small, andthe local velocity model is accurate. Conversely, nonlinearlocation algorithms are found to be a more reliable techniquewhen presented with poorly conditioned location problems(Presti et al., 2008).

Therefore, if a dense seismic network is available, like atEtna, good earthquake location software should be muchmore dependent on the accuracy of the velocity-model esti-mation than on the quality of the arrival-time data, providedthat the noise is not too severe.

To evaluate the accuracy and the robustness of the 3Dvelocity model with the two location techniques, we exam-ined various earthquake location scenarios using synthetictravel times for the 3D model of Patanè et al. (2006). Pre-dicted travel times are computed by the Time2EQ program,with the same Eikonal finite-difference algorithm (Podvinand Lecomte, 1991) used for the NLLoc event locations.These synthetic arrival times are used to relocate the eventsboth with SimulPS and NLLoc, adding variable randomnoise. The synthetic dataset comprises 7 events, each of themconsisting of 49 simulated P and 10 S arrival times atthe network stations. The foci have been positioned in thevolcanic area at different depths, between 0.2 km a.s.l. and10 km b.s.l. (Fig. 6). We generated the set of exact noise-free arrival times; then the 7 events were relocated withSimulPS and NLLoc and the results compared (Fig. 6). Fur-thermore, to analyze the sensitivity of the location procedureto data error, we added noise to the simulated arrival times.Three different levels of random noise were tested, havingvalues of �0:1, �0:2, and �0:5 s.

All three location methods were tested with these threerandom noise levels. The model and source–receiver layoutwere the same as the earlier noise-free cases. For compari-son, the location results in the noise-free case were taken asthe reference.

Figure 5. Differences in (a) depth, (b) latitude, (c) longitude,(d) rms, and (e) GAP between 1D Hypoellipse (black histograms)and 3D SimulPS (gray histograms) versus 3D NLLoc locations, re-spectively.

2452 A. Mostaccio, T. Tuvè, D. Patanè, G. Barberi, and L. Zuccarello

From Figure 6a–c, the following observations can bemade: the NLLoc location approach is less sensitive to therandom noise in the arrival times with respect to the SimulPSapproach, but in both cases hypocentral and origin time es-timation errors will linearly increase with increasing randomnoise. Then, the differences in the focal coordinates betweenthe NLLoc and SimulPS noise-free synthetic locations aretypically the same size as the spatial location uncertainty(Fig. 6a–c). The location error in depth has an accuracy com-parable to that in the horizontal plane, at least for the stationdistribution and earthquake geometry used in the synthetics.The NLLoc algorithm seems more efficient compared toSimulPS, as the results indicate that it is less sensitive to tol-erable levels of random noise in the arrival times.

NLLoc: Event Relocations and Uncertainty

To better define the relocation quality with NLLoc foreach earthquake of the dataset, we classified all earthquakesin four quality classes (QC) A, B, C, and D (Table 1), with Aindicating very good location quality and D the poorest qual-ity (unacceptable locations).

To find the events falling within the QC-D, we definedfor each event a quality factor (Q) of the event calculated byusing a comprehensive set of normalized parameters ofuncertainty:

Q � rmsn ×GAPnNn

0:5 ≤ Q ≤ 4;

where the values of rms, GAP, and N are normalized(n subscript) for each event as follows:

rmsn �rms − rmsmin

rmsmax − rmsmin� 1 1 ≤ rmsn ≤ 2

GAPn �GAP − GAPmin

GAPmax − GAPmin� 1 1 ≤ GAPn ≤ 2

Nn �N − Nmin

Nmax − Nmin� 1 1 ≤ Nn ≤ 2:

The Q includes information on the network geometry(GAP), number of observations (N) and their residuals

Figure 6. (a) Map and (b) west–east cross section showing the comparison between 3D NLLoc (gray circles) and 3D SimulPS (whitediamonds) locations for 7 synthetic events. In the central and right graphs, the differences in the hypocentral coordinates between 3D NLLoc(black squares) and 3D SimulPS (black triangles) locations, obtained by adding random noise to the simulated arrival times for each eventand the noise-free synthetics, are also shown. To perturb synthetic arrival times, three different levels of random noise (�0:1 s [dark graylocations], �0:2 s [light gray locations] and �0:5 s [white locations]) have been considered.

Table 1Definition of Quality Classes (QC) for Earthquake Locations

of Mt. Etna

Quality Class (QC) Selection Criteria

A (very good) Q ≤ 2:0; DIFF < 0:5 km;AE < 1:0 km

B (good) Q ≤ 2:0; DIFF < 0:5 km; AE ≥ 1:0 kmC (reasonable) Q ≤ 2:0; DIFF ≥ 0:5 kmD (unacceptable locations) Q > 2:0

Q corresponds to the quality factor average value; DIFF is the differencebetween maximum likelihood and expectation hypocenter locations; AE isthe average length of the three axes of the 68% error ellipsoid(average error).

Improving Seismic Surveillance at Mt. Etna Volcano by Probabilistic Earthquake Location in a 3D Model 2453

(rms), so that we are able to define the poorest locationsthrough a single value.

The Q threshold equal to 2 (average GAP of 130°� 50,average N of 15� 6, and average rms of 0:15 s� 0:07) isthe factor for which our locations can be considered accept-able. Therefore, we chose to put all the events with Q > 2:0in QC-D.

NLLoc produces two different hypocenter locations, thefirst given by its maximum-likelihood value (probabilisticsolution), while the second is the traditional Gaussian esti-mation of the hypocenter location (expectation-hypocenterlocation).

Large differences between the expected and maximum-likelihood hypocenter location, as estimated from the PDF,can be assumed the result of an ill-conditioned location prob-lem (Lomax et al., 2000). In this case, as already observed byHusen and Smith (2004), Gaussian location estimation doesnot represent adequate uncertainty estimations.

Then, following Husen and Smith (2004), we also basedour classification on two further uncertainty parameters:(1) the difference between the maximum-likelihood and theexpected hypocenter locations (DIFF), and (2) the averagelength of the three axes of the 68% error ellipsoid (averageerror, AE).

Considering also these two uncertainty parameters:(1) for a threshold value of DIFF ≥ 0:5 km, we included thelocated events in QC-C; (2) for a threshold value ofDIFF < 0:5 km and AE ≥ 1 km, the located events areincluded in QC-B; and (3) for a threshold value ofDIFF < 0:5 km and AE < 1 km, the located events areincluded in QC-A.

The results of this classification are shown in Table 2.We note that only 8% of the events may be considered unus-able (QC-D); instead, most events are in QC-A (47%), and theremaining 45% are shared between the other two intermedi-ate classes (28% in QC-C and 17% in QC-B). The averagevalues of the most important location parameters are reportedfor each quality class in Table 2. It is noteworthy that theextreme classes A and D are very well defined, whereas thereis no clear distinction for the intermediate classes. The eventsin the first three classes can be considered well constrained ifwe use the gap between class C and class D.

To show the relationship between our classification anduncertainty, we selected one representative event to plot foreach quality class. In Figure 7a–d, the location uncertaintiesby 3D density-scatter plots, obtained by drawing samplesfrom the PDF with the number of samples proportional tothe probability, are represented (black points). Moreover,for each event we show the 90% confidence region for sourcelocation (gray clouds), the maximum-likelihood hypocenterlocation (filled circles), and the expectation-hypocenter(Gaussian estimates) location (filled squares); an overviewof the network geometry (black triangles) used to locate theevents is also shown. In particular, the QC-A event (Fig. 7a)has a well-defined PDF (volume of about 0:3 km3) with verylittle difference between expectation- and maximum-likelihoodhypocenter location (about 70 m). In this case, the probabilisticlocation is well constrained because stations are distributedover a wide range of distances and azimuths. This exampleallows obtaining good quality parameters and low averageerror (0.4 km). As a whole, all the events classified inQC-A have a very well-defined PDF, small location uncer-tainties, and very good location parameters.

For the events classified in QC-B (Fig. 7b), the differencesbetween expectation and maximum-likelihood hypocenterlocations are also very close (about 250 m on average), eventhough, in this case, the seismic network coverage does notconstrain the probabilistic location well and the PDF cloud islarger than for class A. In this case, two stations, very close tothe epicenter, furnish good quality parameters, even thoughthe location falls in class B.

A typical example of a QC-C event is illustrated inFigure 7c. Here, a wide PDF density cloud around the maxi-mum-likelihood hypocenter location is evident, and, with re-spect to the expected hypocenter, there is a clear shift (morethan 1 km) in depth. Also in this case, the presence of a sta-tion sufficiently close to the epicenter helps better constrainhypocenter parameters. Therefore, we consider these prob-abilistic locations in QC-C as acceptable.

Finally, in Figure 7d, we show the PDF cloud of an eventincluded in QC-D. Here, the two maxima obtained in the PDFlocation are far apart. Lomax et al. (2009) suggest that thisshape of PDF is (1) typical of locations on the border of oroutside the recording network (see Fig. 7d) and (2) near a

Table 2Average Values (AV) of the Most Important Location Parameters for Each Quality Class (QC)

AV (QC) no No GAP (°) rms (s) ERH (km) ERZ (km) Q DIFF (m) AE (m)

A 155 (47%) 18 111 0.17 0.3 0.1 1.3 156 616B 97 (30%) 11 170 0.12 1.1 0.8 1.6 268 1373C 34 (10%) 12 148 0.11 2.4 1.4 1.5 1215 1620D 27 (8%) 9 249 0.18 2.4 2.1 2.2 1266 1894

no represents the total number of earthquakes falling within each quality class (QC); No, GAP, rms, ERH, andERZ represent the location uncertainty; Q corresponds to the quality factor average value; DIFF is the differencebetween maximum likelihood and expectation hypocenter locations; AE is the average length of the three axes ofthe 68% error ellipsoid (average error).

2454 A. Mostaccio, T. Tuvè, D. Patanè, G. Barberi, and L. Zuccarello

sharp horizontal interface in the velocity model across whichthere is a strong velocity contrast. Generally, in these cases,the linearized location procedure never identifies or exploresthe primary maximum of the PDF above the sharp interface,and produces incorrect error information above this interface.A probabilistic direct global-search procedure can determinethe complete location PDF and correctly identify the maxi-mum-likelihood hypocenter located above the sharp interface.The poor representation of the PDF and the bad parameterquality make the locations of the QC-D events unacceptable.

Importance of Source–Receiver Geometry

To verify the effects of the source–receiver geometry, werandomly reduce the number of stations by 20% and 50% foreach set and analyze the obtained results of relocated eventswith NLLoc.

Because we are interested in understanding the variationin earthquake locations by randomly worsening the source–receiver geometry, we keep the relocated events in the pre-viously defined quality classes; then both the average valuesof the location parameters and absolute displacement of thesource (Δ) are recalculated.

In Table 3, the results of a 20% reduction in the numberof seismic stations, chosen at random, are given. The Qaverage value, for every quality class, remains approximatelythe same. The results reveal a moderate deterioration in theaverage values of N and GAP, and a corresponding improve-ment in rms. However, the rms value alone is inadequate tocharacterize the goodness of earthquake locations, as thisparameter provides a measure of the fit of the observedarrival times to the predicted arrival times for the locationand depends on (1) the accuracy of the velocity model;(2) the quality weights assigned to the arrival-time data, and(3) the earthquake location procedure. Therefore, a smaller Ncould give better rms values, as it allows a better fit ofthe data.

In particular, when we apply a small reduction (20%) tothe number of seismic stations, the events belonging toQC-A, -B, and -C suffer only slight changes, and the absolutemovements (Δ) can be considered negligible (on average, nogreater than 1.0 km).

Instead, if we apply a greater reduction to the stationnumber (50%), a marked deterioration in the average valuesof location parameters can be observed (Table 4). Nonethe-less, the obtained final locations in QC-A and -B are onlyshifted relatively slightly with respect to the original earth-quake locations. Obviously, this occurs if we are dealing withdata from a dense seismic network.

Because even today automatic earthquake locations aremainly performed by P-phase arrivals alone, an example ofcomparison between an offline NLLoc earthquake location,in which S arrivals are added by an operator, and the corre-sponding automatic one, is shown in Figure 8.

Compared to the location uncertainties, the correspond-ing shift in epicenter and focal depth is relatively small,thereby indicating that adding S-wave arrivals improves thelocation quality. At Etna, the use of only P arrivals does notpose significant problems for earthquakes located inside thenetwork and recorded at a large number of high-quality sta-tions, at least in terms of epicentral location, but can be prob-lematic for low-magnitude earthquakes recorded at a lessernumber of stations. As even a few S observations can provideimportant constraints on the focal depth of an earthquakelocation, we will consider the automatic determination ofS-arrival times in the future.

Figure 7. Maps of the probabilistic location uncertainties (left)and 3D density-scatter plots of the PDF locations (right) of someexamples of events classified in the quality classes A, B, C, andD (see text for details). The gray clouds mark the 90% confidenceregion for source location (i.e., the region inside which the proba-bility of finding a hypocenter is 90%). The white circles indicate themaximum-likelihood hypocenter location, and the white squares in-dicate the expectation-hypocenter location (Gaussian estimates). Inthe map, the considered stations (black triangles) and the epicenterlocation are also shown.

Improving Seismic Surveillance at Mt. Etna Volcano by Probabilistic Earthquake Location in a 3D Model 2455

Table 3Results Obtained after Reducing the Number of Seismic Stations

by 20%

AV (QC) No GAP (°) rms (s) ERH (km) ERZ (km) Q Δ (m)

A 15 129 0.16 0.3 0.6 1.3 749B 9 178 0.11 1.5 1.8 1.6 908C 10 167 0.10 3.1 3.6 1.5 1124D 7 246 0.11 4.9 3.6 2.1 2505

AV is the average value of the most important location parameters for each qualityclass (QC); No, GAP, rms, ERH, and ERZ represent the location uncertainty; Qcorresponds to the quality factor average value; Δ is the absolute displacementof the source.

Table 4Results Obtained after Reducing the Number of Seismic Stations

by 50%

AV (QC) No GAP (°) rms (s) ERH (km) ERZ (km) Q Δ (m)

A 10 162 0.13 1.6 2.5 1.4 1952B 7 214 0.10 3.7 3.6 1.6 2732C 7 201 0.10 8.0 6.7 1.6 2371D 6 257 0.11 4.2 3.3 2.1 3466

AV is the average value of the most important location parameters for each qualityclass (QC); No, GAP, rms, ERH, and ERZ represent the location uncertainty; Qcorresponds to the quality factor average value; Δ is the absolute displacementof the source.

Figure 8. Epicentral maps of 3D NLLoc earthquake locations obtained by using only P (solid circle) and P- and S- (solid cross) arrivaltimes. Probabilistic location uncertainties, displayed by the 3D density-scatter plots of the location PDF, are shown for both cases. The grayclouds mark the 90% confidence region for source location; the white cross (location with S phases) and the white circle (location withoutS phases) indicate the maximum-likelihood hypocenter location.

2456 A. Mostaccio, T. Tuvè, D. Patanè, G. Barberi, and L. Zuccarello

Discussion and Conclusions

The main focus of this work was to investigate the per-formance of NLLoc in a 3D velocity model, the aim being toimprove earthquake locations for routine surveillance at Mt.Etna volcano, which is currently dependent on a horizontallylayered 1D velocity model.

The analysis of seismic activity is one of the mostimportant components of volcanic surveillance to correctlyevaluate the state of a volcano. Volcanic eruptions are nor-mally preceded by increasing seismicity, and the most reli-able indicators of an impending eruption are the occurrenceof shallow earthquakes and the increase in the tremor ampli-tude (e.g., Chouet, 1996; Patanè et al., 2004; Di Grazia et al.,2009; Alparone et al., 2012). Accurately locating seismicfoci enables identification of an area where an eruptive fis-sure may open and provides both a better understanding ofthe temporal evolution of an eruption and valuable informa-tion on the eruptive dynamic. Well-constrained earthquakelocations, together with other geophysical and geochemicalinformation, are fundamental for scientists and decision mak-ers in order to evaluate the volcano alert level and describe thepotential scenario. Determining the basic parameters of therecorded earthquakes, such as absolute hypocenter locationsand source mechanisms, requires an accurate and realisticvelocity model. In a volcanic area, characterized by irregularsurface topography, and in most tectonically active regionswith complex geological structures, there are marked lateralvariations of near-surface velocities, and consequently largelocation errors could be introduced when using simple 1D lay-ered media for earthquake location. Indeed, though a layeredmodel is always a good first-order approximation, these re-gions can only be represented by fully 3D velocity models.

Although 3D velocity models are available for manytectonic areas and volcanoes, the majority of modern seismicnetworks still depend on 1D velocity models for routineearthquake location.

Greater computing capacity has now made large-scale,grid, and stochastic direct searches feasible for earthquakelocations. Nonlinear location methods such as proposed byLomax et al. (2000) have become increasingly available andimplemented, but these are not yet standard tools in seismicnetworks.

To achieve more reliable earthquake locations at Etnawith a 3D velocity model, a dataset of 328 earthquakes wasused to evaluate the performances of combining nonlinearglobal search algorithms with probabilistic earthquake loca-tions. The PDF of earthquake relocation clearly shows theactive tectonic structures during the 2002–2003 eruption; theseismicity pattern focus is well drawn, sloping southwardsalong the north–south cross section, and eastwards alongthe west–east cross section.

These obtained hypocentral relocations and relateduncertainties are also compared with those computed by a3D linearized inversion (Thurber, 1983; Eberhart-Phillips,1993; Eberhart-Phillips and Reyners, 1997) and with those

using standard traditional and linearized earthquake locationalgorithms such as Hypoellipse.

As expected, the NLLoc, combined with a 3D VP modeland a constant VP=VS ratio, yields more precise and reliablehypocenter locations. The NLLoc also leads to a reliableidentification of poorly constrained hypocenter locations,associated with the larger location uncertainties, and pro-vides stable solutions for these events. Indeed, the complexvelocity structure beneath Etna causes errors in 1D hypocen-ter determination, especially in the focal depth.

Raytracing is a key factor in the performance of an earth-quake-location algorithm. It is well known that incorporatinglaterally heterogeneous Earth models into earthquake loca-tion procedures requires an efficient method for seismicwave travel-time calculations in 3D structures. The resultsof Haslinger (1998), Husen et al. (1999), and Kissling et al.(2001) indicate that the ART-PB raytracing method imple-mented in the local earthquake tomography code SimulPSworks efficiently up to epicentral distances of about 60 km,which are greater than the Etna network aperture. Somenumerical tests with simulated realistic 3D cases show thatthe NLLoc procedure works as well as SimulPS, with rela-tively small location errors for both origin time and sourcehypocenter for tolerable levels of random noise in the arrivaltimes. The location error in depth has comparable accuracyto that in the horizontal plane, at least for the consideredstation geometry and earthquake distribution. This entailshaving depths comparable to the network aperture and notsignificantly smaller than the epicentral distances.

The relocated dataset of 328 earthquakes recorded bythe Etna seismic network by NLLoc, showed similar resultsto those obtained by SimulPS, implying the goodness of thePatanè et al. (2006) 3D velocity model. As suggested byLomax and Michelini (2001), we can associate the minor var-iations in the locations to differences in the model parameter-ization, travel-time calculation, and optimization algorithms ofthe two methods, and may affirm that the two sets of locationsare compatible.

The use of the Oct-Tree algorithm gives rapid and im-proved hypocenter locations, mainly by better constrainingthe depth of the events.

To classify the obtained hypocenter locations withNLLoc, we followed the approach of Husen and Smith(2004), who considered four quality classes (QC) A, B, C,and D, and included a quality factor (Q) for each event, cal-culated on a comprehensive set of normalized parameters ofuncertainty. The importance of the source–receiver geometryand how both the azimuthal GAP and epicentral distance tothe closest station influence the hypocenter location qualityis also confirmed from a number of tests. The nearest stationdistance proves a fundamental constraint on focal depth ac-curacy, especially for shallower earthquakes, even in a densenetwork.

In conclusion, the combination of nonlinear globalsearch algorithms and probabilistic earthquake locationprovides a fast and reliable tool for earthquake location

Improving Seismic Surveillance at Mt. Etna Volcano by Probabilistic Earthquake Location in a 3D Model 2457

(Husen et al., 2003). Among the benefits of using 3D veloc-ity models, probabilistic earthquake location gives a com-plete description of location uncertainties (Lomax et al.,2000). This is crucial, as realtime seismic surveillance mustbe capable of producing estimates and uncertainties on thelocation and size of an earthquake a few seconds after theevent is detected. We may therefore consider the NLLoc pro-gram a reliable tool for routine earthquake locations on Etna.

Data and Resources

All data used in this paper derive from published sourceshttp://www.ct.ingv.it/ufs/analisti/ (last accessed April 2013).The earthquakes dataset was relocated using the NonLinLocalgorithm (NLLoc, http://alomax.net/nlloc; last accessedJanuary 2013).

Acknowledgments

We thank G. Di Grazia for the collaboration, discussion, and criticism;G. Saccorotti for the useful advice on the data analysis; and the SeismicAnalysis Group of INGV-Catania for earthquake data. We also wish to thankM. D’Agostino for his technical support and Stephen Conway and BorisBehncke for correcting and improving the English language of this manu-script. This work has been partially supported by the Spanish projectCGL2011-29499-C02-01. Thoughtful reviews by two anonymous reviewersand Associate Editor Charlotte A. Rowe significantly improved themanuscript.

References

Alparone, S., G. Barberi, O. Cocina, E. Giampiccolo, C. Musumeci, and D.Patanè (2012), Intrusive mechanism of the 2008–2009 Mt. Etna erup-tion: Constraints by tomographic images and stress tensor analysis, J. Vol-canol. Geoth. Res. 229, 50–63, doi: 10.1016/j.jvolgeores.2012.04.001.

Bonaccorso, A., S. D’Amico, M. Mattia, and D. Patanè (2004). Intrusivemechanisms at Mt.Etna forerunning the July–August 2001 eruptionfrom seismic and ground deformation data, Pure Appl. Geophys.161, 1469–1487.

Chiarabba, C., A. Amato, E. Boschi, and F. Barberi (2000). Recent seismic-ity and tomographic modeling of the Mount Etna plumbing system,J. Geophys. Res. 105, 10,923–10,938.

Chiarabba, C., P. De Gori, and D. Patanè (2004). The Mt. Etna plumbingsystem: The contribution of seismic tomography, inMt. Etna Volcano,A. Bonaccorso, S. Calvari, M. Coltelli, C. Del Negro, and S. Falsaperla(Editors), American Geophysical Monograph, Vol. 143, 191–204.

Chouet, B. A. (1996). New methods and future trends in seismological vol-cano monitoring, in Monitoring and Mitigation of Volcanic Hazards,R. Scarpa and R. I. Tilling (Editors), Springer-Verlag, Berlin, 23–88.

De Luca, G., L. Filippi, D. Caccamo, G. Neri, and R. Scarpa (1997). Crustalstructure and seismicity of southern Tyrrhenian basin, Phys. EarthPlanet. Int. 103, 117–133.

Del Pezzo, E., S. De Martino, S. Gresta, M. Martini, G. Milana, D. Patanè,and C. Sabbarese (1993).Velocity and spectral characteristics of thevolcanic tremor at Etna deduced by a small seismometer array,J. Volcanol. Geoth. Res. 56, 369–378.

Di Grazia, G., A. Cannata, P. Montalto, D. Patanè, E. Privitera,L. Zuccarello, and E. Boschi (2009). A multiparameter approach tovolcano monitoring based on 4D analyses of seismo-volcanic andacoustic signals: The 2008 Mt. Etna eruption, Geophys. Res. Lett.36, L18307, doi: 10.1029/2009GL039567.

Eberhart-Phillips, D. (1990). Three-dimensional P and S velocity structure inthe Coalinga Region, California, J. Geophys. Res. 95, 15,343–15,363.

Eberhart-Phillips, D. (1993). Local earthquake tomography: Earthquakesource regions, in Seismic Tomography: Theory and Practice,H. M. Iyer and K. Hirahara (Editors), Chapman & Hall, London,614–643.

Eberhart-Phillips, D., and M. Reyners (1997). Continental subduction andthree-dimensional crustal structure: The northern South Island, NewZealand, J. Geophys. Res. 102, no. B6, 11,843–11,861.

Ferrucci, F., and D. Patanè (1993). Seismic activity accompanying the out-break of the 1991–1993 eruption of Mt. Etna (Italy), J. Volcanol.Geoth. Res. 57, 125–135.

Gambino, S., A. Mostaccio, D. Patanè, L. Scarfì, and A. Ursino (2004).High-precision locations of the microseismicity preceding the 2002–2003 Mt. Etna eruption,Geophys. Res. Lett. 31, L18604. doi: 10.1029/2004GL020499.

Got, J. L., J. Fréchet, and F. W. Klein (1994). Deep fault plane geometryinferred from multiplet relative 488 relocation beneath the south flankof Kilauea, J. Geophys. Res. 99, no. B8, 15,375–15,386.

Haslinger, F. (1998). Velocity structure, seismicity, and seismotectonics ofNorthwestern Greece between the Gulf of Arta and Zakynthos, Ph.D.Thesis, Swiss Federal Institute of Technology of Zürich, Zürich, Swiss.

Hirn, A., A. Nercessian, M. Sapin, F. Ferrucci, and G. Wittlinger (1991).Seismic heterogeneity of Mt Etna: Structure and activity, Geophys.J. Int. 105, 139–153.

Husen, S, and R. B. Smith (2004). Probabilistic earthquake relocation inthree-dimensional velocity models for the Yellow stone National Parkregion, Wyoming, Bull. Seismol. Soc. Am. 94, 880–896.

Husen, S., E. Kissling, N. Deichmann, S. Wiemer, D. Giardini, and M. Baer(2003). Probabilistic earthquake location in complex 3-D velocitymodels: Application to Switzerland, J. Geophys. Res. 108, doi10.1029/2002JB001778.

Husen, S, E. Kissling, E. Flueh, and G. Asch (1999). Accurate hypocentredetermination in the seismogenic zone of the subducting Nazca plate innorthern Chile using a combined on/offshore network, Geophys. J. Int.138, 687–701

Lahr, J. C. (1989). HYPOELLIPSE/VERSION 2.0: A computer programfor determining local earthquake hypocentral parameters, magnitudeand first motion pattern, U.S. Geol. Surv., Open-File Rept. 512, 89-116, 81 pp.

Laigle, M., A. Hirn, M. Sapin, J. C. Lépine, J. Diaz, J. Gallart, and R.Nicolich (2000). Mount Etna dense array local earthquake P and Stomography and implications for volcanic plumbing, J. Geophys.Res. 105, no. B9, 21,633–21,646.

Lee, W. H. K., and J. C. Lahr (1972). HYPO-71 a computer program fordetermining hypocenter, magnitude, and first motion pattern of localearthquakes, U.S. Geol. Surv. Open-File Rept. 72-224, Menlo Park.

Lomax, A., and A. Curtis (2001). Fast, probabilistic earthquake locationin 3D models using Oct-Tree Importance sampling, Geophys. Res.Abstr., 3.

Lomax, A., and A. Michelini (2001). Comparison of NonLinLoc and linearearthquake locations in a 3D model, http://alomax.free.fr/nlloc/alberto_comp/alberto_comp.html.

Lomax, A., A. Michelini, and A. Curtis (2009). Earthquake location, direct,global-search methods, in Encyclopedia of Complexity and SystemScience, Part 5, Springer, New York, 2449–2473, doi: 10.1007/978-0-387-30440-3.

Lomax, A., J. Virieux, P. Volant, and C. Thierry-Berge (2000). Probabilisticearthquake location in 3D and layered models, in Advances in SeismicEvent Location, C. H. Thurber and N. Rabinowitz (Editors), KluwerAcademic Publishers, Dordrecht/Boston/London, 101–134.

Kissling, E., S. Husen, and F. Haslinger (2001). Model parametrization inseismic tomography: A choice of consequence for the solution quality,Phys. Earth Planet. Int. 123, 89–101.

Klein, F. W. (1978). Hypocenter location program HYPOINVERSE, U.S.Geol. Surv. Open-File Rept. 78-694, 113 pp.

Maurer, H., and U. Kradolfer (1996). Hypocentral parameters and velocityestimation in the western Swiss alps by simultaneous inversion ofP- and S-wave data, Bull. Seismol. Soc. Am. 86, 32–42.

2458 A. Mostaccio, T. Tuvè, D. Patanè, G. Barberi, and L. Zuccarello

Moser, T. J., T. van Eck, and G. Nolet (1992). Hypocenter determination instrongly heterogeneous earth models using the shortest path method,J. Geophys. Res. 97, 6563–6572.

Musumeci, C., O. Cocina, P. De Gori, and D. Patanè (2004). Seismologicalevidence of stress induced by dyke injection during the 2011 Mt. Etnaeruption, Geophys. Res. Letters 31, doi: 10.1029/2003GL019367.

Patanè, D., and E. Giampiccolo (2004). Faulting processes and earthquakesource parameters at Mt.Etna: State of the art and perspectives, in Mt.Etna: Volcano Laboratory, A. Bonaccorso, S. Calvari, M. Coltelli,C. Del Negro, and S. Falsaperla (Editors), American GeophysicalMonograph 143, 167–189.

Patanè, D., G. Barberi, O. Cocina, P. De Gori, and C. Chiarabba (2006).Time-resolved seismic tomography detects magma intrusions at MountEtna, Science 313, 821–823.

Patanè, D., C. Chiarabba, O. Cocina, P. De Gori, M. Moretti, and E. Boschi(2002). Tomographic images and 3D earthquakes locations of theseismic swarm preceding the 2001 Mt. Etna eruption: Evidence fora dyke intrusion, Geophys. Res. Lett. 29, 10.

Patanè, D., O. Cocina, S. Falsaperla, E. Privitera, and S. Spampinato (2004).Mt. Etna volcano: A seismological framework, in Mt. Etna VolcanoLaboratory, A. Bonaccorso, S. Calvari, M. Coltelli, C. Del Negro,and S. Falsaperla (Editors), American Geophysical Monograph 143,147–165.

Patanè, D., P. De Gori, C. Chiarabba, and A. Bonaccorso (2003). Magmaascent and the pressurization of Mt. Etna’s volcanic system, Science299, 2061–2063, doi: 10.1126/science.1080653.

Patanè, D., F. Ferrucci, E. Giampiccolo, and L. Scaramuzzino (1997).Source scaling of microearthquakes at Mt. Etna volcano and in theCalabrian Arc (southern Italy), Geophys. Res. Lett. 24, 1879–1882.

Patanè, D., M. Mattia, and M. Aloisi (2005). Shallow intrusive processesduring 2002–2004 and current volcanic activity on Mt Etna, Geophys.Res. Lett. 32, L06302, doi: 10.1029/2004GL021773.

Patanè, D., E. Privitera, S. Gresta, A. Akinci, S. Alparone, G. Barberi,L. Chiaraluce, O. Cocina, S. D’Amico, P. De Gori, G. Di Grazia,S. Falsaperla, F. Ferrari, S. Gambino, E. Giampiccolo, H. Langer,V. Maiolino, M. Moretti, A. Mostaccio, C. Musumeci, D. Piccinini,D. Reitano, L. Scarfi, S. Spampinato, A. Ursino, and L. Zuccarello(2003) Seismological constraints for the dike emplacement of July–August 2001 lateral eruption at Mt. Etna volcano, Italy, Ann. Geophys.46, 599–608.

Podvin, P., and I. Lecomte (1991). Finite-difference computation of travel-times in very contrasted velocity models: A massively parallelapproach and its associated tools, Geophys. J. Int. 105, 271–284.

Poupinet, G., W. L. Ellsworth, and J. Frechet (1984). Monitoring velocityvariations in the crust using earthquake doublets: An application to theCalaveras fault, California, J. Geophys. Res. 89, 5719–5731.

Presti, D., B. Orecchio, G. Falcone, and G. Neri (2008). Linear versusnonlinear earthquake location and seismogenic fault detection in thesouthern Tyrrhenian Sea, Italy, Geophys. J. Int. 172, 607–618.

Richards-Dinger, K. B., and P. M. Shearer (2000). Earthquake locations insouthern California obtained using source-specific station terms, J.Geophys. Res. 105, 10,939–10,960.

Saccorotti, G., F. Bianco, M. Castellano, and E. Del Pezzo (2001). TheJuly–August 2000 seismic swarms at Campi Flegrei volcanic complex,Italy, Geophys. Res. Lett. 28, no. 13, 2525–2528.

Sharp, A. D. L., P. M. Davis, F. Gray, and G. C. P. King (1980). Preliminaryresults from a seismic array study of Mount Etna, in United KingdomResearch on Mkount Etna, 1977–1979, The Royal Society, London,15–18.

Tarantola, A. (1987). Inverse Problem Theory: Methods For Data Fittingand Model Parameter Estimation, Elsevier, Amsterdam, 613 pp.

Tarantola, A., and B. Valette (1982). Inverse problems = Quest forinformation, J. Geophys. Res. 50, 159–170.

Thurber, C. H. (1983). Earthquake locations and three-dimensional crustalstructure in the CoyoteLake area, central California, J. Geophys. Res.88, 8226–8236.

Wittlinger, G., G. Herquel, and T. Nakache (1993). Earthquake location instrongly heterogeneous media, Geophys. J. Int. 115, 759–777.

Istituto Nazionale di Geofisica e VulcanologiaOsservatorio Etneo, Sezione di CataniaCatania, Italymostaccio@ct.ingv.ittuve@ct.ingv.itpatane@ct.ingv.itbarberi@ct.ingv.itzuccarello@ct.ingv.it

Manuscript received 20 July 2011

Improving Seismic Surveillance at Mt. Etna Volcano by Probabilistic Earthquake Location in a 3D Model 2459