Accessibility / Report Error

Travel-time seismic tomography for the seismic stratigraphic interpretation of the crust around the San Rafael knickpoint at Coca River, Ecuador

Abstract

This research aims to obtain tomographic images to interpret the seismic stratigraphic profile around the San Rafael (SR) knickpoint using seismic tomography. The SR waterfall on the Coca River sinkhole in February 2020 caused regressive soil erosion, suggesting that the knickpoint is highly unstable. Seismic tomography for the P-wave velocity model, utilizing data from 14 stations in the permanent Ecuadorian seismic network, produced vertical cross-sections at azimuths of 0°, 45°, 90°, and 80°. These sections reveal geological formations and their transition to volcanic deposits, while a horizontal cross-section illustrates the extent of the volcanic deposits. The tomographic images facilitated the identification and correlation of seismic stratigraphy with the area's geological features. The most unstable materials for erosion lie from 0 km to 1.5 km with velocities of up to 3.5 km/s. The erosive process impacts the population and infrastructure settled on the banks of the Coca River with a high susceptibility to erosion in zone 1, with velocities of 3 km/s and a 10 km extension. The Coca Codo Sinclair Dam (CCSD) is located in zone 2, with velocities of 3.8 km/s extending up to 40 km with medium susceptibility to erosion.

KEYWORDS:
Coca River; San Rafael knickpoint; volcano deposit; erosion; seismic tomography

INTRODUCTION

The modern San Rafael waterfall was formed due to a slide of volcanic debris and lava flow from the adjacent Reventador volcano. The lava flow caused the Coca River to overtop a high-strength lava barrier, causing slowdowns in its upstream retreat (Reyes et al. 2021Reyes P., Procel S., Sevilla J., Cabero A., Orozco A., Córdova J., Lima F., Vasconez F. 2021. Exceptionally uncommon overburden collapse behind a natural lava dam: Abandonment of the San-Rafael Waterfall in northeastern Ecuador. Journal of South American Earth Sciences, 110, 103353. https://doi.org/10.1016/j.jsames.2021.103353
https://doi.org/10.1016/j.jsames.2021.10...
). The connection between the lava flow and the underlying volcanic debris was preserved before the cascade, reflecting the paleolevel. The composition of the lava flow is between andesite and basaltic andesite (Sevilla 1990Sevilla J. 1990. Un exemple d'importants glissements de terrain en Equateur. International Congress International Association of Engineering Geology, 1713-1717.).

Behind the San Rafael lava dam, volcanic debris sediments formed an extensive network that allowed water to seep in, progressively weakening the top load in the 2020s. Gradual loss of mechanical support led to a collapse and vertical descent, creating an oval sinkhole and causing the immediate abandonment of the cascade. A new interior pathway was created between the modern alluvial surface, the bottom of the sinkhole, and the excavated arch beneath the uncollapsed lava dam (Reyes et al. 2021Reyes P., Procel S., Sevilla J., Cabero A., Orozco A., Córdova J., Lima F., Vasconez F. 2021. Exceptionally uncommon overburden collapse behind a natural lava dam: Abandonment of the San-Rafael Waterfall in northeastern Ecuador. Journal of South American Earth Sciences, 110, 103353. https://doi.org/10.1016/j.jsames.2021.103353
https://doi.org/10.1016/j.jsames.2021.10...
).

The rapid retreat of the Coca River escarpment was triggered by the formation of the sinkhole, following a parallel retreat pattern. This parallel retreat pattern requires a resistant upper layer, a non-resistant underlying material, and active removal of the underlying material to allow basal erosion of the scarp front, which coincides with the erosional process of the Coca River (Gardner 1983Gardner T. W. 1983. Experimental study of knickpoint and longitudinal profile evolution in cohesive, homogeneous material. Geological Society of America Bulletin, 94(5), 664-672. https://doi.org/10.1130/0016-7606(1983)94<664:ESOKAL>2.0.CO;2
https://doi.org/10.1130/0016-7606(1983)9...
).

The causes of the regressive erosion process could explain why, from May 6 to 9, 2020, almost 2.5 km of the slopes and cliffs of the Coca River collapsed. The headward erosion has affected essential infrastructure such as oil pipelines, bridges, and the main road between Ecuador's capital and oil fields in the Amazon Basin. Also, the oil transportation system was the most affected, according to the Servicio Nacional de Gestión de Riegos y Emergencias (SNGRE 2020Servicio Nacional de Gestión de Riegos y Emergencias (SNGRE). 2020. Informe de Situación – Socavamiento Napo. SNGRE.). Geological and geophysical investigations are indispensable for reducing the risk of the collapse of these infrastructures.

To understand the stratigraphy using seismic data, we use travel-time seismic tomography. This method provides the velocity models for the P- and S-waves and the relocation of the seismic events in the study area. The Corporación Eléctrica del Ecuador (CELEC EP 2020Corporación Eléctrica del Ecuador (CELEC EP). 2020. Investigación Del Subsuelo Con Métodos De Geofísica, en el tramo entre el salado y la Cascada San Rafael, debido a la erosión regresiva en el Río Coca. CELEC EP.) conducted low-range seismic and geoelectric studies, detecting that most superficial materials correspond to humid and erodible avalanche sediments. In addition, seismic refraction tomography has been performed in the San Luis Zone, which is located approximately 3 km from the San Rafael knickpoint, where a volcanic breccia was detected that has possibly stopped the regressive erosion in this area (Araujo et al. 2023Araujo S., Guzmán O., Guamán A., Espín R., García I., Chulde E. 2023. Seismic refraction tomography in San Luis, headward Coca River erosion zone. Journal of Applied Geophysics, 212, 104981. https://doi.org/10.1016/j.jappgeo.2023.104981
https://doi.org/10.1016/j.jappgeo.2023.1...
).

However, geophysical methods that can be extended to kilometers at the surface and depth have yet to be performed. The main objective of this research is to obtain tomographic images to interpret the seismic stratigraphic profile around the San Rafael knickpoint. Our research aims to understand the backward erosion and geological risk of the dam from Ecuador's largest hydroelectric power plant, Coca Codo Sinclair. In addition, the area's known geology and tectonic structure allow for correlation and image interpretation.

STUDY AREA

The study region dimensions are 51.5 km in the west-east (WE) direction, 49 km in the north-south (NS) direction, and 40 km in depth in the Coca River basin. The Coca River basin is located in the eastern part of Ecuador, with an area of 5283.74 km2 and an elevation range of 26-5.790 masl. The temperature ranges from 18 to 24 °C, and its rainfall can reach 3,000 mm per year. It is known for its rich water resources, fed by the Salado and Quijos rivers (EPN-PNUMA 2016EPN-PNUMA. 2016. Proyecto TEEB – cuenca del rio coca. Available at: https://www.epn.edu.ec/wp-content/uploads/2018/03/Informe-Final-TEEB-Cuenca-Rio-Coca.pdf. Accessed on: Dec. 14, 2023.
https://www.epn.edu.ec/wp-content/upload...
). The Coca River, flowing northeast, is the largest and most dynamic trunk drainage that receives most of the input from adjacent streams draining from the slopes of the Reventador volcanic edifice (Pourrut et al. 1995Pourrut P., Gomez G., Segovia A. 1995. El agua en el Ecuador: clima, precipitaciones, escorrentía. ). The study area then extends around the upper part of the Coca River basin with the Reventador volcano, as shown in Fig. 1.

Figure 1
(A) Geological map of the study area (Tibaldi 2008Tibaldi A. 2008. Contractional tectonics and magma paths in volcanoes. Journal of Volcanology and Geothermal Research, 176(2), 291-301. https://doi.org/10.1016/j.jvolgeores.2008.04.008
https://doi.org/10.1016/j.jvolgeores.200...
, Egüez et al. 2017Egüez A., Gaona M., Albán A. 2017. Mapa geológico de la República del Ecuador. Ministerio de Minería e Instituto Nacional de Investigación Geológico Minero Metalúrgico., Kawsus 2020Kawsus. 2020. Informe final: contrato de consultoría para levantamiento topográfico con ortofotografia –lidar y diagnostico geológico geotécnico de erosión del lecho del río Coca. Kawsus.). The red star is the San Rafael knickpoint, and the yellow star is the Coca Codo Sinclair Dam (CCSD). (B) Map of the tomography box with topography scale. A thin black line shows the analyzed area that matches the geology map. Black triangles are the 14 seismic stations with names. The red star is the knickpoint, CCSD is in a blue star, and orange dots are the a priori seismicity inside the box.

The Reventador volcano is located in the sub-Andean zone of the Napo uplift. The Napo uplift presents some volcanic edifices: Reventador, Sumaco, Pan de Azúcar, and Yanaurco. The Reventador volcano is the most active, and its eruptions dictate the geomorphological structure of the region (Barragán and Baby 2004Barragán R., Baby P. 2004. Evolución magmática actual de la zona subandina: volcanes el Reventador y Sumaco, modelos geodinámicos preliminares. La Cuenca Oriente: Geología y Petróleo, 183-201.). The Reventador is a stratovolcano whose last significant eruption was in 2002. This eruption had a sub-Plinian activity with VEI 4 (Samaniego et al. 2008Samaniego P., Eissen J.-P., Le Pennec J.-L., Robin C., Hall M. L., Mothes P., Chavrit D., Cotten J. 2008. Pre-eruptive physical conditions of El Reventador volcano (Ecuador) inferred from the petrology of the 2002 and 2004-05 eruptions. Journal of Volcanology and Geothermal Research, 176(1), 82-93. https://doi.org/10.1016/j.jvolgeores.2008.03.004
https://doi.org/10.1016/j.jvolgeores.200...
). The volcano has had other eruptions that were less explosive in the years 2004 and 2005 (Samaniego et al. 2008Samaniego P., Eissen J.-P., Le Pennec J.-L., Robin C., Hall M. L., Mothes P., Chavrit D., Cotten J. 2008. Pre-eruptive physical conditions of El Reventador volcano (Ecuador) inferred from the petrology of the 2002 and 2004-05 eruptions. Journal of Volcanology and Geothermal Research, 176(1), 82-93. https://doi.org/10.1016/j.jvolgeores.2008.03.004
https://doi.org/10.1016/j.jvolgeores.200...
), and the activity continues to this day.

Geology of the study area

The upper part of the Coca River basin is located on an avalanche of debris and sediment accumulated from El Reventador Volcano (Schuster et al. 1996Schuster R. L., NietoThomas A. S., D. O’Rourke T., Crespo E., Plaza-Nieto G. 1996. Mass wasting triggered by the 5 March 1987 ecuador earthquakes. Engineering Geology, 42(1), 1-23. https://doi.org/10.1016/0013-7952(95)00024-0
https://doi.org/10.1016/0013-7952(95)000...
). The outcrops are composed of basaltic lavas, andesites, and dacites. These massive volcanic and volcaniclastic deposits formed in the middle-late Pleistocene to the Holocene depict the geology around El Reventador volcano (Sevilla 1990Sevilla J. 1990. Un exemple d'importants glissements de terrain en Equateur. International Congress International Association of Engineering Geology, 1713-1717., Tibaldi 2008Tibaldi A. 2008. Contractional tectonics and magma paths in volcanoes. Journal of Volcanology and Geothermal Research, 176(2), 291-301. https://doi.org/10.1016/j.jvolgeores.2008.04.008
https://doi.org/10.1016/j.jvolgeores.200...
). Figure 1A shows the position of the San Rafael knickpoint (red star) in the local and regional contexts and the Coca Codo Sinclair Dam (CCSD) (yellow star). Colors represent the principal geological formations: Holocene fluvial-lacustrine sediments (magenta), Holocene volcanic deposits (blue), 19 ky BP debris avalanche (yellow), Late Pleistocene debris avalanche (light blue), Misahualli Formation (green), Granodiorite intrusions (orange), Tena Formation (dark orange), Napo Formation (dark green), Volcanic breccias (red), Coca synthem (black), and Malo synthem (pink). The arrow also indicates the direction of water flow in the Coca River (Tibaldi 2008Tibaldi A. 2008. Contractional tectonics and magma paths in volcanoes. Journal of Volcanology and Geothermal Research, 176(2), 291-301. https://doi.org/10.1016/j.jvolgeores.2008.04.008
https://doi.org/10.1016/j.jvolgeores.200...
, Egüez et al. 2017Egüez A., Gaona M., Albán A. 2017. Mapa geológico de la República del Ecuador. Ministerio de Minería e Instituto Nacional de Investigación Geológico Minero Metalúrgico., Kawsus 2020Kawsus. 2020. Informe final: contrato de consultoría para levantamiento topográfico con ortofotografia –lidar y diagnostico geológico geotécnico de erosión del lecho del río Coca. Kawsus.).

In our study region, the Coca River flows through the southern edge of the Reventador volcano cone (IGM 1988IGM. 1988. Carta Topográfica: Atenas, Escala 1:50.000 (CT-OIII-B1). IGM.). The position of the knickpoint in the San Rafael waterfall is -77.58°W, -0.10°S, at 1,189 m above sea level (Fig. 1). The knickpoint was formed when the lava flow from the Reventador volcano blocked the Coca riverbed (Schuster et al. 1996Schuster R. L., NietoThomas A. S., D. O’Rourke T., Crespo E., Plaza-Nieto G. 1996. Mass wasting triggered by the 5 March 1987 ecuador earthquakes. Engineering Geology, 42(1), 1-23. https://doi.org/10.1016/0013-7952(95)00024-0
https://doi.org/10.1016/0013-7952(95)000...
). The Jurassic Misahualli formation composed of welded ash-flow tuffs characterizes the knickpoint zone (INECEL 1992INECEL. 1992. Anexo G (Vulcanología) Estudio de Factibilidad del Proyecto Hidroeléctrico Coca Codo Sinclair - Fase B. INECEL.). On the other hand, the volcanic composition has three differentiated units: the Coca Synthesis of the Pleistocene, the Malo Synthesis of the middle-late Pleistocene, and the Holocene deposits (Tibaldi 2008Tibaldi A. 2008. Contractional tectonics and magma paths in volcanoes. Journal of Volcanology and Geothermal Research, 176(2), 291-301. https://doi.org/10.1016/j.jvolgeores.2008.04.008
https://doi.org/10.1016/j.jvolgeores.200...
). To the east and northeast of the volcano are two granodiorite igneous intrusions (INECEL 1992INECEL. 1992. Anexo G (Vulcanología) Estudio de Factibilidad del Proyecto Hidroeléctrico Coca Codo Sinclair - Fase B. INECEL., Tibaldi 2008Tibaldi A. 2008. Contractional tectonics and magma paths in volcanoes. Journal of Volcanology and Geothermal Research, 176(2), 291-301. https://doi.org/10.1016/j.jvolgeores.2008.04.008
https://doi.org/10.1016/j.jvolgeores.200...
).

DATA AND STATIONS

The study region has 14 seismic stations on the Ecuadorian Seismic Network (RENSIG). Most RENSIG stations are equipped with three-component seismometers and are broadband. However, 9 stations of the 14 we chose for the study lie in this category. Five stations of the Reventador Seismic Network are very close to the crater and, therefore, are more exposed to eruptive processes; these stations are only three-component short period.

The activity survey of the Reventador volcano shows a high number of seismic stations in the zone. The results obtained in the first tomographic inversion over the Ecuadorian region allowed us to relocate the seismicity from 1988 to 2016 in Ecuador (Araujo et al. 2021Araujo 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
https://doi.org/10.1016/j.jsames.2021.10...
). Subsequently, the data were completed up to September 2019 with the corresponding catalog. Figure 1B shows the geographical location of the 14 seismic stations (black triangles).

To select the data for our tomographic survey, we started with the results of tomography over the entire region of Ecuador (Araujo et al. 2021Araujo 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
https://doi.org/10.1016/j.jsames.2021.10...
). For this model, the cell sizes are 5 km × 5 km × 3 km, aligned with latitude, longitude, and depth, respectively. We attach the additional RENSIG records from April 2016 to April 2019 to this data until a set of 55,181 seismic events is completed. We then chose only earthquakes recorded by 3 or more of the 14 seismic stations we used in our study. With this filter, we obtained 3,923 earthquakes to start the tomography. It is unnecessary to impose spatial filtering on the earthquakes because the INSIGHT software automatically eliminates those earthquakes outside the tomography box. Of these events within the tomography box, 1,631 remain, as represented in Figure 1B.

This selection gives a total of 15,487 data points to invert, where 11,973 are P-phases, and 3,514 are S-phases. To test the data coverage at all stations, the number of events for each station is depicted in Table 1. In this box, the cells are cubes with a 0.5 km edge. The number of cells in the box is 103 in X (east-west), 98 in Y (north-south), and 80 in Z (depth).

Table 1
Number of P-wave phases and S-wave phases for each station.

The INSIGHT and TIME3D software work only in a cartesian reference system, so the symmetry of the modeling is only rectangular. There is no other symmetry imposed on the velocity models.

TOMOGRAPHY METHOD

The INSIGHT software is a set of routines written in Fortran 90 that use parallel computing libraries. This software implements the stochastic solution of the inverse problem (Valette 2011Valette B. 2011. Inversion of geophysical data. In EFIDIR: Extraction and Fusion of Information for Displacement measurement from SAR Imagery. Ecole de Physique Des Houches, Chamonix. Available at: http://efidir.poleterresolide.fr/index.php/efidir-seminars/spring-school-2011. Accessed on: Feb. 6, 2023.
http://efidir.poleterresolide.fr/index.p...
) and was programmed by Potin (2016)Potin B. 2016. Les Alpes occidentales: tomographie, localisation de séismes et topographie du Moho. Doctoral dissertation, Grenoble Alpes.. It has several associated tools that allow us to choose the seismicity and the size of the tomography box and build the a priori model. Due to the size of the problems posed in tomography, INSIGHT works on a computer cluster, in our case, at the CIMENT of Grenoble University. The connection to the cluster is made through the SSH protocol that allows loading the seismicity files and the a priori model, launching the tomography jobs, and downloading the result of the inverse problem.

INSIGHT solves the inverse problem that comes from the direct problem equation (Eq. 1):

(1) d o b s = g ( m )

In Eq. 1, the data vector dobs contains the travel times of P- and S-waves, the model vector m has the velocity of the P-wave, the ratio between the velocity of the P-wave and the S-wave, the hypocenters, and the origin time of the earthquakes. Then, INSIGHT solves the problem of seismic event location within the same inversion problem. At the output of each iteration of the software, we obtain a new model of velocities and the new source locations in this model. An additional parameter in the model vector is the site effect of each seismic station. The linear operator g is a functional theoretical relation between the data and the parameters.

The stochastic inversion leads to a Tikhonov problem. This problem involves searching for the minimum cost function in a regularization space. The cost function has the form (Eq. 2):

(2) C d 1 / 2 d o b s g ( m ) D 2 + T m m p r i o r M 2

where T is the regularization operator that, in the case of the Gaussian stochastic approach, is the square root of the covariance operator Cm of the model. Cd is the covariance matrix of the data, and the a priori model is mprior. To solve this problem numerically, we use the following quasi-Newton algorithm (Tarantola and Valette 1982Tarantola 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
https://doi.org/10.1029/RG020i002p00219...
) (Eq. 3):

(3) m k + 1 m k = C m 1 + G k * C d 1 G k 1 G k * C d 1 g m k d o b s + C m 1 m k - m p r i o r

Equation 3 can be written using the matrix decomposition described in Monteiller et al. (2005)Monteiller V., Got J., Virieux J., Okubo P. 2005. An efficient algorithm for double-difference tomography and location in heterogeneous media, with an application to the Kilauea volcano. Journal of Geophysical Research: Solid Earth, 110(B12), 1-22. https://doi.org/10.1029/2004JB003466
https://doi.org/10.1029/2004JB003466...
(Eq. 4):

(4) m k + 1 m k = A k * A k 1 A k * v k

In each iteration of Eq. 4, we must solve the direct problem: find the propagation times of the seismic waves in the new model obtained in the iteration. In addition, we must trace the seismic rays in the a posteriori model. To solve this problem, INSIGHT uses the TIME3D algorithm (Podvin and Lecomtte 1991Podvin P., Lecomte I. 1991. Finite difference computation of traveltimes in very contrasted velocity models: a massively parallel approach and its associated tools. Geophysical Journal International, 105(1), 271-284. https://doi.org/10.1111/j.1365-246X.1991.tb03461.x
https://doi.org/10.1111/j.1365-246X.1991...
) that traces the diving seismic rays by solving the Eikonal equation by finite differences.

Our data have only direct P- and S-phases of the diving seismic rays, and we do not have data on refractions and reflections, so tracing these rays is unnecessary. The starting data in the tomography survey are travel times; hence, we do not consider the amplitude of seismic rays or waveforms.

The accuracy of TIME3D depends on the size of the chosen cells. The dimensions of the inversion mesh are too large for this task. For this reason, each iteration of the a posteriori model is interpolated from cells of 0.5 km on a side to cells of 0.1 km. The algorithm TIME3D traces the rays over these new cells and computes the new propagation times. The minimum size of the forward modeling mesh is established by the dimensions of the specific problem being solved. With a mesh of 0.1 km per side, we have 51.5 * 49 * 40/0.1 ≈ 1 milion, which is the number of knots that can be processed in a reasonable time by the computer cluster.

The decomposition (4) implies the following definitions:

A k = C d 1 / 2 G k C o r 1 / 2 Σ - 1 , v k = C d 1 / 2 g m k d o b s C o r 1 / 2 Σ 1 m k m p r i o r and, C m = Σ C o r Σ

The operators Cor and are important because they let us analytically compute the value of Cm-1/2. This one avoids the numerical work with the matrix of the model covariance Cm, a sparse matrix with many zeros and, therefore, very difficult to invert with traditional methods. The covariance kernel for two points of the inversion grid x and x´´ with a correlation length ξ and standard deviations σ(x) and σ(x´´) is (Eq. 5):

(5) C m x , x = σ ( x ) σ x e x p i = 1 3 x i x i 2 ξ i 2 1 / 2

where the three indexes for I = 1, 2, and 3 represent the three directions of the space; thus, the covariance kernel takes into account the influence between the points of the grid during the inversion process. The correlation kernel is defined as and Corx,x=Cmx,xσ(x)σx Σ the multiplication operator by σ is (ƒ) = σ(x)ƒ(x).

After being solved, the inverse problem must be regularized. That is, we must choose between the infinite solutions and the optimal solution of Eq. 2. The optimal solution expresses the trade-off between the best data fitting dobs = g(m) and the most reasonable stabilization m - mprior (Zhdanov 2015Zhdanov M. 2015. Inverse theory and applications in geophysics (Vol. 36). Elsevier.). In the INSIGHT stochastic solution, the value of (m - mprior)(x) is proportional to the information σ(x)2 ξ1ξ2ξ3 inside a cell of volume ξ1ξ2ξ3. Then, the software introduces a renormalization of the model standard deviation σ to be independent of the correlation lengths (Eq. 6):

(6) σ = ξ 0 3 ξ 1 ξ 2 ξ 3 σ p h y s

where the length ξ0 is the characteristic size of the volume of matter for which it is possible to define the standard deviation σphys. Then we have all the parameters to regularize the tomography problem: the model standard deviation σphys, the correlation of smoothing lengths ξs, and the attenuation or damping length ξ0.

We choose the values of σp = 0.750 km/s and of σvp/vs = 0.15 based on 10% of the value of the respective maximum values of the a priori velocity models in both Vp and Vp/Vs ratio. It also set the correlation lengths equal in the three space directions ξs = ξx = ξy = ξz.

A total of 16 tests were performed, varying the smoothing parameters length ξs and the damping length ξ0. These tests are represented in Figs. 2A and 2B, with a smoothing range from 2 km to 4 km. In the Vp model (Fig. 2A), the damping ranges from 1 km to 5 km; the Vp/Vs model (Fig. 2B) has a damping range from 0.5 km to 3 km. Each test had an execution time of 10 h using three computer cluster nodes, with 32 processors and 188 GB of memory per node.

Figure 2
L-curves for the tomography tests for (A) vp, and (B) vp/vs. Color circles show the smoothing lengths, and color lines show the damping lengths. The L-curve corner is ξ0 = 3 km and ξ0 = 4 km is marked with a red arrow. The histograms of the residuals for the (C) P-wave and (D) S-wave are in thin lines for the a priori and colored bars for the a posteriori.

Figures 2A and 2B are L-curves. Then, we use the L-curve criteria (Hansen 1992Hansen P. 1992. Analysis of discrete ill-posed problems by means of the L-curve. SIAM Review, 34(4), 561-580. https://doi.org/10.1137/1034115
https://doi.org/10.1137/1034115...
) to find the optimization parameters in the corner of the curve, that is ξ0 = 3 km and ξ0 = 4 km.

The correct convergence of the INSIGHT algorithm is observed by computing the residuals for both the P-wave and the S-wave, as shown in Figs. 2C and 2D. We have calculated the residuals of the two phases before and after tomography and represented them as histograms with a thin line for the a priori residuals and colored bars for the a posteriori residuals. We see how, in both cases, the histograms narrow after the tomography, indicating the improvement in the model fitting to the data.

Regarding the data, one problem that solves the INSIGHT software is the long-tail probability function distribution of the outliers. The solution is to take a hyperbolic secant distribution for the data new variable d (Eq. 7):

(7) ρ o b s ( d ) = 1 2 σ o b s 1 c o s h π 2 d d o b s σ o b s

σobs is the error in the data attributed to their quality picking.

The stochastic approach works only with Gaussian probability distribution functions; thus, it is necessary to define a new variable, which is the inversion (Eq. 8):

(8) y ( d ) = E r f - 1 2 π a r c t a n s i n h π 2 d d o b s σ o b s

To test the resolution level of our tomography experiment, we start with the concept of the resolution operator. This operator is a way to test the resolving power of a tomography study (Vergely et al. 2010Vergely J.-L., Valette B., Lallement R., Raimond S. 2010. Spatial distribution of interstellar dust in the Sun's vicinity. Astronomy and Astrophysics, 518, A31. https://doi.org/10.1051/0004-6361/200913962
https://doi.org/10.1051/0004-6361/200913...
). From the stochastic solution of the inverse problem (Valette 2011Valette B. 2011. Inversion of geophysical data. In EFIDIR: Extraction and Fusion of Information for Displacement measurement from SAR Imagery. Ecole de Physique Des Houches, Chamonix. Available at: http://efidir.poleterresolide.fr/index.php/efidir-seminars/spring-school-2011. Accessed on: Feb. 6, 2023.
http://efidir.poleterresolide.fr/index.p...
), the resolution operator is calculated as (Eq. 9):

(9) R = C m G * C d + G C m G * 1 G = I C m 1 + G * C d 1 G 1 C m 1

With the decomposition of Eq. 4, the resolution operator is (Eq. 10):

(10) R = I A * A 1 C m 1

In tomography problems, it is not practical to calculate the value of the operator at all points of the inversion grid. The proposal then is to replace the value of R for a parameter p at a point xi of the grid with the average of the values of R at the xj points around xi (Eq. 11):

(11) J p x i = j R p x i , p x j

This average value is the restitution index (RI) (Vergely et al. 2010Vergely J.-L., Valette B., Lallement R., Raimond S. 2010. Spatial distribution of interstellar dust in the Sun's vicinity. Astronomy and Astrophysics, 518, A31. https://doi.org/10.1051/0004-6361/200913962
https://doi.org/10.1051/0004-6361/200913...
). The RI takes values close to one when the averaged values are close to the actual value, indicating good resolution. On the contrary, if the average is far from one, the average does not represent the actual value, and the resolution is low.

To calculate the RI with the INSIGHT software, we replace Eq. 11 in Eq. 10 to obtain Eq. 12:

(12) J i = 1 m i

where mi is a solution of the system (4) with the vector vk having only the model component while the data component is zero. This value of the RI is obtained for each point of the three-dimensional (3D) grid, and we can observe it graphically through sections in Fig. 3. The directions for the cross-sections are the same as for the tomography cross-sections represented in Fig. 4.

Figure 3
Cross-sections of the restitution index (RI) for the azimuths of (A) 0°, (B) 45°, and (C) 90°. The values of RI are marked with iso-value curves and a color scale. The transparent coloring in the value 3.5 marks the limits of the resolution.
Figure 4
Topography map with the seismic hypocenters after tomography and topography scale. The thin black line is the analyzed region. The map shows the lines plotting the azimuths of the three tomography cross sections at the San Rafael contact point (red star) with a posteriori seismicity and Coca Codo Sinclair Dam (CCSD) with a blue star. In addition, we performed a cross-section of the Reventador volcano.

Since the INSIGHT software has already been used to obtain some tomography results in different contexts, we can know the value we can expect in our experiment for the RI. In the case of a well-developed seismic network over the Alps, the resolution obtained is higher than the value of 0.7 (Potin 2016Potin B. 2016. Les Alpes occidentales: tomographie, localisation de séismes et topographie du Moho. Doctoral dissertation, Grenoble Alpes.). For experiments in a context of active subduction, such as the Chile Trench, where adequate azimuthal coverage throughout the territory is not possible, the limit for resolution drops to 0.6 (Pastén-Araya et al. 2021Pastén-Araya F., Potin B., Ruiz S., Zerbst L., Aden-Antoniów F., Azúa K., Rivera E., Rietbrock A., Salazar P., Fuenzalida A. 2021. Seismicity in the upper plate of the Northern Chilean offshore forearc: Evidence of splay fault south of the Mejillones Peninsula. Tectonophysics, 800, 228706. https://doi.org/10.1016/j.tecto.2020.228706
https://doi.org/10.1016/j.tecto.2020.228...
, Pastén-Araya et al. 2022Pastén-Araya F., Potin B., Azúa K., Sáez M., Aden-Antoniów F., Ruiz S., Cabrera L., Ampuero J. P., Nocquet J. M., Rivera L., Duputel Z. 2022. Along-dip segmentation of the slip behavior and rheology of the Copiapó ridge subducted in north-central Chile. Geophysical Research Letters, 49(4), e2021GL095471. https://doi.org/10.1029/2021GL095471
https://doi.org/10.1029/2021GL095471...
). In the case of the tomography over Ecuador (Araujo et al. 2021Araujo 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
https://doi.org/10.1016/j.jsames.2021.10...
), the resolution value is reduced to 0.5 due to using a historical database with very heterogeneous seismic networks.

For our experiment, we have limited the interpretation of the results to values of the RI greater than 3.5. This cutoff value is shown as a shaded area in the cross-sections in Fig. 3 and quantifies the scope and limitations of our resulting model. As seen in Fig. 1B, we do not have correct azimuthal coverage of the seismic stations around the San Rafael knickpoint, which explains a lower value in our RI for our study compared to previous research cited in the preceding paragraph.

RESULTS

The tomography cross-sections for the San Rafael knickpoint have three azimuths (Fig. 4): 0°, 45°, and 90°. The azimuth of 45° follows the Coca River course (Fig. 4) and shows the structure outside the Reventador edifice. The seismicity represented in the cross-sections is the result of the relocalizations obtained after the tomography. There are three different models as INSIGHT results: vp absolute velocity in km/s, the relative velocity: Δvp in percentage from the a priori model, and the vp/vs ratio. Each of these three models gives information about the geological structure of the region. In addition, we can add a cut at 80° in the direction of the Reventador volcano.

Cross-section at azimuths of 0°, 45°, and 90° of the model absolute vp velocity

Two velocity changes can be identified in the stratigraphy of the cross-sections in Fig. 5. The first velocity change in stratigraphy is between 1 km and 4 km depth; values oscillate between 3.5 and 4.5 km/s. The thickness of the unit varies from end to end. The first end (D) is thicker than the other (C) in Fig. 5A. Here, we can observe in the central part (from -10 to 5 km) a lateral discontinuity of convex shape upward (anticline) from the second change in velocities. The second change in stratigraphy is located between 4 and 9 km depth; here, we have values between 4.6 and 5.4 km/s.

Figure 5
Cross-section with azimuths of 0°, 45°, and 90° for the local tomography result of the model absolute vp velocity. The altitude of the topography is between 1,200 and 2,480 masl.

We can appreciate a lateral discontinuity convex upward that deforms the material of the stratigraphy in the central part, and a discontinuity can be observed to the right with less visibility. We can observe the Napo uplift with its different geological formations. In the most superficial part, we have volcanic products of the Reventador, such as avalanche deposits between 2 and 3.8 km/s in Fig. 5B. In addition, these deposits rest on the bedrock of the Misahualli Formation with velocities between 3.8 and 4.5 km/s and a depth between 2and 5 km deep (Fig. 5B).

The cross-section of Fig. 5C covers the southern part (left) of the Reventador volcano. Here, we can observe 2.5 km/s up to 3 km depth velocities. These velocities are from volcanic deposits derived from the volcano.

Cross-section at azimuths of 0°, 45°, and 90° of the model relative Δvp velocity

The P-wave's relative velocity model helps clarify the anomaly observed in the absolute vp images in Fig. 5.

We can appreciate a lateral discontinuity convex upward that deforms the material of the stratigraphy in the central part, and a discontinuity can be observed to the right with less visibility.

However, in Fig. 6A, we can see that the lateral discontinuity of the convex (anticline) shape in the central and right parts may be other geological structures differing from the surrounding material. The first (central) structure has relative velocity values of 1−8%, and the second structure has values ranging from 1 to 4%.

Figure 6
Cross-section with azimuths of 0°, 45°, and 90° for the local tomography result of the model relative Δvp velocity. The altitude of the topography is between 1,200 and 2,480 masl.

In addition, the surrounding materials are between -30 and 0% (Fig. 6A). This is a possible interpretation considering the region's strong thrust regime. Nevertheless, the anomaly in P-wave velocity goes until 4 km below the surface, and it is visible in the three cross-sections, which discards the fold interpretation because the anticlines present a preferential folding direction.

In Fig. 6Bvp), -20% is observed in the first 2 km of depth in the central part (knickpoint), while the convex anomaly reaches values of 10%.

The P-wave's relative velocity model Δvp helps clarify the anomaly observed in the absolute vp images. We can observe a volumetric structure defined by Δvp ≥ +10% and highlighted by the red color in Fig. 6B. This structure is approximately 8 km long in the NS direction, 5 km in the WE direction, and 5 km below the surface.

We can observe regions with high uncertainty shaded in white in the three cross-sections. In these places, obtaining a specific result but making a guess according to the results with low uncertainty is impossible. The lack of seismic information causes these uncertainties due to the position of the seismic stations distributed almost all around the volcano.

The horizontal cross-section in the study area

Figure 7 shows the horizontal slice of the volcanic deposits produced by its proximity to the Reventador volcano. We present the vp model and the relative velocity model Δvp. Here we can observe the shape and extent of the volcanic deposits and the possible materials that have eroded and possibly will erode. In the vp in the direction of the knickpoint, we have velocities of 2.5-3 km/s. While in the Δvp, we have an oval-shaped extension ranging from -40 to -12% with a length of approximately 40 km.

Figure 7
Horizontal cross-sections 1 km below (A) the surface of the vp absolute velocity and (B) the Δvp relative velocity. The SR knickpoint and Coca Codo Sinclair Dam (CCSD) (black star). Shows the Coca River, Salado River, and Quijos River with line blue and green circles as the seismic stations. The altitude of the topography is between 1,200 and 2,480 masl.

Cross-section at azimuth 80°

In addition, it was possible to perform the tomographic cross-section at 80° to the east of the Reventador volcano since the relocation of the seismic data obtained by INSIGHT is essential at this location, as shown in Fig. 4. After the relocalization, the earthquakes define a structure almost vertical, 12 km below the surface (Fig. 8). The resulting model B also shows a deeper magmatic chamber 15 km to the east and a 25 km depth defined by a high value of the velocity's ratio vp/vs = 1.7.

Figure 8
Cross-section with an azimuth of 80°. Absolute vp velocity model A and ratio vp/vs velocity model B (km = scale). The altitude of the topography is between 1,300 and 3,360 masl.

DISCUSSION

The seismic stratigraphic interpretation is carried out from the oldest to the most recent geological events (Fig. 9). An igneous intrusion is centered at San Rafael knickpoint below Unit 2 at a depth of 10 km in hard rock. There is evidence of another two intrusions of the same type in the region (Fig. 1A) to 10 km NE. We can mention that the tomographic cross-section of one of these intrusions can be seen in Fig. 5A for the vp model and Fig. 6A for the Δvp model; here, we can see that the power of this intrusion is 8 km.

Figure 9
Seismic stratigraphic interpretation of the 45° cross-section. The seismic stratigraphic units are identified in the central part of the line (-5 km left), where there is convex deformation (anomaly).

The previous geological studies (INECEL 1992INECEL. 1992. Anexo G (Vulcanología) Estudio de Factibilidad del Proyecto Hidroeléctrico Coca Codo Sinclair - Fase B. INECEL., Tibaldi 2008Tibaldi A. 2008. Contractional tectonics and magma paths in volcanoes. Journal of Volcanology and Geothermal Research, 176(2), 291-301. https://doi.org/10.1016/j.jvolgeores.2008.04.008
https://doi.org/10.1016/j.jvolgeores.200...
) do not show evidence of this intrusion, probably because it was recovered with the debris avalanche of Reventador. However, seismic tomography studies show that in many geological contexts (e.g., in the Iberian Massif, Spain), the granitic plutons are detected by 2D reflection seismic tomography with velocities vp ≥ 5.5 km/s (Flecha et al. 2006Flecha I., Carbonell R., Zeyen H., Martí D., Palomeras I., Simancas F., Pérez-Estaún A. 2006. Imaging granitic plutons along the IBERSEIS profile. Tectonophysics, 420(1-2), 37-47. https://doi.org/10.1016/j.tecto.2006.01.019
https://doi.org/10.1016/j.tecto.2006.01....
). In the Ngatamarki geothermal field in New Zealand, the diorite plutons are identified using seismic travel-time tomography with velocities vp = 5.5-6 km/s and relative velocity Δvp = +15%; in this volcanic zone, plutons do not reach the surface and are 5−10 km deep (Sherburn et al. 2003Sherburn S., Bannister S., Bibby H. 2003. Seismic velocity structure of the central Taupo Volcanic Zone, New Zealand, from local earthquake tomography. Journal of Volcanology and Geothermal Research, 122(1-2), 69-88. https://doi.org/10.1016/S0377-0273(02)00470-5
https://doi.org/10.1016/S0377-0273(02)00...
).

Another reason to favor the granodioritic intrusion as an explanation for the high vp anomaly found in our tomography survey is that volcanic activity in a compressional tectonic realm, such as the Reventador volcano (Tibaldi 2005Tibaldi A. 2005. Volcanism in compressional tectonic settings: Is it possible? Geophysical Research Letters, 32(6), L06309. https://doi.org/10.1029/2004GL021798
https://doi.org/10.1029/2004GL021798...
), produces granite magma migration and its emplacement in the form of plutons, which can be demonstrated by the Reventador volcano using an analogical model (Ferré et al. 2012Ferré E. C., Galland O., Montanari D., Kalakay T. J. 2012. Granite magma migration and emplacement along thrusts. International Journal of Earth Sciences, 101(7), 1673-1688. https://doi.org/10.1007/s00531-012-0747-6
https://doi.org/10.1007/s00531-012-0747-...
). This igneous intrusion has a power of 8 km and produces convex deformation in the zone below the knickpoint, allowing the overlying materials to deform and the more recent materials to become unstable as they consolidate.

Due to the few studies carried out in the sub-Andean zone at greater depths, there is no information on the stratigraphy since it does not outcrop in the Amazon basin. However, Unit 2 could be considered the Santiago Formation based on the studies carried out by Tschopp (1953)Tschopp H. 1953. Oil Explorations in the Oriente of Ecuador, 1938-1950. AAPG Bulletin, 37(10), 2303-2347. https://doi.org/10.1306/5CEADD94-16BB-11D7-8645000102C1865D
https://doi.org/10.1306/5CEADD94-16BB-11...
. This formation has a thickness ranging from 1.5 to 2.7 km. This formation consists of thin-bedded black siliceous limestones alternating with calcareous sandstones and sometimes bituminous black shales (Tschopp 1953Tschopp H. 1953. Oil Explorations in the Oriente of Ecuador, 1938-1950. AAPG Bulletin, 37(10), 2303-2347. https://doi.org/10.1306/5CEADD94-16BB-11D7-8645000102C1865D
https://doi.org/10.1306/5CEADD94-16BB-11...
). This formation is compact and not susceptible to erosion, has a convex deformation caused by igneous intrusion beneath it, and is in stratigraphic contact below the Misahualli Formation. This cycle is syn-tectonic and records the opening of a "rift" (Christophoul et al. 1999Christophoul F., Baby P., Dávila C. 1999. Descripción de las influencias eustáticas y tectónicas en la Cuenca Oriente ecuatoriana desde el Aptiano hasta el Oligoceno.), evidenced by continental tholeiitic volcanism (Romeuf et al. 1995Romeuf 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.), and is related to the Tetian opening that influenced the evolution of the Northern Andes (Jaillard et al. 1990Jaillard 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
https://doi.org/10.1144/gsjgs.147.6.1009...
).

In Unit 1, we have the Misahualli Formation in direct contact with the volcanic materials laying above this formation (erodible materials). Being in contact with the possible Santiago Formation, this also has a convex deformation, causing instability of the erodible volcanic materials. This formation outcrops at -15 km (Fig. 9, SSW); this outcrop is the last non-erodible layer. The Misahualli Formation thickness usually ranges from 1 to 4 km (Tschopp 1953Tschopp H. 1953. Oil Explorations in the Oriente of Ecuador, 1938-1950. AAPG Bulletin, 37(10), 2303-2347. https://doi.org/10.1306/5CEADD94-16BB-11D7-8645000102C1865D
https://doi.org/10.1306/5CEADD94-16BB-11...
), and we obtain 4 km in our tomography result. This interpretation matches the geological background of the region where the Jurassic Misahualli Formation is the basement of the Coca River course (INECEL 1992INECEL. 1992. Anexo G (Vulcanología) Estudio de Factibilidad del Proyecto Hidroeléctrico Coca Codo Sinclair - Fase B. INECEL.).

Tschopp (1953)Tschopp H. 1953. Oil Explorations in the Oriente of Ecuador, 1938-1950. AAPG Bulletin, 37(10), 2303-2347. https://doi.org/10.1306/5CEADD94-16BB-11D7-8645000102C1865D
https://doi.org/10.1306/5CEADD94-16BB-11...
defined the Misahualli Formation as an upper member of the Chapiza Formation. The Misahualli Formation is a geological formation of tuffs and basalts, very hard rocks. Therefore, Misahualli and the Chapiza Formations are taken as unit(s) 1 in Fig. 10. The Misahualli Formation consists of volcanic accumulations and intrusions (Aspden et al. 1990Aspden J. A., Rundle C., Viteri F., Bermudez R., Harrison S. 1990. Edades radiométricas del Batolito de Zamora-rio Mayo. Boletín Geológico Ecuatoriano, 1(1), 85-88., Aspden and Litherland 1992Aspden 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
https://doi.org/10.1016/0040-1951(92)904...
, Romeuf et al. 1995Romeuf 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.). The geodynamic event terminates the opening of the "rift" ("aborted rift"), which causes uplift and emersion of the "rift," resulting in the erosional surface of the base of the Chapiza Formation (Christophoul et al. 1999Christophoul F., Baby P., Dávila C. 1999. Descripción de las influencias eustáticas y tectónicas en la Cuenca Oriente ecuatoriana desde el Aptiano hasta el Oligoceno.). The continental Chapiza Formation is contemporaneous with the Misahualli volcanic arc of continental origin (Baby et al. 2004Baby P., Rivadeniera M., Barragán R. 2004. La cuenca Oriente: geología y petróleo (Vol. 144). In P. Baby, M. Rivadeneira, R. Barragán (Eds.), Institut français d’études andines. ).

Figure 10
Horizontal cross-sections 1 km below the surface show the extent of volcanic deposits. The knickpoint position (black star): zone 1 (black lines): high susceptibility to erosion. The Coca Codo Sinclair Dam (CCSD) (black star): zone 2 (magenta lines) medium susceptibility to erosion. White triangles indicate human settlements, and the Coca River, Salado River, and Quijos River (blue lines; black arrow; water flow direction) and green circles for the station seismic. The colors show the velocity vp.

In the study area, we have the deposits of the Reventador volcano that are currently eroding. These deposits are on top of Unit 1, composed of gravitational debris avalanches and water flows that settled in this place. Erodible materials are dacites, andesites, and rhyolites interrupted by basaltic lavas at the knickpoint (Tibaldi 2008Tibaldi A. 2008. Contractional tectonics and magma paths in volcanoes. Journal of Volcanology and Geothermal Research, 176(2), 291-301. https://doi.org/10.1016/j.jvolgeores.2008.04.008
https://doi.org/10.1016/j.jvolgeores.200...
).

The mechanical model of Gardner (1983)Gardner T. W. 1983. Experimental study of knickpoint and longitudinal profile evolution in cohesive, homogeneous material. Geological Society of America Bulletin, 94(5), 664-672. https://doi.org/10.1130/0016-7606(1983)94<664:ESOKAL>2.0.CO;2
https://doi.org/10.1130/0016-7606(1983)9...
states that a parallel retreat model must meet three conditions: (1) a resistant upper layer; (2) an underlying level of non-resistant material; and (3) active removal of the underlying material to allow basal scour of the headwall, which meets the conditions of the Coca River erosional process.

Long periods of slow retreat alternate with short periods of rapid retreat, explaining that the sequence of volcaniclastic sediments is a heterogeneous material with variable resistance to erosion, mainly in the sections of the river upstream of the cascade point (Reyes et al. 2021Reyes P., Procel S., Sevilla J., Cabero A., Orozco A., Córdova J., Lima F., Vasconez F. 2021. Exceptionally uncommon overburden collapse behind a natural lava dam: Abandonment of the San-Rafael Waterfall in northeastern Ecuador. Journal of South American Earth Sciences, 110, 103353. https://doi.org/10.1016/j.jsames.2021.103353
https://doi.org/10.1016/j.jsames.2021.10...
).

The collapse of the San Rafael cascade starts with the transport of avalanche deposits below the lava flow. This process relied on the weakly resistant materials of the volcanic deposits, leaving newly abandoned alluvial terrace plains in the downstream reaches. Whereas, upstream of the knickpoint, the erosional process showed a non-constant rate of retreat, suggesting that the volcanic deposit, despite its loose nature, exhibited a strong internal heterogeneity (Reyes et al. 2021Reyes P., Procel S., Sevilla J., Cabero A., Orozco A., Córdova J., Lima F., Vasconez F. 2021. Exceptionally uncommon overburden collapse behind a natural lava dam: Abandonment of the San-Rafael Waterfall in northeastern Ecuador. Journal of South American Earth Sciences, 110, 103353. https://doi.org/10.1016/j.jsames.2021.103353
https://doi.org/10.1016/j.jsames.2021.10...
). Figure 10 of Zone 1 (black lines) shows velocities < 3.5 km/s over a distance of 10 km. This zone is very close to the knickpoint (black star), indicating a high susceptibility to erosion, as low velocities (< 3.5 km/s) indicate looser materials.

Furthermore, it is expected that backward erosion will reach the point of the CCSD. However, predicting when the erosion will reach the dam is uncertain because the erodible material is heterogeneous and the backward velocity varies depending on factors such as geology, stratigraphy, tectonics, and seismicity that can accelerate or slow down the erosion process (Reyes et al. 2021Reyes P., Procel S., Sevilla J., Cabero A., Orozco A., Córdova J., Lima F., Vasconez F. 2021. Exceptionally uncommon overburden collapse behind a natural lava dam: Abandonment of the San-Rafael Waterfall in northeastern Ecuador. Journal of South American Earth Sciences, 110, 103353. https://doi.org/10.1016/j.jsames.2021.103353
https://doi.org/10.1016/j.jsames.2021.10...
). Figure 10 shows that the CCSD (black star) is located where the velocity changes from 3.5 to 3.8 km/s, which indicates a medium susceptibility to erosion. Besides, the CCSD, human settlements, roads, pipelines, and other infrastructure will be damaged due to backward erosion (Reyes et al. 2021Reyes P., Procel S., Sevilla J., Cabero A., Orozco A., Córdova J., Lima F., Vasconez F. 2021. Exceptionally uncommon overburden collapse behind a natural lava dam: Abandonment of the San-Rafael Waterfall in northeastern Ecuador. Journal of South American Earth Sciences, 110, 103353. https://doi.org/10.1016/j.jsames.2021.103353
https://doi.org/10.1016/j.jsames.2021.10...
).

The regressive erosion is enhanced in the main impact zone, located approximately 100 m from the community of San Luis, where the risk of slope instability is high. In addition, the erosion front is located at a distance of 7.9 km from the CCSD and remains in the zone of medium susceptibility to erosion, according to the Dirección de Monitoreo de Eventos Adversos (2022Dirección de Monitoreo de Eventos Adversos. 2022. Informe de Situación No. 68 Erosión Hídrica Regresiva en el tramo fluvial: Hidroeléctrica Coca Codo Sinclair – Río Napo. Available at: https://www.gestionderiesgos.gob.ec/wp-content/uploads/downloads/2022/07/SITREP-No-68-Erosion-Hidrica-NapoOrellana-15072022.pdf. Accessed on: May 29, 2023.
https://www.gestionderiesgos.gob.ec/wp-c...
).

In addition, a cross-section that passes through the Reventador volcano was analyzed (Fig. 8). In this area, it is possible to obtain seismic stratigraphic data. This activity may indicate an active feeding system for the volcano through faults with a depth of 10 km. However, the seismic activity in the Reventador volcano is characterized by a high number of tectonic and volcano-tectonic events that make it difficult to separate these two phenomena (Hall et al. 2004Hall M., Ramón P., Mothes P., LePennec J. L., García A., Samaniego P., Yepes H. 2004. Volcanic eruptions with little warning: the case of Volcán Reventador's Surprise November 3, 2002 Eruption, Ecuador. Revista Geológica de Chile, 31(2), 349-358. https://doi.org/10.4067/S0716-02082004000200010
https://doi.org/10.4067/S0716-0208200400...
). Some precursory earthquakes were localized to the south-west of the volcano before the paroxysmal eruption in November 2002 (Hall et al. 2004Hall M., Ramón P., Mothes P., LePennec J. L., García A., Samaniego P., Yepes H. 2004. Volcanic eruptions with little warning: the case of Volcán Reventador's Surprise November 3, 2002 Eruption, Ecuador. Revista Geológica de Chile, 31(2), 349-358. https://doi.org/10.4067/S0716-02082004000200010
https://doi.org/10.4067/S0716-0208200400...
). Like other active volcanoes, the Reventador seismic recurrence produces tremors and explosions at irregular intervals, favoring landslides in the area (Lees et al. 2008Lees J. M., Johnson J. B., Ruiz M., Troncoso L., Welsh M. 2008. Reventador Volcano 2005: Eruptive activity inferred from seismo-acoustic observation. Journal of Volcanology and Geothermal Research, 176(1), 179-190. https://doi.org/10.1016/j.jvolgeores.2007.10.006
https://doi.org/10.1016/j.jvolgeores.200...
).

CONCLUSION

The tomographic images provided insight into the geology and structures beneath the San Rafael knickpoint, allowing strata identification and correlation with the geology of the area. The tomographic analysis allowed the identification of unstable materials highly susceptible to erosion at a depth of 0-1.5 km with velocities of up to 3.5 km/s.

In addition, the igneous intrusion was detected below the SR knickpoint, which allowed the deformation of the Misahualli and possible Santiago geological formations, causing instability in the area.

The erosion process impacts the population and essential infrastructure settled along the banks of the Coca River with a high susceptibility to erosion in zone 1, with velocities of 3 km/s and an extension of 10 km. It can be expected that the regressive erosion will eventually reach the CCSD, which lies in a region with unconsolidated volcanic deposits with velocities of 3.8 km/s that are susceptible to erosion.

Another factor that allows instability for the regressive erosion process in the area is the high seismicity of the Reventador volcano.

ACKNOWLEDGMENTS

The running of the tomography software INSIGHT for this study is performed in the computer cluster DAHU of the center Calcul Intensif/Modélisation/Expérimentation Numérique et Technologique CIMENT at Grenoble University. The research team thanks the three reviewers for their comments that enhanced the quality of this article.

REFERENCES

  • Araujo S., Guzmán O., Guamán A., Espín R., García I., Chulde E. 2023. Seismic refraction tomography in San Luis, headward Coca River erosion zone. Journal of Applied Geophysics, 212, 104981. https://doi.org/10.1016/j.jappgeo.2023.104981
    » https://doi.org/10.1016/j.jappgeo.2023.104981
  • 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
    » 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
    » https://doi.org/10.1016/0040-1951(92)90426-7
  • Aspden J. A., Rundle C., Viteri F., Bermudez R., Harrison S. 1990. Edades radiométricas del Batolito de Zamora-rio Mayo. Boletín Geológico Ecuatoriano, 1(1), 85-88.
  • Baby P., Rivadeniera M., Barragán R. 2004. La cuenca Oriente: geología y petróleo (Vol. 144). In P. Baby, M. Rivadeneira, R. Barragán (Eds.), Institut français d’études andines.
  • Barragán R., Baby P. 2004. Evolución magmática actual de la zona subandina: volcanes el Reventador y Sumaco, modelos geodinámicos preliminares. La Cuenca Oriente: Geología y Petróleo, 183-201.
  • Christophoul F., Baby P., Dávila C. 1999. Descripción de las influencias eustáticas y tectónicas en la Cuenca Oriente ecuatoriana desde el Aptiano hasta el Oligoceno
  • Corporación Eléctrica del Ecuador (CELEC EP). 2020. Investigación Del Subsuelo Con Métodos De Geofísica, en el tramo entre el salado y la Cascada San Rafael, debido a la erosión regresiva en el Río Coca CELEC EP.
  • Dirección de Monitoreo de Eventos Adversos. 2022. Informe de Situación No. 68 Erosión Hídrica Regresiva en el tramo fluvial: Hidroeléctrica Coca Codo Sinclair – Río Napo Available at: https://www.gestionderiesgos.gob.ec/wp-content/uploads/downloads/2022/07/SITREP-No-68-Erosion-Hidrica-NapoOrellana-15072022.pdf Accessed on: May 29, 2023.
    » https://www.gestionderiesgos.gob.ec/wp-content/uploads/downloads/2022/07/SITREP-No-68-Erosion-Hidrica-NapoOrellana-15072022.pdf
  • Egüez A., Gaona M., Albán A. 2017. Mapa geológico de la República del Ecuador Ministerio de Minería e Instituto Nacional de Investigación Geológico Minero Metalúrgico.
  • EPN-PNUMA. 2016. Proyecto TEEB – cuenca del rio coca Available at: https://www.epn.edu.ec/wp-content/uploads/2018/03/Informe-Final-TEEB-Cuenca-Rio-Coca.pdf Accessed on: Dec. 14, 2023.
    » https://www.epn.edu.ec/wp-content/uploads/2018/03/Informe-Final-TEEB-Cuenca-Rio-Coca.pdf
  • Ferré E. C., Galland O., Montanari D., Kalakay T. J. 2012. Granite magma migration and emplacement along thrusts. International Journal of Earth Sciences, 101(7), 1673-1688. https://doi.org/10.1007/s00531-012-0747-6
    » https://doi.org/10.1007/s00531-012-0747-6
  • Flecha I., Carbonell R., Zeyen H., Martí D., Palomeras I., Simancas F., Pérez-Estaún A. 2006. Imaging granitic plutons along the IBERSEIS profile. Tectonophysics, 420(1-2), 37-47. https://doi.org/10.1016/j.tecto.2006.01.019
    » https://doi.org/10.1016/j.tecto.2006.01.019
  • Gardner T. W. 1983. Experimental study of knickpoint and longitudinal profile evolution in cohesive, homogeneous material. Geological Society of America Bulletin, 94(5), 664-672. https://doi.org/10.1130/0016-7606(1983)94<664:ESOKAL>2.0.CO;2
    » https://doi.org/10.1130/0016-7606(1983)94<664:ESOKAL>2.0.CO;2
  • Hall M., Ramón P., Mothes P., LePennec J. L., García A., Samaniego P., Yepes H. 2004. Volcanic eruptions with little warning: the case of Volcán Reventador's Surprise November 3, 2002 Eruption, Ecuador. Revista Geológica de Chile, 31(2), 349-358. https://doi.org/10.4067/S0716-02082004000200010
    » https://doi.org/10.4067/S0716-02082004000200010
  • Hansen P. 1992. Analysis of discrete ill-posed problems by means of the L-curve. SIAM Review, 34(4), 561-580. https://doi.org/10.1137/1034115
    » https://doi.org/10.1137/1034115
  • IGM. 1988. Carta Topográfica: Atenas, Escala 1:50.000 (CT-OIII-B1) IGM.
  • INECEL. 1992. Anexo G (Vulcanología) Estudio de Factibilidad del Proyecto Hidroeléctrico Coca Codo Sinclair - Fase B. INECEL.
  • 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
    » https://doi.org/10.1144/gsjgs.147.6.1009
  • Kawsus. 2020. Informe final: contrato de consultoría para levantamiento topográfico con ortofotografia –lidar y diagnostico geológico geotécnico de erosión del lecho del río Coca Kawsus.
  • Lees J. M., Johnson J. B., Ruiz M., Troncoso L., Welsh M. 2008. Reventador Volcano 2005: Eruptive activity inferred from seismo-acoustic observation. Journal of Volcanology and Geothermal Research, 176(1), 179-190. https://doi.org/10.1016/j.jvolgeores.2007.10.006
    » https://doi.org/10.1016/j.jvolgeores.2007.10.006
  • Monteiller V., Got J., Virieux J., Okubo P. 2005. An efficient algorithm for double-difference tomography and location in heterogeneous media, with an application to the Kilauea volcano. Journal of Geophysical Research: Solid Earth, 110(B12), 1-22. https://doi.org/10.1029/2004JB003466
    » https://doi.org/10.1029/2004JB003466
  • Pastén-Araya F., Potin B., Azúa K., Sáez M., Aden-Antoniów F., Ruiz S., Cabrera L., Ampuero J. P., Nocquet J. M., Rivera L., Duputel Z. 2022. Along-dip segmentation of the slip behavior and rheology of the Copiapó ridge subducted in north-central Chile. Geophysical Research Letters, 49(4), e2021GL095471. https://doi.org/10.1029/2021GL095471
    » https://doi.org/10.1029/2021GL095471
  • Pastén-Araya F., Potin B., Ruiz S., Zerbst L., Aden-Antoniów F., Azúa K., Rivera E., Rietbrock A., Salazar P., Fuenzalida A. 2021. Seismicity in the upper plate of the Northern Chilean offshore forearc: Evidence of splay fault south of the Mejillones Peninsula. Tectonophysics, 800, 228706. https://doi.org/10.1016/j.tecto.2020.228706
    » https://doi.org/10.1016/j.tecto.2020.228706
  • Podvin P., Lecomte I. 1991. Finite difference computation of traveltimes in very contrasted velocity models: a massively parallel approach and its associated tools. Geophysical Journal International, 105(1), 271-284. https://doi.org/10.1111/j.1365-246X.1991.tb03461.x
    » https://doi.org/10.1111/j.1365-246X.1991.tb03461.x
  • Potin B. 2016. Les Alpes occidentales: tomographie, localisation de séismes et topographie du Moho Doctoral dissertation, Grenoble Alpes.
  • Pourrut P., Gomez G., Segovia A. 1995. El agua en el Ecuador: clima, precipitaciones, escorrentía.
  • Reyes P., Procel S., Sevilla J., Cabero A., Orozco A., Córdova J., Lima F., Vasconez F. 2021. Exceptionally uncommon overburden collapse behind a natural lava dam: Abandonment of the San-Rafael Waterfall in northeastern Ecuador. Journal of South American Earth Sciences, 110, 103353. https://doi.org/10.1016/j.jsames.2021.103353
    » https://doi.org/10.1016/j.jsames.2021.103353
  • 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.
  • Samaniego P., Eissen J.-P., Le Pennec J.-L., Robin C., Hall M. L., Mothes P., Chavrit D., Cotten J. 2008. Pre-eruptive physical conditions of El Reventador volcano (Ecuador) inferred from the petrology of the 2002 and 2004-05 eruptions. Journal of Volcanology and Geothermal Research, 176(1), 82-93. https://doi.org/10.1016/j.jvolgeores.2008.03.004
    » https://doi.org/10.1016/j.jvolgeores.2008.03.004
  • Schuster R. L., NietoThomas A. S., D. O’Rourke T., Crespo E., Plaza-Nieto G. 1996. Mass wasting triggered by the 5 March 1987 ecuador earthquakes. Engineering Geology, 42(1), 1-23. https://doi.org/10.1016/0013-7952(95)00024-0
    » https://doi.org/10.1016/0013-7952(95)00024-0
  • Servicio Nacional de Gestión de Riegos y Emergencias (SNGRE). 2020. Informe de Situación – Socavamiento Napo SNGRE.
  • Sevilla J. 1990. Un exemple d'importants glissements de terrain en Equateur. International Congress International Association of Engineering Geology, 1713-1717.
  • Sherburn S., Bannister S., Bibby H. 2003. Seismic velocity structure of the central Taupo Volcanic Zone, New Zealand, from local earthquake tomography. Journal of Volcanology and Geothermal Research, 122(1-2), 69-88. https://doi.org/10.1016/S0377-0273(02)00470-5
    » https://doi.org/10.1016/S0377-0273(02)00470-5
  • 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
    » https://doi.org/10.1029/RG020i002p00219
  • Tibaldi A. 2005. Volcanism in compressional tectonic settings: Is it possible? Geophysical Research Letters, 32(6), L06309. https://doi.org/10.1029/2004GL021798
    » https://doi.org/10.1029/2004GL021798
  • Tibaldi A. 2008. Contractional tectonics and magma paths in volcanoes. Journal of Volcanology and Geothermal Research, 176(2), 291-301. https://doi.org/10.1016/j.jvolgeores.2008.04.008
    » https://doi.org/10.1016/j.jvolgeores.2008.04.008
  • Tschopp H. 1953. Oil Explorations in the Oriente of Ecuador, 1938-1950. AAPG Bulletin, 37(10), 2303-2347. https://doi.org/10.1306/5CEADD94-16BB-11D7-8645000102C1865D
    » https://doi.org/10.1306/5CEADD94-16BB-11D7-8645000102C1865D
  • Valette B. 2011. Inversion of geophysical data. In EFIDIR: Extraction and Fusion of Information for Displacement measurement from SAR Imagery Ecole de Physique Des Houches, Chamonix. Available at: http://efidir.poleterresolide.fr/index.php/efidir-seminars/spring-school-2011 Accessed on: Feb. 6, 2023.
    » http://efidir.poleterresolide.fr/index.php/efidir-seminars/spring-school-2011
  • Vergely J.-L., Valette B., Lallement R., Raimond S. 2010. Spatial distribution of interstellar dust in the Sun's vicinity. Astronomy and Astrophysics, 518, A31. https://doi.org/10.1051/0004-6361/200913962
    » https://doi.org/10.1051/0004-6361/200913962
  • Zhdanov M. 2015. Inverse theory and applications in geophysics (Vol. 36). Elsevier.

Publication Dates

  • Publication in this collection
    15 Apr 2024
  • Date of issue
    2024

History

  • Received
    05 Sept 2023
  • Accepted
    23 Jan 2024
Sociedade Brasileira de Geologia R. do Lago, 562 - Cidade Universitária, 05466-040 São Paulo SP Brasil, Tel.: (55 11) 3459-5940 - São Paulo - SP - Brazil
E-mail: sbgeol@uol.com.br