New locations of volcano-tectonic earthquakes under Popocatépetl Volcano applying a Genetic Search Algorithm

Pia Berger^1*, F. Alejandro Nava Pichardo², Carlos Valdés González³ and Alicia Martínez Bringas⁴

¹ Instituto de Geofísica, Departamento de Sismología, Universidad Nacional Autónoma de México, Ciudad Universitaria, Delegación Coyoacán, 04510. México D.F., México. ^*Corresponding author: piaahoi@gmail.com

² Centro de Investigacion Científica y Educación Superior de Ensenada 22860, BC, Mexico

³ Instituto de Geofísica, Departamento de Sismología, Universidad Nacional Autónoma de México, Ciudad Universitaria, Delegación Coyoacán, 04510, México D.F., México

Received: August 27, 2010
accepted: December 8, 2010
published on line: June 30, 2011

Resumen

En los inicios de los años noventa, después de ~70 años de quietud, la actividad sísmica bajo el Volcán Popocatépetl en el centro de México se reactivo, comenzando un período de alta actividad volcánica. Esta actividad continúa hasta la fecha y ha provocado emisiones importantes de gases volcánicos, así como formación de domos y erupciones moderadas. Hemos aplicado un algoritmo genético de búsqueda, para ajustar las diferencias de los tiempos de arribo de cada sismo individualmente, para obtener relocalizaciones de alta precisión (~200 m) para 405 sismos volcano-tectónicos (VT) de un grupo inicial de 968 eventos registrados por una red local de 1995 a 2006. Comparamos los resultados de la relocalización con los obtenidos aplicando el método de dobles diferencias. El objetivo de este trabajo es el de caracterizar la sismicidad, determinar parametros hipocentrales y explorar fuentes sísmicas.

Los resultados del presente análisis mejoran la determinación de la distribución de la actividad sísmica, permitiendo observar características que ocultaba la dispersión en las localizaciones. La agrupación difusa de eventos, asociada con una falla previamente identificada en el SE, aparece ahora como un sistema de fallas que consiste de al menos una falla con rumbo NW-SE, atravesada por otra falla con rumbo NE-SW.

Se encontró que la ocurrencia de sismos ha variado localmente con el tiempo. Algunas alineaciones de eventos coinciden con fallas previamente identificadas. Otras agrupaciones lineales de eventos sugieren fallas ocultas reactivadas por actividad volcánica, como intrusiones de diques, o degasificación de un cuerpo magmático bajo el Volcán Popocatépetl.

Abstract

In the early 1990's, after ~70 years of quiescence, seismic activity was renewed beneath Popocatépetl Volcano, central Mexico, and a period of high volcanic activity began. This activity continues today, and has featured high emissions of volcanic gases and fumaroles, dome filling processes, and moderate eruptions. The three largest explosive eruptions occurred in 1997, 2001, and 2003.

We applied a genetic search algorithm to fit arrival-time differences on an individual basis, and we obtain high precision (~200 m) relocations for 405 volcano-tectonic (VT) events from 968 events recorded by a local network from 1995 to 2006. We compare these results to relocations obtained by applying a double-difference algorithm. The objective is to characterize the seismicity, determine hypocentral parameters, and explore seismic sources.

The results shed light on the distribution of seismic activity, revealing features previously hidden by location scatter. A diffuse cluster associated with a previously-identified SE-trending fault now appears to be associated with a fault system consisting of at least one NW-SE trending fault crossed by a NE-SW trending fault.

Event occurrence was found to be locally time-dependent. Some aligned events coincide with previously-proposed faults. Other linear clusters may indicate hidden faults being activated by volcanic activity, such as dike intrusion or degassing of the magma body below Popocatépetl Volcano.

Key words: Popocatépetl Volcano, relocalization, arrival-time difference genetic algorithm, double-difference.

Introduction

Accurate event location is important and necessary in order to understand the behavior and characteristics of active volcanoes like Popocatépetl. Whereas on a global scale earthquake location allows imaging of the major geodynamic features of the Earth, its accuracy remains too low to allow description of many seismo-tectonic features at local scales and in highly heterogeneous media such as volcanoes. At local scales, an unfavorable geometry of either the earthquake spatial distribution and/ or the seismic network stations generates tradeoffis between model parameters and induces uncertainties in earthquake locations. Double-difference relocation methods developed in the last decades (e.g., Jordan and Sverdrup, 1981; Poupinet et al., 1984; Ito, 1985; Got and Frechet, 1994; Slunga et al., 1995; Shearer, 1997; Rubin et al., 1998; Waldhauser and Ellsworth, 2000) have shown that knowledge about a seismogenic region may be significantly changed after accurate relocation of earthquakes (see Wolfe, 2002) for a discussion of double-difference algorithms). The limitations of these methods are not always taken into account, however, so that too large a confidence is sometimes placed on their results.

In this study, we apply both methods: (1) a common double-difference method and (2) a newly developed relocation method, which is based on a genetic algorithm, a global optimization search method.

Geological setting

Popocatépetl Volcano (5,452 m) is one of the most active stratovolcanoes in Mexico (Fig. 1), located in the central part of the Mexican Volcanic Belt, a volcanic arc related to the subduction of the oceanic Cocos and Rivera plates beneath the North American plate. Popocatépetl poses a major geological hazard for Mexico, as a sudden eruption could threaten highly populated areas, including Mexico City (60 km northwest of the crater) and Puebla (40 km east of the crater); in all, more than 10 million people are living within 70 km of the volcano (Siebe et al., 1996; De la Cruz-Reyna and Siebe, 1997; Macías and Siebe, 2005).

The current phase of seismic activity in Popocatépetl began in 1990 (De la Cruz-Reyna et al., 2008), with a large increase occurring in 1993 and explosive reactivation starting in December 1994, after nearly 70 years of dormancy (Siebe et al., 1996; De la Cruz-Reyna and Siebe, 1997). Since 1996, several dome emplacement-destruction processes have taken place (Arciniega Ceballos et al., 2000; Wright et al., 2002). Due to its high potential risk, many geological and geophysical studies have been carried out on the volcano (e.g., Valdés et al., 1995; Campillo et al., 1996; Shapiro et al., 1997; Espíndola et al., 2004; Espinasa-Perena and Martín Del Pozzo, 2006).

López Ramos (1983) proposes the existence of a basement of limestones and granodiorites beneath Popocatépetl volcano, extrapolating geologic data from the surrounding areas. Valdés González and Comité (1994) propose that regional basement structure occurs at a depth of 9 km b.s.l. (below sea level). Geological studies from e.g. Fries (1965), Meritano-Arenas et al. (1998) indicate that the volcano rests upon a stratum of limestone (-1.5 km b.s.l.) and metamorphic rock (~ sea level): limestones and metamorphic rocks are exposed as isolated outcrops related to horst and graben fault structures more than 20 km south of Popocatépetl's summit. Geophysical, geological and geochemical studies on Popocatépetl have led to controversial results, so that until now the crustal seismic structure beneath Popocatépetl is not well understood. Straub and Martín-Del Pozzo (2001) suggest that an andesitic magma rises from the moho and mixes with another dacitic melt at depths of ~4 to 13 km below the crater. Based on gravimetric measurements, Espíndola et al. (2004) modeled density difference below the volcano, and interpreted a 25 km³ negative density contrast at a depth of 7 km b.s.l. as the magma chamber. Studies on fluid inclusions by Atlas et al. (2006) and Roberge et al. (2007) negate the existence of a large magma chamber at depths shallower than 4.5 km b.s.l. and propose a dike and sill system instead.

Popocatépetl is subjected to the regional stress state, which is SE-NW for the minimal stress s₃, SW-NE for the medium stress s₂, and vertical for the principal stress s₁ (Ego and Ansan, 2002; Arámbula-Mendoza et al., 2010). The study by Arámbula-Mendoza et al. (2010), who calculated the focal mechanisms of volcanic events, showed that magma movement and volcanic activity can locally influence these stress axes.

The volcano-tectonic activity of Popocatépetl began in 1990 (De la Cruz-Reyna et al., 2008), while its explosive reactivation started in December 1994. In 1995, seven seismic stations (PPM, PPC PPP, PPX, PPQ, PPS, PPB) from the National Seismological Service (SSN) were operating at distances between 1.5 and 10.5 km from the volcano (Fig. 1b, c, d). Through the efforts of the SSN and the National Center for Disaster Prevention, Mexico (CENAPRED), the seismograph network eventually grew to seven stations. The stations are sited on the volcano's fanks, at altitudes from 2500 to 4450 m. Five of the stations (PPM, PPC, PPX, PPQ, and PPP) are three-component seismometers, and two of them (PPS and PPB) are vertical-component only. All stations have 1 Hz natural frequency sensors. In the following years, the network expanded to include three triaxial stations of 1 Hz natural frequency: PPN (September 1995), PPT (May 1996), and PPJ (November 1997). In March 1998, the network was augmented when the triaxial broadband station PPA was installed and the short period station PPP was upgraded to broadband. In July 1998, the array was completed with an upgrade of station PPX to broad-band capabilities. The Popocatépetl network, shown in Figure 1, operated continuously with 11 stations, until station PPB failed in June 1997, station PPN in November 1997, and station PPA in August 1998. Since 2000, the network has been operating with nine stations.

For the period from November 1995 to December 2006, about 1,800 identifiable volcano-tectonic (VT) earthquakes were recorded, with coda magnitudes ranging from M_C 1.17 to 3.80. VT earthquakes are indistinguishable from common double-couple tectonic earthquakes; they have been interpreted as the brittle response of volcanic rocks fracturing due to fluid pressure superimposed on the regional stress field (e.g., Chouet, 1996), and act as indicators of stress concentrations within the area surrounding magma reservoirs and conduits. Because they originate as abrupt shear motion along faults, the seismograms of VT-events typically show an impulsive onset followed by coda waves, whose spectral peaks are broadly distributed between 3 and 18 Hz.

Waveform data from the Popocatépetl network were telemetered to CENAPRED, and a central GPS-clock was used to control time for the detection system. First-arrival times on the seismograms from VT earthquakes were manually identified, and earthquake locations were determined using the Hypocenter earthquake location algorithm (Lienert and Havskov, 1995; Havskov, 2003); and the 1-D layered P-wave velocity model from Valdes-González and Comité (1994) (Line 3 blue, dash-dotted, Fig. 4). The description of this and other velocity models is given in section 4. Initially, several velocity models were tested for the location of the largest VT-earthquakes, which had clear arrivals and were detected by most of the stations. The model which produced the smallest arrival-time residuals with Hypocenter was the one proposed by Valdes-González and Comité (1994), which resulted in an average root mean square (RMS) residual of 0.12 sec for about 800 located VT events. The maximum adjusted location error was calculated to be 0.25 km, as an average of the horizontal and vertical errors reported by Hypocenter.

Velocity models

For our study, preliminary locations (see § 5.1) were obtained using the single event program Hypocenter (Lienert and Havskov, 1995) and the Valdés-González and Comité (1994) 1-D velocity model.

The model from Valdés-González and Comité (1994) and two other velocity models from Cruz-Atienza et al. (2001) and De Barros et al. (2008), all proposed for Popocatépetl volcano, were used to relocate our data set, using the genetic algorithm program (§ 5.2) and the double-difference algorithm (§ 5.3).

The shear wave velocity model from Cruz-Atienza et al. (2001) was obtained by inversion of receiver functions using four teleseismic events from South America at station PPIG (PPM), located 4 km north of the Popocatépetl crater; this model includes a low velocity zone between 6 and 10 km depth (3-7 km b.s.l.) (Fig. 4).

The model from De Barros et al. (2008) was obtained from analysis of Rayleigh waves and a recalculation of the phase velocities corresponding to several previous models. The shear-wave velocities in this model are similar to the MVB crust, and are close to those of the Valdés et al. (1986) model, which is based on a seismic refraction study in Oaxaca (Line 5, Fig. 4). P-wave velocity models were estimated from S-wave velocities models and vice versa, using a Poisson ration of Vp/Vs = 1.75, a reasonable value for this region.

Location methods

First locations

A first study locating Popocatépetl earthquakes was reported by Valdés et al. (1995). They located 55 earthquakes of the types A, B and AB (VT, low frequency and hybrid events, respectively) using digital seismograms recorded at seven seismic stations. These events, recorded during December 1994 and March 1995, were mainly concentrated below the crater and only three events were located in the southeast region.

In the unpublished work by Valdés-González and Comité (1994), with parts shown in Arámbula-Mendoza et al. (2010), nearly 2,000 VT events recorded from 1994 to 2006 below Popocatépetl were located using the Hypocenter program, within a range of about 20 km. They found the events formed two connected clusters, one below the summit (Cluster A) and the other southeast of it (Cluster B) (Fig. 2). Earthquake sources were distributed at depths between -3 km and 7 km b.s.l. under the crater region, and down to 5 km b.s.l. in the SE-zone. These locations are used in this study for relocation with the genetic search algorithm, called DisLoca, and described in section 5.2. For better comparison, Fig. 3 shows the same 405 events located with Hypocenter, later relocated with DisLoca using the velocity model from De-Barros et al. (2008) (Line 1 dark green, dashed, of Fig. 4 (see also section 4)).

Lermo-Samaniego et al. (2006) analysed volcano-tectonic earthquakes recorded by a minimum of 5 stations from the Popocatépetl Seismic Network during the period 1994-1999.

Their hypocentral locations were calculated using the SEISAN program (Havskov, 2003), with the one-dimensional (1-D) velocity model proposed by Valdés-González and Comité (1994) (see Fig. 4) and a Poisson ratio of Vp/Vs=1.76; coda magnitudes were estimated from parameters, as proposed by Chavacán et al. (2004). Lermo-Samaniego et al. (2006) verify the existence of the two main clusters mentioned above. Their events have coda magnitudes M_c< 3.2 and depths below 12 km, with hypocentral location errors estimated as < 1 km. Based on the asumption that VT-earthquakes under the crater are caused by direct volcanic activity (i.e., magma motion), as proposed by Minakami (1974) and Karpin and Thurber (1987); Lermo-Samaniego et al. (2006) proposed that seismicity under the crater could be related to fracture processes in the volcano fanks, as a result of magmatic upwelling.

According to the studies by De la Cruz-Reyna and Siebe (1997); De Cserna et al. (1988); Meritano-Arenas et al. (1998); Lermo Samaniego et al. (2006); Arámbula-Mendoza et al. (2010), events in the SE zone are probably tectonic events occurring on a NW-SE striking fault, which according to Arámbula-Mendoza et al. (2010) appears to be activated by the movement of rising magma in the Popocatépetl edifice. Lermo-Samaniego et al. (2006) associate events in the SE zone to normal faulting in a NW-SE direction that may have been caused by tensional stresses that activated or reactivated a fault system caused by the persistent activity of Popocatépetl during the last 10 years.

The first hypocentral locations obtained by using SEISAN and Hypocenter do not show any relationship to the faults proposed by De Cserna et al. (1988) and Meritano-Arenas et al. (1998) (except for Tlaltzompa Fissure (TF) (see Fig. 2, where vents are mapped by Espinasa-Perena and Martín-Del-Pozzo 2006) nor do they allow the interpretation of a magma chamber. We will explore the possibility that more sophisticated location algorithms may shed some light on these expected features.

In order to relocate the events of the Popocatépetl volcano-tectonic catalog data recorded between 1995 and 2006, and preliminarily located with the location program Hypocenter, we applied a recently-developed absolute location algorithm, called DisLoca. DisLoca (Nava, 2010) estimates hypocentral coordinates by finding the hypocenter which results in the minimum L1 (mean absolute) residual between calculated and observed arrival time differences between all recorded P_g and S_g phases, by a genetic search in space, followed by the determination of the origin time. Location determination begins by considering hypocenters at a given grid of strategically-distributed points within and around the station array; it can also include a previously-determined hypocenter for relocation. A given number of children are generated by varying the parent (initial) locations with normally distributed pseudo-random numbers, with X, Y, and Z standard deviations proportional to the corresponding standard deviations for the parent hypocenters and a given mutation probability; other children are constructed as averages of pairs of parents. Next, new parents are chosen from among the entire new population, including parents and children, and the process is repeated until residual and spatial convergence criteria are met for an acceptable location, or until a given number of generations is reached, without meeting the residual creteria, in which case the location is considered unacceptable.

Once an acceptable location is found, the origin time is estimated from the travel times to the stations. The resulting hypocentral parameters are reported, together with arrival time residuals. The average difference time residual is reported as a measure of the fit error, and the space around the preferred hypocenter is explored to determine the X, Y, and Z ranges for which the location error increases by less than a given criterion, which gives a good estimate of the location uncertainty.

DisLoca uses a layered model, and can feature a truncated cone having the first layer's velocity, to roughly model the volcanic edifice. The program accepts station corrections, which allows the possibility of using some selected hypocenter as a master event for the location of other events close to it. The use of arrival time differences, besides reducing the dimension of the search space, makes DisLoca locations more stable in cases of less-than-optimal azimuthal coverage versus straight arrival time fitting. A principal difference between DisLoca and many other common location programs such as Hypo71, Hypocenter or HypoDD, is that DisLoca accounts for station elevation and includes an approximation to the shape and altitude of the volcanic edifice. This is why DisLoca may locate events within the volcanic cone, events that are considered "airquakes" by other programs such as HypoDD or Hypocenter.

Relocation method - HypoDD

In order to compare the results obtained from DisLoca to other relocation results, we applied the double-difference relocation algorithm HypoDD (Waldhauser and Ellsworth, 2000; Waldhauser, 2001). HypoDD is a relative location program based on the double-difference algorithm which uses the double-difference equation:

Where dr^ij_k is the residual between observed and calculated differential travel times for the two events i and j, recorded at station k. HypoDD minimizes residual double-differences for pairs of earthquakes by adjusting the vector differences between their hypocenters, and determines the interevent distances between correlated events, without need for station corrections (Waldhauser and Ellsworth, 2000; Waldhauser, 2001).

Data

The original locations of VT events compraise around 1,800 events and form two main clusters of about 5 km diameter, one below the summit (Cluster A) and a second SE of the summit zone (Cluster B). Out of this data set we selected 968 events, including 4,321 P-wave phases and 3,009 S-wave phases, for the relocation process (Fig. 2). This selection included events recorded by at least four stations with a minimum weight of 0.1, where 1 is the best and 0 is the worst. The RMS residual error was not chosen as a selection parameter because it results, among other factors, from uncertainties in the 1-D velocity model, which cannot correctly represent all parts of the heterogeneous volcanic body.

Results

Results from DisLoca

The hypocenter relocations obtained by the time-difference genetic search approach, DisLoca, provide a clear new picture of VT events beneath Popocatépetl (compare Fig. 2, Fig. 3 and Figs. 5a, b, c). Out of 968 events, we relocated 331 to 405 volcano-tectonic events, depending on the velocity model applied. Events show at least four P-phases and weights of 0.1 - 1, (see Tab. 1 for details) and were recorded between Nov. 1995 and Dec. 2006. Earthquakes with a final location error higher than 0.1 s were discarded by the location algorithm.

Although DisLoca minimizes the L1 residuals, relocation yielded a substantial reduction in the RMS of the relocated events. The initial weighted RMS values for catalog data were 0.128 s. After relocation, these errors were reduced to 0.045 s (Tab. 1) for the same group of hypocenters.

The number of relocated events differs for the different velocity models used for the inversion (see Figs. 5a, b, c); we relocated 331 events for the velocity model from Cruz-Atienza et al. (2001), 386 for the velocity model from Valdés-Gonzáléz and Comité (1994) and 405 events for the velocity model from De-Barros et al. (2008). We will base our interpretation on the results using the velocity model from De-Barros et al. (2008) (Fig. 5c), as its application resulted in the highest number of relocated events.

Relocated events below the crater area still accumulate in Cluster A, with depths down to 4  km b.s.l. and 7 to 11.4 km b.s.l., with their main occurrence limited to -1.5 to +1 km b.s.l. Relocated events in the southeast zone (about 5  km southeast of the crater), however, are not relocated in a clear Cluster B (see Fig. 1), unlike the originally-located events and as cited by Arámbula-Mendoza et al. (2010) and Lermo-Samaniego et al. (2006). Events in the southeast zone are relocated in perpendicularly-aligned clusters, with NE-SW and NW-SE strikes, marked in Figs. 5c and 253 6a (see § 8). The maximum depth of events in the southeast zone is 4.3 km, with main event occurrence at -1.2 to +0.9 km b.s.l.

We observe a strong correlation in local and temporal event occurrence. It is notable that events recorded in the time period from 1995 to the eruption on 30th June 1997 are relocated mainly in the southern part of Cluster A and to the southeast. A small number of events of this Period 1 are relocated in the northern part of Cluster A, but not in its center (see Figure 6a). We observe as well a strong correlation between event occurrence and volcanic episodes, which are defined after Arámbula-Mendoza (2007). Table 2 shows the 13 volcanic episodes differentiated in time and their mode of activity. In Figure 6b-f we observe that Cluster A (shown in Figure 6a) consists of relocated events from nearly all thirteen episodes, except from Episode 1. Episode 1 is characterized by gas and ash emission and its events are relocated in a NW-SE aligned cluster of about 10 km length located north and east from the central part of the volcano. Episode 3 is a phase of moderate ash emission, as well as spasmodic tremors. The events of Episode 3 are grouped into two clusters, one located 1.5 km south of the crater and a second aligned in an E-W direction begining at 4 km south of the crater and up to 9 km southeast of it.

Events from Episodes 7 and 12, post-eruptive and relaxation phases, are located 4 to 8 km south of the crater region. They are located at depths from -3.5 to 10.5 km b.s.l., and therefore more than 6 km deeper than events located on the NW-SE and NE-SW striking faults in the SE region. Events located south of the central volcanic region are from episodes of dome construction and eruptions: (4), 5, 6, (9), 10, 11, (13), where episodes marked in parentheses contain only one event, others two or three events. Figure 7 shows the events from Episodes 3, 10 and 12 in 3-D cross-section, in order to display inclined event alignments not visible from the map views shown in Fig. 6. The events are inclined, suggesting normal or inverse faults (Figs. 6 and 7).

Differing station corrections for different ray paths can indicate the presence of low or high velocity regions, such as magmatic chambers or conduits. DisLoca permits the estimation of station corrections from the mean (re)location error (in seconds) and its standard deviation for P- and S-wave rays. Here, we show station correction examples for events located with the velocity model from De-Barros et al. (2008).

We estimated the station correction for each station for the entire dataset (Table 3) and then searched for differences in station corrections as functions of the earthquake's occurrence in time and location (Tables 1, 3 and 4). We observe a high variation in values of station corrections in space and time, which does not allow for a systematic location including station correction.

Results from HypoDD

With HypoDD, depending on the velocity model applied, we relocated between 430 and 561 out of 739 events, which were preselected by a program that searches for possible pairs of earthquakes (Fig. 8a). Events were automatically discarded by the relocation algorithm if they were badly connected or if their depths decreased to less than the subsurface (airquakes, 3 km above sea level in the case of Popocatépetl).

There is a strong dependence on the velocity model for HypoDD, which is clearly observed when comparing the results obtained by applying the velocity model of Cruz-Atienza et al. (2001), which includes a low velocity zone, with the results obtained by applying the velocity models of Valdés-González and Comité (1994) and DeBarros et al. (2008) (Fig. 8).

In the map views of the relocation results, relocation concentrates the events in two clusters (Cluster A and B) for all three velocity models (map views, Fig. 8). We observe no visible alignment in event distribution correlating to the fault zones, as indicated by De Cserna et al. (1988) and Meritano-Arenas et al. (1998). The area of event distribution was slightly reduced, compared to the initially located events in Figure 2. Events are located at depths of +3 and -5 km b.s.l. The value of the maximum depth was diminished by only 1 km (5 km b.s.l. located to 4 km b.s.l. relocated).

The velocity model of Cruz-Atienza et al. (2001) shows a horizontal event accumulation at -1.5 to 0 km depth b.s.l. including events from below and around the crater region and the southeastern region. Below 312 the crater we image an inclined "conduct" of events at depths of -2 to 4.5 km b.s.l. and a horizontal diameter of about 2 km.

The relocated Cluster A, using HypoDD and the models of Valdés-González and Comité (1994) and De-Barros et al. (2008) looks pretty similar to the Cluster A resulting from the relocation using DisLoca.

Discussion and Conclusions

This study applied a genetic algorithm location method to relocate 331, 386 and 405 volcano-tectonic events, depending on the velocity model applied, recorded between Nov. 1995 and Dec. 2006 beneath Popocatépetl Volcano, and initially located using the program Hypocenter. These results are supplemented by applying the double-difference algorithm HypoDD (Waldhauser and Ellsworth, 2000) on the initially event locations.

Different relocation results between HypoDD and DisLoca may arise from differences in the algorithms: (1) DisLoca includes topography and station elevation and may accept events that are discarded as airquakes in HypoDD, (2) the genetic algorithm searches for the global optimum while the least squares regression applied in the double-difference algorithm may fall easily in a local optimum. While the RMS in the relocation with DisLoca gets diminished by nearly one third, the RMS in the HypoDD relocation increased from 0.18 s to ~ 0.3 s. The most probable reason for this increase is the lack of waveforms, which may enhance relocation results significantly, as clearly emphasized by Waldhauser and Ellsworth (2000) and Waldhauser (2001).

Comparing the relocation results obtained by DisLoca (Figs. 5 and 6) with those obtained using HypoDD (Fig. 8) clearly shows that the events relocated using the genetic algorithm are aligned and may image predicted fault zones. Aligned epicenters are important because they may indicate strain localization within the volcanic edice, which in turn may indicate increased local risk of collapse. Having the alignments coincide with known faults makes them doubly important. Relocated events using HypoDD, however, result in two main spherical or plate-shaped clusters without any event-alignments. Nearly-spherical shapes were derived using similar velocity models of Valdés-González and Comité (1994) and De-Barros et al. (2008), while plate-shaped clusters were derived using the velocity model of Cruz-Atienza et al. (2001), which includes a low velocity zone from 3 to 6 km b.s.l. The high dependence on the velocity model using HypoDD is another indication of the instability of the relocation results using the double-difference algorithm with our available data. Hence, we will discuss the results obtained by the genetic algorithm location method DisLoca in more detail.

The number of acceptable relocated events using DisLoca differs depending on the velocity model used during the inversion (see Fig. 4 and Figs. 5a, b, c: 331 events using the Cruz-Atienza et al. (2001) model, 386 events using the Valdés-González and Comité (1994) model and 405 events using the De-Barros et al. (2008) model). The RMS error was considerably improved from 0.128 s to 0.045 s. These relocations provide a picture of the spatial and temporal patterns of this seismicity.

Events from the first locations with Hypocenter were relocated slightly differently in the relocation process using DisLoca. This is clearly visible in the two main clusters determined using the first locations, Cluster A below the crater and Cluster B in the southeast. Cluster A is still distinguishable in the relocation, but the former clear accumulation of Cluster B is relocated in two aligned linear clusters extending SE-NW, NE-S W, revealing probable faults in direction of the axes of the regional stress field, as previously determined by Meritano-Arenas et al. (1998) and De-Cserna et al. (1988).

We also observe a clear correlation between event location and volcanic episodes (Tab. 2 and Figs. 6a-f). Events produced during dome construction or eruptive phases occur at different depths than events caused during relaxation. Typical depths for events produced during dome construction are: -4.5 to +2 km b.s.l., while events during explosion phases are relocated to deeper areas, at -3 to +10 km b.s.l.

The relocated aligned clusters, which we interpret as faults, coincide with a special type of volcanic activity, limited in time. In Figures 7a-c, we show the Episodes 3, 10 and 13 inclined aligned structures, which we interpret as faults. These alignments are only seen from the SE or SW direction, in accordance with the regional stress field.

Faults were confirmed from the studies of Meritano-Arenas et al. (1998) and De-Cserna et al. (1988) are (Figure 6): Tetela Fault and SE-NW fault are both imaged by alignments of events from volcanic Episode 11 (post eruptive and construction of small domes) and Episode 12 (small explosions and relaxation). Tlamacas Fault can be confirmed by aligned clusters from Episodes 7 (post eruptive and relaxation) and 8 (relaxation). Events of Episode 8 image Tlaltzompa Fissure and a W-E striking fault proposed by De-Cserna et al. (1988); Espinasa-Perena and Martín-Del-Pozzo (2006) mapped vents SW and NE of the crater, where Tlaltzampa Fissure and Atexca Fault are proposed by Meritano-Arenas et al. (1998) and De-Cserna et al. (1988); this is why Atexca Fault is better dened as a ssure than a fault. Relocated events from Episodes 2 (dome construction), 3 (ash and explosions), 4 (dome construction, effusive) and 5 (explosions and energy accumulation) are found to form aligned clusters about 5-10 km south of the crater; these assumed faults strike E-W and N-S.

Zones where seismicity is lacking are of importance, as they may indicate areas of hot or molten material, such as small magma reservoirs. At depths shallower than 4 km b.s.l., we observe a strong decrease in seismic events; in agreement with studies from e.g., Espíndola et al. (2003); Atlas et al. (2006) and Roberge et al. (2007), a larger magma body or presence of magma could be interpreted below that area. We observe as well areas of no seismicity, e.g., below the crater at -2 km b.s.l., that may be explained by the fact that several significant changes occur at shallow depths (1-3 km) in the highly fractured upper part of Popocatépetl (Novelo Casanova et al., 2007).

Observations on the agreement of event location and volcanic behavior indicate that the behavior of magmatic material is predictable, after observing several known patterns. Faults often become activated by magma movement, as proposed for the southeast fault zone by Lermo-Samaniego et al. (2006) and Arámbula-Mendoza et al. (2010). But faults may be activated as well in times of volcanic relaxation, as is the case for Episode 8, for fault zones west of the volcano (Fig. 6c). Although we observe some correlation between the location of the VT seismic events and previously proposed faults, we do not find large concentrations of earthquakes associated with the SE faults, which could indicate a weakening or future collapse of the volcano edifice, and that thus would represent a volcanic risk in the near future. Based on the thesis of Ramirez-Olvera (2003) and the work of Novelo-Casanova et al. (2007), we conclude that a shallow (< 3 km) earthquake of magnitude 4.5 to 5 is necessary, to create a risk of collapse or to create a fissure. The earthquake that produced the fank collapse at Mount Saint Helens, for example, had a magnitude of 5.1 (e.g., Mullineaux and Crandell, 1981). The possibility of an earthquake of magnitude 4.5 to 5, while unlikely, cannot be excluded. In the past, the volcano has collapsed to the south and south-east (e.g., Siebe and Macías, 1997). It is possible that this is the mechanism that may occur in the future.

Acknowledgements

We want to thank the volcano monitoring group of CENAPRED (Nacional Center for Disaster Prevention of Mexico), for keeping the seismic monitoring network at Popocatépetl Volcano. Part of this project was funded by the CONACYT (National Council of Science and Technology of Mexico). Thanks also go to Christine Gans and two anonymous reviewers who provided helpful criticisms of the manuscript.

Bibliography

Arámbula-Mendoza R., 2007, Estado de esfuerzos en el volcán Popocatépetl obtenidos con mecanismos focales en el periodo de actividad de 1996 a 2003 . Tesis de maestría Instituto de Geofísica, UNAM Mexico.         [ Links ]

Arámbula-Mendoza R., Valdés-González C., Martínez-Bringas A., 2010, Temporal and spatial variation of the stress state of Popocatépetl volcano, Mexico. J. Volcan. Geoth. Res., 196 , 156-168.         [ Links ]

Arciniega-Ceballos M.A., González C.V., Dawson P., 2000, Temporal and spectral characteristics of seismicity observed at Popocatépetl volcano, central Mexico. J. Volcan. Geoth. Res., 102 , 207-216.         [ Links ]

Atlas Z.D., Dixon J.E., Sen G., Finnya M., Martín del Pozzo A.L., 2006, Melt inclusions from Volcan Popocatepetl and Volcan de Colima, Mexico: Melt evolution due to vapor-saturated crystallization during ascent. J. Volcan. Geoth. Res., 153 , 221-240. Doi:10.1016/j. jvolgeores.2005.06.010.         [ Links ]

Campillo M., Singh S., Shapiro N., Pacheco J., and Herrmann R., 1996, Crustal structure South of the Mexican Volcanic belt, based on group velocity dispersion. Geofísica Internacional, 35 , 361-370.         [ Links ]

Chavacán M., Lermo J., Quintanar L., 2004, Determinación de una escala de magnitud para temblores corticales en la parte central del Cinturón Volcánico Mexicano. In: I Conf. Int. Peligrosidad e Ingeniería Sísmica. Santiago de Cuba, Cuba: Universidad de Oriente. CD.ROM.         [ Links ]

Chouet B.A., 1996, Long-period volcano seismicity: Its source and use in eruption forecasting. Nature, 380 , 309-316.         [ Links ]

Cruz-Atienza V., Pacheco J., Singh S., Shapiro N., Valdés C., Iglesias A., 2001, Size of Popocatépetl volcano explosions (1997-2001) from waveform inversion. Geophy. Res. Lett., 28, 4,027-4,030.         [ Links ]

De la Cruz-Reyna S., Siebe C., 1997, Volcanology: The giant Popocatépetl stirs. Nature, 388, 227.         [ Links ]

De la Cruz-Reyna S., Yokoyama I., Bringas A.M., Ramos E., 2008, Precursory seismicity of the 1994 eruption of Popocatépetl Volcano, Central Mexico. Bull. Volc., 70, 753-767. Doi:10.1007/s00445-008-0195-0.         [ Links ]

De-Barros L., Pederson H., Métaxian J., Valdés-González C., Lesage, P. (2008). Crustal structure below Popocatépetl Volcano (Mexico) from analysis of Rayleigh waves. J. Volcan. Geoth. Res., 170 , 5-11.         [ Links ]

De-Cserna Z., Fuente-Duch M., Palacios-Nieto M., Trial L., Miltre-Salazar L.M., Mota-Palomino R., 1988, Estructura geológica, gravimétrica, sismicidad y relaciones geotectónicas regionales de la cuenca de México. Instituto de Geología, UNAM. Bolletin 104.         [ Links ]

Ego F., Ansan V., 2002, Why is the Central Trans-Mexican Volcanic Belt (102-99ºW) in transtensive deformation?, Tectonophysics, 359, 189-208.         [ Links ]

Espinasa-Perena R., Martín Del Pozzo A.L., 2006, Morphostratigraphic evolution of Popocatépetl volcano, Mexico. GSA Special Papers, 402, 115-137. Doi:10.1130/2006.2402(05).         [ Links ]

Espíndola J.M., Godínez M.D.L., Espíndola V.H., 2004, Models of ground deformation and eruption magnitude from a deep source at Popocatépetl Volcano, Central Mexico. Nat. Haz., 31, 191-207.         [ Links ]

Espíndola J.M., Macías J.L., Godínez M.L., 2003, Un modelo físico de cámara magmática para la estimacion de los cambios premonitorios en sísmicidad y deformación esperados en volcanes activos. Geos, 14 , 356.         [ Links ]

Fries C.Jr., 1965, Geología de la Hoja Cuernavaca, Estados de Morelos, México, Guerrero y Puebla. Instituto de Geofísica, Universidad Nacional Autónoma de México. Carta Geológica de México, scale 1:100 000.         [ Links ]

Geiger L., 1910, Herdbestimmung bei Erdbeben aus den Ankunftszeiten. Koenigliche Gesellschaft der Wissenschaften zu Goettingen, 4, 331-349.         [ Links ]

Got J.-L., Fréchet J., 1994, Deep fault plane geometry inferred from multiplet relative relocation beneath the southfank of Kilauea. J. Geophys. Res., 99 , 15,375-15,386.         [ Links ]

Havskov O.L.J., 2003, SEISAN: The earthquake analysis software, version 8.0 . Department of Earth Science, University of Bergen (ed.), 250 p. Bergen, Norway.         [ Links ]

Ito A., 1985, High resolution relative hypocenters of similar earthquakes by cross-spectral analysis method. J. Phys. Earth, 33 , 279-294.         [ Links ]

Jordan T.H., Sverdrup K.A., 1981, Teleseismic location techniques and their application to earthquake clusters in the south-central Pacic. Bull. Seis. Soc. Am., 71, 1,105-1,130.         [ Links ]

Karpin T., Thurber C., 1987, The relationship between earthquake swarms and magma transport: Kilauea Volcano, Hawaii. P. App. Geophys., 125 , 971-991.         [ Links ]

Klein F.W., 1985, User's Guide to HYPOINVERSE, a Program for VAX and Professional 350 Computers to Solve for Earthquake Locations. Open-File Report 85-515 U.S. Geological Survey.         [ Links ]

Lee W., Lahr J., 1972, HYPO71: A computer program for determining hypocenter, magnitude, and first motion pattern of local earthquakes. U.S. Geological Survey. Open File Report.         [ Links ]

Lermo Samaniego J., Antayhua Vera Y., Chaván Avila M., 2006, Análisis de la actividad en el volcan Popocatépetl (México) durante el periodo 1994-1997. Bol. Soc. Geol. México, pp. 253-257.         [ Links ]

Lienert B., 1994, HYPOCENTER 3.2 Manual. Honolulu: HIGP.         [ Links ]

Lienert B., Berg E., Frazer L.N., 1986, HYPOCENTER: An earthquake location method using centered, scaled and adaptively damped least squares. Bull. Seis. Soc. Am., 76, 771-783.         [ Links ]

Lienert B., Havskov J., 1995, Computer program for locating earthquakes both locally and globally. Seis. Res. Lett., 66, 26-36.         [ Links ]

Lomnitz C., 2006, Three Theorems of Earthquake Location. Bull. Seis. Soc. Am., 96, 306-312. Doi:10.1785/0120050039.         [ Links ]

López Ramos E., 1983, Geología de México. (3rd ed.). Mexico, D.F., Editorial México.         [ Links ]

Macías J., Siebe C., 2005, Popocatépetl crater fillled to the brim: significance for hazard evaluation. J. Volc. Geoth. Res., 141, 327-330.         [ Links ]

Meritano-Arenas J.D., Rosales-Gómez J., Conde-Asiain A., 1998, Carta Geológica y Carta Morfoestructural del Popocatépetl. Consejo de Recursos Minerales Ministry of Comerce and Industry (Secretaría de Comercio y Fomento Industrial) (SECOFI).         [ Links ]

Minakami T., 1974, Seismology of Volcanoes in Japan: Physical Volcanology. Amsterdam: Elsevier.         [ Links ]

Mullineaux D., Crandell D., 1981, The 1980 eruptions of mount St. Helens, Washington. chapter The Eruptive History of Mount St. Helens. (p. 844). U.S. Geological Survey Professional Paper 1250.         [ Links ]

Nava A.F., 2010, A Local Earthquake Location Method, based on time difference fitting through Genetic Search, for Tectonic and Volcanological Applications. Comunicaciones Academicas (Open-file report) CICESE, 8, 195-209.         [ Links ]

Novelo-Casanova D.A., Valdés-González C., Ramírez-Olvera G., 2007, A numerical model for the mechanical behavior of Popocatépetl volcano (Central Mexico). Phys. Earth Plan. Int., 162, 99-110. Doi:10.1016/j.jvolgeores.2007.03.001.         [ Links ]

Poupinet G., Ellsworth W., Fréchet J., 1984, Monitoring velocity variations in the crust using earthquake doublets: An application to the Calaveras fault, California. J. Phys. Earth, 89, 719-731.         [ Links ]

Ramírez-Olvera, G., 2003, Mallado tridimensional para elemento finito del Volcán Popocatépetl, correlacionado con la sismicidad. Bachelor's thesis Universidad Autónoma Nacional de Mexico.         [ Links ]

Roberge J., Delgado-Granados H., Wallace P.J., 2007, Determination of the degassing depth at Popocatépetl volcano (Mexico) using data from olivine-hosted melt inclusions. In American Geophysical Union, Fall Meeting 2007. San Francisco. Abstract No. V42B-05.         [ Links ]

Rubin A., Gillard D., Got J.-L., 1998, A re-examination of seismicity associated with the January 1983 dike intrusion at Kilauea volcano, Hawaii. J. Geoph. Res., 103, 10,003-10,015.         [ Links ]

Shapiro N.M., Campillo M., Paul A., Singh S., Jongmans D., Sanchez-Sesma F., 1997, Surface wave propagation across the Mexican volcanic belt and the origin of the long period seismic wave amplication in the Valley of Mexico. Geophy. J. Int., 128, 151-166.         [ Links ]

Shearer P., 1997, Improving local earthquake locations using the l1 norm and waveform cross-correlation: Application to the Whittier Narrows, California, aftershock sequence. J. Geophy. Res., 102, 8,269-8,283.         [ Links ]

Siebe C., Abrams M., Macías J.L., Obenholzner J., 1996, Repeated volcanic disaster in prehispanic time at Popocatépetl, Central Mexico: Past key to the future? Geology, 24, 399-402.         [ Links ]

Siebe C., Macías J., 1997, Volcanic hazards in the Mexico City metropolitan area from eruptions at Popocatépetl, Nevado de Toluca, and Jocotitlán stratovolcanoes and monogenetic scoria cones in the Sierra de Chichinautzin volcanic eld. Field Trip Guidebook, Geol. Soc. Am., 1, 77.         [ Links ]

Slunga R., Roegnvaldsson S., Boedvarsson R., 1995, Absolute and relative locations of similar events with application to microearthquakes in southern Iceland. Geophy. J. Int., 123, 409-419.         [ Links ]

Straub S.M., Martín-Del-Pozzo A.L., 2001, The significance of phenocryst diversity in tephra from recent eruptions of Popocatépetl Stratovolcano (central Mexico): Contributions to Mineralogy and Petrology. Geophy. J. Int., 140, 487-510.         [ Links ]

Valdés C., González G., Arcieniega A., Guzmán M., Nava E., Gutiérrez C., Santoyo M., 1995, Volcan Popocatépetl, Estudios realizados durante la crisis de 1994-1995. chapter Sismicidad del Volcán Popocatépetl a partir del 21 de Diciembre de 1994 al 30 de Marzo de 1995. pp. 129-138). México, D.F.: CENAPRED-UNAM (ed.         [ Links ]).

Valdés C.M., Mooney W., Singh S., Meyer R., Lomnitz C., Luetgert J.H., Hesley B., Lewis B., Mena M., 1986, Crustal structure of Oaxaca, Mexico from seismic refraction measurements. Bull. Seis. Soc. Am., 76, 547-564.         [ Links ]

Valdés-Gonzalez C.M., Comité C.d.V.P., 1994, P-wave 1D velocity model, Popocatépetl, Mexico.         [ Links ]

Waldhauser F., 2001, HypoDD - A program to compute Double-Difference hypocenter locations. U.S. Geol. Survey Open File Report, 1-113, 1-25.         [ Links ]

Waldhauser F., Ellsworth W.L., 2000, Double difference earthquake location algorithm: method and application to the Northern Hayward fault, California. Bull. Seis. Soc. Am., 90, 1,353-1,368.         [ Links ]

Wolfe C.J., 2002, On the mathematics of using difference operators to relocate earthquakes. Bull. Seis. Soc. Am., 92, 2,879-2,892.         [ Links ]

Wright R., De la Cruz-Reyna S., Harris A., Flynn L., Gómez-Palacios J., 2002, Infrared satellite monitoring at Popocatépetl: explosions, exhalations, and cycle of dome growth. J. Geophy. Res., 107 B8, 1,029-1,045.         [ Links ]