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. 2021). 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 1990).

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. 2021).

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 1983).

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 2020). 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 2020) 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. 2023).

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 km² 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 2016). 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. 1995). 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 2008, Egüez et al. 2017, Kawsus 2020). 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 2004). 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. 2008). The volcano has had other eruptions that were less explosive in the years 2004 and 2005 (Samaniego et al. 2008), 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. 1996). 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 1990, Tibaldi 2008). 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 2008, Egüez et al. 2017, Kawsus 2020).

In our study region, the Coca River flows through the southern edge of the Reventador volcano cone (IGM 1988). 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. 1996). The Jurassic Misahualli formation composed of welded ash-flow tuffs characterizes the knickpoint zone (INECEL 1992). 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 2008). To the east and northeast of the volcano are two granodiorite igneous intrusions (INECEL 1992, Tibaldi 2008).

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. 2021). 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. 2021). 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.

#	Station Names	#P-phases	#S-phases
1	AZU1	363	182
2	BAMB	33	14
3	CASC	1,249	110
4	CHAR	719	306
5	CONE	2,813	1,342
6	COPE	368	149
7	LAV2	185	60
8	LAV3	1,092	642
9	LAV4	1,494	209
10	PALM	398	255
11	REVE	121	61
12	REVN	1,512	92
13	REVS	1,602	73
14	TRES	24	19

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 2011) and was programmed by Potin (2016). 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) $${\mathit{d}}_{\mathit{o}\mathit{b}\mathit{s}}=\mathit{g}\text{(}\mathit{m}\text{)}$$

In Eq. 1, the data vector d_obs 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) $${‖{\mathit{C}}_{\mathit{d}}^{-\mathbf{1}\mathit{/}\mathbf{2}}\left({\mathit{d}}^{\mathit{o}\mathit{b}\mathit{s}}-\mathit{g}(\mathit{m})\right)‖}_{\mathbb{D}}^{\mathbf{2}}+{‖\mathit{T}\left(\mathit{m}-{\mathit{m}}_{\mathit{p}\mathit{r}\mathit{i}\mathit{o}\mathit{r}}\right)‖}_{\mathbb{M}}^{\mathbf{2}}$$

where T is the regularization operator that, in the case of the Gaussian stochastic approach, is the square root of the covariance operator C_m of the model. C_d is the covariance matrix of the data, and the a priori model is m_prior. To solve this problem numerically, we use the following quasi-Newton algorithm (Tarantola and Valette 1982) (Eq. 3):

(3) $${\mathit{m}}_{\mathit{k}\mathit{+}\mathbf{1}}-{\mathit{m}}_{\mathit{k}}=-{\left({\mathit{C}}_{\mathit{m}}^{-\mathbf{1}}+{\mathit{G}}_{\mathit{k}}^{\mathit{*}}{\mathit{C}}_{\mathit{d}}^{-\mathbf{1}}{\mathit{G}}_{\mathit{k}}\right)}^{-\mathbf{1}}\left({\mathit{G}}_{\mathit{k}}^{\mathit{*}}{\mathit{C}}_{\mathit{d}}^{-\mathbf{1}}\left(\mathit{g}\left({\mathit{m}}_{\mathit{k}}\right)-{\mathit{d}}_{\mathit{o}\mathit{b}\mathit{s}}\right)+{\mathit{C}}_{\mathit{m}}^{-\mathbf{1}}\left({\mathit{m}}_{\mathit{k}}\mathit{-}{\mathit{m}}_{\mathit{p}\mathit{r}\mathit{i}\mathit{o}\mathit{r}}\right)\right)$$

Equation 3 can be written using the matrix decomposition described in Monteiller et al. (2005) (Eq. 4):

(4) $${\mathit{m}}_{\mathit{k}\mathit{+}\mathbf{1}}-{\mathit{m}}_{\mathit{k}}={\left({\mathit{A}}_{\mathit{k}}^{\mathit{*}}{\mathit{A}}_{\mathit{k}}\right)}^{-\mathbf{1}}\left({\mathit{A}}_{\mathit{k}}^{\mathit{*}}{\mathit{v}}_{\mathit{k}}\right)$$

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 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:

$${\mathit{A}}_{\mathit{k}}=\left[\begin{array}{c}{\mathit{C}}_{\mathit{d}}^{-\mathbf{1}\mathit{/}\mathbf{2}}{\mathit{G}}_{\mathit{k}}\\ \mathit{C}\mathit{o}{\mathit{r}}^{-\mathbf{1}\mathit{/}\mathbf{2}}{\mathit{\Sigma }}^{\mathit{-}\mathbf{1}}\end{array}\right],{\mathit{v}}_{\mathit{k}}=\left[\begin{array}{c}{\mathit{C}}_{\mathit{d}}^{-\mathbf{1}\mathit{/}\mathbf{2}}\left(\mathit{g}\left({\mathit{m}}_{\mathit{k}}\right)-{\mathit{d}}_{\mathit{o}\mathit{b}\mathit{s}}\right)\\ \mathit{C}\mathit{o}{\mathit{r}}^{-\mathbf{1}\mathit{/}\mathbf{2}}{\mathit{\Sigma }}^{-\mathbf{1}}\left({\mathit{m}}_{\mathit{k}}-{\mathit{m}}_{\mathit{p}\mathit{r}\mathit{i}\mathit{o}\mathit{r}}\right)\end{array}\right]\text{and,}{\mathit{C}}_{\mathit{m}}=\mathit{\Sigma }\mathit{C}\mathit{o}\mathit{r}\mathit{\Sigma }$$

The operators Cor and ∑ are important because they let us analytically compute the value of C_m^-1/2. This one avoids the numerical work with the matrix of the model covariance C_m, 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) $${\mathit{C}}_{\mathit{m}}\left(\mathit{x},{\mathit{x}}^{\prime }\right)=\mathit{\sigma }(\mathit{x})\mathit{\sigma }\left({\mathit{x}}^{\prime }\right)\mathit{e}\mathit{x}\mathit{p}\left[{\left({\displaystyle \sum _{\mathit{i}\mathit{=}\mathbf{1}}^{\mathbf{3}}\frac{{\left({\mathit{x}}_{\mathit{i}}-{\mathit{x}}_{\mathit{i}}^{\prime }\right)}^{\mathbf{2}}}{{\mathit{\xi }}_{\mathit{i}}^{\mathbf{2}}}}\right)}^{\mathbf{1}/\mathbf{2}}\right]$$

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 $\mathrm{Cor}\left(x,x^{\prime }\right)=\frac{C_m\left(x,x^{\prime }\right)}{\sigma (x)\sigma \left(x^{\prime }\right)}$ Σ 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 d_obs = g(m) and the most reasonable stabilization m - m_prior (Zhdanov 2015). In the INSIGHT stochastic solution, the value of (m - m_prior)(x) is proportional to the information σ(x)² ξ¹ξ²ξ³ inside a cell of volume ξ¹ξ²ξ³. Then, the software introduces a renormalization of the model standard deviation σ to be independent of the correlation lengths (Eq. 6):

(6) $$\mathit{\sigma }=\frac{{\mathit{\xi }}_{\mathbf{0}}^{\mathbf{3}}}{{\mathit{\xi }}_{\mathbf{1}}{\mathit{\xi }}_{\mathbf{2}}{\mathit{\xi }}_{\mathbf{3}}}{\mathit{\sigma }}_{\mathit{p}\mathit{h}\mathit{y}\mathit{s}}$$

where the length ξ₀ 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 ξ₀.

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 V_p and V_p/V_s 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 ξ₀. These tests are represented in Figs. 2A and 2B, with a smoothing range from 2 km to 4 km. In the V_p model (Fig. 2A), the damping ranges from 1 km to 5 km; the V_p/V_s 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) v_p, and (B) v_p/v_s. Color circles show the smoothing lengths, and color lines show the damping lengths. The L-curve corner is ξ₀ = 3 km and ξ₀ = 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 1992) to find the optimization parameters in the corner of the curve, that is ξ₀ = 3 km and ξ₀ = 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) $${\mathit{\rho }}_{\mathit{o}\mathit{b}\mathit{s}}(\mathit{d})=\left(\frac{\mathbf{1}}{\mathbf{2}{\mathit{\sigma }}_{\mathit{o}\mathit{b}\mathit{s}}}\right)\frac{\mathbf{1}}{\mathit{c}\mathit{o}\mathit{s}\mathit{h}\left(\frac{\mathit{\pi }}{\mathbf{2}}\frac{\mathit{d}-{\mathit{d}}_{\mathit{o}\mathit{b}\mathit{s}}}{{\mathit{\sigma }}_{\mathit{o}\mathit{b}\mathit{s}}}\right)}$$

σ_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) $$\mathit{y}(\mathit{d})=\mathit{E}\mathit{r}{\mathit{f}}^{\mathit{-}\mathbf{1}}\left[\frac{\mathbf{2}}{\mathit{\pi }}\mathit{a}\mathit{r}\mathit{c}\mathit{t}\mathit{a}\mathit{n}\left(\mathit{s}\mathit{i}\mathit{n}\mathit{h}\left(\frac{\mathit{\pi }}{\mathbf{2}}\frac{\mathit{d}-{\mathit{d}}_{\mathit{o}\mathit{b}\mathit{s}}}{{\mathit{\sigma }}_{\mathit{o}\mathit{b}\mathit{s}}}\right)\right)\right]$$

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. 2010). From the stochastic solution of the inverse problem (Valette 2011), the resolution operator is calculated as (Eq. 9):

(9) $$\mathit{R}={\mathit{C}}_{\mathit{m}}{\mathit{G}}^{*}{\left({\mathit{C}}_{\mathit{d}}+\mathit{G}{\mathit{C}}_{\mathit{m}}{\mathit{G}}^{\mathit{*}}\right)}^{-\mathbf{1}}\mathit{G}=\mathit{I}-{\left({\mathit{C}}_{\mathit{m}}^{-\mathbf{1}}+{\mathit{G}}^{\mathit{*}}{\mathit{C}}_{\mathit{d}}^{-\mathbf{1}}\mathit{G}\right)}^{-\mathbf{1}}{\mathit{C}}_{\mathit{m}}^{-\mathbf{1}}$$

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

(10) $$\mathit{R}=\mathit{I}-{\left({\mathit{A}}^{\mathit{*}}\mathit{A}\right)}^{-\mathbf{1}}{\mathit{C}}_{\mathit{m}}^{-\mathbf{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 x_i of the grid with the average of the values of R at the x_j points around x_i (Eq. 11):

(11) $${\mathit{J}}_{\mathit{p}}\left({\mathit{x}}_{\mathit{i}}\right)={\displaystyle \sum _{\mathit{j}}\mathit{R}}\left(\mathit{p}\left({\mathit{x}}_{\mathit{i}}\right)\mathit{,}\mathit{p}\left({\mathit{x}}_{\mathit{j}}\right)\right)$$

This average value is the restitution index (RI) (Vergely et al. 2010). 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) $${\mathit{J}}_{\mathit{i}}=\mathbf{1}-{\mathit{m}}_{\mathit{i}}$$

where m_i is a solution of the system (4) with the vector v_k 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 2016). 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. 2021, Pastén-Araya et al. 2022). In the case of the tomography over Ecuador (Araujo et al. 2021), 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: v_p absolute velocity in km/s, the relative velocity: Δv_p in percentage from the a priori model, and the v_p/v_s 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 v_p 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 v_p 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 Δv_p velocity

The P-wave's relative velocity model helps clarify the anomaly observed in the absolute v_p 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 Δv_p 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. 6B (Δv_p), -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 Δv_p helps clarify the anomaly observed in the absolute v_p images. We can observe a volumetric structure defined by Δv_p ≥ +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 v_p model and the relative velocity model Δv_p. 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 v_p in the direction of the knickpoint, we have velocities of 2.5-3 km/s. While in the Δv_p, 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 v_p absolute velocity and (B) the Δv_p 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 v_p/v_s = 1.7.

Figure 8
Cross-section with an azimuth of 80°. Absolute v_p velocity model A and ratio v_p/v_s 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 v_p model and Fig. 6A for the Δv_p 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 1992, Tibaldi 2008) 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 v_p ≥ 5.5 km/s (Flecha et al. 2006). In the Ngatamarki geothermal field in New Zealand, the diorite plutons are identified using seismic travel-time tomography with velocities v_p = 5.5-6 km/s and relative velocity Δv_p = +15%; in this volcanic zone, plutons do not reach the surface and are 5−10 km deep (Sherburn et al. 2003).

Another reason to favor the granodioritic intrusion as an explanation for the high v_p anomaly found in our tomography survey is that volcanic activity in a compressional tectonic realm, such as the Reventador volcano (Tibaldi 2005), 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. 2012). 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). 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 1953). 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. 1999), evidenced by continental tholeiitic volcanism (Romeuf et al. 1995), and is related to the Tetian opening that influenced the evolution of the Northern Andes (Jaillard et al. 1990).

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 1953), 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 1992).

Tschopp (1953) 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. 1990, Aspden and Litherland 1992, Romeuf et al. 1995). 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. 1999). The continental Chapiza Formation is contemporaneous with the Misahualli volcanic arc of continental origin (Baby et al. 2004).

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 v_p.

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 2008).

The mechanical model of Gardner (1983) 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. 2021).

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. 2021). 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. 2021). 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. 2021).

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 (2022).

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. 2004). Some precursory earthquakes were localized to the south-west of the volcano before the paroxysmal eruption in November 2002 (Hall et al. 2004). Like other active volcanoes, the Reventador seismic recurrence produces tremors and explosions at irregular intervals, favoring landslides in the area (Lees et al. 2008).

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