Raphael S M De Plaen

Aurélien Mordret

Raul Arámbula- Mendoza

Dulce Vargas-Bracamontes

Victor Hugo Márquez- Ramírez

Thomas Lecocq

Carlos Ariel Ramírez Vázquez

Miguel González

Raul Arámbula-Mendoza

Dulce Vargas- Bracamontes

Victor Hugo Márquez-Ramírez

: Supervision

Software Carlos

Ariel Ramírez Vázquez

Data Curation Miguel González

The shallow three-dimensional structure of Volcán de Colima revealed by ambient seismic noise tomography

Keywords: ambient noise tomography, volcano, magmatic plumbing system Conceptualization, Investigation, Visualization, Writing -Review & Editing, Formal analysis. Aurélien Mordret: Supervision, Methodology, Software, Investigation. Raul Arámbula-Mendoza: Supervision, Data curation, Writing -original draft

Volcán de Colima is one of the most active volcanoes in North America. Even so, it has a poorly constrained upper crustal structure. Here, we present the highest-resolution three-dimensional shear-wave velocity tomography to date of the volcano.

We measured group velocity dispersion curves of Rayleigh and Love waves extracted from ambient seismic noise recorded on the Colima Volcanic Complex and regionalised them into 2-D velocity maps. Those were then locally inverted using a neighbourhood algorithm to obtain accurate shear-velocity models down to 4 km below sea level. The resulting three-dimensional shear wave velocity model gives us a new insight into the internal structure of the volcano to specifically better understand the shallow magma storage and magmatic plumbing system and, in turn, how it can be placed in a wider geotectonic context.

Our results highlight a network of deeply rooted NE-SW low velocity zones oriented along a local fault system. The southward orientation of this lowvelocity anomaly also roughly aligns with the north-south trend of the volcanoes that compose the quaternary Colima Volcanic Complex and could be associated with the gradual trenchward shift of the magmatic front of the volcanic complex dating from the formation of the Cántaro Volcano. The lowvelocity anomaly overlaps a negative radial anisotropy indicating that magma follows vertically oriented structures, such as interfingered dikes or faults and cracks with a substantial vertical component.

The low velocity anomaly under Volcán de Colima contrasts with a distinct high-velocity anomaly under Nevado de Colima that highlights the difference J o u r n a l P r e -p r o o f Journal Pre-proof between the former active system, filled with solidified dikes and sills, and the current one, associated with a network of fluid-filled dikes.

Introduction

When studying active volcanoes, understanding the mechanisms behind the emplacement and ascent of the magma inside the crust can provide significant insight into the expected behavior of a volcano. The magmatic processes in the volcanic plumbing system directly impact the size and the type of an eruption, motivating the careful determination of the geometry of fluid pathways and magma reservoirs. High resolution imaging of the structure of a volcano offers a way to infer its geological history which, in turn, can be placed in a wider geotectonic context. Volcán de Colima, also referred to as Fuego de Colima, is one of the most active volcanoes in Mexico and in North America with eruptions and sector collapses posing a high risk for more than 500,000 people in surrounding cities [START_REF] Capra | Re-assessing volcanic hazard zonation of Volcán de Colima, México[END_REF][START_REF] Norini | Structure of the Colima Volcanic Complex: Origin and Behaviour of Active Fault Systems in the Edifice[END_REF]. The volcano is located in western Mexico (Figure 1), at the exceptionally complex conjunction of tectonic processes (e.g., [START_REF] Ferrari | The dynamic J o u r n a l P r e -p r o o f Journal Pre-proof history of the Trans-Mexican Volcanic Belt and the Mexico subduction zone[END_REF] .

While it has the highest probability of a VEI ≥ 4 eruption in Mexico within the next 20 years (Mendoza-Rosas and De la [START_REF] Mendoza-Rosas | A statistical method linking geological and historical eruption time series for volcanic hazard estimations: Applications to active polygenetic volcanoes[END_REF] , Volcán de Colima has gone through various styles of activity throughout its eruptive history ranging from lessintense effusive and explosive periods to explosive sub-Plinian to Plinian eruptions. The last two sub-Plinian to Plinian eruptions, both classified as VEI = 4, occurred in 1818 and in 1913. The most recent activity has been characterized by an alternation of multiple vulcanian eruptions, lava flows, dome formations, and dome collapses [START_REF] Reyes-Dávila | Volcán de Colima dome collapse of July, 2015 and associated pyroclastic density currents[END_REF] . However, a large event such as the 1913 eruption would now pose a severe risk for the >300,000 people living in the vicinity of the volcano [START_REF] Saucedo | Eyewitness, stratigraphy, chemistry, and eruptive dynamics of the 1913 Plinian eruption of Volcán de Colima, México[END_REF] .

Fully understanding the structural control on the evolution of the volcanic complex and its activity requires an accurate insight into the relationship between the regional tectonics, the magmatic system and the weak basement. Such information can be extracted by using seismic tomography to analyze the seismic velocity structure of the volcano. Seismic tomography is a powerful tool in volcanic environments that can rely on earthquakes or ambient seismic noise. The crustal structure beneath the Colima volcanic complex (CVC) has been studied and imaged using both source types and different strategies with distinct performances at different depths.

A series of body-wave tomographic studies used teleseisms [START_REF] Yang | Seismic structure beneath the Rivera subduction zone from finite-frequency seismic tomography[END_REF] , regional [START_REF] Ochoa-Chávez | P-Wave Velocity Tomography from Local Earthquakes in Western Mexico[END_REF] , and local earthquakes [START_REF] Sychev | Fault-Associated Magma Conduits Beneath Volcán de Colima Revealed by Seismic Velocity and Attenuation Tomography Studies[END_REF] to image the crust down to the upper mantle. Besides the velocity structure under the CVC, [START_REF] Sychev | Fault-Associated Magma Conduits Beneath Volcán de Colima Revealed by Seismic Velocity and Attenuation Tomography Studies[END_REF] notably estimated the seismic attenuation in the region.

The S-wave velocity structure of the CVC has also been imaged using ambient seismic noise at regional scale by [START_REF] Castellanos | Imaging the Eastern Trans-Mexican Volcanic Belt With Ambient Seismic Noise: Evidence for a Slab Tear[END_REF] and [START_REF] Spica | Anatomy of the Colima volcano magmatic system, Mexico[END_REF] . The main advantage of using ambient noise instead of local and regional earthquakes is that the noise tomography is not limited by the spatial and/or temporal distribution of seismicity. Nevertheless, those regional models offer a low resolution under volcán de Colima. [START_REF] Escudero | Ambient seismic noise tomography of the Colima Volcano Complex[END_REF] worked at local scale instead to produce J o u r n a l P r e -p r o o f ambient noise Rayleigh-wave group velocity maps at 3 and 5s, but did not perform a full depth inversion.

Most seismic tomographies of Volcán de Colima resulted in regional models with poor resolution under the edifice, especially in the upper crust. As a result, the structure and geometry of plumbing system in the upper 10 km remains poorly defined, with some speculation of a shallow magma chamber at ~7km depth that still has to be imaged (e.g., [START_REF] Spica | Anatomy of the Colima volcano magmatic system, Mexico[END_REF].

Here, we address this short-coming and specifically study the upper crust by performing a joint inversion of Rayleigh and Love wave group velocity dispersion curves to assess both the shear wave velocity and radial anisotropy structure under Volcán de Colima at a local scale. J o u r n a l P r e -p r o o f

Geological setting

The volcán de Colima is part of the Colima Volcanic Complex (CVC) consisting of three N-S aligned calc-alkaline andesitic stratovolcanoes, along with a number of monogenetic scoria cones of presumably Pleistocene age [START_REF] Luhr | The Colima Volcanic complex, Mexico -I. Post-caldera andesites from Volcán Colima[END_REF][START_REF] Robin | Eruptive history of the Colima volcanic complex (Mexico)[END_REF] : the now extinct and strongly eroded El Cántaro (3,860 m), the large Nevado de Colima (4,340 m), and Volcán de Colima (3,860 m) (Figure 2). The CVC is located along the Middle America subduction zone and, more notably, near one of the world's most complex convergent margins in the world. The subduction presumably controlled the development of the Trans-Mexican Volcanic Belt (TMVB), a ~1000-kilometer-long continental volcanic arc with a peculiar lack of parallelism with the Middle America trench (Figure 1).

Here the Rivera and the Cocos plates are subducting under the North American plate with different ages, compositions, convergence vectors, and dipping angles while likely being separated by a trench-orthogonal tear starting just north of Volcán de Colima. This gap between the two subducting plates likely played a critical role in the formation of the CVC [START_REF] Ferrari | Slab detachment control on mafic volcanic pulse and mantle heterogeneity in central Mexico[END_REF][START_REF] Yang | Seismic structure beneath the Rivera subduction zone from finite-frequency seismic tomography[END_REF] . Within the study area, the subduction is also responsible for the fragmentation of the Jalisco and the Michoacán blocks from the overriding continental North American plate. Those distinct geological units are separated by extensional corridors called Tepic-Zacoalco rift (TZR) in a NW-SE direction, and Colima rift (CR) in a N-S orientation [START_REF] Ferrari | The dynamic J o u r n a l P r e -p r o o f Journal Pre-proof history of the Trans-Mexican Volcanic Belt and the Mexico subduction zone[END_REF][START_REF] Ferrari | Late Miocene to Quaternary extension at the northern boundary of the Jalisco block, western Mexico: The Tepic-Zacoalco rift revised[END_REF][START_REF] Yang | Seismic structure beneath the Rivera subduction zone from finite-frequency seismic tomography[END_REF] . The TZR and the CR connect to the Chapala-Tula rifts (CTR) to form the Guadalajara triple junction [START_REF] Allan | Geology of the Northern Colima and Zacoalco Grabens, southwest Mexico: Late Cenozoic rifting in the Mexican volcanic belt[END_REF][START_REF] Ferrari | Late Miocene to Quaternary extension at the northern boundary of the Jalisco block, western Mexico: The Tepic-Zacoalco rift revised[END_REF][START_REF] Luhr | Active rifting in southwestern Mexico: Manifestations of an incipient eastward spreading-ridge jump[END_REF] .

Besides those well-known tectonic features, the Tamazula fault (TF) has a more disputed origin and influence that has been described in the literature by two contrasting hypotheses. On the one hand, it been interpreted as a major NE-SWstriking regional fault cutting through the CR and the CVC and extending for more than 160 km from the Michoacán block to the Manzanillo Bay [START_REF] Cortés | Geologic mapping of the Colima volcanic complex (Mexico) and implications for hazard assessment[END_REF][START_REF] Garduño-Monroy | La falla tamazula, límite suroriental del bloque jalisco, y sus relaciones con el complejo volcánico de Colima, México[END_REF][START_REF] Rosas-Elguera | Continental boundaries of the Jalisco block and their influence in the Pliocene-Quaternary kinematics of western Mexico[END_REF] . The system mostly includes normal faults generating a NW-SE-trending extensional activity, cutting through the CVC and potentially at the origin of the many past episodes of collapse with preferred direction toward the S and SE [START_REF] Garduño-Monroy | La falla tamazula, límite suroriental del bloque jalisco, y sus relaciones con el complejo volcánico de Colima, México[END_REF] and toward the SW [START_REF] Cortés | Geologic mapping of the Colima volcanic complex (Mexico) and implications for hazard assessment[END_REF] . The current activity of the TF has been inferred from the alignment of crustal earthquakes in the Manzanillo area and along parallel crustal structures (e.g., [START_REF] Domínguez | The fracturing in volcanic edifice before an eruption: The June-July 1998 high-frequency earthquake swarm at Volcán de Colima, México[END_REF][START_REF] Núñez-Cornú | Stress field estimations for Colima volcano, Mexico, based on seismic data[END_REF][START_REF] Pacheco | The Colima, Mexico, earthquake (Mw 5.3) of 7 March 2000: Seismic activity along the southern Colima rift[END_REF][START_REF] Zobin | Overview of the 1997-2000 activity of Volcán de Colima, México[END_REF] . The intersection of the fault system with the CR is interpreted as a zone of weakened crust which, in turn, allowed the development of the volcanic complex. Hence, the occurrence of magma migration along the TF is reflected in the NE-SW alignment of Volcán de Colima and the Volcancito parasitic cone.

On the other hand, [START_REF] Norini | Structure of the Colima Volcanic Complex: Origin and Behaviour of Active Fault Systems in the Edifice[END_REF][START_REF] Norini | Structural architecture of the Colima Volcanic Complex[END_REF] offer a distinct interpretation and associate several fault scarps previously associated to the TF to a 30 km E-W-trending fault system intersecting the CVC across the eastern and western flanks of Nevado de Colima. This unnamed system is characterized by normal kinematics and N-S extension. Along with the TF, they attribute this local volcano-tectonic structure to J o u r n a l P r e -p r o o f the southward gravitational spreading of the volcanic complex due to the loading of the volcanic edifice over the weak sedimentary fill of the CR gently dipping to the south.

These two distinct explanations on the origin of the TF system eventually impact the overall interpretation of the structural architecture and evolutionary mechanisms of the CVC.

Ambient noise tomography

The ambient noise tomography was implemented in four main steps. The first involved the computation of the cross-correlations functions between pairs of stations for combinations of vertical, radial and transverse components (e.g., ZZ, RR, TT). The frequency dependent group travel times were then measured to later be used to construct the 2D group velocity maps at different periods. Finally, the depth structure was revealed by taking advantage of the regionalized dispersion curves for local 1D shear velocity models in every cell of a grid using a neighborhood algorithm.

Seismic data and computation of cross-correlations

Data

J o u r n a l P r e -p r o o f

This study integrates data from two distinct acquisitions: the Volcán de Colima Deep Seismic Experiment (CODEX), and the Telemetric Seismic Network of Colima (RESCO) (Figure 2).

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. RESCO has stations distributed over the state of Colima and over the volcano. Here we use data recorded between 2012 and 2017 by 10 stations of the network that consist of broadband Guralp CMG-40TD and CMG-6TD that record in three components at 100 Hz. The CODEX experiment was an IRIS/PASSCAL supported seismic array deployed around Volcán de Colima by the Geophysical Institute at the University of Alaska Fairbanks in collaboration with the Universidad de Colima-Observatorio Vulcanológico [START_REF] West | The Colima Deep Seismic Experiment: Imaging the Magmatic Root of Colima Volcano[END_REF]. The array consisted of

Cross-correlation

For all stations, the data is first resampled to 20 Hz and band-pass filtered between 0.01 and 8 Hz. Periods of recordings from the RESCO network that corresponded to eruptions of Volcán de Colima are eliminated. The instrument response for all stations is then removed, along with the mean and the trend.

Earthquakes and local high amplitude spikes that could contaminate the ambient noise wave-field are mitigated by implementing a temporal normalization at 3 times the standard deviation of the daily trace and a spectral whitening between 0.05 and 5 Hz. Station pairs with synchronous recordings are then cross-correlated using MSNoise [START_REF] Lecocq | MSNoise, a Python Package for Monitoring Seismic Velocity Changes Using Ambient Seismic Noise[END_REF] for the vertical-vertical (ZZ), and rotated transversetransverse (TT), and radial-radial (RR). The resulting cross-correlation functions (CCFs) for each pair of stations are then stacked using a phase weighted stacking (PWS), and their symmetric parts are retained by averaging their positive and negative lag times.

The PWS has the advantage of enabling the detection of weak but coherent arrivals by exploiting the phase coherence in individual causal and a-causal correlograms and therefore significantly improves the signal to noise ratio in the stacks [START_REF] Schimmel | Phase cross-correlations: Design, comparisons, and applications[END_REF][START_REF] Schimmel | Frequency-dependent phase coherence for noise suppression in seismic array data[END_REF][START_REF] Schimmel | Noise reduction and detection of weak, coherent signals through phase-weighted stacks[END_REF] . The whole process resulted in 311 CCFs for each component, respectively.

Dispersion measurement

The group velocity dispersion curves are estimated using a frequency-time analysis (FTAN) from 0.5 to 10 s [START_REF] Levshin | On a frequency-time analysis of oscillations[END_REF][START_REF] Ritzwoller | Eurasian surface wave tomography: Group velocities[END_REF] . To prevent the wrong automatic picking of the dispersion curve caused by highamplitude surface wave overtones or scattered waves, we use a Graphical Users Interface that involves analyst validation (Mordret et al., 2014a) . Group velocities from inter-station distances smaller than 1.5 wavelengths and SNR<5 are not considered. The analysis is limited to reasonably covered period ranges with at least 50 measurements, which includes periods between 0.5 and 10 seconds for Rayleigh waves and 0.8 to 8.4 seconds for Love waves.

J o u r n a l P r e -p r o o f

Journal Pre-proof

Rayleigh and Love Wave Group Velocity Maps

The tomographic inversion of the group-velocity measurements is performed using the Cartesian version of algorithms described by [START_REF] Barmin | A fast and reliable method for surface wave tomography[END_REF] and implemented by [START_REF] Mordret | Ambient noise surface wave tomography to determine the shallow shear velocity structure at Valhall: Depth inversion with a Neighbourhood Algorithm[END_REF] . The algorithms are based on ray theory and involve a regularization function composed of a spatial Gaussian smoothing function along with a constraint on the amplitude of the perturbation based on local path density. We use a grid composed of 54× 53 = 2862 square cells with a grid size of 1 km. The inversion is performed in two steps, with an initial model that has a constant velocity taken as the mean group velocity for each period (Mordret et al., 2014a;[START_REF] Moschetti | Surface wave tomography of the western United States from ambient seismic noise: Rayleigh wave group velocity maps[END_REF] . The first inversion produces a very smooth map used to identify and reject measurements with travel time residuals greater than two standard deviations. The remaining measurements are used in the second inversion to produce the final group velocity maps. The procedure does not account for the topography. The error of the flat topography approximation on the velocity was estimated by [START_REF] Brenguier | 3-D surface wave tomography of the Piton de la Fournaise volcano using seismic noise correlations[END_REF] and Mordret et al. (2014b) to be less than 5%. Such an error should therefore be negligible compared to the group velocity variations of over 30% at each period.

J o u r n a l P r e -p r o o f

The group velocity measurements and coverage along with the path density for the respective periods are shown for Rayleigh and Love waves in Figure 4 and5 (d-f andg-i), respectively. Both Volcán de Colima and the Nevado de Colima are covered by a dense coverage of at least 8 ray paths for all 3 periods. The group-velocity maps for Rayleigh and Love waves, where lateral velocity variations can be observed at periods of 2, 4 and 6 s, are shown in Figure 4 and 5 (a-c), respectively. The locations shown and interpreted have at least four rays per cell. A variance reduction of the travel-time residuals exhibits values ranging from 72 to 79% and 88 to 93% for Rayleigh and Love waves respectively, indicating that the retrieved velocity maps fit the data well. The mean velocity of the maps increases slightly with the period for both Rayleigh and Love waves, going from 1.79 to 1.94 km/s and 1.29 to 1.63 km/s respectively, in good agreement with the dispersion curves shown in Figure 3. At all periods, we observe an elongated low velocity zone that extends from the under Volcán de Colima to the SW direction (Figure 4 and5 a-c). The Love wave group velocity maps display a similar global velocity structure pattern as for Rayleigh waves, with the exception of small low velocity zones.

J o u r n a l P r e -p r o o f

The accuracy of the size and location of the retrieved anomalies is assessed through the evaluation of the resolution matrix associated with each inversion [START_REF] Barmin | A fast and reliable method for surface wave tomography[END_REF]Mordret et al., 2014b[START_REF] Mordret | Ambient noise surface wave tomography to determine the shallow shear velocity structure at Valhall: Depth inversion with a Neighbourhood Algorithm[END_REF] . For each cell of the model, the corresponding line of the resolution matrix represents the response of the tomographic process to a Delta function type anomaly with the size, the shape, and the location of the output spot indicating how accurately the anomaly can be retrieved. The precision of the location of the recovered spike varies slightly depending on the station coverage. The entire resolution matrix at 2, 4, and 6 s for Rayleigh and Love waves are shown in Figure 6 and Figure 7, respectively. Here, with a grid cell size of 1 km, a resolution shift <0.5 km means that the center of the retrieved anomaly is in the same cell as the input spike.

J o u r n a l P r e -p r o o f

Journal Pre-proof

Depth inversion

Rayleigh and Love waves group velocity maps are like a gather of local group velocity dispersion curves with group velocity measurements at different periods for every particular location on the map.

A 1-D local layered velocity model can then be assessed by inverting these local dispersion curves for each cell of the grid. For this inversion we use a Monte-Carlo global direct-search technique, the Neighborhood Algorithm developed by [START_REF] Sambridge | Geophysical inversion with a neighbourhood algorithm-I. Searching a parameter space[END_REF] and described in detail in [START_REF] Mordret | Ambient noise surface wave tomography to determine the shallow shear velocity structure at Valhall: Depth inversion with a Neighbourhood Algorithm[END_REF]. The solution space is sequentially and non-uniformly sampled while considering the data fit from old samples.

The parametrization here is similar to Mordret et al. (2014b) and uses a power law backbone velocity profile to take into account the possible slow velocities in the shallow layers due to weathering and alteration overlaid with smoothly varying velocity anomalies and anisotropy. Along with two parameters controlling the power law backbone velocity profile (the surface velocity and the curvature of the velocity profile, respectively) we use five cubic B splines that parametrize the overlaid velocity anomalies. To reduce the number of parameters to be inverted, these continuous functions are discretized into 33 layers. For each layer, we consider a constant shear-wave velocity (Vs) and the parameter ranges allowed are generally J o u r n a l P r e -p r o o f defined based on the plausible crustal models in a volcanic region, allowing a large variability to resolve both low and high velocity body.

For the initial isotropic inversion, a total 25,000 Vsv and Vsh are sampled for each cell of the grid, and the best-fitting 300 models are kept to calculate a mean model resulting in a 3D Vs model of the volcano. The anisotropic and isotropic models are compared to evaluate the necessity of an anisotropic inversion. The anisotropic inversion was parametrized using an additional four dedicated cubic B splines for a complete parametrization of 13 parameters. In this case, a total 19,000 models are samples for each cell of the grid.

The Vs isotropic model or effective S wave velocity at each location is computed from Vsh and Vsv via a Voigt's average which better reflects the isotropic variations in elastic properties [START_REF] Dziewonski | Preliminary reference Earth model[END_REF][START_REF] Ekström | The unique anisotropy of the Pacific upper mantle[END_REF] :

𝑉𝑠 = √ 2𝑉𝑠𝑣² + 𝑉𝑠ℎ² 3
The average model is removed from the resulting 3D Vs model to obtain the 3D Vs anomaly model presented and interpreted below. The difference between Vsv and Vsh, called radial anisotropy (𝜉), is represented as the percentage difference between Vsh and Vsv in the medium [START_REF] Jaxybulatov | A large magmatic sill complex beneath the Toba caldera[END_REF]Mordret et al., 2014b) :

𝜉 = 𝑉𝑠ℎ -𝑉𝑠𝑣 𝑉𝑠

Figure 8 shows the best fit to isotropic and anisotropic dispersion curves at two different locations. On average over the whole studied area, anisotropic parameterization reduces the misfit by ∼ 40%. The isotropic inversion has particularly poor results below ~4 s period, after which a better fit can be observed between the measured and the synthetic dispersion curves. This indicated the need for an anisotropic inversion with the additional parametrization focused on the upper 5 km depth. This could indicate that the radial anisotropy is limited to shallow depths, although this is probably caused by a lack of dispersion measurements for Love waves at longer periods combined with a shallower depth sensitivity compared to Rayleigh wave at corresponding periods (Figure 9).

J o u r n a l P r e -p r o o f

J o u r n a l P r e -p r o o f

Journal Pre-proof

Interpretation and Discussion

We built a high-resolution shear-wave velocity model of Volcán de Colima down to 6 km b.s.l. The model includes the Nevado de Colima, Volcán de Colima and a good portion of its southern flank. The Cántaro is the only part of the CVC that is notable missing from our observations. Nevertheless, this study includes most structures of the CVC observed in other geophysical and geological studies, with a finer resolution than past ambient seismic noise tomographies. .

J o u r n a l P r e -p r o o f

Journal Pre-proof

Shear wave velocity and anomaly

The group velocity maps show a consistent low velocity zone trending from NE to SW with very little lateral variation from 2s to 6s (Figures 4 and5 a-c). This anomaly has been similarly observed by the local-scale ambient noise tomography by [START_REF] Escudero | Ambient seismic noise tomography of the Colima Volcano Complex[END_REF], who also used CODEX data to retrieve Rayleigh wave velocity maps of the CVC at 3 and 5s. Notably, their mean velocity of 1.8 km/s for Rayleigh waves is relatively close to the ones observed in this study ranging from 1.79 to 1.94 between 2 and 6 s.

The ambient seismic noise based regional Vs model by [START_REF] Spica | Anatomy of the Colima volcano magmatic system, Mexico[END_REF] displays a similar low velocity anomaly at 4 km depth under Volcán de Colima. In their case, the anomaly appears much larger, which they acknowledge as a possible consequence of the low horizontal resolution of their regional tomography, and that, in turn, could cause spatial smearing when trying to image a body of this size. The low velocity in our 3D Vs model has a distinct NE-SW orientation and is limited at its northern border by the Nevado de Colima (Figure 11). A major portion of the Nevado de Colima is associated with a high-velocity anomaly that extends to the north and at depth.

High velocity zones on volcanoes are not uncommon and are typically associated with solidified dikes or sills that formed during past eruption history of the volcano. The Piton the la Fournaise volcano has such magma path intrusions at depth that were revealed as high velocity anomaly by ambient noise tomography (Mordret et al., 2014b) , gravimetry [START_REF] Gailler | Gravity structure of Piton de la Fournaise volcano and inferred mass transfer during the 2007 crisis[END_REF] , and drilling [START_REF] Lerebour | The Grand Brûlé exploration drilling: New data on the deep framework of the Piton de la Fournaise volcano. Part 2: Secondary minerals[END_REF] . Similar high velocity anomalies have been revealed by local earthquake tomography beneath the Popocatepetl [START_REF] Kuznetsov | The three-dimensional structure beneath the Popocatépetl volcano (Mexico) based on local earthquake seismic tomography[END_REF] and the Asama Volcano [START_REF] Aoki | Pwave velocity structure beneath Asama Volcano, Japan, inferred from active source seismic experiment[END_REF] .

The southward orientation of the low-velocity anomaly is roughly consistent with the north-south trend of the chain of three andesitic stratovolcanoes that compose the quaternary CVC. This configuration has caused speculations about a gradual trenchward shift of the magmatic front of the CVC in the frame of a long-lasting southward migration of the magma pathways dating from the formation of the Cántaro Volcano (e.g., [START_REF] Alvarez | Southward Migration of Magmatic Activity in the Colima Volcanic Complex, Mexico: An Ongoing Process[END_REF][START_REF] Luhr | The Colima Volcanic complex, Mexico -I. Post-caldera andesites from Volcán Colima[END_REF].

The CVC had at least 12 gravitational collapses during the past 45,000 years, nine of which occurred at volcán de Colima during the past 22,000 years [START_REF] Borselli | Flank collapse scenarios at Volcán de Colima, Mexico: A relative instability analysis[END_REF]. [START_REF] Roverato | Stratigraphic reconstruction of two debris avalanche deposits at Colima Volcano (Mexico): Insights into pre-failure conditions and climate influence[END_REF] attribute such a collapse rate to the continuous magmatic recharge in the volcanic system, the regional and the local stress regime by the Colima Graben, and the Tamazula Fault whose extensional NW-SE vector has promoted sector collapses of the CVC to the SW. The location of the low-velocity anomaly to the south of the Nevado de Colima is consistent with the multiple sector collapses of both Nevado de Colima and Paleofuego volcanoes, although it is too deep to be explained by this element alone (e.g., [START_REF] Sychev | Fault-Associated Magma Conduits Beneath Volcán de Colima Revealed by Seismic Velocity and Attenuation Tomography Studies[END_REF] . Such collapses generate stress perturbations due to unloading caused by mass redistribution. At some volcanoes, the lateral flank failures were followed by a modification of the shallow plumbing system (Maccaferri et al. 2017). The low-velocity anomaly could contain a magma reservoir that cannot be easily distinguished from shallow volcanic deposits J o u r n a l P r e -p r o o f and soft fluid saturated rocks. The volcano tectonic hypocenters likely provide evidence that the low velocity zone below Volcán de Colima hosts the magmatic plumbing system feeding the volcano. The location of the hypothetical shallow magma reservoir to the south of the main vent has led to speculations that the parasitic domes of Los Hijos, on the southern flanks of Volcán de Colima, feeds from the same source, further indicating a southward migration of the volcanism of the CVC [START_REF] Alvarez | Southward Migration of Magmatic Activity in the Colima Volcanic Complex, Mexico: An Ongoing Process[END_REF][START_REF] Sychev | Fault-Associated Magma Conduits Beneath Volcán de Colima Revealed by Seismic Velocity and Attenuation Tomography Studies[END_REF].

The presence of Los Hijos, along with the Volcancito parasitic cone on the northeastern flank of Volcán de Colima, led [START_REF] Garduño-Monroy | La falla tamazula, límite suroriental del bloque jalisco, y sus relaciones con el complejo volcánico de Colima, México[END_REF] to associate the aforementioned potential migration of volcanism with the NE-SW Tamazula Fault (TF) system. This system, mostly composed of normal faults, is believed to have weakened the crust and serve as a path for magma migration and, along with the associated stress field, is responsible for the flank instability of the CVC [START_REF] Cortés | Geologic mapping of the Colima volcanic complex (Mexico) and implications for hazard assessment[END_REF][START_REF] Garduño-Monroy | La falla tamazula, límite suroriental del bloque jalisco, y sus relaciones con el complejo volcánico de Colima, México[END_REF] . In contrast, some of the features attributed to the TF have been also interpreted as the result of active gravitational spreading of the volcanic complex by [START_REF] Norini | Structure of the Colima Volcanic Complex: Origin and Behaviour of Active Fault Systems in the Edifice[END_REF][START_REF] Norini | Structural architecture of the Colima Volcanic Complex[END_REF] . With the NE-SW orientation of the low velocity anomaly in our 3D Vs model parallel to the TF system (Figure 11) we can speculate that this is more in agreement with a regional whole-crustal structure that facilitates magma migration beneath CVC. [START_REF] Sychev | Fault-Associated Magma Conduits Beneath Volcán de Colima Revealed by Seismic Velocity and Attenuation Tomography Studies[END_REF] built 3D models of the P and S wave velocities and attenuation for the CVC from the body waves of local earthquakes recorded by CODEX. They observed the same distinction between a high velocity zone under the Nevado de Colima and a low velocity anomaly to the south of Volcán de Colima. By extrapolating it at depth as a south dipping structure that reaches the lower crust, they interpret the TF as a mechanically weakened zone that forms a pathway to feed CVC after a transition at ~15 km depth.

The alignment of our low velocity zone along the TF does not unequivocally prove that the fault system is the deeply rooted structure guiding magma to the surface. We therefore cannot exclude the interpretation by [START_REF] Norini | Structure of the Colima Volcanic Complex: Origin and Behaviour of Active Fault Systems in the Edifice[END_REF][START_REF] Norini | Structural architecture of the Colima Volcanic Complex[END_REF]) that the fault system has been originated by gravitational spreading rather than by regional active tectonics. Nevertheless, there is increasing geophysical evidence showing a relationship between the TF and a regional whole-crustal structure that contributes to magma migration towards the surface.

The shallow extensional activity of the segment of TF within the Colima Graben could be the result of a deeper regional active tectonics combined with topographic configuration of the volcanic complex, the weakness of the underlying sedimentary graben fill, and of the geometry of the regional basement. Such a concept would also help explain the lack of seismicity specifically associated with this segment of the fault system, as opposed to the crustal earthquakes observed in its Manzanillo section.

Unlike [START_REF] Sychev | Fault-Associated Magma Conduits Beneath Volcán de Colima Revealed by Seismic Velocity and Attenuation Tomography Studies[END_REF], our model does not reach 15 km depth. However, the low velocity anomaly appears as a consolidated structure down to 6 km b.s.l. that is ~15 km wide in the B-B'' profile and ~10 in the C-C'' profile. The A-A'' profile, however, displays two distinct vertical structures that connect at ~2 km b.s.l (L1 and L2, Figure 12a). Previous estimates place a shallow magma chamber at 6.3 km below J o u r n a l P r e -p r o o f the summit based on the depth of the hornblende stability region for the 1998-1999 erupted lavas [START_REF] Luhr | Petrology and geochemistry of the 1991 and 1998-1999 lava flows from Volcán de Colima, México: Implications for the end of the current eruptive cycle[END_REF] and the top at about 1.5 km depth below sea level based on a low-density body in the gravimetric model of [START_REF] Medina | A gravity model of the Colima, Mexico region[END_REF].

The observation of an aseismic zone between 4 and 7 km depth during the 1991 and 1997-2000 activity of the volcano is believed to further indicate the existence of this shallow magma reservoir below the summit [START_REF] Núñez-Cornú | Seismic activity related to the 1991 eruption of Colima Volcano, Mexico[END_REF][START_REF] Zobin | Overview of the 1997-2000 activity of Volcán de Colima, México[END_REF] . All those estimates place a shallow magma chamber at a depth where our model provides reliable observation of a low velocity anomaly.

Radial anisotropy

The combination of the shallower maximum sensitivity depth of Love waves compared to Rayleigh waves and the lack of long period measurements for Love wave dispersion curves caused the radial anisotropy to be apparently limited to above 5 km depth in our results.

Nevertheless, the association of a fault system with the magma transport under the CVC is supported by the radial anisotropy results. We focus this part on the interpretation on the sign of the radial anisotropy because the simplifications in the inversion procedure, ignoring the topography and the P wave anisotropy, likely have biased its absolute value. In a similar context, however, [START_REF] Xie | Crustal radial anisotropy across eastern Tibet and the Western Yangtze Craton[END_REF] and Mordret et al. (2014b) showed that the sign should be robust.

In a volcanic environment and at crustal depths, radial anisotropy can be a very good indicator of the dominant fabric within the magmatic system and, consequently, of the dominant type of intrusions. Partial melting can produce radial anisotropy with oriented volcanic fluid-filled cracks strongly affecting the macro-mechanical properties of the material (e.g., [START_REF] Gudmundsson | Deflection of dykes into sills at discontinuities and magma-chamber formation[END_REF] , with common interpretations pointing towards the presence of sill layering, lava flows, dike networks or lenses of partial melt. For example, sill layering tends to produce positive radial anisotropy because Vsv is globally slowed down due to the presence of molten material piled horizontally [START_REF] Jaxybulatov | A large magmatic sill complex beneath the Toba caldera[END_REF] . Alternatively, active dike networks would yield a negative anisotropy because their vertically oriented structure would slow down Vsh (e.g., [START_REF] Lee | Upper crustal shear wave velocity and radial anisotropy beneath Jeju Island volcanoes from ambient noise tomography[END_REF]Mordret et al., 2014b) .

A negative radial anisotropy dominates the edifice above 5 km depth indicating dominant vertical structures such as steep normal faults, jointing, and magmatic dike intrusions. Some of these structures are exposed at the surface, such as the three subvertical N-S-trending dikes crop out identified by [START_REF] Norini | Structural architecture of the Colima Volcanic Complex[END_REF] in the inner wall of the Paleofuego collapse scar. The N-S orientation of these dikes parallels the alignment of the Nevado de Colima, Volcán de Colima, and Los Hijos domes and are interpreted as evidence that the magma storage of the active CVC is drained by a network of N-S-trending dikes.

Magmatic fluids are typically expected to accumulate in shallow magma-chambers in horizontally extensive, overlapping sill-like structures (e.g., [START_REF] Jaxybulatov | A large magmatic sill complex beneath the Toba caldera[END_REF] J o u r n a l P r e -p r o o f or horizontal lenses [START_REF] Fournier | Hydrothermal processes related to movement of fluid from plastic into brittle rock in the magmatic-epithermal environment[END_REF] . Nevertheless, the alignment of the low velocity anomaly in our model with the TF suggests a development of magma ascent along fault lines and shear zones that is further confirmed by the negative radial anisotropy. This association indicates that faults and vertical magma-filled dikes are the likely dominant structure in the edifice. Furthermore, the negative radial anisotropy identified by [START_REF] Spica | Anatomy of the Colima volcano magmatic system, Mexico[END_REF] under the CVC seems to extend deeper through the crust, giving additional evidence of vertically extensive fluid-filled fractures in the crust.

A preexisting fault zone such as the CR included potentially deep reaching structures. This provides a natural pathway for fluids, with the high-permeability damage zone and the low-permeability fault core, to allow magma to ascent to the surface (e.g., [START_REF] Shelly | A fluid-driven earthquake swarm on the margin of the Yellowstone caldera[END_REF] .

Another interesting feature comes from the concurrent modeling of the gravimetric and magnetic anomalies by [START_REF] Alvarez | Southward Migration of Magmatic Activity in the Colima Volcanic Complex, Mexico: An Ongoing Process[END_REF] . They identified two low density/ high magnetic susceptibility anomalies, interpreted as magma bodies, up to ~1.5 km b.s.l. that correspond with the first and the second negative anomaly along the A-A'' profile of our model (L1 and L2 respectively, Figure 12a). Between them they identify a higher density anomaly at a depth of ~1.2 km b.s.l. below the Yerbabuena depression that corresponds with the middle vertical anomaly along the A-A'' profile of our model. They interpret this structure as a now solidified intrusion that induced a local gravity collapse following its cooling. These features are additional evidence pointing towards a southward migration of the volcanism in the CVC, although they don't provide clear insights on the mechanism behind that change. [START_REF] Alvarez | Southward Migration of Magmatic Activity in the Colima Volcanic Complex, Mexico: An Ongoing Process[END_REF] identified another low-density region south of the Yerbabuena depression that they called the "Southern Magma Body" and interpreted as molten material. Although our results cannot confirm the presence of magma in that area of the model, that low density body corresponds to the southernmost vertical low velocity anomaly along the A-A'' profile (L2, Figure 12a) and should be further investigated.

In the absence of supporting evidence for a magma reservoir at such shallow depths, the hydrothermal alteration could offer an alternative explanation for the observed velocity anomaly. Low-velocity anomalies in the shallow crust have been found on other volcanoes such as Uturuncu volcano in the Andes [START_REF] Jay | Shallow seismicity, triggered seismicity, and ambient noise tomography at the long-dormant Uturuncu Volcano, Bolivia[END_REF], and Mount St-Helens [START_REF] Wang | Ambient noise tomography across Mount St. Helens using a dense seismic array[END_REF] where they interpreted the low-velocity anomaly as a highly fractured volume containing the magma conduit and hydrothermal systems, but with most of the pore space in the volume likely filled with hydrothermal fluids rather than melt. [START_REF] Heap | Hydrothermal alteration can result in pore pressurization and volcano instability[END_REF] looked at laboratory data for altered rhyodacites from Chaos Crags (Lassen volcanic center) to ultimately emphasize that both the porositydecreasing alteration they explored and porosity-increasing alteration can promote volcano instability and collapse, although by different mechanisms. This, in turn, could increase the likelihood of volcano spreading, flank collapses, and phreatic/phreatomagmatic explosions. This could be linked with the elevated collapse rate of the CVC with the progressive weakening of the volcanic edifice by hydrothermal alteration as a main indirect endogenous factor [START_REF] Roverato | Stratigraphic reconstruction of two debris avalanche deposits at Colima Volcano (Mexico): Insights into pre-failure conditions and climate influence[END_REF].

J o u r n a l P r e -p r o o f

The location of the low velocity zones relative to the TF establishes a strong apparent relationship between the fault system and intrusive bodies that use such a natural path of weakness to ascend through networks of interfingered dikes. Basement fault activity beneath the volcanic edifice could have contributed to both the formation of the TF and the frequent flank collapses at the CVC (e.g., Vidal and Merle, 2000). A more comprehensive study of the seismicity using the decades of existing seismic data would help to better investigate the tectonic structure in the CVC and how they relate to active magma pathways in the crust.

Because the radial anisotropy alone does not fully support a fault-assisted magma migration, further attention should be paid to possible processes behind it. Theoretical and field studies have overwhelmingly indicated that magma pathways are determined by host rock stresses, with magma moving along the path of minimum energy within a volcanic edifice (e.g., Rivalta et al., 2019). Magma ascends to the surface through dikes that align themselves with the most energy-efficient orientation, roughly perpendicular to the least compressive principal stress axis σ3 .

Along with structural features, stress is attributed to generating the anisotropy observed in the crust (e.g., Johnson et al., 2011). In turn, both the local stress field and inherent structures (topography, fault, collapsed caldera structure, etc.) also impact the development of complex magma plumbing systems (Chaussard and Amelung, 2014;Tibaldi, 2015). This allows for distinguishing between regions of stress-induced and structural anisotropy around volcanoes (e.g., Johnson et al., 2011).

A better understanding of the stress conditions in the crust of the CVC could help put in perspective their impact on volcanic structures and processes, the observed anisotropy, and potentially the TF. Massaro et al ( 2020) analyzed the stress-filled of the CVC by integrating published geophysical, volcanological, and petrological data into a 2D finite-element model. Their modeling indicates that the majority of stress variations are located at the tips of the magma chambers, implying a close-toequilibrium state of the volcano associated with a stress distribution induced by a feeding system directly connected to the surface. This is in agreement with the longlasting open conduit dynamics that have lasted since 1913. Nevertheless, the lowvelocity anomaly could come in part from the contribution of a highly fractured region beneath the volcanic edifice generated by stress reorganization above the magma storage, providing magma pathways to the surface (e.g., [START_REF] Wang | Ambient noise tomography across Mount St. Helens using a dense seismic array[END_REF] .

When it comes to the negative radial anisotropy, [START_REF] Spica | Anatomy of the Colima volcano magmatic system, Mexico[END_REF] emphasized that the shear strains extensional contexts preferentially orient the seismic slow axes of the minerals along the vertical axis which is expected to induce a positive radial anisotropy (e.g., Mainprice et al., 2000) . The negative radial anisotropy is, therefore, evidence that the shape-preferred orientation of vertical structures in the crust, likely guiding fluids toward the surface, dominates over the crystallographic preferred orientation at the CVC.

Insights from petrology

Petrological and geochemical studies focused on the CVC suggest the existence of a complex underlying magmatic plumbing system, involving multiple magma storage J o u r n a l P r e -p r o o f levels from the base of the crust to a shallow-level upper crustal magma storage region, with deposits such as those of the major 1913 Plinian eruption of volcán de Colima revealing a particular process of magma mixing or mingling (Crummy et al., 2014) . Pressure estimates from two-pyroxene barometry by Hughes et al. (2021) suggest a magma mixing, crystallization, and remobilization at 15 km depth, in the deep mushy magma storage identified by other seismic tomographies [START_REF] Spica | Anatomy of the Colima volcano magmatic system, Mexico[END_REF][START_REF] Sychev | Fault-Associated Magma Conduits Beneath Volcán de Colima Revealed by Seismic Velocity and Attenuation Tomography Studies[END_REF], prior to migration to a 'spongy' shallow magma reservoir.

Studies based on melt inclusions consistently show melt inclusion entrapment depths in vapor saturated conditions of typically <8 km for interplinian activity [START_REF] Atlas | Melt inclusions from Volcán Popocatépetl and Volcán de Colima, Mexico: Melt evolution due to vapor-saturated crystallization during ascent[END_REF][START_REF] Reubi | Volatiles contents, degassing and crystallisation of intermediate magmas at Volcan de Colima, Mexico, inferred from melt inclusions[END_REF]Reubi and Blundy, 2008). While the melt inclusions from the major 1913 Plinian eruption reveal depths of entrapment ranging from 8.8 to 11.1 km, depths recorded by the 1998-2005 melt inclusions indicate that the magma from those recent eruptions was shallowly stored [START_REF] Reubi | Volatiles contents, degassing and crystallisation of intermediate magmas at Volcan de Colima, Mexico, inferred from melt inclusions[END_REF]. This Likely matched the seismically active low-velocity anomaly identified in our results. [START_REF] Atlas | Melt inclusions from Volcán Popocatépetl and Volcán de Colima, Mexico: Melt evolution due to vapor-saturated crystallization during ascent[END_REF] also analyzed melt inclusions of Volcán de Colima that they determined were trapped at depths of >12 km, and suggest that vapor-saturated crystallization most likely occurred during the ascent in a series of interconnected conduits in which mixing is allowed rather than a large stratified magma chamber. This corresponds to the overwhelming negative radial anisotropy we observed pointing toward vertically oriented structures, such as dikes and cracks, as the source of the negative anomaly.

Conclusion

We used the group dispersion curves of Rayleigh and Love waves extracted from ambient seismic noise recorded on Volcán de Colima Complex to build a highresolution 3D model of the shear wave velocity and the radial anisotropy in the uppercrust.

The results indicate a network of deeply rooted NE-SW oriented low velocity zones, parallel to the Tamazula fault, that connect at ~2 km b.s.l. This structure overlaps a region of negative radial anisotropy that indicates that magma under volcán de Colima follows vertically oriented structures. In contrast, a distinct high-velocity anomaly under Nevado de Colima highlights the distinction between the former active system, filled with solidified dikes and sills, and the current one, associated with a network of fluid-filled dikes.

This also indicate that the current shallow magma storage is dominated by vertically oriented intrusion, rather than a large single volume of melt. As other seismic tomographies (e.g., [START_REF] Escudero | Ambient seismic noise tomography of the Colima Volcano Complex[END_REF][START_REF] Spica | Anatomy of the Colima volcano magmatic system, Mexico[END_REF][START_REF] Sychev | Fault-Associated Magma Conduits Beneath Volcán de Colima Revealed by Seismic Velocity and Attenuation Tomography Studies[END_REF] , we speculate that the regional tectonics, specifically the Colima Rift combined with a deeply rooted fault system, play a critical role in the transport of magma from a deeper system through the weakened crust of the volcanic complex. Our results offer the clearest insight to date on the internal structure of Volcán de Colima, including the southward migration of the intrusions. This evolution is already evidenced at the surface with structures Los Hijos parasitic domes and the J o u r n a l P r e -p r o o f

Figure 1 .

 1 Figure 1. Tectonic setting of Western Mexico: the Jalisco block (JB), the Michoacán block (MB), the Tepic-Zacoalco Rift (TZR), the Chapala Rift (ChR), the Colima Rift (CR), El Gordo Graben (EGG), and the Colima Graben (CG). The Pliocene-Quaternary Trans-Mexican Volcanic Belt (TMVB) comprises the area within the orange hatched area. The gray triangles are volcanoes, the red triangle shows the location of Volcán de Colima. The inset shows the regional location of the map in Mexico.

Figure 2 .

 2 Figure 2. Location map for the stations from CODEX (blue triangles) and RESCO (orange triangle) used in this study. The Cántaro, Nevado de Colima, and Volcán de Colima form the Colima volcanic complex. The red line corresponds to Tamazula fault. YD is the Yerbabuena depression, CCG is the Central Colima Graben and NCG is the North Colima Graben (modified from Norini et al., 2019). J o u r n a l P r e -p r o o f

Figure 3 .

 3 Figure 3. Group velocity measurements for Rayleigh (a) and Love waves (b) obtained using a FTAN analysis. Number of group-velocity measurements as a function of period for Rayleigh (c) and Love waves (d). We restrict our analysis to periods with at least 50 measurements (dashed vertical line).

Figure 4 .

 4 Figure 4. Rayleigh group velocity maps at 2, 4 and 6 seconds (a-c respectively). Rayleigh-wave group velocities associated with each path at 2, 4 and 6 seconds (d-f). The red contours show the limit of the study area with 4 paths per cell. Ray-path densities per cell at 2, 4 and 6 seconds (g-i).

Figure 5 .

 5 Figure 5. Love-wave group velocity maps at 2, 4 and 6 seconds (a-c respectively). Love-wave group velocities associated with each path at 2, 4 and 6 seconds (d-f). The red contours show the limit of the study area with 4 paths per cell. Ray-path densities per cell at 2, 4 and 6 seconds (g-i).J o u r n a l P r e -p r o o f

Figure 7 .

 7 Figure 7. Resolution shift for Love waves at 2,4 and 6s (a-c). A resolution shift <0.5 km indicates that the center of the retrieved anomaly is in the same cell as the input spike. Spatial resolution maps (d-f).

Figure 8 .Figure 9 .

 89 Figure 8. Best synthetic Rayleigh and Love group velocity dispersion curves at points A and B for an isotropic and anisotropic inversion. The blue and green curves are the best fitting isotropic dispersion curves for Rayleigh and Love waves, respectively. The blue and green dashed-curves are the best fitting anisotropic dispersion curves for Rayleigh and Love waves, respectively. The blue and green dots with the blue and green areas are the measured Rayleigh and Love wave dispersion curves, respectively, with their uncertainties. The map shows the misfit of the final model and the location of points A and B (yellow and red markers, respectively). J o u r n a l P r e -p r o o f

Figure 11 .

 11 Figure 10. Uncertainties of the final 3D Vs model at 1 km a.s.l. (a), sea level (b), 4 km b.s.l. (c) and the corresponding vertical slices (d). The red line is the Tamazula fault system and the white dashed line outlines the Yerbabuena depression. J o u r n a l P r e -p r o o f

Figure 12 .

 12 Figure 12. Vertical slices through the Vs anomaly model (a-c), the Vs velocity model (d-f), and the radial anisotropy (g-i). The location and orientation of the slices are shown in Figure 11a. The distance along the profile is calculated from the crater of Volcán de Colima. The white dashed lines highlight the -10 and +10% contours. The location of the Nevado de Colima along profile B-B'' is identified as NC and YD identifies the location of the Yerbabuena depression along profile A-A''. The black dots indicate the projection of the hypocenters within 500 m of the profile of volcano tectonic events from 2019-2021 calculated using the velocity model of Núñez-Cornú et al. (1994).

J o u r n a l P r e -p r o o fJournal Pre-proof

Acknowledgements

This research was supported by the CONACYT-299766 and CONACYT CF 2019 -58547 grants and Raphael S. M. De Plaen acknowledges support from a UNAM-DGAPA postdoctoral scholarship. D. M. Vargas-Bracamontes is grateful to CONACYT's project PDCAPN 2579. We thank Luca De Siena and an anonymous reviewer for their in-depth review and comments that helped improve the manuscript.

Data Availability

CODEX data are available through IRIS (network ZA, dio: 10.7914/SN/ZA_2006) and the RESCO data are available through the University of Colima. The Velocity model has been deposited at Zenodo (10.5281/zenodo.6475834).

Declaration of interests

☐ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

☐The authors declare the following financial interests/personal relationships which may be considered as potential competing interests:

J o u r n a l P r e -p r o o f

Highlights

We built the highest-resolution anisotropic shear wave velocity models to date of Volcán de Colima A distinct high-velocity anomaly under Nevado de Colima highlights the distinction between the former, now solidified, system

We observe a network of deeply rooted NE-SW oriented low velocity zones, parallel to the Tamazula fault The negative radial anisotropy indicates that magma under volcán de Colima follows vertically oriented structures such as dikes and cracks J o u r n a l P r e -p r o o f