Introduction
The origin and distribution of felsic volcanism along the Michoacán-Guanajuato Volcanic Field (MGVF) is not completely understood. The MGVF is recognized worldwide as one of the biggest monogenetic volcanic fields on Earth and as homeland of Parícutin volcano, the only scoria cone studied since its inception and through its entire evolution (Foshang and González, 1956). Probably less famous abroad, but remarkably historic within the MGVF, is Jorullo volcano, one of the volcanoes climbed by A. von Humbold during his expedition to America. A high percentage of volcanoes within the MGVF resemble Parícutin’s volcanic style, and most analyzed rocks (88 %, see references through text) fall within a compositional cluster characterized by mafic to intermediate compositions.
Given the overwhelming volume of mafic-intermediate rocks, felsic volcanic products (dacite-rhyolite) have been ignored in terms of their genesis and the magmatic processes that integrated their final composition. As opposed to intermediate rocks, silicic rocks are not ubiquitous in the MGVF, they are scarce, clustered in few locations and erupted in a wide span of time (see below). Understanding their distribution and the processes that differentiate their chemical composition can contribute to shed light on the genesis of felsic magmas in other volcanic monogenetic fields in Mexico and around the world.
In this contribution in honor to our friend Víctor Hugo Garduño Monroy, we present the results of a bibliographic research on rock chemistry, structural geology and experimental petrology of the MGVF. The main aim for using available published data is to establish the relative abundance of silicic rocks within the MGVF, and the relation between structural systems and the magmatic evolution in central Mexico. In addition, we discuss the results of two geophysical surveys; a seismic tomography to visualize at depth the plumbing system that bring magmas to the upper-crust, and a gravimetric-magnetometric study to investigate the depth at which the magmas might acquire their final silicic composition.
Geologic background
The MGVF is located in the central portion of the Trans-Mexican volcanic belt (TMVB), a continental arc built on the southern edge of the North America plate, which overrides the Rivera microplate and the northern part of the Cocos plate (Figure 1). The MGVF was originally described by Hasenaka and Carmichael (1985) as a monogenetic volcanic field formed by more than 1200 volcanoes, most of them monogenetic, covering an area of approximately 40000 km2. Scoria cones are ubiquitous along the volcanic field and are coeval in some locations with stratovolcanoes (Ownby et al., 2011), medium-sized shield volcanoes (Hasenaka, 1994; Osorio-Ocampo et al., 2018), fissural lava flows (Hasenaka and Carmichael, 1985a), maar-craters (e.g., Valle de Santiago; Zacapu), and felsic domes (e.g., Pérez-Orozco et al., 2018; Osorio-Ocampo et al., 2018). Inception of volcanism in the MGVF is dated at ~7 Ma (Avellán et al., 2020); although, the MGVF is worldwide known because of its recent eruptions: Parícutin in 1943-1952 and Jorullo in 1759-1774.

Figure 1 Regional tectonic regime of Central Mexico that shows the active subduction of the Rivera and Cocos Plates beneath the North American Plate at the Middle-American Trench. The dotted white line depicts the boundary of the Trans-Mexican Volcanic Belt (TMVB). The Michoacán-Guanajuato volcanic field (MGVF) and major volcanoes are shown for reference: La Pimavera caldera (LP), Colima (Co), Tancítaro (Ta) and Jorullo (Jo), the latter two located in the MGVF, Amealco caldera (Am), Huichapan caldera (Hc), Acoculco caldera (Ac), Los Humeros caldera (Hm), Nevado de Toluca (Nt), Popocatépetl (P), La Malinche (M), Cofre de Perote (Pe), Pico de Orizaba (O), San Martín (SM). GDAL: Guadalajara city, EPR: East Pacific Rise.
Rocks from the MGVF are mainly calc-alkaline (Hasenaka and Carmichael, 1987; Gómez-Tuena et al., 2005; Gómez-Tuena et al., 2018), although some intraplate and potassic rocks occur (Cavazos-Tovar, 2006; Ortega-Gutierrez et al., 2014b; Losantos et al., 2017). Rocks are predominantly intermediate, 40 % of all known rocks are andesites and 33 % are basaltic-andesites (Figure 2). Previous works proposed that the mafic and intermediate magmas erupted in the MGVF were derived from a combination of partial melting of a heterogeneous mantle doped with fluids from subducted sediments and oceanic crust (Larrea et al., 2019; Gilbaud et al., 2019), partial melting of the deep crust (Ownby et al., 2011), assimilation and fractional crystallization (AFC) processes (McBirney et al., 1987; Ceibrá et al., 2011; Losantos et al., 2017), and fractional crystallization (Luhr and Carmichael, 1985; Johnson et al., 2008). The origin of more felsic magmas in the MGVF has been related to the crystallization of mafic magmas and the assimilation of middle-crust lithologies (Ownby et al., 2011), and to the assimilation of upper-crust granites favored by extensional-transtensional tectonics (Pérez-Orozco et al., 2018).

Figure 2 Total alkalis vs. SiO2 diagram for rocks of the Michoacán-Guanajuato volcanic field (MGVF). Lavas= blue circles, ignimbrites= red circles. The whole rock composition of silicic rocks, granite xenoliths and more mafic rocks are presented in Supplementary File S1. Data were extracted from http://georoc.mpch-mainz.gwdg.de
The MGVF occurs in an extensional tectonic environment associated with plate boundary processes such as the rollback of the Cocos Plate (Singh and Pardo, 1993) and the sinistral rotation of the Michoacán block due to the oblique convergence of the Cocos plate relative to the North American plate. The crustal configuration underneath Michoacán-Guanajuato is defined by a complex arrangement of fault zones that control the kinematics of different regional blocks (e.g., Johnson and Harrison, 1989, 1990; Pasquaré et al., 1991; Garduño-Monroy et al., 2009; Guilbaud et al., 2012; Kshirsagar et al., 2016) and the distribution of magmatism (Avellán et al., 2020; Gómez-Vasconcelos et al. , 2020; Olvera-García et al., 2020). The MGVF is deformed by two main fault systems, the Morelia-Acambay fault system (MAFS) with ENE-WSW and NE trending faults and the Taxco-San Miguel fault system (TSMFS) with NNW-SSE and NW-SE striking faults (Alanís-Álvarez et al., 2002). The MAFS is displaced by reactivated older TSMFS faults oriented to the NNW.
Data collection methods
We obtained structural data from 38 articles and 11 theses published during the 1985-2020 period. We found 466 whole-rock analyses published in 39 articles and 4 theses during the 1954-2020 period. Data were extracted from the GEOROC global compilation. For the sake of space and fluency we report the data and references to articles and theses in Supplementary File S1. We also used experimental data from one thesis (Boijseauneau-López, 2018). We recognize that overall, volcanic products are not widely studied in terms of stratigraphy and distance from the vent; in most works, with only one sample is pretended to characterize a full deposit without acknowledging that lava fronts could be different to the last emissions tapped from the vent. The published works we used for this investigation lack information to ascertain if the deposits sampled are homogeneous or comprise different compositional domains. Thus, our results represent minimum relative percentages.
Geophysical methods
We used one ambient noise tomography from Spica et al. (2016). The tomography was performed around Colima volcano with a horizontal resolution of ~14×16 km2 that results from the inversion of the group and phase dispersion curves of the fundamental Rayleigh and Love surface wave. We superimposed the historical seismicity reported since 1900 and the down-going slab of the Coco's plate identified from the Wadati-Benioff zone (Figure 3).

Figure 3 Shear wave velocity tomography for the Colima-Tancítaro region. Left: location of earthquakes (black dots) along the section represented by the tomography. Top right: topography of the section. Bottom right: vertical section of the shear wave velocity (Vs) tomography. Shear wave velocity is color coded. Historical earthquakes are superimposed (empty circles); the size of the circles is proportional to the earthquake magnitude. The Cocos slab is shown as a red line.
Magnetometry data were acquired with an Overhauser Gem Systems GSM-19WG magnetometer according to the methods of López-Loera et al. (2021). Gravimetry data were acquired with a Scintrex Autograv CG5 gravimeter. Each point was georeferenced with a PENTAX G3100-R2 GPS. We performed a reduction to the pole (Baranov y Naudy, 1964) to all magnetic anomalies acquired in the Parícutin area during the 1943-1952 period. Geomagnetic reductions were performed with data from the Geomagnetic Observatory of Teoloyucan. In order to reduce the short wavelength signals, the magnetic and the Bouguer anomalies were treated with a 50 m ascendant analytic continuum (Henderson, 1970); afterwards, we applied a regional-residual separation (LaFehr y Nabmposighian, 2012).
The 2D model was generated with the GM-SYS software. The anomalies are modeled through polynomial bodies with variable density and magnetization. Input values are as follows: 1) magnetic and density susceptibility from Telford et al. (1990); 2) Foshang and González (1956) observations; 3) magnetic parameters from Urrutia-Fucugauchi et al. (2004); 4) topographic model from Corona-Chávez (2002).
Geophysical results
Figure 3 shows a vertically polarized shear wave (SV) tomography with a nearly horizontal anomaly at 35 km depth. The anomaly begins NE of Colima volcano and goes toward the Tancítaro area. The shear wave anomaly is a zone with lower velocity than the surrounding velocity average; it is more pronounced on the vertically polarized shear wave (SV) tomography than on the horizontally polarized shear wave (SH) tomography (only visible on the anisotropic tomography and not shown here) . The difference in both polarizations indicates that the propagation of the SV wave is inhibited when compared to the SH wave propagation.
The magnetometry-gravimetry conceptual model is formed by two components (Figure 4). The first one comprises an anomaly at 5000 m depth below the sea level and was generated with methods from Bott and Smith (1958) and Blakely (1995). The second component suggests an eruptive sequence (Table 1) according to the descriptions of Foshang and González (1956), Corona-Chávez (2002) and Alcantara-Ayala (2010).

Figure 4 2D model of Parícutin volcano area based on magnetic and density susceptibilities (Table 1). a) Architecture of the Parícutin cone, the differences in the local basement (andesites vs granites) and the anomaly at depth confirmed by density and magnetic susceptibilities. Although the geometry of the anomaly at depth is idealized, we suggest it represents a lateral extension of the Tancítaro magma reservoir; b) Close up of the Parícutin cone; c) and d) adjustments of the observed magnetic anomaly and the observed Bouguer anomaly; anomalies were calculated with the GM-SYS program (Geosoft, 2013).
Table 1 Input parameters for the magnetometry-gravimetry 2D model of the Parícutin area. Parícutin events are named after the eruption date (DDMMAA format, e.g., 120145).
Composition and distribution of silicic volcanism in the MGVF
The compositional range of the volcanic rocks found in the MGVF literature is as follows: 40 % andesite, 33 % basaltic andesite, 15 % basalt, 2 % trachybasalt to trachyandesite and 10 % dacite-rhyolite. Although our study might not account for some unpublished data and the composition of most rocks of the MGVF has not been analyzed, the known volume of mafic-intermediate rocks ensures an underwhelming relative percentage of silicic rocks. Silicic magmatism in the MGVF comprises 45 known rock samples that range in composition from dacite (63-71 SiO2 wt.%) to rhyolite (70.5-76 SiO2 wt.%) (Figure 2).
The silicic samples comprise 37 lavas and 8 ignimbrites. We made sure to include only ignimbrites whose vents are, with high probability, located within the MGVF. Silicic lavas occur in three zones: the first zone (Z1) comprises El Tzirade range, El Águila, and La Muela volcanoes (between Pátzcuaro and Morelia), the second zone (Z2) is around Tancítaro stratovolcano, and the third zone (Z3) is located NW of Zacapu (Figure 5). Silicic rocks were formed between 2.9 and 0.003 My (Cardona-Melchor, 2015; Reyes-Guzmán, 2017; Kshirsagar, et al., 2015; Osorio-Ocampo, et al., 2018; Reyes-Guzmán, 2020). We did not find any relation between age, composition, and location. All silicic rocks yield characteristics of typical arc rocks: high LILE-LREE values, negative anomalies of Nb-Ta, and positive anomalies of Pb (Figure 6). All lavas, mostly dacitic, lack Eu anomalies, whereas ignimbrites, mostly rhyolitic, show a negative Eu anomaly.

Figure 5 Digital elevation model of the Michoacán-Guanajuato volcanic field (MGVF) showing the location of silica-rich andesites (black squares), silicic lavas (blue circles) and ignimbrites (red circles). The rocks located between Morelia and Pátzcuaro define the Z1, the rock near Tancítaro defines the Z2 and the rocks NW of Zacapu define the Z3.

Figure 6 Multielement diagrams normalized to mid-ocean ridge basalt (MORB) for silicic rocks, silica-rich andesites and granite xenoliths found in the Michoacán-Guanajuato volcanic field (MGVF). Note the more pronounced Eu anomaly in some granites, whereas in the silicic rocks is less prominent in the silica-rich andesites is absent in most rocks.
Discussion
The petrogenesis of silicic magmas from the MGVF is a pending task. Reported petrogenetic studies are exclusive of mafic-intermediate magmas in the Tancítaro-Jorullo areas (e.g., Ownby et al., 2011; Larrea et al., 2019, Larrea et al., 2021; Gilbaud et al., 2019; Johnson et al., 2008) and of some isolated locations in the NE border of the volcanic field (Losantos et al., 2017). Although the genesis of silicic rocks in the MGVF has been barely explored (Ownby et al., 2011; Pérez-Orozco et al., 2018), the mechanisms that produce dacites and rhyolites are worldwide recognized. Overall, silicic magmas can be created through partial melting of the low-middle crust + AFC processes (e.g., Petrone et al., 2014), melt extraction from a crystal mush (Lipman and Bachmann, 2015; Bachmann and Huber, 2016), remobilization of granites, near or below the solidus by heating events (Wolff and Ramos, 2013), and rejuvenation of partially crystallized granites produced by hydrous mafic injections (e.g., Zou and Ma, 2020).
Given the overall recognition that magmas from central Mexico are generated in a subduction zone (e.g., Gómez-Tuena et al., 2018), on this work we assume that a series of calc-alkaline basaltic magmas evolved by different mechanisms to andesitic melts and reached the medium-upper crust below the MGVF. We then consider silicic rocks in the MGVF as derived from silica-rich intermediate magmas forced to evolve within diverse structural systems. Magmas were trapped in compression zones and partially melted the granitic hosting rocks to form crystal mushes (e.g., Huber et al. , 2011; Cashman et al., 2017; Zou and Ma 2020). Mixing between the intermediate magmas and the granite partial melts formed the silicic rocks.
Tectonic control on the production of silicic compositions
Some monogenetic volcanic fields around the world (e.g., San Francisco, Higashi-Izu, Kaikohe-Bay and Chichinautzin volcanic fields) contain rocks with chemical compositions similar to those of the MGVF; they all contain calc-alkaline rocks, a small portion of alkaline rocks and some have granitic local basements (e.g., Tanaka et al., 1986; Nichols et al., 2012). We found two common characteristics in the examined monogenetic volcanic fields: 1) the volcanic fields are located at the intersection of two or more fault systems, at least one of which has strike-slip movement (e.g., Fenton and Niederman 2014; Nichols et al., 2012.); 2) silicic magmas are created through partial melting of the mantle and AFC processes in the crust (e.g., Tanaka et al., 1986; Schmidt et al., 2016).
Magma movement and accumulation is considered to be controlled by local and regional stress fields, following trajectories normal to the minimum compressive stress, and parallel to the maximum compressive stress (e.g., Tibaldi and Pasquarè, 2008; Bolós et al., 2015; Martí et al., 2016). Orthogonal fault systems can create local compressive zones parallel to the maximum stress where magma can be trapped and forced to evolve. A reorientation of local stress, or the formation of tensional faults in a transtensional system, would allow the ascent of magma and its eruption would proceed, even for silicic degassed magmas.
However, compressional stress is not the only mechanism to trap magmas, stress barriers could contribute, and in fact, might act alone. During decompression, magma could be laterally diverted along regional tectonic structures or high-density bodies produced by old volcanic roots or batholiths (e.g., Menand, 2008; Martí et al., 2016). Lateral magma entrapments could grow into magma chambers if magma-heat supply reach critical conditions; if magma-heat supply is hindered, magma cools and crystallize. If the rocks hosting magma in sills and dikes cannot hold eruptions, due to overpressure, or the local stress field is modified, eruption proceeds (e.g., Martí et al., 2016; Huber et al., 2019).
Both mechanisms, compression and stress barriers, seem to be acting over the MGVF. The intersection of two orthogonal fault systems in the MGVF (ENE-WSW and NNW-SSE) generates extension in the NW-SE direction (Suter et al., 2001; Ego and Asan, 2002). The kinematics of the ENE -WSW faults is inferred as normal and transtensional with a sinistral component (Suter et al., 2001; Suter, 2016; Garduño-Monroy et al., 2009); this same extension has reactivated NNW-SSE faults (Avellán et al., 2020; Gómez-Vasconcelos et al., 2020). These orthogonal fault systems exert control on the NE-SW alignment of volcaninc cones (Garduño-Monroy et al., 2009; Mendoza-Ponce et al., 2018) and on the final composition of silicic rocks where local stress produce compression (Pérez-Orozco et al., 2018).
Local compression zones have been described in the MGVF. Pérez-Orozco et al. (2018) proposed that felsic volcanism is located in transtensional zones of maximum extension (NNE-SSW), where local compression trap magmas. Thus, the question is if other silicic rocks in the MGVF are created by similar processes. Local compression zones are located in the central-northwestern part of the volcanic field, along the intersection of the Morelia-Acambay fault system (ENE-WSW) and the Taxco-San Miguel de Allende fault system (NNW-SSE), and in the southwestern part of MGVF along the intersection of the Chapala-Oaxaca (NW-SE) fault system and NE-SW trend faults. Such compressional zones can be found along the Z1 zone, between Morelia and Pátzcuaro (e.g., Pérez-Orozco et al., 2018). The Z3 zone, near Zacapu, is structurally similar to Z1. The Zacapu area is affected by the Cuitzeo Fault System (ENE - WSW) and the Querétaro-Taxco Fault System (NNW- SSE) (Kshirsagar et al., 2016). The western boundary of the Zacapu graben contains high-SiO2 andesitic lavas, silicic rocks and monogenetic volcanoes dissected by E-W and NW-SW trending faults (Figure 7). Northwest of the Z3 zone occurs the Penjamillo graben, a NNW-oriented regional structure related to the Basin and Range province (Suter et al., 2001).

Figure 7 Digital elevation model of the surroundings of the Z3. Note the NE-SW faults and the N-S Penjamillo graben. Dotted lines show cone alignments most likely following the trace of ENE-WSW and N-S faults. The intersection of these systems of faults might be the stress barrier that trap intermediate magmas and force them to evolve to silicic melts.
We thus suggest that local compressions act over the Z1 and Z3 zones to trap andesitic magmas and make them evolve to dacites-rhyolites (Figure 8). The only rock found in the Z2 might be reflecting the same mechanism that form the silicic rocks in Z1 and Z3, but in this case, associated to the Chapala-Oaxaca regional fault and NE-SW trend faults

Figure 8 Tectonic evolution of the main structural systems in the Michoacán-Guanajuato volcanic field (MGVF). a) NNW-SSE-oriented normal faults of the Taxco-San Miguel de Allende fault system (TSAFS) developed during the Oligocene; b) The Morelia-Acambay fault system (MAFS) developed ENE-SSW to E-W sinistral strike-slip faults during the Oligocene-Miocene and displaced faults of the TSAFS; c) From middle Miocene to present, extension in the MGVF area developed normal faults with a minor left-lateral component in the MAFS and dextral strike-slip faults in the TSAFS. Since the Quaternary, the kinematics is summarized as a left-lateral transtensional system, developing “en echelon” and “pull apart” structures, displacing and reactivating older NNW faults as transfer zones (Suter et al., 1992, 2001; Ego and Ansan, 2002; Garduño-Monroy et al., 2009; Avellán et al., 2020). d) The intersection of fault systems, which produce transfer zones, generated local compression zones in the crust. The regional extension (σ3; minimum compressional stress) caused normal faults and fractures that facilitated the ascent of magma in areas related to “en echelon” structures.
We recognize that compression might not be acting alone, and the batholith below the MGVF could act as a stress barrier forcing the magmas to accumulate. In the following sections, we show evidence of how granitic bodies below the MGVF might influence the evolution of intermediate magmas to produce the silicic rocks.
Depth of generation of silicic magmas
We used gravimetric-magnetometric data to investigate the rocks conforming the local basement of the MGVF. Classic works on Parícutin and Jorullo volcanoes (Wilcox, 1954; Luhr and Carmichael, 1985; McBirney et al., 1987), and the surroundings of Pátzcuaro (Corona-Chávez et al., 2006) reported the occurrence of upper-crustal granitic xenoliths in lavas. Moreover, intrusive bodies (diorites, granodiorites and granites) of Oligocene age crop out in the Jorullo area (Guilbaud et al., 2011). Although granite xenoliths have not been reported in the Zacapu area, quartz xenocrysts are reported in andesitic lavas (Reyes-Gúzman et al., 2018; Ramírez-Uribe et al., 2019). Thus, the overall thought is that the local basement below the MGVF, at least in the Pátzcuaro-Tancítaro-Zapacu-Jorullo areas, is granitic.
Our gravimetric-aeromagnetic data in the surroundings of Parícutin volcano indicate the granitic basement is located 2.7 km below the surface and can be found at least down to 8 km depth (Figure 4). Thus, if magmas are trapped and assimilate the local basement, we suggest that these magmas are hosted by granitic rocks. If magmas get trapped in upper crust compressional zones and evolve to more silicic compositions, their reservoirs must be able to grow, holding the expulsion of magma, until they become eruptable. The depth of such conditions has been indirectly determined for decades through experimental phase diagrams and melt inclusions analyses (e.g., Gardner et al., 1995; Arce et al., 2006; Sosa-Ceballos et al., 2014). Over the decades, 150-250 MPa (~6-10 km) and 4-6 H2O wt.% have been indicated as common pre-eruptive conditions for intermediate-felsic magmas. More recently, Huber et al. (2019) and Townsend and Huber (2020) proposed that chambers can form and grow at 150-250 MPa as a response to neutral bouyancy, rheological contrast and the reorientation of stress. The main idea is that magma recharge deform the upper-crust in a ductile fashion, so that pressure increments are hindered and eruptions are inhibited. Shallower accumulations of magma occur, but they are generally transient, whereas deeper accumulations are the feeder systems of the 150-250 MPa magma chambers. We used our combined gravimetric-magnetic results in the Parícutin area to provide some insights about magma storage below the MGVF.
An anomaly observed and calculated below Parícutin (~7.5 km depth) is proposed from two geophysical surveys (Figure 4). Using two different geophysical methods reduces interpretation ambiguities and enhance a single solution for physical anomalies. Given the composition of the rocks and the monogenetic nature of Parícutin, we do not expect that the anomaly shown in Figure 4 represents a magma chamber that accumulated and fed the magma of the 1943-1952 eruption. Instead, we think this anomaly might represent a lateral extension of the Tancítaro volcano magma chamber, formed by partially melted sills and dikes. Our methods lack the resolution to resolve the geometry of this reservoir. Although this zone could act as a stress barrier and retain a portion of magma tapped during the Parícutin eruption, we must consider that feeding conduits are not necessarily located straight below the vent and they can propagate oblique to the surface. Regardless of its origin, the anomaly strongly suggests that the observed depth (~7.5 km) is prone to accumulate magma within granitic host rocks.
Partial melting of granitic rocks to produce silicic magmas
Binary models and trace elements
The mixing of magmas, or digestion of rocks by hot magma, can be modeled on the basis of binary linear trends that characterize the amount of mixing between two end members (e.g., Sosa-Ceballos et al., 2015). We performed binary models using silica-rich andesites, dacites-rhyolites and granites, the latter found as xenoliths in lavas, to explain the chemical variability of the silicic rocks in the MGVF.
We selected a series of elements that could be enriched in granites and could be traceable in a series of rocks produced by the assimilation or partial melting of granites (major elements, Sc, Zr, Th, U, Rb, Sr, Ba, Cs, La, Lu, Hf, Y, Nb, Zn, Cu). We recognize that fractional crystallization could enrich these elements in a similar fashion. Nevertheless, the trends in Figure 9 indicate that fractional crystallization cannot explain the compositional variability of the silicic and silica-rich andesitic rocks. Instead, partial melting appears to be a dominant factor (Figure 9).

Figure 9 La/Sm vs. La and Rb/Ba vs. Rb diagrams for rocks of the Michoacán-Guanajuato volcanic field. Assuming plagioclase and clinopyroxene as the main crystallized phases, the La-Sm and Rb-Ba trends suggest partial melting dominate over fractional crystallization. Higher La/Sm and Rb/Ba represents low degree partial melts. Silicic rocks more enriched than granites in La/Sm and Rb/Ba could be the result of pure low-degree partial melts; less enriched silicic rocks could be explained by the mixing of intermediate magmas and granitic partial melts. Lavas: blue circles, ignimbrites: red circles, granites: pentagons, silica-rich andesites: squares.
To approximate how partial melting contribute to the generation of silicic magmas, we modeled the chemical variation of silicic rocks (>63 wt.% SiO2) (Supplementary File S2) using as end memebers silica-rich andesites (61-62 wt.% SiO2) and five different granites (data from Luhr and Carmichael, 1985; Cardona-Melchor, 2015; Kshirsagar, et al., 2015; Rasoazanamparany et al., 2016; Reyes-Guzmán, 2017; Osorio-Ocampo et al., 2018; Pérez-Orozco et al. , 2018; Ramírez-Uribe et al., 2019; Reyes-Guzmán, 2020; Avellán et al., 2020); the location of sampling sites are shown in Figure 5. Silica-rich andesites were chosen given their co-occurrence with silicic rocks (Figure 5) and because their composition would produce more silicic magmas with the minimum amount of assimilation or mixing with partial melts of granitic rocks. We recognize that silica-rich andesites could also be a product of assimilation, but this idea coupled to an experimental petrology study will be presented in a separate contribution. We also recognize that the composition of granites does not represent low degree partial melts, hence our models only represent the mixing between silica-rich andesite and granitic melt produced by advanced degrees of partial melting.
We calculated the correlation factor of the linear models for each element as a function of SiO2 content and use it as a proxy of how the two end members interacted (Figure 10). Then, we calculated the relative percentage of elements with correlation factors >0.7 in the models. Our proposal is that a higher percentage of elements with correlation factors >0.7 yields a higher certainty about the interaction of two end members (granites and silica-rich andesites) to produce a series of rocks (Figure 10). Overall, the five granites explain better the compositional variation of rocks from the Pátzcuaro-Morelia area (Z1) than the variation in the Zacapu area (Z3). The lack of correlation in the Zacapu area could be explained if mixing occurs with multiple end members, if the partial melting occurred within a compositionally zoned intrusive body, if the basement is locally altered by hydrothermalism or if samples contained xenoliths different to the granites. In fact, we did not include the ignimbrites in the models because their composition is quite variable, and we think some of them might include lithics not belonging to the magmatic system.

Figure 10 Diagram that describes the percentage of elements (see text for details) that have an average correlation factor (r) >0.7 for binary models with silica-rich andesites and granite xenoliths as end members (see Supplementary Files S1 and S2 for data and models). All rocks were modeled with five different granites. Z1 samples seems to be better explained by mixing of silica-rich andesites and partial melts produced from granites. Low percentage of elements with correlation factors >0.7 could suggests that partial melts are strongly dependent of the mineral assemblage in each granite, hence in order to explain the mixing relations of each silicic rock, the granitic rocks of the corresponding local basement should be sampled, which is impossible to date.
The normalized trace element compositions of silica-rich andesites and granites suggest that mixing between andesites and partial melts produced in a crystal mush contributed to form silicic rocks. Overall, silica-rich andesites show a wide range of REE and HFSE concentration, whereas silicic rocks have less variable trace element abundances and are more enriched in LILE (Figure 6). Granites can be grouped either with silicic rocks or with the silica-rich andesite with the lowest elemental concentrations (Figure 6). Even though each silicic rock could be produced by variable proportions of granite and silica-rich andesite, we focus on Eu anomalies. Eu+2 is well known to substitute Ca+2 in plagioclase and clinopyroxene (e.g., Bindeman and Davis, 2000); all granites and silicic rocks show variable negative Eu anomalies (overall granites have more negative anomalies), whereas silica-rich andesites lack any Eu anomaly. The intermediate anomalies in the silicic rocks could be explained if the silica-rich andesitic melts interacted with granites exerting a higher degree of melting on plagioclase, which would release Eu and result in less developed Eu negative anomalies in the hybrid melt with respect to the granite. Although fractional crystallization of feldspars could produce the same effect, the variability of some incompatible elements suggests fractional crystallization play a minor role (Figure 9). In order to produce ignimbrites with prominent negative Eu anomalies, a granite near to the solidus could be mobilized by more mafic injections. Although this process has been proposed before (Wolff and Ramos, 2013), we do not have enough information to support this hypothesis.
Experimental petrology
We have an ongoing investigation performing hydrothermal-doped experiments to investigate if silicic rocks could be formed when silica-rich andesitic melts assimilate granitic rocks. Preliminary results suggest that the high volume of peritectic crystals produced during the experiment, near the assimilation zone, could generate silicic residual melts, more silica rich than the doped assimilant (Boijseauneau-López, 2018). Moreover, under the same experimental conditions, granitic quartz, feldspar and pyroxene are digested at different rates (Figure 11). Pyroxene and feldspar are digested more rapidly than quartz. The reason for this could be that pyroxene and plagioclase are more reactive than quartz, because of their prominent cleavage, and because of viscosity-diffusion issues. When plagioclase and quartz are heated, the first portion of melt has the effect of stabilizing the crystal at the crystal-boundary interface layer; if this layer is not removed, dissolution is hindered. Plagioclase and pyroxene release Na, Ca, Al, and Si as melting proceeds. This multi-component liquid is similar to the melt surrounding the crystal-boundary layer and diffuses away. The process is facilitated because the diffusivity of Na-K is decoupled from the viscosity of the melt, hence, alkalis can move through the melt without requiring readjustments on the topology of the tetrahedrally coordinated network structure (Mungall, 2002). Contrary to plagio-clase, the melting of quartz produces a highly viscous crystal-boundary interface layer. Quartz is a single component phase, that only releases SiO2 upon melting. Thus, dissolution will stop when the melt at the crystal surface is saturated in SiO2 and can only resume as SiO2 diffuses away. However, the diffusivity of SiO2 is slow because the diffusion of network-modifying cations with intermediate field strength proceeds only with local rearrangements of the melt structure (Mungall, 2002). This mechanism explains why many andesites in the MGVF contains quartz with reaction rims, but xenocrysts-megacrysts of plagioclase and pyroxene are far less common. We think dacites could be created if intrusive rocks are partially melted (plagioclase+pyroxene preferentially melted over quartz) in short periods of time and at relatively low temperatures; instead, rhyolites could be created if the interaction that produced the partial melting is protracted, the temperature is relatively higher, and quartz is effectively melted.

Figure 11 Backscattered electron image of the products of a hydrothermal experiment doped with granite to investigate the role of assimilation on the generation of dacites in the Michoacán-Guanajuato volcanic field (Boijseauneau-López, 2018). Note that plagioclase-feldspar have dissolution zones (sieved texture), pyroxene dissolved extensively and produced a series of peritectic crystals, and quartz seems to not be greatly affected.
Partial melting triggered by frequent-hot magma injections
Despite the partial melting of cold intrusive rocks is not easily achieved at shallow depths (e.g., Dufek and Bergantz, 2005), we think that the intensely deformed upper-crust in the MGVF and the arrival and accumulation of hot magma in the upper crust could promote partial melting of the local basement (e.g., Huber et al., 2011). It is well recognized that local tectonics and the rheology of the upper-crust allow for the accumulation of more magma at depth than that erupted to the surface (e.g., Martí Molist, 2017; Townsend and Huber, 2020), hence protracted interactions between hot magmas and the local basement are enforced.
The seismogenic zone below Michoacán reaches approximately 15 km depth (Rodríguez-Pérez and Zuñiga, 2017), hence, the local granitic basement below the MGFV is deformed by the faults observed at the surface. The upper crust below the MGVF can be considered as a heat sink. The MGVF has one of the largest concentrations of monogenetic volcanoes in the world. Although the volume of magma that cannot reach the surface is unaccountable, its contribution for heating the crust cannot be neglected. The most recent swarm of earthquakes around the Tancítaro-Parícutin area is a good example. During January 2020 the Servicio Sismológico Nacional (National Seismological Service) reported 1080 earthquakes near the Tancítaro volcano, although this is only a preliminary report and many more earthquakes were registered. This is not the only seismic event that has been detected; during 1997, 230 earthquakes were also studied on the Tancítaro area (Pacheco et al., 1999). The hypocenters for the most recent swarm of earthquakes were calculated at 11-14 km depth (M. Perton, personal communication). This depth matches the bottom of the seismogenic zone proposed for Michoacán (Rodríguez-Pérez and Zuñiga, 2017). We thus propose that if magma migrates by hydraulic fracturing, a series of stress barriers, as the bottom of the granitic intrusion, could hinder the propagation of fractures.
If the upper crust is relatively cold, high fluxes of magma result in dyke opening and eruption proceeds (Jellinek and DePaolo, 2003); otherwise, if magma flux rate is low, the time span between magma fluxes is enough to produce a visco-elastic relaxation of the crust, and magma generally reach its solidus (Annen, 2009). If these accumulated magmas are steadily injected or heated by small volumes of magma, magma chambers can develop and grow with a mixture of eruptable melt and crystals.
The unknown volume of magma stalled beneath the MGVF could be a great heat supplier to the upper-crust via convection of fluids or conduction (e.g., Guerrero et al., 2021). Magma eruption rates calculated for the MGVF are 8×10-4 - 2×10-5 km3/year in areas of 15000 and 560 km2, respectively (Hasenaka and Carmichael, 1985a; Osorio-Ocampo et al., 2018). The magma eruption rate (8×10-4 km3/year) is higher compared to the rate needed for the accretion of plutons (1×10-4 km3/year) (Annen, 2009). Hence, we think that 7 My of volcanic activity, and the contribution of unerupted magmas below the MGVF, could have created a hot middle crust prone to deform in a visco-elastic fashion when magmas are trapped for accumulation.
The most voluminous silicic edifice in the MGVF is El Tzirate dome (~3.5 km3). If intermediate magmas were to accumulate and mix with local granite in a proportion of 3:1, the highest magma flux (8×10-4 km3/year) would last 3280 years to produce the silicic magma. Although this incubation time seems feasible, some particular eruptions indicate that considerably less time would be necessary. The Parícutin eruption tapped at least 1.59-1.68 km3 of magma in nine years, whereas other medium-sized volcanoes erupted up to 10 km3 of magma (Siebe et al., 2014), probably in similar periods of time. If only a portion of such high fluxes are captured by stress barriers and local compressive stresses, the maturation of a magma chamber and mobilization of its eruptable portion of magma (e.g., Zou and Ma 2020) would be guaranteed.
Summary: Crustal plumbing system for the generation of silicic magmas in the MGVF
The schematic model on Figure 12 summarizes our proposal for silicic magma generation in the MGVF. We assume that basaltic magmas are created in a subduction zone. After the basaltic magmas leave the mantle wedge, they get stalled in deep crust reservoirs and evolve through crystallization and partial melting of lower crustal rocks (e.g., Ownby et al., 2011). Our tomography results show a normal wave velocity-pressure gradient at depths above 30 km and ~100 km from the trench, (Figure 3). At 30-35 km depth, the tomography suggests a seismic velocity anomaly and at greater depths, where the pressure-wave velocity gradient resumes, the MOHO can be interpreted. We recognize that the 30-35 km low velocity anomaly could be interpreted as a hot solid zone or as partially melted rocks, although a separate study is needed to resolve the ambiguity. Given its location with respect to the trench and the subducted slab, we interpret this anomaly as a swarm of sills, located in the lower crust, whose emplacement triggered the partial melting of crustal rocks to produce the intermediate magmas observed in the MGVF (e.g., Annen et al., 2006). This depth range matches results from previous geophysical studies suggesting that the crust below the MGVF is about 30-40 km thick (Mazzarini et al., 2010). The compositional range observed for the intermediate rocks of the MGVF could be produced by different degrees of partial melting, assimilation of different crustal rocks en route to the surface and/or heterogeneous mantle sources (e.g., Larrea et al., 2021), however, this is beyond the scope of our investigation.

Figure 12 Schematic model of the plumbing system that feeds silicic magmas to the Michoacán-Guanajuato volcanic field. Basaltic magmas created in a subduction zone get stalled in the lower-crust (low shear velocity zone; Figure 3) and evolve through crystallization and assimilation. Magmas reach the seismogenic zone and get diverged into a deformed upper crust to form monogenetic volcanoes. Where fault systems with different orientation interact, a local compression zone is produced and magma accumulate. The accumulated magma releases heat and the granitic hosting rock is partially melted to produce crystal mushes. If local stresses are reoriented in short periods of time and normal faults are reactivated, plagioclase and pyroxene are digested and melts are mixed to produce and tap dacites; if the interaction between the magmas and the crystal mush is protracted, all minerals, including quartz, can be digested and rhyolites are erupted (see text for details).
Intermediate magmas ascend by positive buoyancy (Wanke et al., 2019) and get diverged through an intense network of faults and fractures in the seismogenic zone to produce numerous monogenetic volcanoes at the surface (e.g., Bolós et al., 2015.). The intersection of ENE-WSW and NNW-SSE faults within the MGVF produce zones where compression acts locally and traps magma at 4-8 km depth. At this depth, a previously heated crust and high magma fluxes enhance the development and growth of magma chambers. Hot intermediate magmas get stalled in upper-crustal reservoirs and interact with granitic host rocks. If the interaction is brief, the digestion of granites yields dacitic hybrid melts due to the preferential melting of plagioclase and pyroxene; if stagnation is protracted, magmas interact more effectively with granitic quartz, crystal mushes develop, and the dissolution of all granite mineral phases yields rhyolitic melts.
Supplementary material
Supporting supplementary files S1. “Whole rock data” and S2. “Mixing models” can be found in the html version of this paper at the RMCG webpage <http://www.rmcg.unam.mx>.










nueva página del texto (beta)











