# Imaging magma storage below Teide volcano (Tenerife)
using scattered seismic wavefields

L. de Barros, F. Martini, C.J. Bean, A. Garcia-Yeguas, J. Ibanez

► To cite this version:

L. de Barros, F. Martini, C.J. Bean, A. Garcia-Yeguas, J. Ibanez. Imaging magma storage below Teide volcano (Tenerife) using scattered seismic wavefields. Geophysical Journal International, 2012, 191 (2), pp.695-706. 10.1111/j.1365-246X.2012.05637.x. hal-00812632

HAL Id: hal-00812632

https://hal.science/hal-00812632

Submitted on 8 Nov 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Ridge (Singh et al. 2006) and along the East Pacific Rise (Kent et al. 2000) but beneath volcanoes, contrary to some petrological expectations, seismic studies find little sign of large magmatic chambers. In practice detecting reflections from chambers is fraught with problems, due to the highly heterogeneous nature of the edifice. The strong layering of the materials with high impedance contrasts (from competent basalt to ash) and irregular layer interfaces, the highly fractured nature of rocks filled by different fluids and the complex topography all severely distort the wavefield. Isolating weak reflected arrivals is consequently a difficult task, as the wavefield may be dominated by scattered waves rather than coherent reflected waves. Beam-forming approaches attempt to address this problem and provide an efficient tool to image the energy that is scattered back by the structures (e.g. Lay 1987; Hedlin et al. 1991; Kruger et al. 1996; Scherbaum et al. 1997; Rietbrock & Scherbaum 1999; Maercklin et al. 2004). The method looks for coherent energy in the seismograms by either stacking or correlating traces. It is usually used in the frequency–wavenumber domain using array data to stack the energy coming from the same azimuth (e.g. Lacoss et al. 1969; Ringdal & Husebye 1982). Here we take advantage of an existing well-resolved tomographic model to look for coherent energy coming from a grid of potential scatter points at the delays computed in the tomographic model. As scatter points can be seen as secondary sources of the seismic energy, this method is somewhat similar to the semblance technique, normally used for seismic source location (Neidell & Taner 1971; Almendros & Chouet 2003). These techniques usually use an array of stations, but by using the Green’s function reciprocity, they can also take advantage of a cluster of sources (Spudich & Bostwick 1987). In our study, we apply the beam-forming approach on both source and receiver arrays, simultaneously, leading to the Double Beam Method (DBM, Kruger et al. 1996). A similar method was previously used by Maercklin (2008) in the Campi Flegrei caldera, in the Bay of Pozzuoli (Italy): the study inferred seismic scattering from a shallow buried caldera rim.

Here, we apply this approach to a data set from Tenerife, the largest and the most inhabited island of the Canary volcanic archipelago (Spain), situated in the North Atlantic Ocean, about 200 km west of Morocco. This volcanic island is dominated by Teide, a 3718 m stratovolcano. In spring 2004 an unusual increase in seismic activity was observed on the island. More than 500 earthquakes were recorded before the activity decreased in the first half of 2005 (Almendros et al. 2007). On the basis of this activity a reawakening of the volcano was suggested and strongly debated at the time (Carracedo et al. 2006). For this reason, a seismic survey was undertaken on Tenerife island in 2007 January through the TOM-TEIDEVS project (Ibáñez et al. 2008), to derive a 3-D seismic velocity tomographic image and gain a better understanding of the volcanic system. The acquired data set comprises more than 6000 off-shore shots recorded by 137 land-based seismic stations. A high-resolution, reliable 3-D velocity tomography of the volcanic island of Tenerife has been obtained by García-Yeguas et al. (2012), which represents the starting model for our analysis. Hence, there is a good general geological framework in place, including a complex velocity model and an excellent data set, to try to image magmatic storage using the scattered wavefield.

In this paper we first describe the geological context (Section 2), the data (Section 3) and the beam-forming method used in this study (Section 4). Synthetic tests are then shown (Section 5) to explore the robustness and the power of this approach to resolve structures. The application of the method to the Tenerife data set is then described in Section 6: two main scattering structures are identified beneath the

island. Their interpretation as magmatic storage areas is discussed: it is found to be consistent with the two magmatic reservoirs inferred by petrology and other geophysical studies.

## 2 TENERIFE: A VOLCANIC ISLAND

Tenerife island, one of the seven Canary islands, has a volcanic history of about 12 Ma. The volcanic processes, still partly unknown, led to a complex geological setting (from basaltic to felsic materials) and a sharp topography. The central part of the island is occupied by a large caldera (Las Cañadas caldera). The origin of this large depression (16 km in diameter, with up to 500 m high walls) is still under debate, as it could be produced by flank failure or vertical collapse (Blanco-Montenegro et al. 2011). Pico Teide and Pico Viejo, two large strato-volcanoes have developed inside the caldera. The Teide peak reaches 3718 m a.s.l., and more than 7000 m above the seafloor.

During the recent volcanic activity, two kinds of volcanism occurred (Martí et al. 2008): (1) a basaltic fissure volcanism, located mainly on the NE and NW ridges and in the south of the island, which produced monogenic cones; and (2) phonolitic eruptions in the Las Cañadas complex, which includes the Pico Viejo-Teide edifice. The basaltic eruptions, fed by the deeper magmatic system, are more frequent (recurrent period of 100–200 yr), the last one occurred in 1909. In 2004, the seismic crisis was interpreted as a possible sign of reawakening of the volcanic complex (Martí et al. 2009).

Martí & Geyer (2009) showed that the eruption vent locations, that is, from the summit or on the flanks, are controlled by the shape of the magmatic chamber and the stress distribution around them. Understanding the geometry of this system is of crucial importance to be able to estimate the size of a possible eruption. Several studies have investigated the volcanic structure below the caldera using geological and geophysical data (Martí et al. 1994; Ablay & Martí 2000; Araña et al. 2000) and they infer a shallow phonolitic magma chamber and a deeper reservoir containing basaltic magma. Petrological evidence, including melt inclusion studies, suggest that the Pico Teide phonolites were stored, prior to eruption, at a shallow depth of about 2 km b.s.l. (Ablay & Martí 2000). The occurrence of summit caldera collapse also supports this evidence of a shallow magma chamber reservoir. Magnetic and gravity studies (Araña et al. 2000; Blanco-Montenegro et al. 2011) infer an anomaly related to magma at a depth between 6 and 10 km. The estimated source depth of Volcano-Tectonic (VT) seismicity recorded in 2004 (Cerdeña et al. 2011) is also compatible with this depth for a basaltic magma supply zone beneath the Las Cañadas caldera.

## 3 THE DATA SET

An active seismic experiment, within the TOM-TEIDEVS project (Ibáñez et al. 2008), was performed in 2007 January to obtain a 3-D seismic velocity tomographic image of Tenerife island (García-Yeguas et al. 2012). Six BOLT 1500LL air-guns, with a maximum capacity of 3520 cubic inches per shot, were used for a total of 6459 shots all around the island (Fig. 1a). Shots were separated by 2 min intervals, giving a spatial interval of around 300 m. The shots were recorded by 137 land-based seismic sensors recording continuously with a sampling rate of 200 Hz. The stations used in this study are broad-band Guralp CMG-6TD, with natural frequency response from 0.33 to 100 Hz. We only use the vertical component, as *P* waves are recorded primarily in this direction. These stations were
<figure><img src="image_2.png" /><figcaption>**Figure 1.** (a) Map view of Tenerife island, including topography, bathymetry and station and shot locations during the TOM-TEIDEVS project. The interval between elevation lines is 1000 m. Red triangles represent the broad-band station locations, and white dots indicate offshore shots. Red triangles and crosses represent stations and shot locations, respectively, in (b) south-north cross-section, at the UTM easting of 314.1 km stations, and (c) west-east profile, at the UTM northing of 3127 km. The analysis described in this paper is performed using these two 2-D profiles.</figcaption></figure>

distributed across the island (Fig. 1a) with the highest density in the region of Teide. Subsets of the stations were aligned in a south-north and a west-east profile crossing the summit of the volcano. Shot lines were also aligned on each side on those profiles, leading to 2-D profiles of shots and receivers. In this study, we focus on those two profiles. The south-north profile comprises 26 stations, with two shot lines in the south (52 shots each) and one shot line (120 shots) in the north (Fig. 1b). The west-east profile has 23 stations with shot lines of 53 shots on both sides (see Fig. 1c). We use the same number of shots for each profile (i.e. 52 and 53, respectively), which leads to a total number of 2704 and 2438 traces for the S-N and W-E profiles, respectively.

The frequency content of the data ranges from 5 to 12 Hz. However, the energy is mainly concentrated between 5 and 7 Hz, as the highest frequencies are highly attenuated. As the data were recorded by instruments with different recording spectra, the data have been

corrected for the relevant instrument response. Two receiver gathers are shown in Fig. 2: the first onsets are very pronounced. Identifiable *P*-wave arrivals are recognized up to 45 km offset, in some cases even up to 60 km. Later in the signals, the wavefield is very complex and includes *S* wave, surface wave and scattered waves. Some secondary waves (reflections, *P*-to-*S* conversion at the ocean bottom, etc.) can be seen, but their energy is very weak and they are partly buried in noise. A specific method has thus to be designed to extract information from these arrivals.

## 4 DOUBLE BEAM IMAGING METHOD

The signal-to-noise ratio can be enhanced, for a particular seismic phase, by stacking the coherent arrivals after correction for the delay times. By this method, it is possible to exploit weak arrivals, such as
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Example of two receiver gather data: (a) a single station in the south of the island, with a shot line in the south and (b) a single station in the west of the island, with a shot line in the east.</figcaption></figure>

reflected waves or scattered waves. Stacking can be performed using either shot (active or passive source), or receiver arrays, or a combination of the two, which leads to the DBM (Kruger et al. 1993). An extensively used approach for locating seismicity, especially in volcanic environments, is searching for the location where the most coherent energy is generated. As scatter points can be considered as secondary energy sources, an image of the scattering structures can be extracted by applying this procedure to active experiment data. As we are dealing with active sources and a sparse station network and because scatter points are close to the receivers, the standard plane wave assumption (which is useful to apply the method in the frequency–wavenumber domain) cannot be used here. Therefore, the method is applied here in the time–space domain following Rietbrock & Scherbaum (1999), and it is based on the traveltimes, which are used to align coherent arrivals. This method is not an inversion process, but rather a signal processing technique that can be seen as a hyperbola summation migration (Yilmaz 2001). As each seismic path is assumed to cross only one scatter point, we are dealing with a single scattering assumption.

First, we build a spatial grid of potential scatter points. We then compute the theoretical propagation time from the source to each scatter point and from the scatter point to each receiver, using the tomographic model of García-Yeguas et al. (2012). This computation is performed using the 2-D ray tracer developed by Zelt & Smith (1992). The 2-D geometry is a crude assumption, which can lead to a simplification of the travel paths. However, as the tomographic model is rather smooth, and to keep the approach consistent with the numerical modelling (see Section 5), only 2-D propagation is assumed in this study. The time delays are then used to shift all traces, to align the waveforms associated to the source-scattered point-receiver path at the origin time. A time mute function (Hanning) is then applied around the origin time, to keep only the waves of interest. We use 0.4 s length windows to include two periods of signal. The next step is to evaluate the coherency between all traces recorded for all source–receiver pairs. A semblance method, based on the normalized stack of all traces (Neidell & Taner 1971), is the most standard approach. In double beam imaging, the general formula for the semblance coherency (double beam–semblance,

DBS) is given by (Rietbrock & Scherbaum 1999):

$$ \text{DBS}(k) = \frac{\sum_{t=-T/2}^{T/2} \left[ \sum_{i=1}^{N} S_{ik} |x_i(t - t_{ik})|^{1/n} \right]^{2n}}{N \sum_{t=-T/2}^{T/2} \sum_{i=1}^{N} x_i^2(t - t_{ik})}, \\ \text{with } S_{ik} = \text{sign}[x_i(t - t_{ik})], \qquad (1) $$

where $x_i$ is the trace for the ith receiver–station pair, $t_{ik}$ is the traveltime in the ith pair for the wave passing by the scatter point k, N represents the number of source–receiver pairs and T is the length of the Hanning windows. If n is equal to 1, eq. (1) exactly reduces to the semblance (Neidell & Taner 1971). For n greater than 1, this expression is the nth root semblance (Kanasewich et al. 1973). The latter measure of the coherency is more dependent on the phase of the signals (Schimmel & Paulssen 1997). Alternatively, another estimation of the coherency between traces can be obtained by a normalized cross-correlation function (double beam–correlation, DBC):

$$ \text{DBC}(k) = \frac{2}{N(N-1)} \sum_{i}^{N} \sum_{j}^{i-1} \frac{f(\gamma_{ij,k}(t))}{\sqrt{\gamma_{ii,k}(0) \gamma_{jj,k}(0)}}, \\ \text{with } \gamma_{ij,k}(t) = \text{corr}[x_i(t - t_{ik}); x_j(t - t_{jk})]. \qquad (2) $$

The function f can be: (1) the value at t = 0 ($f(\gamma_{ij,k}(t)) = \gamma_{ij,k}(0)$), with the assumption that the traveltime used to shift the data is perfectly exact (hereafter referred to as correlation), and (2) the maxima of the correlation function around t = 0 using $f(\gamma_{ij,k}(t)) = \max(\gamma_{ij,k}(t))$ if the traveltimes used to shift the signals are considered imprecise (hereafter referred as corrected correlation). The latter is computationally more expensive, as the correlation functions have to be fully computed for all receiver–station pairs and every scattering point. While other measurements of the coherency can be used (e.g. the phase-weighted stack, Schimmel & Paulssen 1997), they are not investigated herein. The values obtained using eqs (1) and (2) are referred hereatfer as beampower.

Eqs (1) or (2) are computed for each scatter point of the grid, using the traveltimes computed in the ray tracer of Zelt & Smith
(1992). Assuming that the noise is not coherent, a high coherency value means that coherent energy is being scattered back, leading to an image of the medium in terms of scattering structures. The method is sensitive to the geometry of the experiment, as shown by Maercklin *et al.* (2004). Moreover, the presence of undesirable waves (such as direct waves, *S* and surface waves) can contaminate the beampower image. To remove these, the beam method is applied to the data as well as to the synthetic full wavefield data computed in the reference model (i.e. the relatively smooth tomography model used for traveltime computation). The resulting image is constructed by taking the difference between the data and synthetic images. This is an extension of the beam correction used by Hedlin *et al.* (1991) and Maercklin (2008). As the process is linear and the coherency measurement is normalized, this allows an easy suppression of the unwanted effects associated to the reference model.

# 5 BEAM FORMING ON SYNTHETIC DATA: TESTING THE METHOD

## 5.1 Numerical test data

Synthetic data are computed to (1) check the accuracy of the method; (2) look for the best measure of the wave coherency and (3) estimate which type of structure can be reconstructed by such a method. Data are computed using a 2-D rotated staggered grid finite difference

scheme, with a fourth order in space and second order in time. We initially focus on the south-north profile: the model size is 105.6 km long and 20.9 km deep, which also includes 2.8 km wide absorbing boundaries on both sides and at the base of the medium. The west-east profile, used in the data section, is 91.2 × 21.8 km. To accurately represent the topography, a minimum of 20 gridpoints per minimum wavelength are required for code stability. To comply with this, the medium is discretized with a 10 m grid spacing. The 2-D velocity model contains up to 26 million gridpoints. Because of the size of the models, the length of the signals (40 s, with a time sample of $dt = 0.0005$ s) and the number of stations, the computation of synthetic data in a full 3-D model is not computationally feasible. The test model has a vertical *P*-wave velocity ($Vp$) gradient which follows the topography and bathymetry, increasing from 3.5 km s⁻¹ at the surface to 6 km s⁻¹ at 10 km below the free surface. The *S*-wave velocity is defined as $Vp/\sqrt{3}$. Model perturbations are then applied to the starting model: (1) model A, a reflecting structure at 10 km depth (elliptical body, Fig. 3) and (2) model B, a scattering structure (Julia fractal). To produce scattered and reflected waves without modifying the traveltimes in the medium, the *P*-wave velocity is increased by 20 per cent only in the 200-m-thick edge of both of those structures. In a volcanic medium, both the *P*- and the *S*-wave velocities are modified by the presence of magma and volatiles. However, as we only have a *P*-wave tomography model for Tenerife island, an analysis of the scattered *S* wave and converted

<figure><img src="image_4.png" /><figcaption>**Figure 3.** (a) Test velocity model (model A) used to compute synthetic data. The acquisition geometry is given in Fig. 1, that is, with offshore shots on both sides of the island and receivers on land. A 200-m-thick-edged elliptical body (black line) with a perturbation of +20 per cent in the *P*- and *S*-wave velocities is inserted into the velocity model. This ellipse is 2 km wide and 30 km long, to mimic a large anomaly related to magma. The white rectangle marks the area where the beam-forming analysis is performed. (b) Receiver gather of the synthetic data computed in the reference model (without the elliptical perturbation) for the same configuration as in Fig. 2. (c) As (b) using the model with the perturbation (model A), and (d) differential data (i.e. data of panel c minus data of panel b). Amplitudes of the differential data are an order of magnitude smaller than the amplitude of the data in panels b and c. Note that the amplitudes in this figure are clipped to highlight later arrivals.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Results of beam-forming approach to reconstruct the elliptical perturbation in model A, Fig. 3. Coherency among traces is measured by (a) semblance, (b) second root semblance, (c) correlation and (d) corrected correlation. The colour scale gives the beampower, with dark colours indicating locations with highly scattered energy. The dashed line indicates the location of the elliptical perturbation inserted into the velocity model.</figcaption></figure>

*P*-to-*S* wave is not feasible. The synthetic models are thus designed for the analysis of the *P* waves only.

Examples of synthetic data generated in model A are shown in Fig. 3. With a simple velocity gradient model, the wavefield appears very complex. This is due mainly to the topography scattering, which produces a very long coda (as shown by O’Brien & Bean 2009; Kumagai *et al.* 2011), as well as to the water/elastic transition. As expected, the data computed from the reference model and from model A are very similar: the response to the perturbation (i.e. the difference between seismograms) is very weak (an order of magnitude smaller than the total response of the medium) and cannot be seen directly in the data. In the differential data (Fig. 3d), we can see that the major difference is in the *P*-reflected wave, and later, with smaller amplitude, we can also identify changes in the converted *P*-to-*S* wave, in the (converted at the sea bottom) reflected *S*-wave and in the surface waves.

## 5.2 Synthetic tests results

The methodology described in Section 4 is first applied to model A (data shown in Fig. 3), to investigate the most efficient measure of coherency among the four proposed in Section 4. The full data set is used, i.e. no time mute is applied to the data before the analysis. The grid of scattering points is 40 km long, from the sea level to 14 km depth, with a spacing of 200 m in both the horizontal and vertical direction. It is restricted to the part of the profile beneath the island.

The reconstruction of the elliptical structure is shown in Fig. 4, where coherency is measured with semblance, second root semblance, correlation and time-corrected correlation. The geometry of the experiment produces a sensitivity to the scattered points which varies along the profile. In all cases, the structure is well reconstructed on the southern side of the profile, less in the northern side and poorly in the middle of the profile. The results obtained using the second root semblance (Fig. 4b) are better than the ones obtained using the semblance (Fig. 4a), as the beampower maxima are more focused on the elliptical structure. As shown by Schimmel & Paulssen (1997), the root semblance gives more weight to the phase coherency and leads to stable results. Correlation-based coherency (Figs 4c and d) also gives good results, with all the scattered energy focused on the elliptical body. As expected, time-corrected correlation gives a smoother image, as this measure finds coherency for the neighbouring points of the point that actually produces scattered energy. It is consequently more stable, but suffers from a smearing effect at the edge of the structures as well as highly expensive computational costs.

The deciding factor in applying a coherency measure method is whether or not it reconstructs mislocated energy that can be misinterpreted as geological structure. The semblance measures lead to the strongest artefacts, while the correlation and the second root semblance do not produce such patterns. This means that the waves involved in the artefacts have a high amplitude but do not share perfectly similar waveforms (e.g. Schimmel & Paulssen 1997). When looking at the data, the artefacts are mainly produced by the mixed
<figure><img src="image_6.png" /><figcaption>**Figure 5.** (a) Results of beam-forming approach applied on the data computed in: (top panel) the reference model (without perturbation) and middle) the model A (with elliptical perturbation). The bottom panel is the difference between the two upper panels (i.e. same result as given in Fig. 4). Coherency between traces is measured by correlation. The colour scale gives the beampower, with dark colours indicating locations with highly scattered energy. The dashed line indicates the elliptical perturbation inserted into the velocity model. (b) Receiver gather with highlighted waveforms (direct *P* wave) used in the reconstruction of the scatter point at the position (3123 km, 9 km) and differential receiver gather with highlighted portion of the data used for the scattering point at the position (3105 km, 11.5 km), that is, the elliptical structure.</figcaption></figure>

coherency analysis of the direct *P* waves together with the reflected *P* waves. A coherency measure that is sensitive to the phase of signals rather than to the amplitude appears to be a better candidate for imaging scattering structures. These tests show that both the second root semblance and the correlation measures lead to an accurate location of the scattering structure, which is also quite well defined. Therefore, we will use the correlation method (as described by eq. 2) to measure the coherency between traces.

Fig. 5 shows the beampower obtained using data computed in the test reference model and in the model with perturbations. Both images show a high-energy patch in the middle of the profile. It is produced by the coherency analysis of the direct *P* waves, as the wave generated on both sides of the island are crossing this area. This is confirmed on the data (Fig. 5b): the waveforms used to reconstruct this area correspond to the direct *P* waves. When taking the differences between the two images, this high-energy area is efficiently suppressed, revealing the perturbation features.

The same method is applied to the second synthetic data set, computed from the model B (fractal perturbation, Fig. 6). While the elliptical body produces reflected waves, the rough fractal structure will mainly generate scattered waves. The beampower result is shown in Fig. 6. While the ellipse could be identified (see Fig. 4), the geometry of the fractal perturbation is badly reconstructed. The scattered wavefield is weaker than the reflected one, and constructive and destructive interferences occur and lead to a very different amplitude response among the scatter points. However, scattered energy is correctly located (Fig. 6). We conclude that, even if it is not possible with this method to accurately reconstruct the geometry of strongly scattering structures, the location of the coherent energy is reliable.

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Results of beam forming approach to reconstruct the fractal perturbation in model B (grey line). The velocity model is the same as in Fig. 3, with a fractal structures instead of the elliptical feature. Coherency among traces is measured by correlation. The colour scale gives the beampower, with dark colours indicating locations with highly scattered energy.</figcaption></figure>

For the application to real data, we will use a tomographic model, which can be slightly inaccurate. Hence, we first test the sensitivity of this method to imperfect velocity models by modifying the traveltimes used to shift the traces in our synthetic tests. Random modifications of ±0.2 per cent, ±2 per cent and ±5 per cent are applied to the traveltimes computed for each source–scatter–receiver paths. The set-up is the same as in Fig. 3, and the analysis is performed using model A. The modified times are then used in the beam method, using a correlation measure, as well as the corrected correlation measure, as it should be less sensitive to velocity errors. Results are shown in Fig. 7. With ±0.2 per cent time difference, results are as good as without errors, while with ±2 per cent and ±5 per cent the structure is not reconstructed at all. In this case, the resulting beampower image does not show clear coherent structures, but a random distribution of energy throughout the profile. This means that the method is very sensitive to the errors in the
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Beam forming results with traveltimes randomly modified by (a) ±0.2 per cent, (b) ±2 per cent and (c) ±5 per cent, to test the influence of a wrong velocity model. The underlying model we are trying to reconstruct is model A (elliptical perturbation, Fig. 3), also used in Fig. 4. The coherency between traces is measured using the correlation. (d) same as (b) (±2 per cent modification) with the corrected correlation. The grey dashed line indicates the elliptical perturbation inserted into the velocity model.</figcaption></figure>

velocity model, which shift waves incoherently in the beam imaging. This has the positive effect of mitigating against the construction of spurious structures in the image. For the reference scattering image, traveltimes are kept not perturbed. We check if the inability to reconstruct the model is not due to the correction by the reference image by using modified traveltimes to build it. The results confirm that the imperfect reconstruction is due to the traveltime perturbation, and not to the reference model. The same tests are carried out using the corrected correlation for ±2 per cent (Fig. 7d). The result does not improve compared to those obtained using the correlation method. Therefore, and because the corrected correlation is computationally expensive, we only use the correlation method for our real data analysis.

In summary, from these numerical tests, we show that (1) semblance should be avoided and replaced by the nth root semblance or a correlation measure; (2) it is hazardous to interpret the geometry of the structures, as these cannot always be properly reconstructed. However, the average location is accurate, as all the scattering energy focuses on the exact structure location; and (3) an incorrect velocity model does not lead to the reconstruction of spurious structures.

# 6 RESULTS AND DISCUSSION: MAGMATIC STORAGE BENEATH TENERIFE ISLAND?

## 6.1 Results

The process described in Section 4 is applied to the data set recorded along the south-north and west-east profiles on Tenerife island. The

grids of scattering points have a spacing of 200 m in both horizontal and vertical directions and range from 0 to 14 km b.s.l. They are restricted to the cross-sections that are located on the island, that is, from 3100 to 3140 km UTM and 321 to 361 km UTM for the S-N and W-E profiles, respectively (each profile is 40 km long). The traveltimes and the synthetic data are computed using a smooth version of the tomographic model derived by García-Yeguas et al. (2012). The beam-forming method is applied to both recorded and synthetic data, the results from the latter are used as reference. As modelling is performed in 2-D, the seismic amplitudes are multiplied by $\sqrt{t}$ to correct for the geometrical effect differences between a 2-D and 3-D media (Operto et al. 2004). As there are no direct comparisons of waveforms propagating in 2-D and 3-D media, and as all the coherency measures are normalized, this assumption is not a critical issue. The beampower results obtained from the data are then corrected by the reference values, as explained in Section 4. All source-receiver pairs (2422 and 2551 pairs for the S-N and W-E profiles, respectively) are used.

The tomographic models and beampower images are shown in Fig. 8, for both profiles. Small patches of coherent scattered energy appears at different locations in the figures. By comparing to the synthetic tests, it appears that these effects can be due to the fact that we are using an imperfect velocity model. However, two main structures can be clearly identified. On the south-north profile, most of the scattered energy seems to be located in the north of the island (between 3115 and 3140 km UTM), at a depth between 5 and 9 km b.s.l. On the west-east profile, most of the reconstructed energy is located in the centre of the profile (335-345 km UTM), at shallower depth (2-4 km b.s.l.). This suggests the presence of two scattering structures, one shallow, less than 10 km wide, centred beneath Teide summit, and a deeper (6-10 km b.s.l.), larger one (up to 25 km wide) below the north of the island. While we are confident about the locations of the structures, as demonstrated by the synthetic tests, the geometries cannot be interpreted with confidence.

The tomography model was obtained using the same data set, including all shots around the island and all the available stations. This model is thus the best possible, though imperfect, image of the subsurface properties. However, when the model contains significant errors, the synthetic tests showed us that the scattered waves are not coherent anymore, and, as a result, the beam image does not show any high scattering energy area. The fact that we obtained well-defined structures in the reconstructed image means that the tomographic model is good enough, which gives us confidence in the results. Moreover, more than 2000 traces are used for each profile: the use of this amount of data stabilizes the results. Another possible source of error is noise in the data. To test its influence on the results and verify if contamination has occurred on the final image, we apply the beam-forming analysis to the noise recorded on Tenerife island in the quiet time interval between airgun shots. The resulting images are shown in Fig. 9, for both profiles. First, the beampower reconstructed from the noise is approximately five times smaller than the one obtained using the active source data. A few higher energy patches appear in the resulting images but these areas do not share a common location with the structures observed in Fig. 8. Therefore we can conclude that the latter are not produced by the recorded noise, which mainly leads to broadly distributed small energy along the profiles.

## 6.2 Possible magmatic storage

The presence of a magmatic storage area could be associated with either positive or negative anomalies in the P wave tomography
<figure><img src="image_9.png" /><figcaption>**Figure 8.** (a) Tomographic model obtained by García-Yeguas *et al.* (2012) for the south-north profile. The white rectangle marks the area where the beam-forming analysis is performed. (b) Beam-forming results. The colour scale gives the beampower, with the dark colours indicating locations with highly scattered energy. (c) and (d) same as (a) and (b), respectively, for the west-east profile. A 2-D median filter (third order) has been applied to the scattering images (b) and (d).</figcaption></figure>

velocity model (Lees 2007). In either case, scattered waves would be generated by the magmatic structures. The difficulty in identifying a magmatic body below volcanoes may lie in the imaging methodology used (complicated by a strongly heterogeneous medium and sharp topography), but may also be caused by the fact that a possible magma chamber is a complex set of dykes, sills and fractures rather than a fully filled continuous body. In this case, the velocity contrast may not be sufficient to be seen by tomography, which has a resolution limited of the grid size (700 m in the Tenerife model, García-Yeguas *et al.* 2012). However, mesoscale structures (dykes and fractures from a few hundreds to thousands of metres long) scatter and reflect waves. The results shown in Fig. 8 are in agreement with this fact: no clear indication of magmatic storage appears in the tomographic model, but coherent scattering is seen in localized areas in the beampower plots. In other words, it appears that those structures are scattering energy without generating any detectable velocity anomaly, which implies strong contrasts and limited spatial scale. A possibility is that those structures are fractures filled by volatiles, which lead to high-velocity contrasts and high-amplitude reflected waves. This is similar to the bright spot effects in exploration seismology, such as observed in CO₂ geological storage monitoring (Arts *et al.* 2004). Hence, the recon-

structed structures could possibly be associated with the gas-filled upper part of the magma storage system. Another possibility is that the scattered waves are produced by a complex set of dykes and sills filled with magma. The magma can indeed produces velocity perturbations of up to 30 per cent (Lees 2007). In both cases, the observed anomalies can be linked to magmatic storage areas.

Petrological studies carried out on Tenerife island infer the possibility of two magma chambers, one deep basaltic and one shallow phonolitic magma chamber. The latter is associated with the Teide-Pico Viejo system (Ablay *et al.* 1998; Martí *et al.* 2008). Araña *et al.* (2000) detected a magnetic anomaly, in agreement with gravity measures, which is interpreted as the top of magmatic body at 5.7 km b.s.l. (and extending down to 12 km b.s.l.), 40 km wide in the W-E direction and narrower (10 km) in the S-N direction. This is supported by analysis from Blanco-Montenegro *et al.* (2011) who detected magnetic anomalies extending to 9 km b.s.l. beneath the northern part of the island. The recent activity in 2004 also occurred in the northern part of the Cañadas caldera. VT events are located by Cerdeña *et al.* (2011) in two main clusters, with one at the northwest of the Teide-Pico Viejo system mainly between 5 and 10 km b.s.l. deep. Almendros *et al.* (2007) found VT and Long-Period (LP) locations at similar, but deeper location (14 km b.s.l.). The depth and
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Beam power response of the noise recorded on Tenerife island for the S–N (top panel) and the W–E profile (bottom panel). The colour scale gives the beampower, with the dark colours indicating locations with highly scattered energy. A 2-D median filter (third order) has been applied. No clear structures can be identified.</figcaption></figure>

the location of this anomaly are consistent with our observations of an anomaly between 6 and 10 km b.s.l. depth. However, we find that this structure is elongated in the north–south direction. As it is not clearly visible in the west–east profile, the structure is either rather small in this direction, or the profile we used is located to the south of the main scattering body.

Using gravity methods, Ablay & Kearey (2000) identified a high-density body, extending from 328 to 342 km UTM easting. It is interpreted as the upper part of the mafic plutonic complex, the Boca Tauce volcano. This structure forms the high-density core of the island (Martí *et al.* 2008), and can be clearly seen in the tomographic image (Fig. 8), and consequently is not likely to produce the scattered energy observed in Fig. 8. Gottsmann *et al.* (2008) used gravity methods to infer a 10-km-wide cylindrical body below Teide-Pico Viejo, extending from the sea level to *ca.* 3 km b.s.l. They interpret this structure as the current plumbing system of the volcanic complex. The presence of a shallow magmatic chamber is confirmed by the occurrence of several caldera collapses (Ablay & Martí 2000). Martí *et al.* (2008), using the petrology and the geophysical information available, infer that phonolitic magma has to be stored in a small chamber at 1–2 km depth b.s.l. The structure we identified is in agreement with this shallow phonolitic magma chamber theory.

Hence, the two structures we identified through the beam-forming method (one shallow beneath Las Cañadas, and one deeper below the northern part of the island) seem to be consistent with the magmatic bodies inferred by petrological and other non-seismic geophysical studies.

## 7 CONCLUSION

Tenerife island is a volcanic island dominated by the Teide-Pico Viejo system. An active seismic experiment was performed in 2007 to obtain a traveltime tomography (García-Yeguas *et al.* 2012) of the island to better understand the Teide volcanic system, which had shown signs of reawakening in 2004, and identify the presence of possible magma chamber. Using the obtained tomographic model as a starting point, herein we perform a beam-forming analysis on two

2-D orthogonal profiles recorded across the island during the same experiment. The aim is to image the scattering structures beneath the volcano. Through numerical tests, we show that (1) nth root semblance (rather than semblance) or correlation should be used to measure coherency between scattered waves, (2) an incorrect velocity model destroys the coherency among the wavefields, but importantly does not lead to spurious structures in the image and (3) even if the full geometry of the structures cannot be retrieved, the location of the scattering structures is reliable. The beam-forming method thus proves to be a rather simple, but powerful method to image such structures beneath volcanoes. As the magmatic storage can produce strong S-wave anomalies (Lees 2007), a similar study for scattered S waves is of great importance to accurately describe and interpret with less ambiguity the structures present beneath volcanoes. Where S-wave tomography models are available (they were not in this study), it would be straightforward to extend the application of this method to the analysis of reflected S and converted P-to-S waves.

Applying the beam-forming method on Tenerife data set, two main structures can be identified: (1) one shallow (2–4 km b.s.l., beneath Las Cañadas caldera), and a deeper (6–10 km b.s.l.), wider structure, in the northern part of the island. These locations are consistent with possible magmatic chambers identified by other geophysical and petrography studies. We can hypothesize that the shallow structure corresponds to the phonolitic plumbing system, which is feeding the Las Cañadas and Teide-Pico Viejo complex. The deeper structure, in the northern part of the island, is linked to the mafic reservoir which is responsible for the basaltic eruptions. These scattering structures could be associated with fractures filled with volatiles or with complex sets of dykes and sills, which produce high-amplitude reflected waves from localized structures. It is important to note that neither of these structures can be seen in the velocity tomography model. Information on possible magma storage can be found in the scattered wavefield rather than in the traveltimes. This suggests that these structures have strong velocity contrasts and, individually, are spatially limited, suggestive of a population of fractures and/or dyke-like structures.

## ACKNOWLEDGMENTS

This work was partially funded by the Spanish Projects CGL2008-01660 and CGL2011-29499-C02-01. The seismic instruments and data management facilities were provided under loan no. 812 by SEIS-UK at the University of Leicester to the University of Liverpool. Data collected will be available through the IRIS Data Management Center. The facilities of SEIS-UK are supported by the Natural Environment Research Council under Agreement No. R8/H10/64. LDEB and FM were funded by Science Foundation Ireland (SFI grant 08/RFP/GEO1486), with additional support from the Department of Communications, Energy and Natural Resources (Ireland) under the National Geosciences programme 2007-2013. Computational facilities from SFI/HEA ICHEC (Irish Centre for High End Computing) are acknowledged. Comments made by the editor J. Wassermann, by A. Rietbrock and an anonymous reviewer greatly contributed to improve the manuscript.

## REFERENCES

* Ablay, G.J. & Kearey, P., 2000. Gravity constraints on the structure and volcanic evolution of Tenerife, Canary Islands, *J. geophys. Res.*, **105**(B3), 5783–5796.
Ablay, G.J. & Martí, J., 2000. Stratigraphy, structure, and volcanic evolution of the Pico Teide Pico Viejo formation, Tenerife, Canary Islands, *J. Volc. Geotherm. Res.*, **103**, 175–208.

Ablay, G.J., Carroll, M.R., Palmer, M.R., Martí, J. & Sparks, R.S.J., 1998. Basanite-phonolite lineages of the Teide Pico Viejo volcanic complex, Tenerife, Canary Islands, *J. Petrol.*, **39**(5), 905–936.

Almendros, J. & Chouet, B., 2003. Performance of the radial semblance method for the location of very long period volcanic signals, *Bull. seism. Soc. Am.*, **93**(5), 1890–1903.

Almendros, J., Ibáñez, J.M., Carmona, E. & Zandomeneghi, D., 2007. Array analyses of volcanic earthquakes and tremor recorded for Las Cañadas caldera (Tenerife Island, Spain), during the May 2004 seismic activation of Teide volcano, *J. Volc. Geotherm. Res.*, **160**, 285–299.

Araña, V., Camacho, A.G., Garcia, A., Montesinos, F.G., Blanco, I., Vieira, R. & Felpeto, A., 2000. Internal structure of Tenerife (Canary Islands) based on gravity, aeromagnetic and volcanological data, *J. Volc. Geotherm. Res.*, **103**, 43–64.

Arts, R., Eiken, O., Chadwick, A., Zweigel, P., des Meer, L.V. & Zinszner, B., 2004. Monitoring of CO₂ injected at Sleipner using time-laps seismic data, *Energy*, **29**, 1323–1392.

Berger, P., Got, J.-L., Valdés González, C. & Monteiller, V., 2011. Seismic tomography at Popocatépetl volcano, Mexico, *J. Volc. Geotherm. Res.*, **200**(3–4), 234–244.

Blanco-Montenegro, I., Nicolosi, I., Pignatelli, A., García, A. & Chiappini, M., 2011. New evidence about the structure and growth of ocean island volcanoes from aeromagnetic data: the case of Tenerife, Canary Islands, *J. geophys. Res.*, **116**(B3), B03102, doi:10.1029/2010JB007646.

Carracedo, J.C., Troll, V.R., Torrado, F. J.P., Badiola, E.R., Machín, A.H., Paris, R., Guillou, H. & Scaillet, S., 2006. Recent unrest at Canary islands? Teide volcano?, *EOS, Trans. Am. geophys. Un.*, **87**(43), 462–465.

Cerdeña, I.D., del Fresno, C. & Rivera, L., 2011. New insight on the increasing seismicity during tenerife’s 2004 volcanic reactivation, *J. Volc. Geotherm. Res.*, **206**(1-2), 15–29.

Chouet, B., 2003. Volcano seismology, *Pure appl. geophys.*, **160**, 739–768.

De Barros, L., Pedersen, H.A., Métaxian, J.-P., Valdès-Gonzalez, C. & Lesage, P., 2008. Crustal structure below Popocatépetl Volcano (Mexico) from analysis of Rayleigh waves, *J. Volc. Geotherm. Res.*, **170**(1-2), 1–7.

De Natale, G., Troise, C., Trigila, R., Dolfi, D. & Chiarabba, C., 2004. Seismicity and 3-D substructure at Somma-Vesuvius volcano: evidence for magma quenching, *Earth planet. Sci. Lett.*, **221**(1–4), 181–196.

Del Pezzo, E., Bianco, F., De Siena, L. & Zollo, A., 2006. Small scale shallow attenuation structure at Mt. Vesuvius, Italy, *Phys. Earth planet. Inter.*, **157**(3–4), 257–268.

García-Yeguas, A., Koulakov, I., Ibáñez, J. & Rietbrock, A., 2012. First high resolution P wave velocity structure beneath Tenerife Island (Canary Islands, Spain), *J. geophys. Res.*, in press.

Gottsmann, J., Camacho, A.G., Martí, J., Wooller, L., Fernández, J., García, A. & Rymer, H., 2008. Shallow structure beneath the central volcanic complex of Tenerife from new gravity data: implications for its evolution and recent reactivation, *Phys. Earth planet. Inter.*, **(3–4)**, 212–230.

Hedlin, M.A.H., Minster, J.B. & Orcutt, J.A., 1991. Beam-stack imaging using a small aperture array, *Geophys. Res. Lett.*, **18**(9), 1771–1774.

Ibáñez, J.M., Rietbrock, A. & García-Yeguas, A., 2008. Imaging an active volcano edifice at Tenerife island, Spain, *EOS, Trans. Am. geophys. Un.*, **89**(32), doi:10.1029/2008EO320001.

Kanasewich, E.R., Hemmings, C.D. & Alpaslan, T., 1973. Nth-root stack nonlinear multichannel filter, *Geophysics*, **38**(2), 327–338.

Kent, G. et al., 2000. Evidence from three-dimensional reflectivity images for enhanced melt supply beneath mid-ocean-ridge discontinuities, *Nature*, **406**, 614–618.

Kruger, F., Weber, M., Scherbaum, F. & Schlittenhardt, J., 1993. Double beam analysis of anomalies in the core-mantle boundary region, *Geophys. Res. Lett.*, **20**(14), 1475–1478.

Kruger, F., Scherbaum, F., Weber, M. & Schlittenhardt, J., 1996. Analysis of asymmetric multipathing with a generalization of the double-beam method, *Bull. seism. Soc. Am.*, **86**(3), 737–749.

Kumagai, H., Saito, T., O'Brien, G. & Yamashina, T., 2011. Characterization of scattered seismic wavefields simulated in heterogeneous media with topography, *J. geophys. Res.*, **116**(B3), B03308, doi:10.1029/2010JB007718.

Lacoss, R., Kelly, E. & Toksoz, M., 1969. Estimation of seismic noise structure using arrays, *Geophysics*, **34**, 21–38.

Lay, T., 1987. Analysis of near-source contributions to early P-wave coda for underground explosions. III. Inversion for isotropic scatterers, *Bull. seism. Soc. Am.*, **77**(5), 1767–1783.

Lees, J., 2007. Seismic tomography of magmatic systems, *J. Volc. Geotherm. Res.*, **167**(1–4), 37–56.

Maercklin, N., 2008. Conception, verification, and application of innovative techniques to study active volcanoes, in *Seismic Scatterer Imaging Using Shot Array Beamforming: Method and Application to the Campi Flegrei Caldera*, pp. 261–268, eds, Marzocchi, W. & Zollo, A., Doppiavoce, Naples.

Maercklin, N., Haberland, C., Ryberg, T., Weber, M., Bartov, Y. & DESERT Group, 2004. Imaging the Dead Sea transform with scattered seismic waves, *Geophys. J. Int.*, **158**(1), 179–186.

Marsh, B., 2000. Magma chambers, in *Encyclopedia of Volcanoes*, pp. 191–206, Academic Press, San Diego, CA.

Martí, J. & Geyer, A., 2009. Central vs flank eruptions at Teide Pico Viejo twin stratovolcanoes (Tenerife, Canary Islands), *J. Volc. Geotherm. Res.*, **181**(1–2), 47–60.

Martí, J., Mitjavila, J. & Araña, V., 1994. Stratigraphy, structure and geochronology of the Las Cañadas caldera (Tenerife, Canary Islands), *Geol. Mag.*, **131**, 715–727.

Martí, J., Geyer, A., Andujar, J., Teixid, F. & Costa, F., 2008. Assessing the potential for future explosive activity from Teide Pico Viejo stratovolcanoes (Tenerife, Canary Islands), *J. Volc. Geotherm. Res.*, **178**, 529–542.

Martí, J., Ortiz, R., Gottsmann, J., Garcia, A. & De La Cruz-Reyna, S., 2009. Characterising unrest during the reawakening of the central volcanic complex on Tenerife, Canary Islands, 2004-2005, and implications for assessing hazards and risk mitigation, *J. Volc. Geotherm. Res.*, **182**(1-2), 23–33.

Martínez-Arévalo, C., Patanè, D., Rietbrock, A. & Ibáñez, J.M., 2005. The intrusive process leading to the Mt. Etna 2001 flank eruption: constraints from 3-D attenuation tomography, *Geophys. Res. Lett.*, **32**(21), L21309, doi:10.1029/2005GL023736.

Monteiller, V., Got, J.-L., Virieux, J. & Okubo, P., 2005. An efficient algorithm for double-difference tomography and location in heterogeneous media, with an application to the Kilauea volcano, *J. geophys. Res.*, **110**, B12306, doi:10.1029/2004JB003466.

Neidell, N. & Taner, M., 1971. Semblance and other coherency measures for multichannel data, *Geophysics*, **36**, 482–497.

O'Brien, G.S. & Bean, C.J., 2009. Volcano topography, structure and intrinsic attenuation: their relative influences on a simulated 3D visco-elastic wavefield, *J. Volc. Geotherm. Res.*, **183**(1–2), 122–136.

Operto, S., Ravaut, C., Improta, L., Virieux, J., Herrero, A. & Dell'Aversana, P., 2004. Quantitative imaging of complex structures from multi-fold wide aperture seismic data, *Geophys. Prospect.*, **5**, 625–651.

Park, J., Morgan, J., Zelt, C., Okubo, P., Peters, L. & Benesh, N., 2007. Comparative velocity structure of active Hawaiian volcanoes from 3-D onshore-offshore seismic tomography, *Earth planet. Sci. Lett.*, **259**(3–4), 500–516.

Patanè, D., Barberi, G., Cocina, O., De Gori, P. & Chiarabba, C., 2006. Time-resolved seismic tomography detects magma intrusions at Mount Etna, *Science*, **313**(5788), 821–823.

Prôno, E., Battaglia, J., Monteiller, V., Got, J.-L. & Ferrazzini, V., 2009. P-wave velocity structure of Piton de la Fournaise volcano deduced from seismic data recorded between 1996 and 1999, *J. Volc. Geotherm. Res.*, **184**(1–2), 49–62.

Rietbrock, A. & Scherbaum, F., 1999. Crustal scattering at the KTB from combined microearthquake and receiver analysis, *Geophy. J. Int.*, **136**(1), 57–67.

Ringdal, F. & Husebye, E.S., 1982. Application of arrays in the detection, location, and identification of seismic events, *Bull. seism. Soc. Am.*, **72**(B6), S201–S224.
Scherbaum, F., Krüger, F. & Weber, M., 1997. Double beam imaging: mapping lower mantle heterogeneities using combinations of source and receiver arrays, *J. geophys. Res.*, **102**(B1), 507–522.

Schimmel, M. & Paulssen, H., 1997. Noise reduction and detection of weak, coherent signals through phase-weighted stacks, *Geophy. J. Int.*, **130**(2), 497–505.

Singh, S. et al., 2006. Discovery of a magma chamber and faults beneath a mid-Atlantic ridge hydrothermal field, *Nature*, **442**, 1029–1032.

Spudich, P. & Bostwick, T., 1987. Studies of the seismic coda using an earthquake cluster as a deeply buried seismograph array, *J. geophys. Res.*, **92**(B10), 10526–10546.

Vanorio, T., Virieux, J., Capuano, P. & Russo, G., 2005. Three-dimensional seismic tomography from P wave and S wave microearthquake travel

times and rock physics characterization of the Campi Flegrei Caldera, *J. geophys. Res.*, **110**(B3), B03201, doi:10.1029/2004JB003102.

Yilmaz, O., 2001. *Seismic Data Analysis: Processing, Inversion and Interpretation of Seismic Data*, Society of Exploration Geophysicists, Tulsa, OK.

Zandomeneghi, D., Barclay, A., Almendros, J., Ibaez, J.M., Wilcock, W.S.D. & Ben-Zvi, T., 2009. Crustal structure of deception island volcano from p wave seismic tomography: tectonic and volcanic implications, *J. geophys. Res.*, **114**(B6), B06310, doi:10.1029/2008JB006119.

Zelt, C.A. & Smith, R.B., 1992. Seismic traveltime inversion for 2-d crustal velocity structure, *Geophy. J. Int.*, **108**, 16–34.

Zollo, A. et al., 2003. Evidence for buried rim of Campi Flegrei caldera from 3-D active seismic imaging, *Geophys. Res. Lett.*, **30**(19), doi:10.1029/2003GL018173.