1. Introduction

Geophysical methods have been a useful tool for monitoring the activity of different volcanoes worldwide. Vulcanologists who use these methods have pursued three main objectives: i) characterization of magma movement through the crust and its location within a magmatic chamber, ii) observation of the hydrothermal activity associated with a volcano, and iii) monitoring of the stability of the volcanic edifice (Wynn et al., 2006). Gravimetry has been essential for the study of changes in the stress tensor and the temperature gradient at the subsurface due to intrusions of magmatic bodies. This results in density variations and surface deformations where the body intrudes, which can be detected by fluctuations in the vertical vector of gravitational attraction (e.g., Vajda, 2016). The monitoring of such dynamics is especially relevant since it has been demonstrated that these deformations on the volcanic edifice are precursors of volcanic eruptions. Thus, their monitoring is fundamental in the evaluation of volcanic hazards (e.g., Vajda, 2016).

Our main objective in this contribution is to present a gravimetric study of the Cerro Machin volcano and to obtain (using data acquired from fieldwork between March 16th and March 19th, 2020) a representative model of the density distribution from its crater area at depth. Based on the total Bouguer anomaly obtained from the field data, this model was evaluated and compared with the geology and seismicity of the area, aiming to contribute to the general characterization of the processes associated with the activity of the Cerro Machin volcano.

2. Tectonic-structural setting

The Cerro Machin volcano (2750 m MSL) is located on the eastern flank of the Colombian Central Cordillera, 17 km to the northwest of the city of Ibagué (Murcia et al., 2008). According to the stratigraphic record, it has produced at least six major dacitic eruptions, the last one dated at 900-year B.P. (Laeger et al., 2013). This volcano is part of the Northern Volcanic Zone (NVZ) of the Andean Volcanic Belt, which stretches up to Ecuador. This belt is the result of the subduction of the Nazca Plate beneath the South American Plate (Laeger et al., 2013). The volcanic arc in Colombia runs parallel to the Colombia-Ecuador Trench, located 300 km west of the arc (Rueda, 2005). The arc subdivisions are usually named considering the proximity of the volcanoes to each other (Monsalve and Pulgarín, 1993). Nonetheless, other authors have established the existence of a sole segment (named the Cauca segment), which includes the totality of the Colombian volcanic arc (Rueda, 2005). This segment is subdivided into three subsegments: North, Central, and South. The Cerro Machín volcano is located in the southernmost area of the northern subsegment, along with other volcanoes such as Nevado del Ruiz, Cerro Bravo, Nevado del Tolima, and Nevado del Quindío. This subsegment extends for 90 km from the SW to the NE (Hall and Wood, 1985) and is defined by two main fault systems in the region: the Cauca-Patía (east of the Colombian West Cordillera) fault, the Romeral-Dolores fault (west of the Colombian Central Cordillera) and the secondary faults associated with both main faults (Rueda, 2005).

Cerro Machin is located at the intersection of two faults: the Cajamarca fault, with a N 20° E strike and vertical dip, and the Machin fault, with a N 42° W strike and normal movement with a dextral component (Rueda, 2005). Both faults have shown activity during the Holocene since their movements affect only volcanic deposits and paleosols associated with the activity of the volcano during this period (Rueda, 2005). The formation of Cerro Machin could be attributed to the existence

Figure 1.
a) Geological map of the study area extracted from Geology Map 244 of the Servicio Geológico Colombiano (Mosquera et al., 1982), Machin fault from Méndez Fajury et al. (1996). b) Global and regional location of the Cerro Machin volcano in Colombia. c) N-S Panoramic view of the Cerro Machín volcano

of a pull-apart structure (with undefined dimensions) generated by the lateral movement of the previously mentioned faults of the region. This rifting structure could have facilitated the ascent of magmatic material that led to the formation of the volcano (Rueda, 2005).

2.1. Geological setting

The volcano formed on top a set of metamorphic rocks known as the Cajamarca Complex (Figure 1), which is the central core of the Colombian Central Cordillera (Piedrahita et al., 2018). These rocks are characterized by low-pressure regional metamorphism with a predominance of green schist, quartz-sericitic phyllites and quartzites generated from mafic lava flows, pyroclastic material, and sedimentary rocks of different compositions and grain sizes. These rocks were deposited in the Triassic during rifting between North and South America and subsequently underwent different anatexis events due to the continental arc magmatism that occurred in the region during the Jurassic (Laeger et al., 2013).

The actual volcanic edifice is formed by a dacitic dome and various pyroclastic flow deposits resulting from the last and previous eruptive activities of the volcano (Laeger et al., 2013). These deposits are grouped in a unit called “Anillo”, which is divided into five lithofacies types generated from different pyroclastic flow density currents of high turbulence and low density currents whose formation could involve water (Piedrahita et al., 2018). Furthermore, this and other differentiated units from previous recognizable eruptions have been grouped into a single stratigraphic unit named the Machin Formation, divided into five subunits, each associated with eruptive volcanic activity (Rueda, 2005). Their composition consists mainly of tephra and pumice, as well as pyroclasts and distal deposits related to the lahars formed during the eruptive events (Méndez, 2002).

3. Methods

The study was divided into three sections:

3.1 Data collection campaign at the Cerro Machin volcano

Gravimetric and GPS data were obtained during a field campaign conducted from March 16th to March 20th, 2020. Forty measurements were taken, mostly inside the crater of the volcano around the Machin dome (Figure 2), with two measurements made on the two summits of the dome, using a Scintrex Autograv CG-5 gravimeter. For the base station of the gravity survey, the absolute gravity value (gA) and its coordinates were provided by the Servicio Geológico Colombiano (SGC) (gA: 977480.179 mGal; latitude: 4.48537; longitude: -75.38008; height: 2526.2 m). It was used to calculate the absolute gravity at each station.

3.2 Data processing and Bouguer correction

For the processing of the data obtained during the campaign, we prepared an Excel spreadsheet. Here, the appropriate gravimetric correction equations were applied to the data as follows:

where L1B1 and L2B1 are the gravity measurements taken at the base station at the beginning and at the end of the day, respectively; t1B1 and t2B1 are the times at which these measurements were performed, respectively; Δt is the difference between the time at which a measurement was taken at the station and at the base on a given day; λ is the latitude coordinate of the measurement station; l is the linear distance between the measurement station and the base; h is the height of the measurement station relative to the datum; G is the gravitational constant (6.672 * 10−11 N m2/kg2); and ρ is the plate density (in kg/m3).

To calculate the instrumental drift correction (ΔgD), we applied Equation 1. Using the data collected at the base at the beginning and at the end of each day, we obtained a drift rate, which was then multiplied by the time interval between the

Figure 2.
Map showing the contour level map showing the measurement stations throughout the volcano crater and the transects used on the 2D density models on Oasis

initial and final daily measurements at the base. Subsequently, we either added this value, if the measurement at a station was lower than that at the base, or subtracted this value, if the measurement at a station was larger than that at the base, from the relative gravity data obtained. This value was then added to the absolute gravity at the base (provided by the Servicio Geológico Colombiano) and repeated to tie each station by adding the value of the relative gravity measured at each station with the value obtained at the base. This total variation in gravity was finally added to the absolute gravity at the base to obtain the measured gravity at each station.

The latitude correction (Δglat) was calculated using Equation 2, for which we considered the latitude coordinate of every station and its distance to the base. The free air correction (ΔgFA) was obtained using Equation 3, taking the difference in height between the station and the sea level, considered as the reference. We then applied Equation 4 to calculate the Bouguer plate correction (ΔgBP), for which we needed the difference in height between the station and the reference level, as well as the existing mass between both points; thus, an optimal density for this mass was obtained. An average crustal density of 2670 kg/m3 (Hinze, 2003) was implemented. Finally, for each point, we calculated the free air correction to obtain the simple Bouguer anomaly, which can be expressed by the following equation:

(5)

where gs is the relative gravity at the station in relation to the reference point, ℎ is the height of the measurement station (in meters), and ρ is the density used for the calculation of the Bouguer plate correction (in kg/m3).

Finally, a topographic correction (ΔgT) was applied to the data collected at all stations. This correction considers the irregularities of the surface in the vicinity of the stations. These irregularities may exert an upward pull on the gravimeter for an excess of mass or a downward pull for a lack of mass in the area, which may affect the gravimetric measurements. The topographic correction at each measurement point was calculated by using a computational method in MATLAB provided by Palacios (2017). This implementation required the use of a digital elevation model (DEM) extracted from the Copernicus DEM (OpenTopography, 2021), which was used as an input for the calculation. Another script from Palacios (2017) transformed the topography provided by the DEM into polyhedra of different sizes and two distinct shapes (squared and triangular prisms) and divided the area covered by the DEM into six concentric, square-shaped sections centered on the base station. Over 5 million polyhedra were created with this method. Thus, the correction values for each station depended on the polyhedra where they were located.

All aforementioned calculated values were combined to obtain the total Bouguer anomaly for each station using the following equation:

(6)

where is the absolute gravity, is the gravity value obtained with the gravimeter, and is the theoretical gravity for each location, which was obtained using the International Gravitational Formula (Dentith and Mudge, 2014):

(7)

where is the latitude coordinate of the station.

3.3 Generation of the 2D forward models using Oasis Montaj

The processing of the gravity and topographic data (gridding and imaging) and the forward modeling of the gravity data were carried out using Oasis Montaj software and its GM-SYS extension. To do so, the process consisted of i) the interpolation and generation of the topographic and total Bouguer anomaly maps and ii) the forward modeling of a test configuration of the density distribution of the volcano at depth using GM-SYS.

The forward modeling was carried out with the GM-SYS extension of Oasis Montaj. This modeling type involves the generation of a hypothetical density model where the physical parameters (such as the density and geometry of the model) are adjusted by the interpreter until a match is obtained between the computed response of the model and the observed data. This is an iterative process that requires the model and both the observed and computed responses to be displayed graphically so that the results can be assessed by the interpreter (Dentith and Mudge, 2014, p. 74). GM-SYS is a modeling program that allows an interactive manipulation of the geological model and real-time calculation of the gravity response. The methods used in this program to calculate the gravity and magnetic model responses are based on the methods developed by Talwani et al. (1959) and Talwani and Heirtzler (1964) and make use of the algorithms described in Won and Bevis (1987) (NGA, 2004).

GM-SYS uses a two-dimensional, flat-Earth model for gravity and magnetic calculations; that is, each structural unit or block extends to plus and minus infinity in the direction perpendicular to the profile. The Earth is assumed to have topography but no curvature. The models also extend by default plus and minus 3 × 104 kilometers along the profile to eliminate edge effects (NGA, 2004).

4. Results

4.1 Field-obtained Bouguer anomaly

Figure 3 shows the results after the application of the gravimetric corrections to the field-obtained data (see Table 1 in Supporting Information). The resulting Excel spreadsheet was imported into Oasis Montaj, which was used to generate the Bouguer anomaly map (Figure 3a) of our study area. The locations of the recorded earthquakes in the study area (León, 2019) were included in the Bouguer anomaly map, as one of the main objectives of our work was to identify the possible origins of the observed anomaly. The transects where the 2D depth profiles were obtained are also included in Figure 3b.

Figure 3.
a) Bouguer map calculated using Equation 6 in the study area. b) Bouguer anomaly map, overlapped with the earthquakes shown by black dots, and the N-S and W-E transects for the depth models shown in Sections 3.3 and 4.2

The field-obtained anomaly showed a range of values ranging between approximately −87 and −29 mGal (see values in the Supplementary Data). The lowest anomaly values were found in the area above the topographic highs of the dome (Figure 2). Figure 3b also shows that the earthquakes were mostly located southwest of the volcano. The locations of the earthquakes, all of which were associated with activity below the Machin dome, support the idea of the existence of an anomaly in this area. In addition, the data processing designated a region where the anomaly was less negative (−29 mGal) over the entire study area. This measurement was located at the southernmost station at a low elevation near a hot spring at 2412 meters.

4.2 Oasis Montaj 2D forward modeling

To better understand the way the density distribution of the geological bodies affected the observed gravity, N-S and W-E depth profiles, shown in Figure 2, were analyzed. The resulting models suggest the following distribution: a dacitic dome complex at the top and a weakened section below (simplified as a single unit that comprises a complex system of dikes and volcanic conduits) embedded over the edifice, which comprises volcanoclastic sequences from the previous eruptive activity, whose density increases with depth.

For this distribution, we assumed the volcanic edifice to be surrounded by host rock with an average crustal density of 2670 kg/m3, which is represented by the green layers in Figure 4. Within the study area, we set 3 layers named volcanoclastic sequences that sit over the metamorphic basement, cut by the dacitic complex that connects with deeper magmatic chambers. The density of the fragmented part of the dome, which we set up between the topographic relief and a depth of 0 km, is 2300 kg/m3, while volcaniclastic sequence layers vary between 2500 kg/m3 (uppermost volcaniclastic sequence layer) and 2670 kg/m3 (deeper volcaniclastic sequence layer). Ultimately, the metamorphic complex density was established at 2900 kg/m3, close to the dacitic complex in the northwestern regions of the volcano and tilting downward toward the southeast. It is important to note the presence, below the fragmented dome complex, of a weakened section, modeled as a simplified unit that includes a complex system of dikes and volcanic conduits, with a density of 2400 kg/m3. This distribution allowed us to obtain the best adjustment to the observed anomaly. We associate this fragmentation with the activity of the faults in the area (especially the Cajamarca fault, trending NE-SW) and with the presence of the hydrothermal system of the volcano that may include the infiltration of meteoric waters, the ascent of waters stored in the metamorphic complex or from the magmatic reservoir of the volcano, which is much deeper and outside the study area.

Figure 4.
Tested 2D models for a density configuration of a) the N-S transect and b) the W-E transect. See Figure 4b for profile location over the Bouguer Anomaly

Red dots: Block model editing vertices.

5. Discussion

5.1 Geological interpretation of the Bouguer anomaly

Based on the analysis of the Bouguer anomaly and the 2D forward models, we determined that the dome sector is less dense than its surrounding subsoil. This can be explained by alterations caused by the hydrothermal system associated with the Cerro Machin volcano and the interactions between its different parts. Piedrahita et al. (2018) suggested the presence of aquifers in the Cajamarca Complex, which, given the levels of shearing and fracturing of these rocks, would allow the movement and storage of confined water. Likewise, Cerpa (2018) stated that the hydrothermal system of the Cerro Machin volcano is characterized by the interaction of meteoric water that infiltrates from the surface to the basement, as well as by endogenous waters coming from a deep magmatic reservoir. The transport of these waters is favored by the permeability of the existing rocks (especially recent lava flows) and the fractures associated with the activity of the Cajamarca fault, which suggests the presence of different “mixing zones” between all these types of water. One of these zones stands out, located just below the area of the Cajamarca fault (see Figure 33 of Cerpa, 2018). This suggests that fracturing and infiltration of hydrothermal fluids in the lower part of the Machin dome considerably affected its density, which is reflected in the observed Bouguer anomaly.

Last, the higher values of the Bouguer anomaly in the southern part of the study area (Figure 3b) can also be explained by the hydrothermal processes associated with the volcano. This point is close to two thermal springs, named “Piscinas” and “Estatuas” by Inguagiatto et al. (2016). The first of these sources is characterized by relatively high concentrations of chlorides and silica, while the second is mainly distinguished by the presence of large amorphous silica concretions along a zone of structural discontinuity through which hydrothermal fluids ascend (Inguagiatto et al., 2016). Based on the presence of these concretions, we infer that due to the large accumulation of silica precipitated from the infiltrated fluids and considering that this new material fills the spaces left by the original porosity of the deposited rocks, there might have been a slight increase in the average density of the subsoil in this sector. The effects of these concretions, nonetheless, should be further investigated due to the lack of evidence or studies regarding this matter both in the study area and in other similar tectonic environments.

5.2 Geological interpretation of the Oasis Montaj models

Figure 5a and 5b show that the earthquakes located within our study area overlapped with those over the Oasis Montaj models. The earthquake locations were provided by the Servicio Geológico Colombiano and correspond to the period from 2004 to 2019. The obtained densities for all modeled units coincide with what has been reported in the literature for similar types of volcanoes (dacitic to andesitic stratovolcanoes with domes) outside Colombia and the Cerro Machin. For example, Kueppers et al. (2005) made density measurements for several dry rock samples from the deposits that were generated during the 1990-1995 eruptive cycle of Mount Unzen in Japan. The authors found a bimodal density distribution with peaks of 2000 ± 100 kg/m3 and 2300 ± 100 kg/m3 and reported that the density values could vary between 1400 and 2400 kg/m3. Similarly, Hoblitt and Harmon (1993) found a bimodal density distribution for the dacitic rocks from the 1980 eruption of Mount St. Helens, with peaks of 1600 and 2300 kg/m3. Furthermore, Acosta (2019) calculated the bulk density for several rock samples associated with the activity of three volcanoes in Colombia (Azufral, Doña Juana and Cerro Machin). For the Cerro Machin samples, the author found bulk density values between 2000 and 2300 kg/m3.

The above studies show that the chosen range of density values is plausible with the common geology that is associated with the activity of dacitic volcanoes similar to Cerro Machín. It also opens the opportunity for future investigations and modeling regarding the density of rock samples within the study area to obtain more complete data and more detailed constraints for the models. This is important considering the following: 1) the nonunique solution principle, which indicates that different geological distributions can be the best fit for a set of gravimetric data, and 2) even though there is a coincidence between what is reported here and what was reported in previous studies, there are some other studies where much lower densities were found for the volcanoclastic sequences. For example, Rueda et al. (2013) reported density values that ranged from 920 to 960 kg/m3 for the volcanoclastic products of a 31 ka Plinian–sub-Plinian eruption at Tláloc Volcano in Mexico. Similarly, Arce et al. (2005) reported samples from a volcanoclastic deposit associated with a 12.1 ka Plinian–sub-Plinian eruption of Volcán Nevado de Toluca (Mexico), whose density variation was established from 700 to 820 kg/m3.

These results coincide with the model proposed by Londoño (2011). By performing seismic P-wave velocity tomography to obtain the internal structure of Cerro Machin down to a depth of 4 to 5 km, this author determined that the dacitic dome reached a depth up to 3 km (from topographic relief), which explained why the seismicity was generally located at these depths. Furthermore, Londoño found that although the inner parts of the dome were compact enough to avoid fracturing and thus the circulation of fluid material to the shallowest parts of the system, two zones extending E and SW of the dome featured lower P-wave velocities. This behavior was associated with the genesis of weakening zones in the deepest part of the

Figure 5
Tested 2D profiles, including the seismicity provided by the Servicio Geológico Colombiano. a) Model for the N-S transect and b) model for the W-E transect

Black dots: Projection of the location of seismic events recorded by the Servicio Geológico Colombiano during 2004-2019. Red dots: Block model editing vertices.

dome, which were produced by the hot, fluid material coming from the magmatic reservoir attempting to reach the surface. This would explain the fracturing in these sectors and, therefore, the general seismic environment.

Taking these results into account, we built a general 2D W-E geological schematic model (Figure 6). For this, we extrapolated the density distribution that we obtained with the gravity forward models. The depth in which this extrapolation reaches was decided based on the resolution depth that we calculated using the equation given by Musset and Khan (2000) for obtaining the depth of an anomalous irregular body. This equation relates the maximum gravity anomaly and the maximum slope of the gravimetric profile; therefore, we used the one that we obtained when modeling the W-E forward model for our calculation. Taking this into account, we obtained a depth of approximately 6.2 km, which is the limit that we defined for our forward models. The schematic model reaches 10 km, which we inferred based on the seismic and geological data available for the Cerro Machin volcano, indicating that there may be a magmatic reservoir outside our study area at the selected depth (León, 2019; Londoño, 2011).

A dacitic dome complex occupies the uppermost part of our model. This dome is located over a series of layers associated with the volcaniclastic sequences that were formed during the reported volcanic activity of the Cerro Machin volcano and whose density may have been diminished due to weathering produced by infiltration of meteoric waters. At a depth of 3 km

Figure 6.
Final proposed test model of the geological distribution based on the GM-SYS modeling results

The red dashed line represents the vertical limits of the study. Seismic hypocenters were furnished by the Servicio Geológico Colombiano (Manizales).

from the surface, we infer a transition to a more compact, less fragmented, slightly denser zone of the dacitic complex, just above the Cajamarca metamorphic complex, where the seismic activity associated with the volcanic activity is located. This transition corresponds to altered materials affected by the different geological processes described in this study. Thus, we associate the existence of this sector with the existence of a weakened zone provoked by the general fault activity associated with the Cajamarca fault, the interaction with the hydrothermal system of the volcano and the ascent of hot fluid material that may have risen to the surface through less compact material (as Londoño, 2011). This ascent would have been through a system of dikes and fractures that connects this sector with a magmatic reservoir located to the SE of our study area and at a much deeper depth (Londoño, 2011). The contrast between the volcaniclastic sequences and this weakened zone is what we infer to be the origin of the gravity anomaly.

This general test configuration model offers valuable results as a starting point for future investigation projects in the area. The obtained results may have been slightly affected by the limited quantity of acquired data, which negatively impacts the resolution of the models explicitly defining the origin, spatial distribution and magnitude of the gravimetric anomaly related to the Cerro Machin volcano. We suggest, for future research, to focus on improving the data resolution of the study area, minimizing the distance between stations/points. Further measurements, such as seismic tomography and rock sampling, would provide precise density value information over the study area.

6. Conclusions

Based on observed gravity data acquired between March 16th and March 19th, 2020, we studied the local gravimetric anomalies in the Cerro Machin volcano. The minimum observations ranged from -87 mGal in the Machin dome area and reached up to −29 mGal in a measurement zone near hydrothermal sources in the southern area. These results allowed us to obtain a geological configuration with volcaniclastic sequences whose density values varied between 2500 kg/m3 and 2600 kg/m3 and a dacitic complex with an average density of 2300 kg/m3 over a less altered volcanic conduit with a suggested density of 2400 kg/m3 between the Machin dome and the metamorphic basement (Cajamarca Complex). These density values coincide with what has been reported for similar volcanoes in other parts of the world and for Cerro Machin. We interpret our results with a general weakening of the dacitic dome complex due to weathering caused by infiltration of meteoric waters in the uppermost section, the ascent of hot fluid material, as well as the activity of the faults located in the study area (particularly the Cajamarca fault), and the hydrothermal system associated with the Cerro Machin volcano. All these factors contributed to the generation of a weakened zone in the deepest part of the Machin dome at the interface between the volcanic edifice and the metamorphic basement. Through a system of dikes and fractures, this weakened zone may be connected to another reservoir (the Cerro Machin magmatic reservoir) located to the SE of our study area and at deeper depths (Londoño, 2019).

For future studies, we would like to extend our study area beyond the Machin dome zone and have a broader record of the gravimetric anomaly in the influence zone of the Cerro Machin volcano and a better understanding of the associated geology. This study allowed us to obtain an initial model providing a first glimpse at the geological distribution associated with this volcanic system. This represents a starting point for further understanding this volcano and its associated processes.

7. Acknowledgment

This work was supported by the FAPA project (Fondo de Apoyo para Profesores Asistentes) No. INV-2019-63-1702 from the Universidad de los Andes, Colombia

8. Supplementary data

Supplementary data for this article can be found online at https://doi.org/10.32685/0120-1425/bol.geol.49.1.2022.619

References

Acosta Agudelo, L. (2019). Contrastes texturales en fragmentos juveniles de depósitos de flujos piroclásticos de erupciones vulcanianas en los andes del norte: casos de estudio Cerro Machín, Doña Juana y Azufral, Colombia [Undergraduate thesis]. Universidad de los Andes. https://repositorio.uniandes.edu.co/handle/1992/45549

Arce, J. L., Cervantes, K. E., Macías, J. L., & Mora, J. C. (2005). The 12.1 ka Middle Toluca Pumice: A dacitic Plinian–subplinian eruption of Nevado de Toluca in Central Mexico. Journal of Volcanology and Geothermal Research, 147(1-2), 125-143. https://doi.org/10.1016/j.jvolgeores.2005.03.010

Cerpa, A. C. (2018). Caracterización del potencial geotérmico a partir de análisis geoquímicos de fuentes termales, en el volcán Cerro Machín, Colombia [Undergraduate thesis]. Universidad EIA.

Dentith, M., & Mudge, S. T. (2014). Geophysics for the Mineral Exploration Geoscientist. Cambridge University Press.

Hall, M. L., & Wood, C. A. (1985). Volcano-tectonic segmentation of the northern Andes. Geology, 13(3), 203-207. https://doi.org/10.1130/0091-7613(1985)13<203:VSOTNA>2.0.CO;2

Hinze, W. J. (2003). Bouguer reduction density, why 2.67? Geophysics, 68(5), 1559-1560. https://doi.org/10.1190/1.1620629

Hoblitt, R. P., & Harmon, R. S. (1993). Bimodal Density Distribution of Cryptodome Dacite from the 1980 Eruption of Mount St. Helens, Washington. Bulletin of Volcanology, 55, 421-437. https://doi.org/10.1007/BF00302002

Inguagiatto, S., Londoño, J. M., Chacón, Z., Liotta, M., Gil, E., & Alzate, D. (2016). The hydrothermal system of Cerro Machín volcano (Colombia): New magmatic signals observed during 2011–2013. Chemical Geology, 469, 60-68. http://dx.doi.org/10.1016/j.chemgeo.2016.12.020

Kueppers, U., Scheu, B., Spieler, O., & Dingwell, D. B. (2005). Field-based density measurements as tool to identify pre-eruption dome structure: set-up and first results from Unzen volcano, Japan. Journal of Volcanology and Geothermal Research, 141(1-2), 65-75. https://doi.org/10.1016/j.jvolgeores.2004.09.005

Laeger, K., Halama, R., Hansteen, T., Savov, I. P., Murcia, H. F., Cortés, G. P., & Garbe-Schönberg, D. (2013). Crystallization conditions and petrogenesis of the lava dome from the ~900 years BP eruption of Cerro Machín Volcano, Colombia. Journal of South American Earth Sciences, 48, 193-208. https://doi.org/10.1016/j.jsames.2013.09.009

León, A. (2019). Estudio de Coda Q para sismos volcano-tectónicos en el volcán Cerro Machín (Colombia) [Bachelor thesis]. Universidad de los Andes.

Londoño, J. M. (2011). 3D seismic tomography of Cerro Machín Volcano, Colombia. In XIV Latinamerican Congress of Geology. Medellín, Colombia.

Londoño, J. M., & Kumagai, H. (2018). 4D seismic tomography of Nevado del Ruiz Volcano, Colombia, 2000–2016. Journal of Volcanology and Geothermal Research, 358, 105-123. https://doi.org/10.1016/j.jvolgeores.2018.02.015

Méndez Fajury, R. A., Cepeda Vanegas, H., Monsalve Bustamante, M. L., Murcia Leal, L. A., & Núñez Tello, A. (1996). Volcán Cerro Machín, departamento del Tolima, Colombia: Pasado, presente y futuro. Instituto Nacional de Investigaciones Geológico Mineras (Ingeominas). https://miig.sgc.gov.co/Paginas/Resultados.aspx?k=110020101010002305000000000

Méndez, R. A. (2002). Formación Machín (Qm), cordillera Central, departamento del Tolima. In Catálogo de las Unidades Litoestratigráficas de Colombia. Instituto Colombiano de Geología y Minería (Ingeominas).

Monsalve, M. L., & Pulgarín, B. (1993). Mapa preliminar de amenaza volcánica potencial del volcán Puracé: Memoria explicativa. Revista Ingeominas, 2, 3-27.

Mosquera, D., Núñez, A., & Vesga-Ordóñez, C. J. (1982). Geología de la Plancha 244 Ibagué. Instituto Colombiano de Geología y Minería (Ingeominas). https://miig.sgc.gov.co/Paginas/Resultados.aspx?k=130100101010024375000000000

Musset, A. E., & Khan, M. A. (2000). Looking Into the Earth: An Introduction to Geological Geophysics. Cambridge University Press.

Northwest Geophysical Associates, Inc. (2004). Modeling Concepts. In GM-SYS Gravity/Magnetic Modeling Software: Vol. 4.9.

OpenTopography. (2021). Copernicus GLO-90 Digital Surface Model. https://doi.org/10.5069/G9028PQB

Palacios, N. (2017). A Matlab code to compute gravity measurements’ terrain corrections in Puracé (Cauca), Colombia [Bachelor thesis]. Universidad de los Andes, Bogotá, Colombia.

Piedrahita, D. A., Aguilar-Casallas, C., Arango-Palacio, E., Murcia, H., & Gómez-Arango, J. (2018). Estratigrafía del cráter y morfología del volcán Cerro Machín, Colombia. Boletín de Geología, 40(3), 29-48. http://dx.doi.org/10.18273/revbol.v40n3-2018002

Rueda, H. (2005). Erupciones plinianas del Holoceno en el volcán Cerro Machín, Colombia: Estratigrafía, Petrografía y Dinámica Eruptiva [Master Thesis]. Universidad Nacional Autónoma de México.

Rueda, H., Macías, J. L., Arce, J. L., Gardner, J. E., & Layer, P. W. (2013). The ~ 31 ka rhyolitic Plinian to sub-Plinian eruption of Tlaloc Volcano, Sierra Nevada, central Mexico. Journal of Volcanology and Geothermal Research, 252, 73-91. https://doi.org/10.1016/j.jvolgeores.2012.12.001

Vajda, P. (2016). Recent Developments and Trends in Volcano Gravimetry. In K. Nemeth, Updates in Volcanology – From Volcano Modelling to Volcano Geology. IntechOpen. http://dx.doi.org/10.5772/63420

Wynn, J., Dzurisin, D., Finn, C. A., Kauahikaua, J. P., & LaHusen, R. G. (2006). Applications of Geophysical Methods to Volcano Monitoring. 19th Symposium on the Application of Geophysics to Engineering and Environmental Problems. Seattle, United States.