Imaging the Two Magma Chambers of Tequila Volcano, Mexico, by Gravity Analysis ()
1. Introduction
Tequila volcano is located within the Tepic-Zacoalco rift in central-western Mexico; the rift is a tectonic, extensional region containing several volcanic structures, Tequila volcano among them, whose youngest activity occurred at ~90 ka erupting through the flanks of the main cone [1]. In Figure 1, a group of volcanic edifices is displayed, aligned with the Tequila volcano (TQ). This alignment is identified with the north limit of the Jalisco Block [2]. The distance between TQ and its nearest volcanic neighbor, La Primavera domes, is ~30 km; it is worth noticing the closeness of the times of the latest eruptive activity of both volcanic centers: Tequila volcano developed its latest activity at 90 ka, whilst in La Primavera it occurred at 95 ka [3]-[5], and their locations along the north limit of the Jalisco Block.
Tequila volcano reaches an altitude of 2840 m, rising from a plateau of 1000 m and is composed of ~25 km3 of pyroxene andesite lavas and pyroclastic deposits [6] [7]. When the volcano’s activity ceased, abundant eruptions emerged on its flanks [8]. TQ shows a peculiar summit, since it shows a depression resembling a caldera, which is common in other volcanoes of the region (e.g., [9] [10]), or induced by powerful explosions [11] [12]. However, previous reports found that features associated with those phenomena are not present at TQ [13]; the summit depression probably originates in erosional processes [8]. The latter authors mention that the most prominent feature of TQ is a 300 m spine occupying a summit depression that is breached on the eastern side. Here we will look at this problem, analyzing density cross-sections across the summit.
The latest activity at TQ occurred ~90 ka through the flanks of the volcano; it has been concluded that TQ is an unusual stratovolcano that seems to evolve in pulses of petrologically unrelated magmas [1]. La Primavera caldera was formed at ~95 ka inducing a collapse of ~11 km diameter, closely coinciding with TQ latest activity [3]-[5].
For TQ, a magma chamber at 2 - 3 km depth is inferred from phenocryst assemblages [14]. After a lapse of ~110 kyrs, eruptions followed on the NW and SE flanks of TQ, probably originating in a second magma chamber at ~6 km depth [14]. Cerro Tomasillo (~60 ka) is a small andesitic volcano SE of Tequila Volcano.
Multiple magma chambers feeding volcanic systems are not unusual. We have reported on similar situations in the Sierra Nevada [15] and in the Iztaccihuatl volcano [16]. Comparable dual magma-chamber systems have also been discussed theoretically for active volcanoes [17] and described for Yellowstone, where a shallow rhyolitic chamber overlies a deeper basaltic reservoir [18].
Figure 1. Location map of Tequila volcano (TQ) showing the alignment of five volcanic structures (dashed line). Arrows indicate the exact location of each volcanic edifice. SJ, San Juan. TE, Tepeltitic. CE, Ceboruco. LP, La Primavera. GC, Gulf of California. The alignment defines the inferred northern boundary of the Jalisco Block [2]. Digital Elevation Model from GeomapApp [19].
A digital elevation model (DEM) of the Tequila Volcanic Region appears in Figure 2, where topographic contours help identify some of its surface characteristics. Of particular interest is the NW-SE alignment of the summit of the main structure (TQ), Cerro Tomasillo (CT), and five volcanic cones (VC) that extend along 14 km in the SE direction. Regions designated Younger Flows (YF) and Younger Andesitic Flows (YAF) west of the summit of TQ are reported [14]; they suggest that a prominent NW-SE lineament of TQ and andesite flank lavas is an indication that they overlie a major passageway for the ascent of magmas from the lower to the middle crust, and that the two chambers under TQ were probably the site of mingling between ascending magma batches. Subsequently, we will show the exact location of this passageway.
2. Methods
2.1. Data Acquisition
High-resolution gravity data for the La Primavera-Tequila volcanic area were obtained from the satellite-derived GGMplus model. This data set provides near-global coverage within the latitudinal band between 60˚N and 60˚S, although some offshore regions remain less well constrained because of limitations in the underlying source data.
The Bouguer anomaly (BA) used in this study was calculated following standard gravity reduction procedures commonly applied to satellite-derived gravity data and consistent with recent regional studies in Mexico [20]-[26].
Figure 2. Digital Elevation Model (DEM) of Tequila volcano and immediate surroundings; 25 m contours enhance topographic features. Cerro Tomasillo (CT) is located SE of the summit, continuing with an alignment of volcanic cones (VC) traversing the Younger Flanks (YF) [14]. YAF, Younger Andesitic Flanks.
2.2. Full Gravity-Processing Workflow
To derive the residual Bouguer anomaly, the gravity workflow begins by establishing the geographic limits of the analysis. A typical regional gravity study defines the latitude and longitude ranges. In the present case they are: 20.0˚N - 21.0˚N and 103.0˚W - 104.0˚W, respectively (Figure 3). A grid resolution is required, matched to the digital elevation model and the gravity model (e.g., 250 m for GEBCO/GGMplus). These bounds ensure consistent sampling of DEM, terrain corrections, and gravity grids (e.g. GEBCO 250 [27]), which is consistent with global Bouguer workflows described in [28].
Apply ellipsoidal normal gravity, computing normal gravity (γ) at each grid node using the reference ellipsoid (GRS80 or WGS84). This removes the predictable component of Earth’s gravity field.
(1)
This step is part of the “complete spherical Bouguer anomaly” workflow described in [28]. Next, apply freeair correction, correcting for the elevation above the reference ellipsoid:
(2)
This yields the freeair anomaly.
Apply Bouguer slab correction, removing the gravitational effect of the rock between the station and sea level:
(3)
Typical density: 2.67 g/cm3 for crustal rock. This produces the simple Bouguer anomaly.
Apply terrain correction. Terrain correction is essential for accurate Bouguer anomalies, especially in volcanic or mountainous regions. Terrain corrections account for the gravitational attraction of nearby topography, methods include [29] [30], and highresolution DEMbased terrain correction (GEBCO 250 m); terrain corrections are critical for obtaining physically meaningful Bouguer anomalies and avoiding systematic biases.
This step yields the complete Bouguer anomaly. Apply atmospheric correction (optional but recommended). Correct for the gravitational effect of the air column:
(4)
Small but included in modern workflows.
Finally, compute the Bouguer anomaly. The full Bouguer anomaly is:
(5)
This matches the complete Bouguer anomaly formulation used in global gravity forward models [28].
2.3. Regional-Residual Separation
To isolate shallow crustal features, the Bouguer anomaly is decomposed into a regional field (longwavelength) and a residual field (shortwavelength). Methods include polynomial surface fitting, fitting the 1st - 4th order polynomial to BA, subtracting to obtain residuals, spectral filtering, low-pass filter (regional), high-pass filter (residual), and upward continuation, among others. Subsequently, generate the residual Bouguer anomaly map, which highlights shallow density contrasts, faults, intrusions, volcanic conduits, and sedimentary thickness variations. Thus, residual anomalies are more sensitive to crustal heterogeneities than complete Bouguer anomalies.
To emphasize the shallow-to-intermediate density contrasts associated with the volcanic structures, the gravity field was further processed using a Gaussian filter to separate the residual component. This residual BA dataset was then used as the basis for the structural interpretation presented below (Figure 3). In volcanic regions, negative gravity anomalies commonly reflect low-density zones related to volcanic edifices, magma reservoirs, and associated plumbing systems, e.g., [15] [23] [26] [31] [32].
2.4. The Regional Bouguer Anomaly
Regional gravity anomalies describe the behavior of gravity in extended regions; in the present case it describes the regional behavior of gravity in a region comprising La Primavera volcanic field and the Tequila volcano area. Residual gravity values in the map display a range of 37 mGal; a high gradient is observed on the cone of Tequila volcano, evolving into a NW-SE elongation across it.
![]()
Figure 3. Residual Bouguer Anomaly map of the north portion of the Jalisco Block; the NW corner exhibits a sharp, low-gravity anomaly associated with Tequila volcano. A black and a green rectangle show the regions in which 3D gravity inversions are performed. The SE corner of the black rectangle is located close to the central portion of La Primavera volcanic field. Lines 1 - 4 intersect each other at the location of Tequila volcano and indicate the extent and orientation of the density cross-sections that will be obtained after performing the corresponding gravity inversions. TQ, Tequila. CC, Cerro Cuauhtépetl. LP, La Primavera. GDL, Guadalajara. CHA, Chapala. ZA, Zacoalco de Torres. Contours are topographic.
2.5. 3D Gravimetric Inversion Workflow
Following the calculation of the Bouguer anomaly (BA), a 3D gravity inversion was performed to obtain volumetric density-contrast models and thereby constrain the subsurface structure associated with the main tectonic features of the region. The inversion was carried out in Oasis Montaj (Seequent, version 2025.2) using the workflow described in [33], based on the theoretical framework in [34]. The inversion workflow implemented in Oasis Montaj follows the procedure described in [33], which is based on the general potentialfield inversion framework developed in [34]. Although these authors present the method in the context of magnetic data, the underlying theory is not specific to magnetism. Instead, it relies on the mathematical properties common to all potential fields, namely that both magnetic and gravitational anomalies can be expressed as linear combinations of volumetric source parameters through a potential-field sensitivity kernel. In the magnetic case, the model parameter is the magnetic susceptibility, and the forward operator relates susceptibility contrasts to the observed magnetic anomaly. In the case of gravity, the model parameter is the density contrast, and the forward operator relates density variations to the Bouguer or residual gravity anomaly. Because both problems share the same linear structure
, the same discretization strategy, and the same regularization and misfit-minimization scheme, the inversion algorithm described in [34] can be applied directly to gravity data by substituting the magnetic kernel with the gravitational kernel. This equivalence justifies the use of the workflow in [33] for 3D gravity inversion in Oasis Montaj.
This approach discretizes the subsurface using a Cartesian cut-cell (CCC) representation and estimates the 3D distribution of density contrast by minimizing the misfit between observed and forward-calculated gravity data through an iterative reweighting inversion (IRI) procedure [35]. In the present implementation, iterations were continued until the data misfit was reduced to within 5% of the standard deviation of the observed BA values. This threshold was adopted as a practical convergence criterion to achieve an adequate balance between data fit and model stability under regularization.
The misfit criterion was based on the RMS difference between the observed BA and the model-predicted gravity. A 5% tolerance was evaluated as follows: First, the standard deviation (σ) of the Bouguer residual field was computed, second, the objective was to reach an RMS misfit ≤ 0.05σ, third, iterations continued until this threshold was reached and further adjustments produced no geologically meaningful improvements. This makes certain that the model fits the data within the uncertainty inherent to the gravity field and the DEM-based corrections. What is warranted is that iterations continue as long as the misfit is above the imposed value, which in turn is established before the calculation is executed, according to the problem to be solved. We usually assign a 5% misfit to problems in which the resolution involves 1 km3 prisms, where the averaging of the rock density over such volumes may introduce larger errors.
The inversion input consisted of gridded BA data (mGal), which were transformed into density-contrast models (g/cm3) through the 3D inversion process. Two inversions were computed using cell sizes of 1000 m and 250 m, respectively, to evaluate the density structure at different spatial resolutions; the first cell size corresponds to 1 km3 prisms, whilst the latter increases that resolution 64 times. We used the GGMplus data set, available from http://ddfe.curtin.edu.au/gravitymodels/GGMplus/ with 7.5 arc-sec (~220 m) resolution [36]. The maximum depth of each model was selected according to the adopted discretization, yielding model depths of approximately 9 km for the 1000 m grid and 5 km for the 250 m grid. Thus, the depth extent and lateral dimensions of the recovered density structure were controlled by the selected computational domain and the model resolution. The resolution of the GGMplus gravity model, being smaller than the sizes used in the inversion processes, of 1000 and 250 m, assures that no computational artifacts are introduced in the inversion process. The former resolution is used to reach depths of ~8 km from the surface, to explore such depths at a low-resolution, whilst the latter provides larger detail at depths of ~4 km from the surface.
The resulting density distributions were subsequently analyzed and interpreted in terms of their possible relationship with subsurface geological materials and structural features. Figure 4 presents the inversion result obtained with the 1000 m resolution model, which extends to approximately −6 km below sea level, equivalent to a total depth of about 9 km from the average surface elevation. In the 250 m resolution cross-sections, the corresponding DEM-derived topographic profile was included to facilitate comparison between density variations and surface morphology.
Model performance was assessed by comparing observed and predicted gravity responses, including residual anomaly maps, and by confirming that the final solution satisfied the convergence criterion defined above.
3. Density Cross-Sections
Two resolutions are used in these analyses: 1000 and 250 m. The finer-resolution model provides greater detail in the recovered density structure, although at the expense of a shallower depth extent. The locations of the density cross-sections are defined by Lines 1 - 4 in Figure 3, which traverse the black and green rectangles, respectively.
3.1. 1000 m Resolution
Two orthogonal, density cross-sections through the summit of Tequila volcano are shown in Figure 4. In both, and directly under the summit of TQ, a low-density region shows a minimum at an elevation of +500 m, as shown by the contours; we frequently associate this type of minimum with the location of magmatic chambers, which in the present case coincides with the depth of the shallow magma chamber proposed in [14]; gravimetrically confirming its presence. On the eastern portion of the line, separated by a high-density region, there is a large negative anomaly that displays a minimum between −4000 and −6000 mbsl, coinciding with other minima at similar depths to be discussed ahead. The N-S cross-section shows a high-density region at the summit; this type of response is obtained when domes are extruded close to the summit, although in the present case it appears to be associated with the spine reported in [8]. We shall further comment on this when discussing the 250-m resolution profiles.
In the N-end, a low-density region is associated with Santiago River; the topographic profile corresponds to the presence of a canyon, or graben-like depression, with its associated shallow, low-density region, which may be ascribed to sedimentary materials close to the surface and a possible cortical fault extending to the bottom of the section.
Figure 4. Density cross-sections at 1000-m resolution (uniform density values in 1 km3 prisms) through the summit of Tequila volcano, reaching 6 km bsl. (a) E-W cross-section, L1 in Figure 3. (b) N-S cross-section, L2 in Figure 3. TQ, Tequila volcano. SR, Santiago River. The color scale represents density values relative to +2.67 g/cm3.
Figure 5(a) corresponds to L3 in Figure 3, oriented in the SW-NE direction; the high-density region observed at the top of TQ in Figure 4(b) is not intersected in this orientation. In this projection, the width of the low-density region associated with TQ is ~5 km and is flanked by high-density regions, which appear to be divided by the fault plane of the prominent NW-SE lineament reported in [14]. The location of Santiago River is also intersected as a low-density region of shallow depth. Figure 5(b) shows the NW-SE density cross-section through the summit of TQ. This orientation reveals that feeding of the TQ magma chamber at 500 masl occurs laterally, rather than vertically, as shown by the red arrows bifurcating at 5000 mbsl. This is the only cross-section showing a widespread distribution of low-density materials, drastically contrasting with the extent of the low-density region intersected in the SW-NE density cross-section (Figure 5(a)), implying that it coincides with a fault plane that reliefs magma ascent, confirming the account in [14]. It also intersects the position of Cerro Tomasillo, a small andesite volcano SE of TQ and the regions to the NW and SE, the flows called the Younger Flanks; all appear to be fed by a source deeper than 6 km bsl. Continuing to the SE, the volcanic region ends in a high-density region; at the end of the line is La Primavera caldera, where lower density regions replace portions of the high-density region.
![]()
Figure 5. Density cross-sections through the summit of Tequila volcano; both reach 6 km bsl. (a) SW-NE cross-section, L3 in Figure 3. (b) NW-SE cross-section, L4 in Figure 3; red dashed lines indicate potential magmatic-flow trajectories. The question mark refers to the location of the feeding trajectory of Cerro Tomasillo (CT), which probably occurs perpendicular to this plane. TQ, Tequila volcano. SR, Santiago River. YF, Younger Flank. The color scales represent density values +2.67 g/cm3.
At this resolution, regarding the magmatic chambers proposed in [14]; we can establish the existence of a surficial magma deposit centered at an elevation of +500 m, or 1.5 km depth, corresponding to the magma chamber suggested by them between 2 - 3 km depth. We find a feeding channel displaced ⁓8 km SE of the summit of TQ at a depth of −3300 mbsl that connects to a deeper magma source. Since the NW-SE density cross-section (Figure 5(b)) is the only one showing considerable dispersion of low-density materials, we conclude that this is probably a fault plane that enables magma ascent, particularly since the main edifice (TQ) is not directly above the main feeding channel but displaced ⁓8 km to the NW. A similar lateral feeding was registered for Popocatépetl volcano in the Sierra Nevada [15].
To elucidate the connection between the shallow and the deep magma chambers we obtained horizontal density sections at various depths (Figure 6), where a dashed circle marks the surface location of Tequila volcano. At sea level (0 m) the low-density anomaly appears elongated in the NW-SE direction with its thickest portion to the SE. This anomaly is flanked by high-density regions to the NE and the SW, creating a sort of low-density channel; this arrangement persists down to −5900 m. The encircled, low-density anomaly persists to −3300 m, disappearing at −4500 m, being substituted by a larger, deeper anomaly to the SE that reaches the S limit of the sections, and begins to appear at −2100 m. At −3300 m, within the circle remains a thin, low-density passage that we interpret as the fault that connects the main magma chamber with the one located at +500 m (Figure 5). At −5900 m the NW-SE lowest-density region is clearly defined, which we interpret as the magma chamber at 6 km depth proposed in [14].
![]()
Figure 6. Sequence of horizontal sections of the inverted volume corresponding to the black rectangle in Figure 3, at six depths below sea level (0 m) showing how the low-density region (blue) under TQ varies with depth (red circle), disappearing at 4500 m bsl, the location of the feeding channel, where only remains the low-density region associated with the deeper magma chamber at 6500 m depth. On the N portion of the sea-level section, the trace of Santiago River is clearly discerned as a concave, thin low-density region, whose finer trace begins to disappear at −2100 m. The density scale is the same as that in Figure 5.
3.2. 250 m Resolution
The 250 m resolution sections correspond to the lines inside the green rectangle in Figure 3, for a more detailed description of the upper structure of Tequila volcano.
The E-W cross-section at this resolution (Figure 7) reveals a peculiar bifurcation of the contour lines at the summit, which we interpret as probably originating in the occlusion of the original discharge conduit by the spine described in [8], which forced a new, adjacent discharge trajectory.
Figure 7. E-W density cross-section corresponding to L1 inside the green rectangle in Figure 3 with a resolution of 250 m and a 12 km extent. The two red, vertical lines enhance the correlation between the low-density region and the breached summit of TQ. The red, dashed line bifurcates near the summit, as indicated by the contours in the low-density region. The color scales represent density values +2.67 g/cm3.
The bifurcation of the contour lines is also observed in the N-S projection (Figure 8); additionally, a small, high-density region appears at the summit, which we identify with the spine (SP), confirming the assumption of occlusion of the original discharge conduit and subsequent deviation.
Figure 8. N-S density cross-section corresponding to L2 inside the green rectangle in Figure 3 with a resolution of 250 m. In the upper section appears the corresponding topographic profile with the same horizontal scale as the cross-section. The red dashed line traverses the lowest density region, representing a likely feeding trajectory. The lowest density region reaches the summit where it bifurcates owing to a high-density region that we recognize as the 300 m spine (SP) described in [8], corresponding to the high-density region at the summit. The color scales represent density values +2.67 g/cm3.
The bifurcation of the contour lines is also observed in the NW-SE projection (Figure 9); the small, high-density region also appears at the summit (SP). In this projection mid-density regions (brown arrows) become high-density regions as they approach the surface on the volcano’s flanks. They appear to belong to abandoned flow trajectories that fed flank eruptions.
![]()
Figure 9. SW-NE density cross-section corresponding to L3 inside the green rectangle in Figure 3 with a resolution of 250 m. In the upper section appears the corresponding topographic profile with the same horizontal scale as the cross-section. The central, red dashed line traverses the lowest density region, representing a likely feeding trajectory to the summit of TQ. The two red, vertical lines enhance the correlation between the low-density region and the breached summit of TQ. Similarly to the observation made in Figure 8, the lowest density region reaches the summit also bifurcates owing to the high-density region that we identify with the 300 m spine (SP). The gray dashed lines represent possible, abandoned trajectories of volcanic materials that reached the flanks of TQ, inducing the flank flows. The color scales represent density values +2.67 g/cm3.
Figure 10 presents similar features to those in Figure 8 and Figure 9, strengthening the proposition that the trajectory bifurcation of the exhaust path is responsible for the caldera-type summit in Tequila volcano.
![]()
Figure 10. NW-SE density cross-section corresponding to L4 inside the green rectangle in Figure 3 with a resolution of 250 m. In the upper section appears the corresponding topographic profile with the same horizontal scale as the cross-section. As observed in Figure 9, the central, red dashed line traverses the lowest density region, representing a likely feeding trajectory to the summit of TQ. The two red, vertical lines enhance the correlation between the low-density region and the breached summit of TQ. Similarly to the observation made in Figure 8, the lowest density region that reaches the summit also bifurcates owing to the high-density region that we identify with the 300 m spine (SP). YAF, Younger Andesitic Flank. CT, Cerro Tomasillo. The latteral dashed lines represent possible, abandoned trajectories of volcanic materials that reached the flanks of TQ, inducing the flank flows. The color scales represent density values +2.67 g/cm3.
4. Discussion
Regarding the depth and location of the deeper magma chamber we will consider Figure 5 and Figure 6. The NW-SE cross-section in Figure 5(b) shows that the low-density anomaly feeding TQ and the volcanic formations around it, reaches the bottom of the cross-section (6 km bsl), revealing its downward continuation. The horizontal section at −4500 m in Figure 6 shows the top of the deeper magma chamber as a faint blue, whereas at −5900 m darker blue is established in half the area of the anomaly; the darker blue corresponds to the region where the density is lower, or where magmatic products are concentrating. We conclude that the center of the deeper magma chamber is located between 6200 and 6500 mbsl. From this location magmatic materials are sporadically ejected upwards through the feeding channel.
A peculiar bifurcation of the low-density region systematically appears close to the summit in the 250 m resolution cross-sections (Figures 7-10), which originates in the upper magma chamber located at +500 m. The mechanism we prefer for inducing this bifurcation is the blockage of the initial chimney, most likely by the injection of the spine, potentially followed by an explosion that partially destroyed the summit, creating a new discharge trajectory. The spine appears in the cross-sections of Figures 8-10, as a small high-density region. We have reported similar chimney blockages and subsequent explosions at Nevado de Toluca [23] and Popocatépetl volcanoes [15].
The latest activity in TQ occurred ~90 ka through the flanks of the volcano, probably originating in the second magma chamber at ~6 km depth, as inferred in [14]. In Figure 10, there are discharge trajectories ending at the flanks of TQ, which we interpret as the conduits of the flanks’ lava flows that have now solidified, presenting greater densities.
The path of Santiago River is intersected in Figure 4 and Figure 5 associated with a topographic depression and a shallow, low-density region that we attribute to sediment accumulation. In Figure 6, it displays a thin, semi-circular trajectory that persists from sea level to −2100 m, indicating a crustal incision of over 4 km from the surface. This is the section of SR closest to Tequila volcano and its trajectory appears to be influenced by the dominating NW-SE orientation of volcanic structures in this region.
5. Conclusions
Applying gravimetric treatments to Tequila volcano and its surroundings, we confirmed the existence of the two magma chambers proposed in [14]; we found that only the shallow chamber is directly under the volcanic cone, with the deeper one displaced 8 km to the SE.
The connection has been established between the two magma chambers feeding Tequila volcano and the surrounding areas designated as the Younger Flanks, between elevations of −6 km and +500 m, both connecting through a narrow passage located at −3300 m elevation.
The caldera-type summit is explained based on the bifurcation of the low-density region close to the volcano’s summit, as the result of an obstruction of a former chimney by an intrusive body, called the spine, inducing an explosion, and deviating the exhaust to its new trajectory.
Acknowledgements
This work was supported by Universidad Nacional Autónoma de México Postdoctoral Program. M. Camacho-Ascanio performed his academic work at Instituto de Geofísica, UNAM. This study has been supported by the Institute of Applied Mathematics and Systems (IIMAS), of the National University of Mexico (UNAM); we acknowledge material support from both institutions.
Author Contributions
Conceptualization, R.A. and M.C.-A.; methodology, R.A. and M.C.-A.; validation, R.A. and M.C.-A.; writing—original draft preparation, R.A. and M.C.-A. All authors have read and agreed to the published version of the manuscript.