Repeated seismic swarms near Paricutin volcano: precursors to the birth of a new monogenetic volcano in the Michoacán-Guanajuato volcanic field, México?

The birth of a new monogenetic volcano is difficult to forecast with precision, both in space and time. Nevertheless, seismic activity can alert of the imminence of such an eruption because it usually occurs as small-magnitude earthquake swarms that can last for a few weeks to months prior to an eruption. These swarms are usually related to magma that becomes stalled in the Earth’s crust for variable periods of time before its eventual eruption at the surface. For several reasons, volcanic seismic swarms have rarely been recorded with seismometers before the birth of a new monogenetic volcano. Over the past 25 years, six distinct seismic swarms (in 1997, 1999, 2000, 2006, 2020, and 2021) were detected between Tancítaro and Paricutin volcanoes, in the southwestern part of México’s Michoacán-Guanajuato volcanic field. They are believed to represent repeated attempts of magma to reach the surface hinting that in this region magma might become stalled for some time, so as to not reach the surface in a single ascent event from its source in the mantle. To better understand the magma’s migration path through the crust, we re-located with greater precision some of these seismic swarms by using the same methodology and velocity model to the entire data set. Our results show that these swarms originated within a small area beneath the NE flank of Tancítaro at depths of between 15 and 8 km below sea level (bsl). Apparently, magma is trying to reach the surface within the same conduit network at these crustal depths, but stalls when reaching a depth of ~ 8 km bsl. It is crucial to study these swarms because they might be precursors to a new eruption in this part of the Michoacán-Guanajuato volcanic field. This monogenetic field has been very active, producing several dozen eruptions during the Holocene, the last two Jorullo (1759–1774) and Paricutin (1943–1952). Furthermore, the Tancítaro area displays one of the highest densities of Holocene volcanoes within the entire field, making it a probable candidate location for the birth of a future monogenetic volcano. For these reasons, a permanent seismic network should be installed as soon as possible.


Introduction
When magma erupts at the surface it can either give birth to a new monogenetic volcano or produce an eruption at a long-lived polygenetic volcano (e.g., Gudmundsson et al. 2001;Schmincke 2004;Zellmer and Annen 2008;Cashman and Sparks 2013;Martí et al. 2016;Galland et al. 2018). However, in many cases ascending magma does not reach the surface, and there may be many magma intrusions (with associated seismic swarms) before the magma finally reaches the surface (e.g., Hill et al. 1990;Roman et al. 2004;Legrand et al. 2011). During intrusive events, the magma will stall at depth, sometimes being there many years or even decades, as shown by the occurrence of repeated seismic swarms. Generally, eruptions are preceded by several seismic swarms, often repeating over years, without concomitant eruption, until the eruption finally occurs (Hill et al. 1990;Roman et al. 2004). Some famous cases of repetitive seismic swarms preceding magmatic eruption occurred in La Palma on the Canary Islands (Torres-González et al. 2020), Izu Peninsula and Unzen in Japan (Aoki et al. 1999;Nakada et al. 1999), Chaitén in Chile (Lange et al. 2008), Chichón in México (Yokoyama et al. 1992), and Soufriere Hills in the British West Indies (White and McCausland 2016). At La Palma, seismicity was felt years before the three 1949 eruptive vents of San Juan opened, and months to weeks before the 1971 Teneguía eruption, both located inside the monogenetic field at the southern part of the island (Torres-González et al. 2020). Before the September 2021 Tajogaite (Cumbre Vieja) eruption, several seismic swarms were recorded, for example in October 2017 and February 2018, characterized by high b-values (> > 1) of the Gutenberg-Richter law, suggesting the presence of magmatic fluids, and an impending eruption (Torres-González et al. 2020). The earthquake-depths of the 2017 and 2018 swarms varied between 15 and 25 km, and between 25 and 35 km respectively. During these two swarms, no permanent surface ground-deformation was observed, as can be expected for deep intrusion (Torres-González et al. 2020). However, small transitory elastic deformation events, mainly recorded on the vertical component of a Global Positioning System (GPS), were intermittently observed over short periods during the 4 months between the two seismic swarms and continued for at-least 3 months after the second seismic swarm. Transient surface ground-deformation (of up to 2 cm) were observed (Torres-González et al. 2020), which is high for such a deep intrusion. The seismicity migrated over time, from west to east, with the swarms being interpreted as failed eruptions or ascending magmatic intrusion that stalled (Torres-González et al. 2020). The events, though, were argued to culminate with the 2021 Tajogaite eruption. At the Izu Peninsula, more than 12 seismic swarms occurred from 1978 to 1997. The 1989 swarm terminated with a submarine eruption (Aoki et al. 1999). Several sporadic distal VT swarms were recorded between 1984 and 1990, before the 20 May 1991 Unzen eruption (Nakada et al. 1999). Three years before the 2008 Chaitén eruption (Chile), a seismic swarm was recorded in 2005 by a local temporal seismic network (Lange et al. 2008). Several seismic swarms occurred tens of years before the 29 March 1982 Chichón VEI = 5 eruption (México), in the 1920s, 1930, 1964, 1967-1968, 1980-1981(Yokoyama et al. 1992. Before the 18 July 1995 Soufriere Hills (Montserrat, British West Indies) eruption, intense VT swarms occurred in January 1992, January 1993, 1994, early 1995, and June 1995(White and McCausland 2016. In many cases, the distal VT swarms occurred months or years before the eruption and were not located directly below the eruption site, but were displaced laterally by several to tens of kilometers (White and McCausland 2016).
Such intrusion sequences without eruption, also called "failed magmatic eruptions" (e.g., Hill et al. 1990;Roman et al. 2004;Moran et al. 2011;Martí et al. 2016) have been detected beneath both polygenetic volcanoes and monogenetic (distributed) volcanic fields (Hill et al. 1990;Roman et al. 2004;Gardine et al. 2011;Legrand et al. 2011;Moran et al. 2011;Kundu et al. 2012;Torres-González et al. 2020). How and at which velocity the magma migrates upwards (mainly vertically, but also laterally) before reaching the surface, particularly within the last few kilometers of the upper crust, is still an open question, especially in the case of monogenetic volcanoes. The analysis of earthquake locations during the repeated seismic swarms associated with each magma intrusion allows us to detect and follow such magma migrations in both space and time.
The easiest way to detect magma migration is by recording seismicity induced by these movements (e.g., Toda et al. 2002;López et al. 2012;Becerril et al. 2013;Woods et al. 2018). Seismicity associated with such migration is generally composed of volcano-tectonic earthquakes (VTs) and/or (very) long-period seismic events (VLPs or LPs), hybrid events and/ or tremor (Moran et al. 2011). Volcano-tectonic earthquakes are commonly recorded before and during eruptions, due to stress-tensor changes inside the volcano as, for example, generated by the propagation of magma under pressure. In space and time, VTs can be well separated when magma is stalled, or can be clustered as swarms when magma is migrating. The vertical or horizontal migration of seismicity has been used to forecast eruptions (e.g., Matoza et al. 2014;Michon et al. 2015;Lengliné et al. 2016;White and McCausland 2016;Roman and Cashman 2018;Einarsson and Brandsdóttir 2021). VTs can occur on shear faults connecting the edges of offset opening dikes, without clear temporal evolution (Hill 1977). VTs and long-period earthquakes can also occur at the tip of a propagating dike, with a clear temporal evolution in their hypocenters following the migration of the tip of the dike (Ukawa and Tsukahara 1996). VTs can alternatively occur at the walls (away from the tips and edges) of an inflating dike, with a random temporal distribution (Roman and Cashman 2006). Therefore, the distribution of the seismicity does not directly reflect the shape of the active conduit. Instead the zone of seismicity is generally larger than the conduit, a difference which is exaggerated by the unprecise locations of the VTs.
Different kinds of focal mechanisms are expected to occur in this magma migration process. Earthquakes located at the tip of the opening dike and along the dike are directly related to the overpressure of the magma inside the opening tensile crack. They are expected to be of small magnitude , with a focal mechanism of a rupture of mode I (tensile crack, Rubin 1993). The polarities of the first P-wave arrivals are then pure compression when the crack opens or pure dilatation when the crack closes. But other kinds of earthquakes can occur, and these may not be located in the dike itself but around it and are due to stress changes in the intruded rock generated by the opening of the dike itself. Such seismicity may have more classical mode II (double-couple) focal mechanisms , with a combination of compression and dilatation of the first P-wave polarities. Such events have a focal mechanism that is not related to the regional tectonic stress-tensor, but instead to the local volcanic stress-tensor as modified by the magma migration itself (Legrand et al. 2002;Roman and Cashman 2006).
Seismic swarms of VTs can be tectonic or volcanic in origin. In some cases, a "purely" tectonic origin can be established when solely tectonic faults are involved, without any volcano existing near to the swarm. When a seismic swarm occurs near to a stratovolcano or within a monogenetic volcanic field, the concept of a "purely" volcanic origin is misleading because faults are always present in such a geological setting. In this case, swarms generally have mixed, volcanic and tectonic components, and it is difficult to determine their relative contributions. For example, when a seismic swarm is very shallow and occurs just prior to an eruption, it can be associated to magma movement. But at the same time, this magma usually uses a pre-existing shallow fault to reach the surface (Gudmundsson et al. 2001;Schmincke 2004). This is the main reason why the associated earthquakes are called "volcano-tectonic" referring to both possible origins. In contrast, if the depths of these earthquakes are greater, with a diffuse spatial distribution not related to any fault, the origin may be "purely" volcanic (or magmatic). Other criteria used to distinguish a tectonic from a volcanic origin derive from the Gutenberg-Richter and Omori laws. For tectonic earthquakes, the b-value of the Gutenberg-Richter law is ≤ 1, whereas for a volcanic origin or when fluids are involved the b-value is > 1. The Omori law usually applies to a tectonic origin, such as for aftershock sequences, but does not apply for volcanic origins (e.g., Havskov et al. 1983;Medina et al. 1990;Yokoyama et al. 1992;Legrand et al. 2015).
Six distinct seismic swarms have occurred over the past 25 years, between Tancítaro and Paricutin volcanoes within the Michoacán-Guanajuato monogenetic volcanic field of México. The region around the polygenetic Tancítaro volcano has one of the highest concentrations of Holocene monogenetic volcanoes in the world, including the field's youngest, Paricutin, born in 1943. We show that these six seismic swarms are associated to the generation of six small magma batches. We suspect that these repeated six seismic swarms with increasing frequency the last years could be precursors to a new eruption, and, thus, that the next monogenetic volcano could be born in this region. This makes it necessary to define, study, interpret, and track seismic swarms in regions prone to monogenetic eruptions, as we do here.

Geological setting
The Michoacán-Guanajuato volcanic field (MGVF) occupies an area of ~ 40,000 km 2 within the Trans-Mexican volcanic belt (TMVB) and is one of the largest monogenetic volcanic fields in a subduction context on Earth (Hasenaka and Carmichael 1985a;Valentine and Connor 2015). The MGVF corresponds to the widest segment of the TMVB, associated to a relatively flat subduction with associated depths of ~ 80-120 km for the oceanic Cocos plate (Pardo and Suárez 1995, Fig. 1). Partial hydrous melting of the mantle wedge may be favored by flat subduction over a wide area underneath a ∼40-km-thick continental crust (Chevrel et al. 2016). The MGVF is composed of more than 1,100 scoria cones and associated lava flows, ~ 400 small-to-mediumsized shield volcanoes, ~ 22 phreatomagmatic vents (maars and tuff rings), ~ 43 lava domes and isolated lava flows (e.g., Hasenaka and Carmichael 1985a, b;1987;Mahgoub et al. 2017), and two extinct stratovolcanoes (Ownby et al. 2007). The volumes of the MGVF shield volcanoes vary between 0.5 and 10 km 3 and are much larger than those of scoria cones (average of 0.021 km 3 ), but much smaller than the ∼50-100 km 3 volumes of the two stratovolcanoes (Tancítaro and Patamban) in the MGVF (Chevrel et al. 2016). The spatial distribution of monogenetic volcanoes is generally controlled by regional and local tectonics (Valentine and Perry 2007). Within volcanic fields, monogenetic volcanoes are often not only clustered in space (Hasenaka and Carmichael 1985b;Pérez-López et al. 2011), but also in time. This is also the case in the MGVF where several Holocene clusters have been identified (e.g., Hasenaka and Carmichael 1985a, b;1987;Mahgoub et al. 2017). This information is essential for volcanic hazards because it helps in identifying areas where the next monogenetic eruption is most likely to occur. The exact time and place of birth for a monogenetic volcano are difficult to forecast, but continuous monitoring of seismic activity and other parameters, including ground deformation, hydrothermal activity, and gas emissions, can all help constrain the final phase of magma ascent prior to eruption (e.g., Bell et al. 2018;Peltier et al. 2018;Sigmindsson et al. 2022). The persistent occurrence of seismic swarms during the last 25 years beneath the NE flank of Tancítaro stratovolcano are signs of deep magma intrusions that may eventually reach again the surface in this area. Tancítaro is extinct, but is surrounded by numerous young monogenetic volcanoes. Thus, this area represents an ideal case to study the potential precursory phenomena to the birth of a new monogenetic volcano. In this work, we thus compare the 1997,1999,2006, and 2020 seismic swarms to better assess monogenetic eruption hazards. Connor (1987) and Kurokawa et al. (1995) showed that the margins of the MGVF are characterized by normal faulting (horst-and-graben structures), revealing an overall extensional stress-regime. In this regime, the maximum compressive stress σ 1 is vertical, the maximum horizontal compressive stress σ 2 is oriented in a NE-SW direction, and the smallest stress σ 3 has a NW-SE orientation. In such a configuration, we can expect that a shallow vertical dike would be oriented NE-SW and open in an NW-SE direction (Gardine et al. 2011).

Regional tectonics
The Tancítaro andesitic stratovolcano (3860 m asl) is the highest volcano of the MGVF, with a volume of ~ 100 km 3 (Ownby et al. 2007). It is considered extinct since its last effusive activity occurred well before the Last Glacial Maximum (Vázquez-Selem and Bocco-Verdinelli 2016) and has been dated at 237 ± 34 ka (Ownby et al. 2007). In contrast, the density of young monogenetic volcanoes around Tancítaro is among the highest in the entire MGVF and attests to a high frequency of recent volcanic activity (Connor 1987).
In addition, several of the youngest (Late Holocene) monogenetic volcanoes of the MGVF (Fig. 1) occur either on Tancítaros's lower flanks (e.g., El Astillero, El Pedregal, and Paricutin), or in close vicinity to the north (e.g., Janamo) and east (e.g., El Jabalí and El Metate). Among these monogenetic volcanoes, Paricutin is the most recent and is located just ~ 10 km from Tancítaro's summit on its lower northeastern flank.
Paricutin was famously born in a cornfield on 20 February 1943, and its activity lasted for more than nine years until 4 March 1952 (e.g., Foshag and Gonzalez 1956;Luhr and Simkin 1993;Hasenaka and Carmichael 1987). We here use the term "birth" because it is often used for monogenetic volcanoes. In our case, there was no volcano in the cornfield before the sudden appearance of the new Paricutin volcano in 1943, which is a birth. This taxonomy is also used in the three books: "Birth and Development of Paricutin volcano, México (Foshag and Gonzalez 1956), "Paricutin: The volcano Born in a Mexican Cornfield" (Luhr and Simkin 1993) and "Volcanoes" (Decker and Decker 2005). Paricutin's cone grew mainly during the first year of its activity, and two adventive vents (Sapichu and Taqui) became aligned with the cone in a NE-SW direction (Fig. 2, UNAM-IGEOL 1945).
The El Astillero and El Pedregal monogenetic volcanoes are ~ 3.5 km apart and were formed during the same eruption around AD 500-700 (Larrea et al. 2019b). They are located on the lower SW slope of Tancítaro, ~ 15 km from its summit, and ~ 25 km to the SW of Paricutin. Paricutin, El Astillero and El Pedregal, as well as Tancítaro's summit, are also aligned in a NE-SW direction, and are sited along a major fault that is marked by a scarp (Fig. 2). This NE-SW oriented fault and its subsidiary system was probably used by dikes that ascended through the last kilometers of the crust to feed these eruptions (cf., Gudmundsson et al. 2001;Schmincke 2004;Gudmundsson 2022).
Since El Astillero, El Pedregal, and Paricutin are quite close to each other, are related to the same regional fault, and were all born in recent times (AD 500-700 and AD 1943, respectively), it is not surprising that the area still displays signs of activity where the seismic swarms examined here are indicative of a high potential for the birth of a new monogenetic volcano in this region.
Lava and tephra compositions of volcanoes in the MGVF display a great diversity ranging from basalts to rhyolites (e.g., Carmichael 2002;Chevrel et al. 2016;Larrea et al. 2019a, b). Most products are, however, andesitic and follow a calc-alkaline trend (e.g., Carmichael 2002;Chevrel et al. 2016;Larrea et al. 2019a, b). Rocks vary significantly in their major-and-trace element, and isotopic compositions revealing a complex origin from a heterogenous mantle wedge affected by significant hydrous subduction contributions and only a small degree of crustal contamination (Rasoazanamparany et al. 2016;Larrea et al. 2019a). These compositions are believed to result from ascent of several individual small magma batches that become stalled in the middle/upper crust before a new fresh batch from mantle depths intrudes them and triggers their final, rapid ascent and subsequent eruption (Larrea et al. 2017(Larrea et al. , 2019b(Larrea et al. , 2021. The process is also consistent with their differing crystal cargos, as well as by changes in composition during the course of an eruption. This is especially noticeable when eruptions last for several years (Chevrel et al. 2016).

Methodology
Six seismic swarms were detected between Tancítaro and Paricutin in the past 25 years, in 199725 years, in , 199925 years, in , 200025 years, in , 200625 years, in (Pacheco et al. 1999Gutiérrez and Valdés-González 2000;Gardine et al. 2011;Chacón 2013;Pinzón et al. 2017), and in 2020 and 2021 (Jácome-Paz et al. 2022). However, the geometries of the seismic networks that recorded these swarms were all different. In addition, the velocity models and methodologies used to locate the earthquakes were also different leading to several difficulties when comparing the results. We thus reprocessed the original seismic data and relocated the earthquakes of the 1997, 1999, 2006, and 2020 swarms using the same methodology. The 2000 and 2021 swarms were not recorded by a local temporary seismic network and thus are not reprocessed here. For the 1997 and 1999 swarms, the events were automatically detected using the STA/LTA technique (ratio between the averages of the absolute amplitudes of the signal calculated on a short-time window and on a long-time window, Withers et al. 1998), before they were recorded, whereas for the 2006 and 2020 swarms, detections were made manually from the continuously recording data. The P-and S-wave arrival times of the largest earthquakes were manually picked by the same person (L.D.) for all the studied swarms. This was necessary because many of the original pickings from earlier swarms had been lost or became unavailable, and because homogeneous picking criteria were required for comparing the resulting locations of all the swarms. We checked our readings for possible errors by applying the generalized Wadati-diagram test (double differences between the P-and S-wave arrival times for each earthquake at all the stations, see Legrand et al. 2021 for details). This generalized Wadatidiagram test also allows calculation of the Vp/Vs = 1.71 ratio for all the data sets by assuming that this ratio remained constant over the past 25 years, as confirmed in Fig. 3. Independently, we extracted a local velocity model (see Table 1) from the regional seismic noise tomography of Spica et al. (2016). We found that the modeled mean Vp/Vs = 1.70 between the surface and the depth of 10 km is almost identical to the derived 1.71 ratio. Note that the Vp was not directly measured in the tomography but was deduced from the Vs value measured in the tomography and the Brocher's polynomial relation between Vp and Vs (Brocher 2005). For Vs values between 2.2 and 4 km/s, Brocher's relation gives a Vp/Vs ~ 1.70, corresponding to expected values at depths between 3 and 30 km. We relocated all the earthquakes of the 1997, 1999, 2006, and 2020 swarms with this velocity model, using the HYPOINVERSE-2000 code (Klein 2014). The data were unfiltered because of the good quality of the recorded data (which has high signal to noise ratios). The zero-level used for earthquake depth is sea level. The small epicentral distances derived from the close proximity of seismic stations and their small number (≤ 5) generated a small root-mean-square (RMS) error, generally < 0.1 s. This same methodology was applied to the 1997, 1999, 2006, and 2020 swarms, allowing us to compare the latitude/longitude of the hypocenters, and their depths.

The 1997, 1999, 2000, and 2006 swarms
The 1997 seismic swarm was first detected by the seismic network of the Servicio Sismológico Nacional (SSN), the Mexican national seismological service, in February. There was a peak in activity on 21 February, with five large earthquakes with local magnitudes between 3.9 and 4.0. This caused the SSN to install a local seismic network of five broad-band seismometers from the 1 March (blue inverted triangles in Fig. 2) which was operational until 18 March when it was removed following the cessation of seismic activity (Pacheco et al. 1999). We relocated 252 earthquakes of Pacheco et al. (1999) using at least four seismometers (blue dots in Figs. 2 and 4). These were detected between 2 and 15 March 1997.
The 1999 seismic network was installed by the SSN and comprised four broad-band seismometers (six are shown in Fig. 2 because two were moved to other locations during the swarm). We located 89 earthquakes from 4 to 7 June 1999, as shown by green dots in Figs. 2 and 4. These data have not yet been published.
The 2000 swarm was first noticed by local inhabitants around Tancítaro and Paricutin volcanoes and later officially reported by the Centro Nacional de Prevención de Desastres  (CENAPRED), the Mexican institution in charge of monitoring volcanoes. Unfortunately, this swarm was not recorded by a local temporary seismic network. Hence, we cannot show the locations of this swarm here, but it seems that they were located approximately near the 1997 swarm (Gutiérrez and Valdés-González 2000). For the 2006 swarm, we used the closest ten broadband seismometers of the MARS (Mapping the Rivera Subduction Zone experiment of Yang et al. 2009). These are the yellow inverted triangles in Fig. 2. The MARS data are available on the USGS Earthquake Hazards Program website: https:// earth quake. usgs. gov/. The closest station was located 15-20 km from the swarm, and four stations were located at < 60 km. We processed the largest and bestlocated 143 earthquakes that occurred between 28 May and 29 June 2006. The 2019-2020 seismic experiment and 2020 swarm We installed a 120 s QA-Nanometrics broad-band seismometer (which has a bandwidth between 120 s and 145 Hz) and a Centaur acquisition system on 10 December 2019 at a location near the epicenter of the 1997, 1999, 2000, and 2006 swarms (Figs. 2 and 4). The purpose of this installation was to determine the frequency of small magnitude earthquake swarms in this area. Only 25 days after installation, the 2020 seismic swarm started. Once the first earthquakes were detected by the SSN, four additional broad-band seismometers were installed around Tancítaro (Fig. 2). Two were Nanometrics Trillium 120 s QA-instruments, with a Centaur Nanometrics 24-bit acquisition system, one was a Nanometrics Meridian compact Posthole of 120 s, and two were Guralp CMG-40 T of 60 s, with a CMG-DAS Guralp 24-bit acquisition system. The sampling rate was set at 200 Hz at the beginning of the experiment and switched to 500 Hz after 2 months. The stations had GPS time synchronization, but were without real-time data transmission due to the lack of equipment and difficulty of data transmission in the area. Here, we show the bestlocated 475 earthquakes of 2 February 2020 as detected by the three first seismometers installed, which have very clear P-and S-wave recordings (see one example in Fig. 5). In this data set, the P-waves were very weak on the north-south and east-west horizontal components, because the earthquakes originated almost directly below the station. This makes it difficult to determine the azimuths of the earthquakes by using P-waves from a single station (and hence, to locate the earthquakes). For this reason, we only show here locations obtained using data from these three local stations.

The 2021 swarm
From 30 May to 3 September 2021, an additional swarm was recorded by the SSN and located within the red contour of Fig. 2 (Jácome -Paz et al. 2022). However, these data are not processed here due to the lack of nearby seismometers, the associated large errors on location, and the difficulty of using the same methodology. This meant that results would not be comparable with the locations of the 1997,1999,2006, and 2020 swarms. We summarize the main characteristics of all six swarms in Table 2. The number of earthquakes for each swarm is not given because this strongly depends on the number of seismometers, their locations and the method used to detect them. The depths for all the swarms are mainly between 8 and 17 km bsl (Fig. 4), with largest magnitudes of 3.9-4.2, and typical durations of  1 month, although being longer for the last swarms of 2020 and 2021 which also had the highest magnitudes and the smallest time interval (Table 2).

Discussion
Similar location of the swarms The 1997The , 2000, and 2020 swarms were recorded by different seismic networks, of different geometries and with a different number of stations, and were initially processed using different velocity models and by different people. To compare the earthquake locations of the different swarms, it is not just important to apply a standard methodology, as done here, but to perform error estimates. Such errors depend mainly on: • the number of the seismometers (indicated in Fig. 6); • their spatial distribution with respect to the source of the seismicity (for example there is a station-gap given in Fig. 5); • the quality of the P-and S-wave arrival pickings, which depends mainly on the quality (signal to noise ratio) of the data; • the skills of the person reading arrival times; • the velocity model; • the method used to locate the hypocenters.  Fig. 6 Histograms of the errors in latitude, longitude (converted to km), and depth for the 1997, 1999, 2006, and 2020 swarms. Also shown is the number of seismometers used to locate the earthquakes, the gap (in degrees) of the stations for each earthquake location, and the RMS (root mean square) error (in seconds)

Page 11 of 15 30
In this study, the number of stations and their distributions (issues 1 and 2) are different for each case (swarm). However, a standard methodology applied in reprocessing the data makes error associated with issues 3 through 6 common to all cases. In general, errors automatically provided by a computer program do not consider all six error sources considered here, hence errors are often under-estimated. The main source of error generally stems from the velocity model, plus the number and spatial distribution of the seismometers. Earthquake location errors due to the P-and S-wave arrival time pickings are generally smaller, especially for local earthquakes. The velocity model in our case (Table 1) is well constrained by an independent method (seismic noise tomography, Spica et al. 2016). In our case, the main errors thus derive from the number and spatial distribution of seismometers. This is apparent for the 2006 swarm. Although more stations (10) were available than for the 1997, 1999, and 2020 swarms, in 2006 the seismometers were located beyond the source of the swarm and at greater distances. Therefore, absolute errors were larger. This shows the importance of having a permanent local seismic network to monitor future swarms. Despite these issues, the overall errors for the 1997, 1999, 2006, and 2020 swarms are relatively small (only a few km, see Fig. 6), so that we were able to confirm that these swarms all occurred within a small region beneath the northern flank of Tancítaro volcano.

Magmatic origin of the swarms
The volcano-tectonic earthquakes of the 1997, 1999, 2006, and 2020 swarms have similar locations (Figs. 2 and 4). Although we could not check the locations of the 2000 swarm, a report by CENAPRED (Gutiérrez and Valdés-González 2000) mentions that the earthquakes had similar epicenters to those of the 1997 swarm. An important question is whether these swarms are magmatic or tectonic in origin. Since two prominent fault systems (the San Juanico-Buenavista fault striking NW-SE and another fault striking NE-SW, Fig. 2) intersect in the Tancítaro-Paricutin area, the 1997 swarm was initially thought to be purely tectonic in origin (Pacheco et al. 1999). Nevertheless, the b-value calculated for the 2006 swarm is high (1.47, Pinzón et al. 2017), much larger than values typically obtained for tectonic swarms, and interpreted by Pinzón et al. (2017) as a deep magma intrusion, as also suggested by Gardine et al. (2011). We agree with this conclusion and suggest that the hypocenter depths of 8-15 km below Tancítaro's NE flank reveal a deep magma intrusion. Furthermore, this seismicity cannot relate to the reactivation of Tancítaro volcano because it is considered extinct. Instead, it may represent an initial phase of activity prior to a new monogenetic eruption, that could occur in the vicinity of the present swarm's epicenter. The next monogenetic volcano will not necessarily be born at the seismic swarm epicenter because the magma still has to migrate another ~ 8 km through the crust, and could follow a complicated path due to vertical stress change, or by intersecting pre-existing faults, dikes, and sills. An eruption associated with these epicenters will unlikely, though, occur beyond a radius of a few tens of kilometers from the swarms. Magma will most likely use the pre-existing shallow faults during its ascent through the last km of the upper crust, for example the NE-SW regional fault of Fig. 2, as did the magmas of the youngest monogenetic volcanoes in this area. The same location of the different swarms suggests that the magma is following the same vertical dike network (below ~ 8 km depth), but is not necessarily using the same dikes.

Do the VT locations reflect the magma storage region?
For the 2020 swarm, the magma stalled at a depth of ~ 8 km, which is similar to the mineral-equilibration depth determined in petrological studies of volcanoes in the MGVF, indicating this is a zone of temporary magma storage and degassing (Chevrel et al. 2016). At El Metate, Chevrel et al. (2016) found a magma storage region at depths between 7 and 10 km below the average altitude of El Metate (~ 2 km asl), corresponding to depths between 5 and 8 km bsl and being comparable to our depths of 8-15 km bsl. At other nearby monogenetic volcanoes (i.e., the Zacapu cluster), petrological estimates yield a magma storage region at depths of between 6-8 and 13 km bsl after a 2 km average altitude correction (Reyes-Guzmán et al. 2021), which is also comparable with the region here.
The deeper seismicity observed at 15 km bsl corresponds to magma rising from greater depths. An open question is the relationship that exists between the initiation of the 8-15 km depth seismicity and the 5-8 km depth of magma storage. We consider that the bottom of the batch may be located at 15 km depth, from which magma starts to migrate vertically up dikes and/or horizontally along sills, forming a magma reservoir agglomeration (Fig. 7). Magma usually stalls at a depth corresponding to the level of neutral buoyancy, with a small overshoot of a few kilometers (Lister and Kerr 1991). Magma accumulated at this level is in a gravitational equilibrium but can be intruded by dikes due to buoyancy forces (Lister and Kerr 1991). The magma storage region is identified by the locations of our seismicity. As earthquakes do not occur in magma, this region cannot be a chamber completely filled with magma but rather a heterogeneous region composed of aseismic magma-filled dikes surrounded by colder and more fragile seismogenetic zones. This storage zone architecture was proposed by Hill (1977) with clusters of magma-filled dikes within brittle volumes of the crust, with sequences of earthquake swarms occurring along a system of conjugate fault planes joining in echelon offset dike tips at oblique angles. We illustrate this in Fig. 7 where earthquakes (stars) are plotted beyond the zone of dikes represented by red « fingers».
Our earthquake depths are less well constrained than the epicenters, as is often the case for earthquake locations due to the spatial distribution of the seismometers on a horizontal surface compared with the absence of seismometers located at depth. Figures 2 and 4 show a slight upward temporal migration of hypocenters from deeper in 1997 to shallower in 2020, although errors are admittedly large enough to cast doubt. Nevertheless, the vertical errors are much smaller (Fig. 6) than the depths (8-15 km bsl) of the hypocenters. These depths are away from the shallow faults observable in this area. Hence, they may not be related to these shallow tectonic faults but rather to filled-magma dikes as suggested by the horizontal epicenter distribution. The overall North to South "migration" observed in Fig. 4 of ~ 8 km seems more real. This N-S direction may be the direction of a magma-filled dike or a series of sub-parallel dikes. The errors of a few km on our locations (Fig. 6) are too large to distinguish closely separated dikes. The apparent temporal ~ 8 km-North-South horizontal earthquake migration registered from 1997 to 2020 (Fig. 4) does not necessarily imply horizontal magma migration over such a distance. Instead, it might correspond to several independent small vertical magma intrusions, each with its own concomitant dike network and associated seismic swarm. This motion is underscored by the fact that the seismicity did not move in a systematic fashion from N to S. The 1999 swarm is located to the N of the 1997 swarm, whereas the 2006 swarm is located to the S of the 1997 swarm (Fig. 4), and the 2020 swarm is located to the S of the 2006 swarm. Thus, the seismicity is first moving toward the N from 1997 to 1999 and afterwards to the S from 1997 to 2020. These seismic swarms may reflect repeated vertical magma migrations in the form of independent small magma batches, as proposed by petrological studies at nearby volcanoes (Larrea et al. 2017(Larrea et al. , 2019b(Larrea et al. , 2021.
Further monitoring of these swarms should detect future vertical and lateral migrations, which is fundamental for assessing volcanic hazard in the Tancítaro region. For these reasons, a permanent seismic network should continuously monitor the region because the final magma ascent could occur rapidly, leaving little time to install a new seismic network.

Conclusions
Although past seismic networks used in this study had very different spatial geometries, we show how an integration of such data can be used to confidently constrain storage zones below monogenetic systems. At Tancítaro, the 1997Tancítaro, the , 1999Tancítaro, the , 2000Tancítaro, the , 2006, and 2020 seismic swarms were located within a small region 8-15 km beneath the northeast lower flank of the Tancítaro stratovolcano. Swarms were probably caused by repeated attempts by magma to migrate upwards to the surface. The similar epicenters and depths of all four swarms analyzed suggest the repeated use of the same or nearby magma pathways, to form a sub-parallel dike network. The seismicity stops at a depth of ~ 8 km bsl, where the rising magma presumably stalls and accumulates in a small but growing reservoir. The apparent temporal North-South horizontal earthquake migration registered from 1997 to 2020 might correspond to several independent small vertical magma intrusions, each with its own concomitant dike network and associated seismic swarm.
Monogenetic volcanoes of the Michoacán-Guanajuato volcanic field are not distributed randomly, and several areas have a high frequency of occurrence. One of them is the Tancítaro area, where several late-Holocene volcanoes have been identified, including Paricutin. Hence, these seismic swarms may represent precursors to the birth of a new monogenetic volcano in this region.
The earthquakes forming these swarms have small magnitudes (typically < 4), last for short periods (~ 1 to 2 months), and have occurred repeatedly over a period of more than 25 years. Since monogenetic volcanoes occur at different locations in the Tancítaro-Paricutin area, the paths that their magmas take to ascend the last 8 km of the crust are different for each volcano, although some likely stem from a single deeper magma reservoir system at a depth of between ~ 8 km to at least 15 km. Our experience of monitoring and analyzing the seismic swarms of the Tancítaro-Paricutin area may allow us to place some bounds in the precursors to the birth of a new monogenetic volcano. Hence, it is important to continuously monitor the seismic activity of such monogenetic provinces with a local and continuous seismic network so as to study, understand, and track the spatio-temporal evolution of hypocenters. Given that an eruptive event may start within hours to days of the onset of a new swarm, in high-risk areas such networks should be installed well in advance. It is also important to monitor rates of soil degassing, changes in chemical composition of local spring waters, soil temperature, and ground deformation. A rapid final ascent of magma to the surface, as been observed following other "failed" eruptions, would provide only a short time to react.