INTRODUCTION

The Colima Volcanic Complex (CVC), Mexico, has been extensively studied in the fields of petrology (^{Luhr and Carmichael, 1980}, ²⁰⁰²; ^{Sarocchi et al., 2011}; ^{Crummy, 2013}; ^{Crummy et al., 2014}; ^{Reubi et al., 2017}; ^{Massaro et al., 2020a}, ^2020b), geochemistry (^{Ostrooumov and Taran, 2016}), volcanology (^{Cortés et al., 2010}; ^{Norini et al., 2010}; ^{Roverato et al., 2011}; ^{Vázquez et al., 2016}; ^{Massaro et al., 2018}, ²⁰¹⁹; ^{Crummy et al., 2019}; ^{Macías et al., 2017}), and geophysics (^{Lopez-Loera and Urrutia-Fuccugauchi, 1998}; ^{Núñez-Cornú et al., 2010}; ^{Sychev et al., 2019}; ^{Zamora-Camacho, 2020}; ^{Núñez-Cornú, et al., 1994}; ^{Spica et al., 2017}, ^{Escudero and Bandy, 2017}, ^{Alvarez and Yutsis, 2015}). The CVC is an active and migrating volcanic system located at the western edge of the Mexican Volcanic Belt (Figure 1), within the nearly NS rift structure in the central part of what is known as the Colima Graben (CG). Its origin is related to a well-known mantle window at the common boundary of the Rivera micro-plate and the northwestern most part of the Cocos plate due to differences in subduction velocity and the oceanic plates’ thermophysical properties (^{Yang et al., 2009}; ^{León-Soto et al., 2009}; ^{Manea et al., 2013}, ^{Corbo-Camargo, 2013}). Such mantle window allowed upwelling hydrated mantle materials and magmatism along permeable paths leading to the formation of the CVC (^{Omori et al., 2002}; ^{Manea and Manea; 2008}; ^{Manea et al, 2013}). The older volcanoes of the CVC, Cántaro and Nevado located to the north, started activity about 1.7 Ma ago (^{Allan et al., 1991}) while the southern Colima volcano (CV) is currently active. Presently, the Colima Volcano is identified by the International Association of Volcanology and Chemistry of the Earth Interior (IAVCEI) as one of the most potentially destructive volcanoes in North America (^{Saucedo et al., 2005}).

Figure 1 (a) The studied area is located in the westernmost part of the Trans-Mexican Volcanic Belt (TMVB), within the easterly rift boundary of the Jalisco Block (JB), shown in (b) with the residual Free Air anomaly of the region (^{Basset y Watts, 2015}), which broadly mimics the terrain rugosity of (c). Black stars show the main cities in the studied area. The networks of magnetotelluric and seismic stations used in this work are located in the central part of the Colima Graben (CG), covering most of the Colima Volcanic Complex (d) comprising the Cántaro (EC), Nevado (NC) and Colima volcanos (VC), shown with yellow triangles. 

The development of crustal imaging methods based on passive seismology and magnetotellurics, have provided new means of exploring upper and lower crustal magmatic systems under volcanic zones by taking advantage of wide spectra measurements and data redundancy provided by both geophysical methods (^{Hill et al., 2009}; ^{Bertrand et al., 2012}; ^{Koulakov et al., 2017}; ^{Escudero, 2022}; ^{Foix et al., 2021}; ^{De Plaen, et al., 2019}; ^{Piña-Varas et al., 2018}; ^{Samrock et al., 2018}; ^{Cordell et al., 2022}; ^{Comeau et al., 2016}; ^{Aizawa et al., 2014}; ^{Arzate et al., 2017}). This development has been accompanied with new findings regarding how melt is stored underneath volcanic areas. New conceptual models of trans-crustal magmatic systems have raised questions about the magmatic systems and the processes that control their chemical and physical evolution (^{Sparks et al., 2019}; ^{Cashman et al., 2017}; ^{Jackson et al., 2018}), which are expected to be confirmed by geophysical data. Ambient seismic noise tomography at a regional scale has recently provided a new insight into the magmatic system of the CVC (^{Spica et al., 2017}). These results support the existence of a low shear-wave velocity zone at 15 km below sea level, interpreted by the authors as the magma reservoir of the CVC. The low velocity anomaly under the CVC depicts a flat lying zone consistent with a predicted deeper magma storage and chemical differentiation mechanisms suggested by ^{Jackson et al. (2018)}.

Although the knowledge of the Colima volcano has been advanced by multiple chemical and volcanologic studies, the relic pathway of magmatic fluids and actual hydrothermal connections to the surface remain poorly understood. In this work, we aimed at imaging the anomalous conductivity zones in the shallow and middle crust under the CVC as well as the more probable subsurface structure. We measured, processed and inverted a set of 21 magnetotelluric (MT) soundings acquired in the area (Figure 1). Our results are based on horizontal and vertical sections of the constrained 3D model results and on the comparison with shear-wave velocity models of the CVC based on ambient seismic noise (^{De Plaen et al., 2022}) and surface structural geology to provide a robust magmatic model system of the shallow and deep structure of the CVC volcanic complex.

MT DATA SET AND THE 3D MODEL GRID

The magnetotelluric data were collected at different times from 2015-2017, covering mostly the central portion of the Colima graben. Twenty of the stations were acquired using the MTU-A 2000 Phoenix MT systems, and one station was acquired with a long period LEMI-423 station. Data were processed following standard procedures (e.g., ^{Simpson and Bahr, 2005}; ^{Wight and Bostick, 1980}) to obtain the ground impedance in the frequency range of 10⁴-10^-1 Hz.

The 3D data inversion of the impedance was performed using the academic version of ModEM code (^{Egbert and Kelbert, 2012}; ^{Kelbert et al., 2014}). The final 3D resistivity model of the CVC was run using the four elements of the impedance amplitude and phase tensors at six frequencies per decade in the range of 10⁴-10^-3 Hz. Details of data field acquisition, data quality, and data processing can be found in ^{Romo-Lozano (2021}), and the two-dimensional analysis and inversion results are found in ^{Romo Lozano and Arzate-Flores (2020a)}.

The model mesh for the 3D inversion was generated using the academic version of the 3DGRID code (Meqbel and Egbert, pers. comm., 2015), which consisted of 79 cells in N-S direction (X), 82 cells in E-W direction (Y), plus 14 padding cells at factor of 1.3 in each direction, including those corresponding to the ocean to the south. The horizontal cell size was 1800 × 1800 m, with a total mesh size covering an area of 1462 (X) km × 1470 (Y) km including the padding cells. The vertical nodes (Z) of the mesh started at a constant cell size of 50 m thickness down to 1.5 km depth to consider topography, after which the thickness increases by a factor of 1.2 (Figure 2).

Figure 2 Model mesh configuration for the 3D inversion made with the software 3DGRID (Meqbel and Egbert, 2015, pers. comm.). The mesh consists of 79 cells in N-S direction (X), 82 cells in E-W (Y), plus 14 padding cells at factor of 1.3 in each direction and vertical nodes (Z) of the mesh started at a constant cell size of 50 m thickness down to 1.5 km depth to include topography. 

The preferred inversion run was made using an initial resistivity value of 100 Ωm and a low covariance (0.2) value for the fields after testing a wide range of both parameters (^{Romo Lozano and Arzate-Flores, 2020b}). The chosen resistivity value is consistent with the background resistivity obtained from the 2D model results (^{Romo Lozano and Arzate-Flores, 2020a}). The inversion runs were done in the computer cluster facilities of the Laboratorio Nacional de Visualización Científica Avanzada (LAVIS), using up to 50 core processors working in parallel. The overall RMS misfit was 2.05 after 64 iterations.

Figure 3a shows the location of the MT stations in the studied area with colored circles where the rms misfit is indicated in the color scale, being orange-red the larger misfit. The average rms misfit shown for individual sites is below 3.0, except for site NV5 located in the NE of the covered area, away from the CVC N-S axis. Figure 3b shows the rms variation as a function of period for individual sites (grey lines) and the average rms misfit of all sites (red line). The larger values are concentrated at periods higher than 10-500 s. Figure 4 shows the frequency range of soundings and the skin depths estimated using different ground average resistivities. The maximum depth range reached at 70 % of the sites is around 20 km. At the rest of the sites the maximum depth is larger.

Figure 3 a) The plan view on the left side of the figure shows the location of the MT stations in the studied area with colored circles (with the color scale to the right of the plan view), which indicate the average rms misfit is below 3.0 after 64 iterations for all the sites except for one site NV5. The red triangle indicates the location of the Colima volcano. b) The graph lines in the right panel shows the rms of individual stations (dashed black lines) and the mean rms of all the stations in red. 

Figure 4 Frequency range of MT stations (a) and skin depths as function of period using different half space resistivities. The maximum depth range reached by the data is around 20 km (b). About a quarter of the soundings penetrated deeper than 15-18 km. 

FULL IMPEDANCE MODEL RESULTS

Model results along five EW vertical sections A to E extracted from the 3D inverted impedance data are shown in Figure 5a. Clear changes in the background resistivity are observed from north to south in the studied area. The northern sections A and B are more resistive (> 1000 Ωm) and compact than the southern sections, where a background resistivity of 100-200 Ωm dominates over confined more resistive zones. Although all sections show an anomalous conductor (< ~10 Ωm), their depths and extension changes along a NS axis. The higher conductor observed (C1 in sections A and B) is a relatively shallow anomaly (< 4 km bmsl), which is mainly distributed under the CVC peaks, and it is apparently confined underneath by the resistive (> 1000 ohm·m) and compact upper crust. Southwards, it becomes more conductive and apparently less dense. If this is to be true, the gravity low (Figure 5c) shown in the north of the studied area is responding to the observed shallower conductor in A and B, whereas the southern low density anomalous zone is responding to a narrower conductor (~10-20 Ωm) influenced by a more permeable crust. The anomalous conductive zones C2 and C3 observed in the resistivity sections C and D, extend to depths of ~15 km bmsl and are consistent with the continuation of the gravity low south the CVC, which in turn is consistent with a narrower gravity low and a shallow and a smaller conductor (see section E in figure 5a).

Figure 5 a) Model results along five EW vertical sections (A to E) extracted from the 3D inverted full impedance data. Their location is shown in (b). The northern sections A and B show a more continuous and resistive (~ 1000 Ωm) upper crust (R) than the southern sections C, D and E. The higher conductor observed (C1 in sections A and B) is a relatively shallow anomaly (< 4 km bmsl), which is mainly distributed under the CVC peaks, and it is confined underneath by a resistive (~ 1000 ohm·m) and compact lithology. Southwards, the upper crust (R) becomes more conductive and altered. The anomalous conductive zones C2 and C3 observed in the resistivity sections C and D, extent from surface to depths of ~15 km bmsl and are consistent with the continuation of the gravity low (^{Alvarez and Yutsis, 2015}) southwards the CVC (c). The southernmost vertical section suggest a thinner crust at the center of the profile, coinciding with the gravity low axis, shown as a dashed white line in (c).  

The horizontal slices extracted from the 3D resistivity model at 5, 10, 15 and 25 km depth bmsl are shown in Figure 6. At 5 km depth (Figure 6a) high conductivity zones (< 10 Ωm) approximately follow the ~NS axis of the CVC, including a confined area north of the Nevado and Colima (red triangles), related with remnants of the extinct El Cántaro volcano located further north. Some of the surface faults of the NE-SW fault system that strikes across the CVC seems to have a structural control down to this depth. However, at deeper depths (Figure 6b and 6c) the northerly fault system (thicker lines in the figure) seems to be preponderant. The higher conductivity zone around the CVC prevails confined or nearly confined by resistive rocks (~ 1000 Ωm) down to depths of at least 15 km bmsl (Figure 6c). The axis of the peanut-shaped conductivity anomaly (<20 Ωm) is parallel to the strike of the main faults deriving from the northern system. The enhanced conductivity zone is apparently surrounded by a resistive crust down to depths of at least 15 km, which resembles a channel that forms part of the plumbing system of the volcano.

Figure 6 Horizontal slices extracted from the 3D inverted ground impedance at 5, 10, 15 and 20 km depth bmsl (a to d). From south to north, red triangles show the location of the Colima, Nevado and extinct El Cántaro volcanoes (CVC). The NE-SW fault system that strikes across the CVC is shown in black lines (except in d), whereas the fault system parallel to Colima and Nevado volcanoes is shown with thick black lines. Stars indicate the location of MT stations. Dashed red line surrounds the enhanced conductivity zones at different depths, suggesting a peanut shaped vertical pathway of conductive fluids to the surface at least along 10 km vertically. 

At 25 km depth (Figure 6d) the background of the studied area becomes widespread conductive to the west (~50 ohm·m) and to the south (< 100 ohm·m) except to the northeast where a resistive zone becomes apparent just NE of the CVC. Within the broad anomalous conductivity zone, it is possible to observe resistivity contrasts although within a short range of variation (10-150 ohm·m), which are perhaps marking structural differences. At this depth and southwest of the active CVC, the strike of the NE-SW fault system (Figure 6d) coincides with the trend of the enhanced conductive lower crust, suggesting some sort of structural control to these depths.

AMBIENT NOISE TOMOGRAPHY GRID

Seismic tomography is a modern tool to determine the seismic velocity structure of volcanic areas, which in combination with electric tomographies from MT measurements become a powerful strategy to study the magmatic systems of volcanos. Using the energy from ambient seismic noise instead of local and regional earthquakes has the advantage that such a tomography is not limited by the spatial and/or temporal distribution of seismicity.

Ambient seismic noise tomography (ANT) is based on the experimental and theoretical validation that the cross-correlation function of a pair of continuous passive recordings can be used to retrieve the Green&apos;s function between two receivers (^{Claerbout, 1968}; Rickett and Claerbout, 1996; ^{Lobkis and Weaver, 2001}; ^{Campillo and Paul, 2003}; ^{Shapiro and Campillo, 2004}; ^{Wapenaar, 2004}). This, in turn, allows inverting group velocity under the assumption that the Green&apos;s function for stations at the surface is dominated by the fundamental mode of the surface waves (e.g., ^{Shapiro et al., 2005}). Ambient noise tomography has been increasingly used in the last decade to image the lithosphere from local to continental scales. Recently, ^{De Plaen et al. (2022}) published results of an ANT study at the scale of the Volcán de Colima edifice. As a result, the structure and geometry of the plumbing system in the upper crust was better defined. Details of the methodology and hypothesis can be found in ^{De Plaen et al. (2022)}.

The data base comes from two sources: the Colima Volcano Deep Seismic Experiment (CODEX), and the Telemetric Seismic Network of Colima (RESCO, Figure 1). RESCO is part of the Centro Universitario de Estudios Vulcanológicos (CUEV) of the University of Colima and is in charge of the seismic monitoring of Volcán de Colima. The ANT study used data recorded between 2012 and 2017 by 10 stations of the network that record in three components at 100 Hz. The CODEX experiment was an IRIS/PASSCAL supported seismic array deployed around Volcan de Colima by the Geophysical Institute at the University of Alaska Fairbanks in collaboration with the Universidad de Colima-Observatorio Vulcanológico. The array recorded continuously from 2006 to 2008 in three components at 100 Hz.

RELATION OF VELOCITY ANOMALY AND THE SHALLOW CONDUCTOR

Figure 7 shows the comparison of resistivity slices extracted from the 3D inversion results (Figure 7a) with slices of the shear velocity anomaly (Figure 7b) extracted from the ambient noise tomography. The electric resistivity structure of the more active part of the CVC shows a NE-SW trend, which has an approximate azimuth of 15°. This apparent electric strike in Figure 7a (white dashed line) clearly differs from the geologic strike given by the azimuth of the Tamazula Fault System (TFS, black dashed line). Both azimuths differ in ~45°, and this difference remains even at deeper depths (Figure 6), which suggests that the Tamazula Fault System does not play a role in fluid conduction in the active CVC. However, regionally the distribution of shear-wave velocity (Vs) lows appears to be bounded by the southern trace of the Tamazula Fault System, at least down to a depth of 4 km (Figure 7b). We deduce that, although both fault systems could be tectonically related, the TFS predates the active and conducting fault system (white dashed line in Figure 7b) and it behaves as a hydrological barrier rather than a zone of fluid concentration. At local scale, the trend of the Vs minima (red dashed line in Figure 7b) approaches the electric strike azimuth, which approximately follows the surface trace of a northerly structure that marks the more conductive and permeable active zone interpreted as the remnant of the magmatic fluids pathway of the latest activity of the CVC.

Figure 7 Electric resistivity slices under the CVC (a) showing a nearly NS trend (white dashed line) of a shallow enhanced conductor. For comparison, the black dashed line shows the regional trend of the NE-SW Tamazula fault system of the region. The dotted red line oulines the enhanced conductor at 2 and 4 km bmsl. The red squares correspond to the area covered by the ambient noise tomography (b). High conductivity and negative Vs anomaly trends are consistent in that the NE-SW fault system may be unrelated to the shallow anomalies (down to 4 km bmsl). High conductivity and negative shear velocity are assumed to lie along magmatic fluids pathways because higher permeability in the ground with zones of flowing fluids would decrease the resistivity and lower the shear velocity at the same time.  

VERTICALLY CONDUCTING PATHWAY

A nearly NS resistivity section was extracted from the 3D model results and compared with an almost coincident shear-wave velocity vertical section model. The location and orientation of the profiles are shown in Figure 8. The seismic and MT station networks are shown as green dots and red stars, respectively. In Figure 8, the resistivity profile is combined with surface faults (sf), which were extrapolated downwards (dashed red lines) following the areas of maximum resistivity gradients suggested by the MT-3D image.

Figure 8 (a) NS resistivity section showing surface crossing faults (sf) extrapolated downwards with dashed red lines, following resistivity interfaces separating an extended anomalous conductivity zone (C1, C2, C3) from a resistive lithology (R). Red triangles show the approximate location of the Colima and Nevado volcanos. Dotted lines over the resistivity section show the anomalous component of the horizontal shear wave velocity (VSH) shown in (c). The location of the seismic wave velocity anomaly profile in shown in (b) with dashed thick line, and the electric resistivity profile with a red line. The seismic and MT network of stations are shown as green dots and red stars respectively. The resistivity profile in (e) shows the same fault structure as in (a), however the dotted lines overlying the enhanced conductivity zone corresponds to the vertical component of the shear wave velocity anomaly (VSV) from (d). 

The contrasting resistivity along the profile suggests a sub horizontal shallow (~2 km bmsl) conductor (C1) under the Nevado and Colima volcanoes that matches well with the low shear- wave velocity anomaly (Figure 8c, 8d). Dotted lines overlying the resistivity section in Figures 8a show the anomalous component of the horizontal shear wave velocity (VSH) that emphasize the horizontal structures, namely, the upper layer of volcanoclastic deposits. The lower limit of the shallow conductor C1 lies on top of a resistive body (R1), apparently dipping southwards. North and south surface faults (sf) crossing the CVC, in the resistivity profile can be extrapolated to deeper depths to follow interfaces of contrasting conductivity zones. These boundaries are interpreted as the continuation of the surface faults at depth, that are confining an enhanced conductivity zone (< 5 ohm·m) (C1, C2). The limiting fault structures extend down to depths of 15 km or more. Around crustal depths of 15-20 km the fault plane appears to connect with a more conductive layer to the south and SW of the Colima volcano. The vertical conductor is limited by a more conductive layer from above, which is associated to the mineralized clay layer sealing the hydrothermal system of the volcanic complex (e.g.^{Aizawa, et al., 2022}). Vertical conducting pathways have been reported associated to magmatic fluids pathways in volcanic active zones from middle to lower crustal depths (e.g., Aisawa et al., 2022; ^{Bedrosian et al., 2018}; ^{Arzate et al., 2017}; ^{Comeau et al., 2016}).

The high conductor changes to become more resistive (~50 ohm·m) at deeper depths (> 10 km), namely zones C3, C4, being anomalous with respect to the yet more resistive surroundings (R1 and R2). The resistivity profile in Figure 8e shows the same fault structure as in Figure 8a; however, the dotted lines overlying the enhanced conductivity zone corresponds to the vertical shear wave velocity anomaly (VSV) from Figure 8d. Red triangles show the approximate location of the craters of the Colima and Nevado volcanoes. Here, the anomalous VSV zones are compatible with a vertical pathway of conductive fluids from deeper depths confined in between two fault zones. The shallow traces of the low velocity anomaly combined with the anomalous high conductivity confined zone support a vertical fluid pathway about 20 km wide along a ~NS direction, which connects with a flat lying moderated conductor associated with the magma source of the CVC.

THE MUSH MAGMA RESERVOIR OF THE COLIMA VOLCANO

Based upon petrochemical measurements on samples from recent eruptive events of the Colima volcano, ^{Hughes et al. (2021}) interpret a partially crystallized-partially melted mush reservoir at the depth range of 12-18 km, or equivalently to 4-6 kbar under the Colima volcano. The chemical and textural information in their study support the interpretation of two magmatic end members within the mush reservoir; an evolved member and a more primitive mafic end-member, which are not individually identified with seismic nor the electric modeling. According to ^{Jackson et al. (2018}), a reactive melt flow in mush reservoirs produces the low-crystallinity, chemically differentiated (silicic) magmas that ascent to form shallower intrusions or erupt to the surface. These magmas can host much older crystals, stored at low and even sub-solidus temperatures, which is consistent with crystal chemistry data (e.g.^{Cooper and Kent, 2014}).

Several studies (e.g.^{Leon-Soto et al., 2009}) suggest the existence of sharp changes in the plate subduction angles that takes place beyond 50 km inland from the coast at ~35 km depth, where the Rivera plate becomes steeper. The fluids released from this mantle window along the Graben de Colima Rift are thought to induce faulting (^{Shelly et al., 2013}) of the overlying crust by lowering the effective normal stress, which justifies the existence of non-volcanic deep tremors in the forearc (^{Brudzinski et al., 2016}). The low shear wave velocity layer under the Colima volcano identified by ^{Spica et al. (2017}) is shown in Figure 9a. This flat lying body is interpreted by the authors as a partial melt reservoir trapped beneath the granitic crust of the Jalisco block, fed with upcoming fluids from the mantle window. From the edges to the center, the Vs velocity decreases ~0.4 km/s, which was associated with an increase in sub-solidus temperature of as much as 700 °C. The

Figure 9 (a) The low shear wave velocity layer under the Colima Volcano identified by ^{Spica et al. (2017)} consists of a flat lying body interpreted by the authors as a partial melt reservoir trapped beneath the granitic crust of the Jalisco block. The anomalous Vs zone is centred at a depth of about ~15 km bmsl, it has an approximate thickness of 10 km, and is extended horizontally in a diameter of at least 40 km. (b) Horizontal slice of the electrical model under the Colima volcano at 20 km depth, showing that the low shear wave velocity layer is also more electrically conductive southwards. Red lines represent the Tamazula fault system and the black lines the fault system of the Colima graben. Red triangles are volcanos and grey lines are the traces of the seismic profiles. The green square in (c) indicate the extent of the MT of stations, and the red lines show the location of the seismic Vs profiles. 

3.4 km/s velocity contour is regarded the anatectic carapace of the partially melted mush reservoir. The low shear wave velocity layer is centered at a depth of about ~15 km bmsl, it has an approximate thickness of 10 km, and is extended horizontally in a diameter of at least 40 km. The horizontal slice of the electrical model under the Colima volcano at 20 km depth (Figure 9b) shows that the low shear wave velocity layer is also more electrically conductive, most probably at the anatectic caparace interface of the evolved magma layer.

Changes in local bulk composition caused by reactive melt flow, rather than large increases in temperature, produce the rapid increase in melt fraction that remobilizes the cold- stored crystals (^{Jackson et al., 2018}). The reactive melt flow can be thought as a bouyant melt percolating upwards through and reacting with older and younger crystals (Figure 10a). Bouyancy of undifferentiated mantle materials is promoted by the crustal discontinuity due to differences in subduction angles and properties of the Cocos and the Rivera plates. However, the upflow of a rich mix of unbonded elements would flow at different speeds during the ascent through the crust. At sub-solidus conditions (35-40 km), elements commonly forming hornblenda and feldespar minerals would start to differentiate. Larger elements forming these minerals, such as K, Ca and Na are segregated by their size as their atomic radius are 220, 180 and 180 pm, respectively. Atoms like Al, Fe, Mg and Si of smaller radii (125, 126, 150, and 110 pm, respectively), continue the upward flow together with smaller bounding elements such as oxygen (60 mp) and even smaller and volatile hydrogen (25 mp) atoms. Still at sub-solidus conditions, but at mid crustal depths silicic melts horizons start to differentiate to form sill intrusions. If the mantle plume contains abundance of H, O and Si, its migration continues upwards until the flow is interrupted by older crystal rich layering interfaces that remains temporary stable hosting high-melt fraction evolved magma (Figure 10b).

Figure 10 (a) A thick magma reservoir of differentiated and undifferentiated sill intrusions (modified from ^{Jackson et al., 2018}) during an active volcanic phase; the upper part of the reservoir comprises transient layers containing mafic melt (on top) and remobilisation of crystal mush (layers below). (b) Schematics of the internal structure of the magmatic system of the Colima volcano sketched from the resistivity structure and the low velocity body of Spica et al. (2017). The structure of the magmatic system of the CVC, formed by sequences of evolved and primitive melts, behave as a capacitor at the interfaces, which drives the ionic migration. In a simplified manner, (c) shows an open circuit sketch that represents the electric Earth. 

In this process, the interfaces of the impermeable granitic upper crust and the low melt fraction mush, and the layer hosting evolved magma and the silicic intrusion layer, build up naturally electrical potentials at the charged interfaces (Figure 10b). The electric charges build up at the interfaces are due to free electrons and high mobility of hydrogen ions (H⁺) flow. It is expected that low viscosity, low density, and low resistivity H₂O-rich melts are a source of large amounts of electric charges available relative to the electrically neutral upper crust (< 15 km). Free electrons, protons, ions, and impurities, facilitate the electric polarization and mass transport through the mush to create charged interfaces that store electric potential in the way an electric capacitor do (Figure 10c). The crystallographic and shape-preferred orientations of minerals (^{Spica et al., 2017}) are expected to contribute to polarize further the reservoir. Theoretically, mush crystals are also polarized piezoelectrically, particularly along crystallographic and shape-preferred orientations of minerals (^{Parkhomenko (1971}), which adds to the polarizing charge of the mush.

There remain several puzzling questions to be answered, such as what is the role of this stored electric energy in the reactivation of melt flow. Other open question are if these large capacitors are isolated from the earth´s electric circuit or are somehow connected to the surface and atmosphere, and how if so. We hypothesize that they must be part of an interconnected regional electrical circuitry as illustrated in Figure 10c. However, it is clear that further research in electric and mechanical anisotropic behavior of the polarized mush is required to advance in the knowledge and the nature of the polarity of melt intrusion sills and layers, and what causes the reactivation of melt flow and what triggers the melt effusion.

CONCLUSIONS

3D inversion of the full impedance of 21 mainly audio magnetotelluric (AMT) soundings in the range of 1000-0.05 Hz confirms the existence of a flat lying conductor at a depth of 2-3 km bmsl under the Colima volcano acting as a natural cap seal of the CVC hydrothermal system. Model results along five EW vertical sections extracted from the 3D inverted impedance show that the northern sector of the CVC upper crust is far more resistive (> 1000 ohm·m) and compact than in the southern region, where it becomes more conductive and fragmented.

The combined results of horizontal and vertical resistivity slices and sections reveals a vertical conductor (10-20 ohm·m) of ellipsoidal shape under the CVC, which extend from about 5 km bmsl down to at least 15 km bmsl. The approximate dimensions of the major and minor axes of the elliptical conduit along the fault are about 20 km and 5 km respectively. The major axis of the ellipsoidal enhanced conducting pathway is parallel to the main NS fault system, which suggests that the fault is the one controlling the magmatic fluids ascent.

A flat lying magma reservoir at a depth under 15 km bmsl is interpreted as a result of this study. However, the variable frequency range of the soundings and limited skin depths under 10 s for most of the sites prevented to establish its thickness. At nearly the same depth a flat-lying low shear-wave velocity layer (< 3.2 km/s) under the CVC from a previous study found an approximate thickness of 10 km for this anomalous body. We interpret the flat lying low viscosity, low density, low resistivity body as a mush reservoir containing H₂O-rich melts. The electrically anisotropic behavior of the mush due to conducting flowing fluids (mainly H₂O) “permeating” upwards is promoted by ionic affinity. The physical evidence suggest that the recent and near-future activity of the CV is closely related to a deep flat magmatic source rather than to a shallow magmatic chamber, and that the triggering mechanism of recurrent magmatic activity of the CV may have an electrical component.