1. Introduction

Ecuador is a country located in a convergence zone between the Nazca oceanic plate and the South American continental plate. This subduction process began in the Early Jurassic, resulting in continental-scale deformation, as well as increased seismic and volcanic activity (e.g., Michaud et al., 2009; Bilek, 2010; Schütte et al., 2012). Subsequently, the extensive calc-alkaline continental arc of Misahuallí developed in the Middle Jurassic, accompanied by the intrusion of type I batholiths (e.g., Aspden et al., 1992).

The Abitagua Batholith (Fig. 1) was first described by Colony and Sinclair (1932) while studying igneous and metamorphic rocks in eastern Ecuador. It represents a relatively understudied area due to its difficult access and location along protected zones: the Colonso-Chalupas Reserve and the Llanganates National Park. Consequently, scientific investigations of the batholith have relied on samples taken from riverbanks and areas where access is relatively easier (Aspden and Litherland, 1992; Litherland et al., 1994; Ruiz, 2002). In the early 1990s, gold was discovered in samples taken from rivers sourced from the Abitagua Batholith and potentially associated with some younger intrusions that are part of the Abitagua Batholith itself (Aspden and Litherland, 1992). Due to the field limitations, the use of indirect, geophysical methods was encouraged.

Fig. 1. Regional map: location of the study area within the Sub-Andean Zone of Ecuador. Local map: morphotectonic scheme adapted from Litherland et al. (1994), showing the Loja and Salado terrains and the Amazonian Craton, with the Abitagua Batholith in between. The tomography sections performed in this study are marked with dashed lines. Sections AA’, BB’, CC’, DD’, EE’, FF’ have an azimuth of  N90E, while section GG’ has an azimuth of  N25E. The purple rectangle depicts the area of horizontal cuts.

Seismic tomography is a geophysical technique useful for imaging the subsurface of the Earth in three dimensions using seismic waves from natural earthquakes (e.g., Thurber and Ritsema, 2007; Romanowicz, 2021). Some notable examples where this technique was used included: the description of the upper crust beneath the Chalupas caldera (Paredes and Araujo, 2021), the geodynamic impact due to the subduction of the Nazca plate beneath Ecuador (Araujo et al., 2021), and the location of a weathering front around a granite (Wang et al., 2019).

The objective of this study is to provide a regional characterization of the Abitagua Batholith by applying three different seismic velocity models (Vp, ΔVp, Vp /Vs) derived from the inversion of travel time data provided by Ecuadorian survey networks RENSIG (Red Nacional de Sismógrafos del Instituto Geofísico), RENAC (Red Nacional de Acelerógrafos), and ROVIC (Red de Observatorios Vulcanológicos), as well as some stations of the Colombian Geological Service near the border with Ecuador. In doing so, structures indicative of potential fluid escape pathways during magmatic intrusion can be recognized at depth. These structures can be related to the gold-bearing intrusions theorized to be the source of alluvial gold currently exploited on the banks of the rivers that drain the batholith.

2. Geological context

The Sub-Andean Zone of Ecuador is part of the Amazon Basin (Fig. 1). It is structured by reverse faults with a ~N-S orientation (Bès de Berc et al., 2004) consequence of tectonic inversion related to the Late Triassic-Early Jurassic Tetian rift system (Baby et al., 2004). The Sub-Andean is divided into three morphotectonic units: the Napo Uplift, the Pastaza Depression, and the Cutucú Range (Bès de Berc et al., 2004; Fig. 1).

The Abitagua Batholith is part of the Sub-Andean Zone (Fig. 1). It lies to the east of the Cosanga Fault (Litherland et al., 1994). Petrographically, it has been described as a rose-colored biotite monzogranite, ranging from medium- to coarse-grained to megacrystic K-feldspar-rich (Aspden et al., 1992; Litherland et al., 1994). It constitutes part of the basement of the Sub-Andean Zone and represents the highest non-volcanic relief in the region, with an elevation of ~2,700 m a.s.l. The batholith measures ~120 km long and ~12-15 km wide (Litherland et al., 1994; Fig. 1).

The age of the Abitagua Batholith is reported between 178±7 Ma (K-Ar-biotite; Herbert and Pichler, 1983) and 173±1.3 Ma (U-Pb-zircon; Spikings et al., 2015), so it is well constrained to the Middle-Late Jurassic transition, being therefore part of the intrusive suite of the Misahuallí arc (Romeuf et al., 1995). The geochemical composition of the batholith indicates it is a type-I granite, related to other intrusive bodies nearby such as the Zamora and Rosa Florida batholiths (Hall and Calle, 1982; Aspden et al., 1992; Litherland et al., 1994), and consequence of a change in the subduction dynamics since the Middle Jurassic in the region (Jaillard et al., 1990; Litherland et al., 1994).

3. Methodology

3.1. Stations and data

The tomography model used in this study was obtained from Araujo et al. (2021). These authors used data from 66 stations from the RENSIG seismometer network, 100 stations from the RENAC accelerometer network, and 85 stations from the monitoring networks installed around the Ecuadorian volcanoes (ROVIG), besides shared seismic data from the Colombian Geological Service (14 stations). These networks have operated digitally since 1988 and provide an extensive seismic database on Ecuador. The data period ended in April 2016, before the Mw 7.8 earthquake in the subduction zone north of Ecuador. From the seismic tomography study by Araujo et al. (2021), no new models of seismic velocities have been developed in our area of ​​interest.

The 265 stations used to obtain the velocity models from are shown in figure 2. Araujo et al. (2021) obtained a velocity model for a ~6° x 5° area (large rectangle in Fig. 2), with depths of up to 240 km. The area of the Abitagua Batholith is much smaller (inner rectangle in Fig. 2). Although there is not a high density of stations in the study area, the resolution of the model is relatively high due to the intense seismic activity recorded between 1988 and 2016 (Fig. 2). In consequence, 58,060 seismic events were thus manually retrieved (Fig. 2), with 641,036 and 215,134 arrival times calculated for P and S waves, respectively (Araujo et al., 2021).

Fig. 2. Seismic data from Araujo et al. (2021) and used in this study. The gray triangles represent the seismic stations. The seismic events are on a color scale for different depths: red, <35 km; yellow, 35-75 km; blue, 75-150 km; and black, 150-250 km.  The outer rectangle shows the region where the model was solved. The inner rectangle shows the area of ​​the Abitagua Batholith.

The a priori model from which the tomography process was based considered a one-dimensional model for the P-wave velocity with no prior geometry imposed on the crust or the Nazca slab (Araujo et al., 2021). The a priori model for the S wave was then obtained by using the relation Vp /Vs= √3.

3.2. Tomography software

The velocity models of Araujo et al. (2021), and used in this study, were obtained through the INSIGHT software (Potin, 2016). This software was developed to solve seismic tomography problems in the crust and upper mantle and is based on a stochastic resolution of the inverse problem (e.g., Tarantola and Valette, 1982; Valette, 2012).

In the stochastic resolution of inverse problems, the data d and the model parameters m are random vectors that maintain a functional relationship g: d=g(m). This approach quantifies the parameter information with a measure expressed through probability distribution functions. The data vector d contains information on the P-wave arrival times and the S-to-P wave travel time difference as d^P=t_P  and d^S-P=t_S-t_p, respectively. On the other hand, the model vector m is composed of parameters obtained after the inversion algorithm, namely: P-wave velocities (v_p ), the ratio between P- to S-wave velocities (Vp /Vs), seismicity hypocenters (x), earthquake origin times (t₀ ), and seismic station time delays (Δ^P, Δ^S-P ).

The errors in the data are expressed in the covariance matrix C_d , and those of the model parameters are described in the matrix C_m. The functional relationship takes the form of the matrix G. Then, the seismic tomography software INSIGHT solves this linear system at each iteration k (Eq. 1):

(1) 

Where:       

Araujo et al. (2021) obtained the velocity model by using a grid of 5 km in the horizontal direction and 2 km in the vertical direction. This model was resolved over the entire region of Ecuador and gave a total of 2,342,515 inversion nodes. For the errors of the inversion parameters, they took initial values ​​of 0.75 km/s for Vp, 0.15 for Vp/Vs, 30 km for the event locations (in any direction), and 1,000 s for the origin time of the earthquakes. The parameters governing the inverse problem’s regularization define a damping length and a smoothing length. The optimal values ​​of these parameters were found by using the L-curve method, resulting in 6 and 35 km for the damping and smoothing lengths, respectively.

After resolving the inverse problem with INSIGHT, the final data set comprised 335,498 P- and 111,457 S-wave arrival times. Thus, the number of seismic events was reduced from 58,060 to 25,410 (Araujo et al., 2021).

3.3. Resolution tests

In this study, we used the resolution operator criterion (Vergely et al., 2010) to determine the spatial resolution of the velocity models of Araujo et al. (2021) in the Abitagua Batholith area. This technique quantifies how well the inversion solution fits the model features and allows regularization control to obtain more physically stable solutions.

For the stochastic solution (Valette, 2012), the resolution operator R is calculated as (Eq. 2):

(2) 

If the grid size of the inverse problem is vast, it is not practical to calculate the operator R at all grid points. Instead, the restitution index (RI) is far preferable. RI is defined as an average of R for the parameter p at a given point  (Vergely et al., 2010). The parameter p, in our case study, is the value of the seismic wave velocity at a point of the resolution grid. The average of R is calculated over the points  around  (Eq. 3):

(3) 

RI approaches unity when the parameter p reaches the actual value of velocity, which indicates the good resolution of the model. On the contrary, when p is far from the actual velocity value, RI leans towards zero, indicating poor resolution (Vergely et al., 2010; Valette, 2012).

The model resolution results are shown in figure 3 along some of the profiles from Figure 1. In general, the resolution is high in the Abitagua Batholith area. The minimum RI is 0.84 and the maximum even reaches the value of unity (Fig. 3D). This high resolution is explained due to the significant seismicity recorded in the area for the 1988-2016 period (Fig. 2).

Fig. 3. Cross-sections of the RI value in the study area. The sections correspond to those indicated in figure 1. The RI value is represented with a color scale and isocontour lines. The vertical and horizontal scales are in kilometers.

3.4. Tomography cross-sections

The cross-sectional cuts were generated in Matlab. The sections use riverheads as references, and from north to south they are (Fig. 1): Mulatos River- Jatunyaku River (FF’), Illocullín River (EE’), Piatúa Blanco River (DD’), Piatúa River (CC’), Chontayacu River (BB’), and Yurak Yaku Grande River (AA’). All sections have a length of 40 km and a depth of 70 km with an azimuth of N90E, aiming to distinguish the Abitagua Batholith from surrounding lithologies (Fig. 1). An additional, oblique profile was created to observe changes along the internal structure of the batholith. This last profile has a length of 120 km and a depth of 70 km with an azimuth of N25E (Fig. 1). The horizontal cuts were performed within the quadrant [78.26° W 1.59° S 77.60° W 0.55° S] in order to observe the velocity changes related to different terrains (Fig. 1).

The design software must define the interpolation length to draw the tomography images. For this smoothing length, we chose 0.3 km in vertical and horizontal directions. A file containing the seismic event hypocenters resulting from the tomography inversion was also provided.

4. Results and discussion

The tomography results are presented in figures 4-7. Vertical sections AA’, BB’, and CC’ are shown in figure 4. Sections DD’, EE’, and FF’ are in figure 5, and section GG’ is shown in figure 6. The horizontal section is depicted in figure 7. Each section presents results from the relative P-wave velocity model (ΔVp), absolute P-wave velocity model (Vp), and the ratio between P-wave and S-wave velocities (Vp/Vs). Additionally, all figures include the seismic events occurring within 12 km perpendicular to the tomography cuts.

Fig. 4. Seismic model results: ΔVp (in percentage form), Vp (km/s), and Vp/Vs (dimensionless) along sections A: AA’, B: BB’, and C: CC’ (see section location in Fig. 1). Along the surface, the main morphotectonic units (after Litherland et al., 1994) are indicated. Some relevant anomalies (dotted contours) and seismic event hypocenters (blue dots) are also shown.

Fig. 5. Seismic model results: ΔVp (in percentage form), Vp (km/s), and Vp/Vs (dimensionless) along sections A: DD’, B: EE’, and    C: FF’ (see section location in Fig. 1). Along the surface, the main morphotectonic units (after Litherland et al., 1994) are indicated. Some relevant anomalies (dotted contours) and seismic event hypocenters (blue dots) are also shown.

Fig. 6. Seismic model results. A: ΔVp (in percentage form). B: Vp (km/s), C: Vp/Vs (dimensionless). The intersection of the profiles AA’ to FF’ is marked with black arrows. In A, the regions of lower and higher rock compaction are shown by blue segmented curves. In B, black segmented curves refer to two seismic anomalies, with potential fluid migration pathways depicted as black arrows. Seismic event hypocenters (blue dots) are also shown.

Fig. 7. Map view tomography cross-sections. A and B show the absolute Vp model. C and D show the relative P wave velocity model. Sections A and C were obtained 0.5 km below the geoid while sections B and D 3 km below the geoid. The boundaries of the Salado Terrain, the Abitagua Batholith, and the Amazonian craton are shown as black lines.

The ΔVp model illustrates variations in P-wave velocity across the tomography cross-sections. Seismic anomalies (positive and negative) in this model are presented in percentage form and may be indicative of crustal structures, such as lithological discontinuities (Koulakov, 2012, 2013; Pávez et al., 2016; Vargas et al., 2017). In consequence, all tomography profiles were overlaid with the morphotectonic segmentation of Litherland et al. (1994) to correlate the Abitagua Batholith with the metamorphosed Salado Terrain to the west and the Amazonian Craton to the east.

The ΔVp model results for sections AA’, BB’, CC’, DD’ (Figs. 4 and 5) show a surface variation of ΔVp ranging from -5 to -15% and lower. Neither of these anomalies allow the Abitagua Batholith to be recognized. However, by overlying the geological map by Egüez et al. (2017), it is possible to define its boundaries within the  models. Therefore, we propose the region enclosed between -10% ≥ ΔVp ≥ -15% delineates the Abitagua Batholith with a horizontal extension of approximately 12 km, which is the average size reported in Litherland et al. (1994). The -15% isoline marks the contact with the Amazonian Craton, while the -10% isoline indicates the contact with the metamorphosed Salado Terrain.

In sections EE’ and FF’ (Fig. 5A, B) the surface values are distributed from -10 to 0%. It is worth to notice that  velocities within the area where the Abitagua Batholith is located according to profiles AA’ to DD’ vary when compared to those shown in Figure 4. In fact, after loading the Egüez et al. (2017) map, the Abitagua Batholith would be represented by -5% ≥ ΔVp ≥ -7% in section EE’ (Fig. 5B). Likewise, in the FF’ section (Fig. 5C), the 0% curve would delineate the contact with the Salado Terrain, from which a wedge with respect to the Abitagua Batholith along the Cosanga Fault can be interpreted, similar to the geological section shown in Litherland et al. (1994).

There is an increase in the ΔVp along the surface of the batholith from section DD’ (ΔVp < -10%) to section FF’ (-5% ≥ ΔVp ≥ 0%) (Fig. 5). This longitudinal velocity variation is evident in the  model of section GG’ (Fig. 6) and in the horizontal cuts (Fig. 7). The negative ΔVp anomalies in section GG’ represent areas where the structure of the batholith is less compact, indicating a zone where the rock is fractured (Martí et al., 2002; D’Auria et al., 2022). In contrast, the positive anomalies in the same section are due to the presence of more consolidated rock, which causes an increase in the absolute and relative values of P-wave velocity (Fig. 5A, B).

The absolute P-wave velocity values increase with depth in all profiles (Figs. 4, 5, 7B). The depth of the Abitagua Batholith can be estimated based on the  values for granites (4.6-4.8 km/s for weathered or fractured granites and 5.3-5.5 km/s for fresh, unfractured granites; Gelbke et al., 1989; Martí et al., 2002). From sections AA’ to EE’, the average depth of the batholith is at around 10 km (Figs. 4, 5A, 5B), while in section FF’ it decreases to about 2 km (Fig. 5C). This thinning is evident in Section GG’ north from section FF’ (Fig. 6B). At greater depths (higher Vp values), a structure that could represent the root of the intrusive is inferred at around the 6.8 km/s isoline at around 35 km depth.

In the Vp model result profiles, north from section CC’ there is a sudden increase in the number of seismic events recorded (Figs. 4C and 5). This rise in seismicity would be related to the seismic cluster of Pisayambo described in Araujo et al. (2009), located within the quadrant [78.5° W 1.3° S     78.0° W 0.9° S] and classified as the most active of the four seismic nests detected in Ecuador (Araujo et al., 2009; Araujo, 2012). The Pisayambo cluster is associated with an unknown fault segment of the Chingual-Cosanga-Puná-Pallatanga fault system (CCPP), later named as the Pisayambo Lagoon Fault (PLF) by Champenois et al. (2017).

In the map view sections, the 0.5 and 3 km depth slices show the ΔVp transition zone that separates weathered from fresh rock (Fig. 7C, D). To the west, the Salado Terrain is characterized by positive ΔVp  anomalies, consistent with those reported in Paredes and Araujo (2021). The Amazonian Craton, on the other hand, exhibits a pattern of positive ΔVp values in the north and negative in the south, similar to the batholith.

The Vp/Vs model results can be used to determine the crustal architecture associated with volcanoes and shallow intrusions (e.g., Hua et al., 2019). At greater depths, on the other hand, the S-wave velocity values are typically low in the presence of partial melts, so the Vp/Vs ratio is directly affected, resulting in higher ratios sometimes associated with the presence of magmatic reservoirs (Koulakov, 2012, 2013; Ferri et al., 2016; Vargas et al., 2017; Yang et al., 2023). In the profiles obtained in the present study, it is possible to observe a Vp/Vs anomaly at around 1.77-1.79, at ~45 km depth and below (Figs. 4 and 5). This distribution is more evident in section GG’ (Fig. 6C), in which even a second anomaly can be observed. The first one locates below sections AA’ to FF’ between ~40-60 km depth, while ~80 km to the north, a second, much larger anomaly encloses  Vp/Vs ratios above 1.8.

Model results suggest that the high Vp/Vs anomalies are correlated with basaltic melts generated in the mantle, which can ascend and stagnate at the boundary between the mantle and the crust (Koulakov, 2012, 2013; Ferri et al., 2016). In this location, magmatic processes such as fractional crystallization, crustal assimilation, and magma mixing can take place as well (e.g., Zhao et al., 2019; Bugueño et al., 2022). Additionally, in volcanic terrains, the seismicity around these anomalies is generally attributed to magma movement and hydrothermal fluid ascent mechanisms (Poli and Schmidt, 1995; Koulakov, 2013; D’Auria et al., 2022; Yang et al., 2023). The migration of hydrothermal fluids containing metals results in the deposition of minerals through effective precipitation mechanisms (e.g., Stoffell et al., 2004; Zhao et al., 2019; Zhang et al., 2021; Bugueño et al., 2022). The interaction between intrusions and hydrothermal fluids results in porphyries rich in Au and Cu (Zhao et al., 2019; Bugueño et al., 2022). After shallow emplacement followed by weathering and erosion, fragments of these mineralized bodies are transported and deposited in nearby rivers, generating alluvial gold deposits.

Based on the seismicity observed between the AA’ and BB’ sections (Fig. 4A, B), the ascent of fluids from the ~1.8 Vp/Vs anomaly is suggested. Similar patterns can be observed in the CC’ and DD’ sections (Figs. 4C and 5A). These fluids sometimes result in mineral-laden intrusions, and would explain the presence of gold in rivers and streams around the Yurak Yaku Grande River (Aspden and Litherland, 1992) and the Napo River (Barragan et al., 1991). On other hand, in the GG’ section, the transition zone between positive and negative  values falls ~5 km north from the Jatunyaku River (FF’ profile) (Fig. 6A). In this same section, the Vp/Vs model combined with hypocenter locations (Fig. 6C) may be indicative of potential fluid escape pathways from two magmatic reservoirs. At shallow depths, both pathways seem to converge around the mentioned transition zone. In this place, alluvial gold deposits have been reported (Barragan et al., 1991; Aspden and Litherland, 1992), with their origin probably from still unexplored, post-Abitagua porphyry-type intrusions.

5. Conclusions

Seismic velocity models in the Abitagua Batholith area, central Ecuador, report a transitional zone from positive to negative ΔVp values ~5 km north from the Jatunyaku River. This zone would be the place where potential post-Abitagua mineralized intrusions are more likely.

According to the Vp/Vs model results, two seismic anomalies were determined at depths of >40 km beneath the Abitagua Batholith, associated in this study with magmatic reservoirs at the mantle-crust interface. When combined with seismic hypocenter locations, magma ascent pathways are inferred, which at shallow depths seem to converge at the  ΔVp transitional zone.

If confirmed, these potential intrusions could be the source of the alluvial gold documented in several rivers that drain from the batholith. A more robust seismic network around the batholith area would therefore be necessary to improve model resolutions and test the presence of these gold-bearing intrusions at local scales.

Acknowledgments
The authors acknowledge the review and corrections of three anonymous reviewers who helped improve the original manuscript. The editor also provided valuable observations and corrections. All tomography results in this study were performed using the GRICAD infrastructure (https://gricad.univ-grenoble-alpes.fr), supported by Grenoble research communities.

References

Araujo, S. 2012. Resolución del problema directo de tomografía sísmica: trazado de rayos sísmicos en el clúster de Pisayambo. La Granja 15 (1): 48-54. http://dx.doi.org/10.17163/lgr.n15.2012.04

Araujo, S.; Troncoso, L.; Ruiz, M. 2009. Relocalización por dobles diferencias del cluster sísmico de Pisayambo. La Granja 10 (2): 27-34. https://doi.org/10.17163/lgr.n10.2009.03

Araujo, S.; Valette, B.; Potin, B.; Ruiz, M. 2021. A preliminary seismic travel time tomography beneath Ecuador from data of the national network. Journal of South American Earth Sciences 111: 103486. https://doi.org/10.1016/j.jsames.2021.103486

Aspden, J.A.; Litherland, M. 1992. The geology and Mesozoic collisional history of the Cordillera Real, Ecuador. Tectonophysics 205 (1-3): 187-204. https://doi.org/10.1016/0040-1951(92)90426-7

Baby, P.; Rivadeneira, M.; Barragán, R. 2004. La Cuenca Oriente: Geología y petróleo. Institut Français d’études Andines. IRD Éditions: 296 p. https://doi.org/10.4000/books.ifea.2971

Aspden, J.A.; Fortey, N.; Litherland, M.; Viteri, F.; Harrison, S.M. 1992. Regional S-type granites in the Ecuadorian Andes: Possible remnants of the breakup of western Gondwana. Journal of South American Earth Sciences 6 (3): 123-132. https://doi.org/10.1016/0895-9811(92)90002-G

Barragán, G.; Ortiz, C.; Merlyn, M. 1991. Placeres auríferos en el Ecuador. In Simposio internacional sobre yacimientos aluviales de oro (Hérail, G.; editor). Actas: 23-37. La Paz.

Bès de Berc, S.; Baby, P.; Soula, J.-C.; Rosero, J.; Souris, M.; Christophoul, F.; Vega, J. 2004. La superficie mera-upano: Marcador geomorfológico de la incisión fluviatil y del levantamiento tectónico de la zona subandina ecuatoriana. In La Cuenca Oriente: Geología y petróleo (Baby, P.; Rivadeneira, M.; Barragán, R.; editors). IRD Éditions: 153-167. https://doi.org/10.4000/books.ifea.3009

Bilek, S.L. 2010. Invited review paper: Seismicity along the South American subduction zone: Review of large earthquakes, tsunamis, and subduction zone complexity. Tectonophysics 495 (1-2): 2-14. https://doi.org/10.1016/j.tecto.2009.02.037

Bugueño, F.; Calle-Gardella, D.; Comte, D.; Reyes-Wagner, V.; Ojeda, M.; Rietbrock, A.; Roecker, S. 2022. Subsurface Insights of the Maricunga Gold Belt through Local Earthquake Tomography. Minerals 12 (11): 1437. https://doi.org/10.3390/min12111437

Champenois, J.; Baize, S.; Vallee, M.; Jomard, H.; Alvarado, A.; Espin, P.; Ekström, G.; Audin, L. 2017. Evidences of surface rupture associated with a low‐magnitude (M_w 5.0) shallow earthquake in the Ecuadorian Andes. Journal of Geophysical Research, Solid Earth 122 (10): 8446-8458. https://doi.org/10.1002/2017JB013928

Colony, R.J.; Sinclair, J.H. 1932. Metamorphic and igneous rocks of Eastern Ecuador. Annals of the New York Academy of Sciences 34 (1): 1-53. https://doi.org/10.1111/j.1749-6632.1932.tb56973.x

D’Auria, L.; Koulakov, I.; Prudencio, J.; Cabrera-Pérez, I.; Ibáñez, J.M.; Barrancos, J.; García-Hernández, R.; Martínez van Dorth, D.; Padilla, G.D.; Przeor, M.; Ortega, V.; Hernández, P.; Peréz, N.M. 2022. Rapid magma ascent beneath La Palma revealed by seismic tomography. Scientific Reports 12 (1): 17654. https://doi.org/10.1038/s41598-022-21818-9

Egüez, A.; Gaona, M.; Albán, A. 2017. Mapa geológico de la República del Ecuador. Instituto Nacional de Investigación Geológico Minero Metalúrgico. Escala 1:1.000.000. Quito.

Ferri, F.; Burlini, L.; Cesare, B. 2016. Effect of partial melting on Vp and Vs in crustal enclaves from Mazarrón (SE Spain). Tectonophysics 671 : 139-150. https://doi.org/10.1016/j.tecto.2015.12.030

Gelbke, C.; Miranda, F.; Sattel, G. 1989. Results of a seismic transmission tomography survey at the Grimsel Rock Laboratory. The Log Analyst 30 (4): SPWLA-1989-v30n4a3

Hall, M.L.; Calle, J. 1982. Geochronological control for the main tectonic-magmatic events of Ecuador. Earth-Science Reviews 18 (3-4): 215-239. https://doi.org/10.1016/0012-8252(82)90038-1

Herbert, H.-J.; Pichler, H. 1983. K-Ar ages of rocks from the Eastern Cordillera of Ecuador. Zeitschrift der Deutschen Geologischen Gesellschaft 134 (2): 483-493. https://doi.org/10.1127/zdgg/134/1983/483

Hua, Y.; Zhang, S.; Li, M.; Wu, T.; Zou, C.; Liu, L. 2019. Magma system beneath Tengchong volcanic zone inferred from local earthquake seismic tomography. Journal of Volcanology and Geothermal Research 377: 1-16. https://doi.org/10.1016/j.jvolgeores.2019.04.002

Jaillard, E.; Soler, P.; Carlier, G.; Mourier, T. 1990. Geodynamic evolution of the northern and central Andes during early to middle Mesozoic times: a Tethyan model. Journal of the Geological Society 147 (6): 1009-1022. https://doi.org/10.1144/gsjgs.147.6.1009

Koulakov, I. 2012. Multiscale seismic tomography imaging of volcanic complexes. In Updates in Volcanology -A comprehensive approach to volcanological problems (Stoppa, F.; editor). InTech: 207-242. https://doi.org/10.5772/24653

Koulakov, I. 2013. Studying deep sources of volcanism using multiscale seismic tomography. Journal of Volcanology and Geothermal Research 257: 205-226. https://doi.org/10.1016/j.jvolgeores.2013.03.012

Litherland, M.; Aspden, J.A.; Jemielita, R.A. 1994. The metamorphic belts of Ecuador. British Geological Survey, Overseas Memoir 11. Keyworth: 147 p.

Martí, D.; Carbonell, R.; Tryggvason, A.; Escuder, J.; Pérez-Estaún, A. 2002. Mapping brittle fracture zones in three dimensions: high resolution traveltime seismic tomography in a granitic pluton. Geophysical Journal International 149 (1): 95-105. https://doi.org/10.1046/j.1365-246X.2002.01615.x

Michaud, F.; Witt, C.; Royer, J.-Y. 2009. Influence of the subduction of the Carnegie volcanic ridge on Ecuadorian geology: Reality and fiction. In Backbone of the Americas: Shallow Subduction, Plateau Uplift, and Ridge and Terrane Collision (Kay, S.M.; Ramos, V.A.; Dickinson, W.R.; Editors). Geological Society of America Memoir 204: 217-228. https://doi.org/10.1130/2009.1204(10)

Paredes, C.D.; Araujo, S. 2021. The use of seismic tomography to describe the upper crustal structure beneath the Chalupas Caldera, Ecuador. REM- International Engineering Journal 74 (2): 189-197. http://dx.doi.org/10.1590/0370-44672020740079

Pávez, C.; Tapia, F.; Comte, D.; Gutiérrez, F.; Lira, E.; Charrier, R.; Benavente, O. 2016. Characterization of the hydrothermal system of the Tinguiririca Volcanic Complex, Central Chile, using structural geology and passive seismic tomography. Journal of Volcanology and Geothermal Research 310: 107-117. https://doi.org/10.1016/j.jvolgeores.2015.11.018

Poli, S.; Schmidt, M. 1995. H₂O transport and release in subduction zones: Experimental constraints on basaltic and andesitic systems. Journal of Geophysical Research, Solid Earth: 100 (B11): 22299-22314. https://doi.org/10.1029/95JB01570

Potin, B. 2016. Les Alpes occidentales: tomographie, localisation de séismes et topographie du Moho. Ph.D. Thesis (Unpublished). Université Grenoble Alpes: 245 p.

Romanowicz, B. 2021. Seismic tomography of the Earth’s Mantle. In Encyclopedia of Geology (Second Edition) (Alderton, D.; Elias, S.A.; editors). Elsevier: 587-609. https://doi.org/10.1016/B978-0-08-102908-4.00169-7

Romeuf, N.; Aguirre, L.; Soler, P.; Feraud, G.; Jaillard, E.; Ruffet, G. 1995, Middle Jurassic volcanism in the Northern and Central Andes. Revista Geológica de Chile 22 (2): 245-259.

Ruiz, G.M.H. 2002. Exhumation of the northern Sub-Andean Zone of Ecuador and its source regions: A combined thermochronological and heavy mineral approach. Ph.D. Thesis (Unpublished). Swiss Federal Institute of Technology Zurich: 260 p.

Schütte, P.; Chiaradia, M.; Barra, F.; Villagómez, D.; Beate, B. 2012. Metallogenic features of Miocene porphyry Cu and porphyry-related mineral deposits in Ecuador revealed by Re-Os, ⁴⁰Ar/³⁹Ar, and U-Pb geochronology. Mineralium Deposita 47 (4): 383-410. https://doi.org/10.1007/s00126-011-0378-z

Spikings, R.; Cochrane, R.; Villagómez, D.; Van Der Lelij, R.; Vallejo, C.; Winkler, W.; Beate, B. 2015. The geological history of northwestern South America: from Pangaea to the early collision of the Caribbean Large Igneous Province (290-75 Ma). Gondwana Research 27 (1): 95-139. https://doi.org/10.1016/j.gr.2014.06.004

Stoffell, B. 2004. Metal transport and deposition in hydrothermal veins revealed by 213 nm UV laser ablation microanalysis of single fluid inclusions. American Journal of Science 304 (6): 533-557. https://doi.org/10.2475/ajs.304.6.533

Tarantola, A.; Valette, B. 1982. Generalized nonlinear inverse problems solved using the least squares criterion. Reviews of Geophysics 20 (2): 219-232. https://doi.org/10.1029/rg020i002p00219

Thurber, C.; Ritsema, J. 2007. Theory and Observations - Seismic Tomography and Inverse Methods. Treatise on Geophysics 1: 323-360. https://doi.org/10.1016/B978-044452748-6.00009-2

Valette, B. 2012. Non-parametric regularization of tomographic problems. In International Workshop on “Muon and Neutrino Radiography 2012”. CNRS Clermont-Ferrand. France.

Vargas, C.; Koulakov, I.; Jaupart, C.; Gladkov, V.; Gómez, E.; El Khrepy, S.; Al-Arifi, N. 2017. Breathing of the Nevado del Ruiz volcano reservoir, Colombia, inferred from repeated seismic tomography. Scientific Reports 7 (1): 46094. https://doi.org/10.1038/srep46094

Vergely, J.-L.; Valette, B.; Lallement, R.; Raimond, S. 2010. Spatial distribution of interstellar dust in the Sun’s vicinity, comparison with neutral sodium-bearing gas. Astronomy & Astrophysics 518: A31. https://doi.org/10.1051/0004-6361/200913962

Wang, W.; Chen, P.; Keifer, I.; Dueker, K.; Lee, E.-J.; Mu, D.; Jiao, J.; Zhang, Y.; Carr, B. 2019. Weathering front under a granite ridge revealed through full-3D seismic ambient-noise tomography. Earth and Planetary Science Letters 509: 66-77. https://doi.org/10.1016/j.epsl.2018.12.038

Yang, X.; Roman, D.C.; Haney, M.; Kupres, C.A. 2023. Double reservoirs imaged below Great Sitkin Volcano, Alaska, explain the migration of volcanic seismicity. Geophysical Research Letters 50 (11): e2022GL102438. https://doi.org/10.1029/2022GL102438

Zhang, Y.; Han, R.; Ding, X.; Wang, Y.; Wei, P. 2021. Precipitation reaction mechanisms of mineral deposits simulated with a fluid mixing model. Geofluids 2021: 8881677. https://doi.org/10.1155/2021/8881677

Zhao, Y.; Liu, H.; Qin, K. 2019. Discussion on Au transportation mechanism in melt-magma-fluid in porphyry Cu deposit-A case study from Jinchang porphyry Au (Cu) deposit, Heilongjiang province, China. Ore Geology Reviews 111: 102968. https://doi.org/10.1016/j.oregeorev.2019.102968