The use of seismic tomography to describe the upper crustal structure beneath the Chalupas Caldera, Ecuador

The Chalupas Caldera is a rhyolitic volcano located in the Eastern Cordillera to the Southeast of the Cotopaxi Volcano in the Ecuadorian Andes. It is supposedly fed by a huge magmatic chamber below the caldera. For studying the Chalupas magmatic chamber, the tomographic models obtained by a Tarantola-Valette inversion and studied previously in the entire Ecuadorian region were used. The theoretical basis of the tomography method based on a Bayesian solution of the inverse problem is introduced. The process of the inverse problem regularization is obtained by an L-curve scheme to achieve an optimal solution. The seismic tomography images zoomed in the Chalupas region are well correlated with the main geological terranes and faults documented for the basement of Eastern Cordillera. Based on the tomographic images obtained, it can be assumed that there is not a magmatic chamber of considerable magnitude placed close to the surface underneath the Chalupas Caldera, since it does not present seismic and tomographic evidence to be able to assume a possible volcanic eruption before long. Furthermore, the images allow identifying that the seismic and thermal activity is located beneath the Cotopaxi Volcano as a vertical anomaly beneath the Northeast of Cotopaxi Volcano that corresponds to the phreatomagmatic activity observed in 2015. Andes of Ecuador, Bayesian inversion.


Methodology
A Bayesian approach to solve the inverse problem of the tomography was used (Tarantola, 2005). It starts by defining a theoretical model with blocks in which m are the parameters and d the data. The relationship between them is the direct problem: d = g (m).
The data vector d contains the travel times of P waves and the differ-ences between the S and P wave travel times. The model vector m contains the parameters to be obtained after the inversion: P wave velocity, the ratio of P wave velocity to S wave velocity, the hypocenters of the seismicity, the origin times of the earthquakes, and the delay time in the seismic stations. The linear operator g is a functional theoretical relationship between the data and the parameters.
The Bayesian inversion (Tarantola and Valette, 1982) leads, as stated above, to a Tikhonov problem. This problem consists in searching for the minimum of a cost function in a regularization space. The cost function has the form: In Equation (1), T is the regularization operator that in the case of the Gaussian stochastic approach is the square root of the covariance operator C m -1/2 of the model. C d is the covariance matrix of the data.
To solve this problem numerically, we used the following quasi-Newton algorithm: The INSIGHT tomography software uses the LSQR algorithm (Paige and Saunders, 1982) and the matrix decomposition described in Monteiller et al., 2005: The operator Σ applied to any function f and a vector x with standard devia-tion, σ is defined as ( The correlation kernel has an exponential form (Vergely et al., 2010): With: to 8 km from the surface.
Also, based on geochemical and isotope data collected from the surface water, it cannot construct any reasonable theory that indicates the presence of a deep high-temperature aquifer or magmatic intrusions. For these reasons, geophysical techniques were recommended, especially more profound magnetotelluric (MT), and seismic of high resolution (Beate and Salgado, 2010).
One of the most developed geophysical techniques in recent years is seismic tomography. At present, it is a potent and essential tool to analyze 2D and 3D speed models in volcanic and seismically active areas (Thurber, 1992). The seismic tomography applied to magmatic systems gives vast information concerning one of the most famous volcanic hot spot calderas of the world: Yellowstone. The tomography images provided evidence for gas and magmatic sources beneath Yellowstone (Husen et al, 2004); the slab-plume interaction that generates the Yellowstone complex (Obrebski et al, 2010) and the spatial extent of the Yellowstone crustal magma (Farrell et al., 2014).
The advantage of this technique lies in the geodynamic context because the subduction systems provide excellent coverage of rays generated by deep seismicity related to the slab. Besides, a sufficient quantity and proper distribution of receivers on the surface allows studying volcanism processes in detail (Koulakov, 2012).
One of the most studied types of seismic tomography is travel time tomography, which considers the P and S wave travel time between the source and the receiver. After generating an extensive data set with the travel times for various natural events, the next step is to choose the most suitable inversion procedure. Inversion tries to calculate the velocity of the medium by knowing the propagation times of the seismic rays and is one of the most studied geophysics topics. For this study, the results of seismic tomography developed by Araujo (2016) using the INSIGHT software were applied. This method is based on a stochastic inversion (Tarantola and Valette, 1982) and leads to a Tikhonov problem, which consists of searching for the minimum of a cost function in a regularization space (Araujo et al., 2019). The INSIGHT tomography code uses the LSQR algorithm (Paige and Saunders, 1982) following the matrix decomposition proposed in Monteiller et al., 2005. INSIGHT software allows the use of P and S wave travel times of all seismic events available over the study zone, providing the maximal ray coverage of the study domain and stabilizing it by the "large" data set. The results obtained using INSIGHT software in the previous study (Araujo 2016) are analyzed in the present research to obtain an image of the Chalupas Caldera deep structure.
To obtain the P wave velocity model and the ratio between the P wave and S wave model, Araujo (2016) uses the data furnished by IGEPN (Instituto Geofísico de la Escuela Politécnica Nacional, Quito-Ecuador). There are 45941 seismic events registered in 256 seismic stations throughout Ecuadorian territory. The data cover the period from 1988 to April 2016 before the Pedernales subduction earthquake.
The physical standard deviation for the parameters is 750 for v P and 0.15 for v P ⁄ v S and we set the vertical smoothing length ξ V at 15 km. The two horizontal lengths are considered equal. To evaluate the horizontal smoothing lengths and the damping length, the inverse problem must be regularized. To achieve the regularization, we used the L-curve strategy, where the optimization value is the corner of the curve (Hansen, 1992). Different values of ξ H and ξ 0 were plotted, the data fitting RMS versus the norm for the v P , v P ⁄ v S models, and the events hypocenters displacement. In this case, the L-curve corner is ξ 0 = 6 km and ξ H = 35 km, as seen in Figure 1.
To test the solution's stability, the tomography results obtained in Araujo et al. (2019) were also analyzed. In this study, the database is complemented with 14917 seismic events from April 2016 to May 2019 (IGEPN). The resulting models in the region of interest do not have noticeable changes, and for this reason, the models used in the next sections are reported in Araujo (2016) seismic tomography.

Study site: Chalupas Caldera
The INSIGHT software solutions obtained in (Araujo, 2016) and in (Araujo et al., 2019) provide three text files. These files are the v P , ∆v P (percentage variation) and v P ⁄ v S models with data in each node of the model space. The models have a grid with 119 nodes for x (longitude), 155 nodes for y (latitude), and 127 nodes for z (depth).
The cell size of the model is 5km in the horizontal direction and 2km in the vertical direction.
For generating a cross-section, it is necessary to supply the geographic coordinates of the location (longitude and latitude), the azimuth, and length. The referred to coordinates are presented as the 0 in the cross-section, the azi-muth gives the direction, and the length represents the total distance covered by the cross-section. For establishing the seismic events to show up in the crosssection, it is necessary to provide a distance in km. This distance performs the limits of a seismic box parallel to the azimuth of the cross-section to get the events registered within.
The Chalupas Caldera is located in the Eastern Cordillera, southeast of the Cotopaxi Volcano. The simplified stratigraphic sequence of metamorphic (Paleozoic-Triassic-Jurassic) basement of the Eastern Cordillera consists of four well-differentiated terranes (Litherland et al., 1994). Figure 2 shows these terranes traversing (westward) The operator represents the correlation between two points x i and x i´ in the grid model. There are three cor-relation, or smoothing lengths ξ i , for each space direction. The algorithm also renormalizes the standard deviation using a damping length ξ 0 to maintain constant the information inside the inversion volume ξ 1 ξ 2 ξ 3 : Figure 2 -The left map shows the position of the Chalupas Caldera. The Geological map to the right shows the basement of Eastern Cordillera and the Chalupas Caldera. The map is modified from Litherland (1994).
There are three phases well-differentiated in Chalupas's eruptive record in the Pleistocene period: the pre-collapse, the dominant eruptive phase, and the post-collapse. The pre-collapse eruptions consist of a large amount of mainly an-from the Jurassic batholith (Abitagua) and crossing the Cosanga Fault, the deformed Salado Basin (Salado Terrane) consisting of turbiditic metasedimentary rocks and mafic lavas (Upano Unit), in contact with the I-type calc-alkaline Azafrán Batholith (145-140Ma) (Spikings et al.,2015). Continuing the boundary of Salado Terrane westwards is the Llanganates Fault and further west a poorly differentiated zone between metamorphosed Paleozoic metasedimentary rocks of Chiguinda and Agoyán Units with the Triassic Tres Lagunas Unit Trl (240 -250Ma); this is an S-type granite characterized by the presence of "blue quartz" (Paul et al., 2018). The association of these Paleozoic and Triassic rocks constitutes the Loja Terrane (Litherland et al., 1994;Spikings et al., 2019).
Further west, the Baños Fault and submarine meta-andesites and meta-agglomerates of the Alao-Paute, El Pan y Maguazo Units conform to the Alao Arc. These sequences western margin is the Peltetec Unit, an ophiolitic assemblage constituted by basalt and volcanoclastic rocks, and the Peltetec Fault (Litherland et al., 1994). This fault is the east boundary of a series of slates and quartzites belonging to the Guamote Terrane (Spikings et al., 2019).
The caldera has approximately 17 km of diameter, and the dominant erupted phase is dated as 40Ar/39Ar to 211 ± 14 ky (Hammersley, 2003) and by K-Ar to of 216 ± 5 ky (Bablon et al., 2020) generates massive rhyolitic pyroclastic flow described as ignimbrite. The Chalupas Volcanic materials of the dominant eruptive phase are considered the primary temporary marker of Holocene because fine ashes of this event have been found 900 km away from the emission center, and the estimated deposit volume reaches up to 265 km 3 (Bablon et al., 2020). In the first report provided by OLADE/ INECEL, it recognizes the Chalupas prospect as the highest geothermal project with 283 MWe, based on the volumetric heat-in-place method (INE-CEL, 1983).
Chalupas is the only volcano in Ecuador with a VEI (Volcanic Explosive Index) of 7, placing it as a volcano with ultraplinian eruptions (Bablon et al., 2020). According to Hall and Mothes (2008), the pumice of the rhyolitic eruptions of the Cotopaxi during the Holocene are chemically similar to those of the oldest rhyolites of the Barrancas series, except for the initial eruptive products of the Colorado Canyon series, whose chemistry is similar to ignimbrite (216 ± 5 ky) of the neighboring Chalupas,. resulting in a hypothetical kind of connection or cross-feeding between the magmatic chambers of these two volcanoes. The change in the medium's elastic properties due to a change in lithology is the key factor that induces variations in seismic velocities on rocks. The models v P and ∆v P are commonly used to infer lithological structures and delimited formations. Generally, the rocks with high porosity, such as tuff and pyroclastic deposits, display low v P values, and rocks highly consolidated like granite, or metamorphic rocks, display high values of v P (Lira, 2011).
From the v P cross-section, Figure 4a, it is possible to set the boundaries of the terranes that conforms to the metamorphic basement of the Eastern Cordillera as described in Litherland et al., 1994;Spikings et al., 2019; based on the change in the v P model, the Peltetec, Baños and Llanganates faults (NNE-SSW) can be redrawn. Furthermore, a zone with relatively high v P located 20 km east of Quilindaña at 12 km depth below the surface is identified. The Chalupas does not present significant seismicity below the caldera, where the few events occurring belong to Cotopaxi Volcano.
Based on Figure 4b and 4a, we delimit the terranes that conform to the basement of Chalupas Caldera and northern of Eastern Cordillera, concordant to schemes presented in (Litherland et al, 1994;Spikings et al., 2019). Traversing westward, the Salado oceanic terrane is presented as a zone with relatively high values of v P and ∆v P given by the presence of I-type calc-alkaline Azafrán Batholith (145-140Ma) (Spikings et al., 2019). The Llanganates Fault is the west boundary of the Salado Terrane. To the west side is immediately located the Loja Terrane (continental), particularly Tres Lagunas Unit Trl (240-250 Ma), an S-type granite similar in velocity to Azafrán batholith (Jsaz). There is a hypothesized igneous contact in the location of Llanganates Fault between the Trl Unit and Azafrán batholith of the Salado Terrane (Litherland et al., 1994;Spikings et al., 2019).
The results of tomography are presented on vertical cross sections of the velocity models (v P , ∆v P and v P ⁄ v S ).
Additionally, the tomography shows the hypocenter of seismic events corresponding to 10 km parallel to the cross-sections drawing the seismic box. The azimuths of tomography cross-sections used in this work are presented in Figure 3.

Results of tomography and discussion
desitic periclinal lava flows (Hammersley, 2003). The dominant eruptive phase dated 216 ± 5 ky (Bablon et al., 2020) generates a large pyroclastic flow of ignimbrite and destroys the principal volcano structure.
The Chalupas ignimbrite is formed by an unconsolidated pumitic tuff of light gray ash that presents microcrystals of biotite, plagioclase, and quartz with fragments of pumice (Almeida et al., 1992). The post-collapse eruptions described in Córdova (2018) are less acidic and conform to the successive volcanic events of Quilindaña volcano and Buenavista dome in the last resurgence period.
In the middle of Loja Terrane, Figure 4b, an igneous contact is inferred between the Trl Unit and the metasedimentary rocks of Agoyán Unit (Pzla) (Litherland et al., 1994;Cochrane et al., 2014;Spikings et al., 2019). Pratt et al. (2005) interpret this contact as either a metamorphosed intrusive contact or an unconformity between Paleozoic and Triassic units represented in Figure 2. The west boundary of Loja Terrane is the inferred Baños Fault that crosses close to the Quilindaña summit, covered by Quaternary volcanic deposits, Figure 4a. West to the Baños Fault with a negative value of ∆v P and a relatively low value of v P is located the Alao Terrane (Middle Jurassic), with no reliable age dates. In a consensus, the formation of Alao Arc occurs during the Early Cretaceous within a marine environment (Spikings et al., 2019 (Spikings et al., 2019).
The v P ⁄ v S model is used to determine magmatic and phreatic contexts related to volcanoes on different scales. Low v P ⁄ v S values infer consolidated materials, and high v P ⁄ v S values infer fluid magmatic margins. (Koulakov, 2012).
The v P ⁄ v S model ( Figure 5) shows low values (less than 1.56) directly below the Chalupas Caldera, representing its fill that matches with the caldera edges on the surface (Figure 2), and is limited with the igneous contact within the Loja Terrane (Figure 4a). These values allow assuming that there is no shallow magmatic chamber of considerable magni-tude below the Chalupas Caldera. They neither present seismic and tomographic evidence to assume a possible volcanic eruption before long, nor vertical fluxes of magma feeding any chamber below the Chalupas Caldera. On the other hand, we identified a zone of high values of v P and v P ⁄ v S located 20 to 22 km east of Quilindaña at 12 km depth below the surface that could be associated with a cooled magmatic body, based on parameters describes in Lira (2011). This body is not the main magmatic chamber of Chalupas because it is not placed directly under the caldera and is too far from the surface. Based on Aloisi et al., 2002, this zone could be related to young intrusive bodies and would represent a possible hypothetical heat source of an eventual geothermal resource, which could be heating the water infiltrated by the inferred Llanganates fault, if the fault is permeable.   To compare the Chalupas caldera's seismic and tomographic context with an active volcano, a cross-section centered on Quilindaña with an azimuth of N 34° E addressed to Cotopaxi is presented in Figures 4c and 4d. The cross-sections show a vertical velocity anomaly behind the northeast of the Cotopaxi Volcano. This anomaly is consistent with that inferred by Calahorrano et al. (2019) for the Cotopaxi volcanic complex during the 2015 eruption. This gravimetric study estimates that the deep source of Cotopaxi magma is about 6 km depth and suggests a shallower source of magma 3000 m depth from the summit. This magma interacts with an aquifer and generates infiltration of hydrothermal water on the surface to the northeast of the Cotopaxi volcano. The value of the v P ⁄ v S anomaly beneath Cotopaxi is 1.66, and this is not high enough to represent a magma body rising (Koulakov, 2012). According to Lira (2011), when gas or vapor phases intrude the rock, it produces a decrease of v P ⁄ v S , so that a phreatomagmatic eruption with a high percentage of SO 2 as that produced in 2015 by the Cotopaxi Volcano (Mothes et al., 2017) should present relative smaller values of v P ⁄ v S , agreeing to the model presented in Figure 4d. Figure 6 resumes the 3-dimensional structure of the Chalupas-Quilindaña-Cotopaxi complex. It shows the v P ⁄ v S model beneath the volcanoes, with the Cotopaxi volcano's vertical anomaly in agreement to the eruption of 2015 and the low values below the Chalupas Caldera that infers consolidated rocks. Figure 6 sketches the anomaly below Cotopaxi as an iso-surface v P ⁄ v S =1.66 and the anomaly below the Chalupas-Quilindaña as an iso-volume with v P ⁄ v S ≤1.56 . The topography in Figure 6 is depicted on a real scale.