Seismic tomography is a technique that uses seismic data to build 2-D and 3-D models of the Earth’s interior. The most common application measures travel times and amplitudes of P- and S-waves produced by earthquakes or explosions and recorded at seismic stations. Tomography provides a forward model of how these waves have propagated through the Earth (generally, these waves follow “rays”) and uses data and models to invert, i.e., reconstruct in space for quantities like seismic velocity and attenuation on seismic wave-packets (Rawlinson et al., 2010; Koulakov, 2013; Rawlinson and Spakman, 2016). As the only technique available to scan the Earth at all scales, seismic tomography can be combined with other disciplines providing images and indirect measurements of Earth materials, from the Earth core to the shallow surface (Rawlinson et al., 2010). These techniques work when the amount of coherent waves preserved through propagation is sufficient across the seismic network and have found extensive success when applied to magmatic systems (Lees, 2007; Zandomeneghi et al., 2009; Koulakov, 2013; De Siena et al., 2017b). However, heterogeneity in these media produces multiple-scattering and diffusive wavefields that progressively deteriorate coherent-wave detection with distance. Researchers have thus developed tomographic techniques that in their forward model the sensitivity of coda waves (the waves scattered by this same heterogeneity) to the Earth structures (Del Pezzo et al., 2016, 2018), using them to image volcanoes (Prudencio, 2013; Prudencio et al., 2013; De Siena et al., 2017a).

Active seismic experiments offer ideal coverage to obtain highly-resolved tomographic images of volcanoes. The first applications date back to the 1980s, when travel-time tomography was applied to active data recorded across the Newberry and Medice Lake, in the Oregon and California Cascade Range, respectively (Achauer et al., 1988; Evans and Zucca, 1988). In the 1990s, active seismic experiments funded by national governments and international organisations have increased resolution on crustal volcanic structures. The feeding systems of Mt. Vesuvius [TOMOVES and MUREVES, Gasparini and Group (1998)] and Campi Flegrei [SERAPIS, Zollo et al. (2008)] are largely interpreted thanks to these surveys. In the last two decades active experiments have resolved the interior of Deception Island (TOMODEC 2005, Prudencio et al. (2013; 2015a)), Tenerife [TOMTEIDEVS, García-Yeguas et al. (2012); Prudencio (2013); Prudencio et al. (2015b)], Etna [TOMOETNA, Ibáñez et al. (2016)] and Santorini [PROTEUS, Hooft et al. (2019)]. Onshore, the most prominent active seismic experiments are iMuSH, at Mount St. Helens volcano, and MED-SUV, across several volcanic areas in the Mediterranean Sea. These are the only massive surveys where geophysical responses from both seismic and potential fields were collected (De Landro et al., 2017; Bedrosian et al., 2018; Siniscalchi et al., 2019; Ulberg et al., 2020). The corresponding images provide unique constraints on the shallow feeding systems and processes at these volcanoes, especially if used jointly with passive data (Battaglia et al., 2008).

Both coherent and incoherent waves have been extensively used at Deception Island (South Shetland Islands, Antarctica) primarily thanks to the TOMODEC2005 experiment, carried out during the austral summer of 2005 at Deception Island volcano (Ibáñez et al., 2017) (Figure 1). Over half a million seismic waveforms were produced by shots from a boat circling the island, and, especially, entering the central caldera, recorded by seismometers both on land and at the bottom of the sea (Prudencio et al., 2013). Deception Island is an ideal target for such an experiment for three primary reasons: 1) its morphology, with the navigable inner bay that allows ideal ray coverage on the central caldera; 2) the absence of anthropic noise (Lecocq et al., 2020); 3) a 360° spatial ray-coverage (Ibáñez et al., 2017), a unique feature when compared to other off-shore experiments, which generally provide directional ray coverage (Zollo et al., 2008). For these reasons, the corresponding travel-time, direct- and coda-wave attenuation show unprecedented resolution relative to other volcanic areas (Ben-Zvi et al., 2009; Zandomeneghi et al., 2009; Prudencio et al., 2013, 2015a; Del Pezzo et al., 2016, 2018). The dataset has a final, unique advantage: it is entirely available online (Ibáñez et al., 2017) and thus ideal for developing and testing automated and semi-automated data selection strategies.

There is uncertainty, or lack of benchmark strategies, about the best methods for the selection of waveforms necessary to perform seismic tomography using volcanic datasets; this is contrary to what happens at the global scale (Valentine and Woodhouse, 2010). Data selection strategies have generally been applied to passive seismo-volcanic datasets, and, when using active seismic data, they focus on travel times more than on amplitudes (De Landro et al., 2017). Strategies focused on amplitudes and, particularly, coda waves become increasingly necessary, as attenuation tomography is today a standard for imaging volcanoes (Prudencio et al., 2013, 2015a; Serlenga et al., 2016; De Siena et al., 2017a). Given the quantity of seismic traces produced by these experiments, an automatic or quasi-automatic system preliminary to inversion is necessary.

Automatic data processing strategies that focus on both travel-time and amplitude detection are generally applied to global seismic records (Valentine and Woodhouse, 2010). At smaller scales, they concentrate on phases and earthquake detection (Romero et al., 2016; Pang et al., 2019; Cortés et al., 2020). In volcanic media and especially in active seismic experiments, scattered waves dominate, with chaotic wave packets soon dominating seismic traces (Ibáñez et al., 2017). Semi-automated strategies provide better results in this case, especially if the expected medium heterogeneity is high (Bianco and Zaccarelli, 2009). These approaches combine two steps: dataset selection and application to imaging as a benchmark. The selection generally comprises several iterations, each depending on a single selection parameter (Cortés et al., 2015). It is necessary to set a-priori conditions specific to the inversion using standard practice from neural networks (Valentine and Woodhouse, 2010). Following each iteration, the obtained dataset must be applied to the corresponding imaging code to compare the results quantitatively with the same standard.

Given the increasing importance of coda techniques when imaging volcanoes, in this work, we present a quasi-automatic method for the selection of the best traces for coda-wave tomography using the TOMODEC2005 dataset. The open-access code MuRAT 1, previously applied in volcanic (Prudencio et al., 2015a; De Siena et al., 2016, 2017b; Gabrielli et al., 2020; Sketsiou et al., 2020) and tectonic setting (Borleanu et al., 2017; Napolitano et al., 2020), is used to perform data processing and inversion. The majority of the analysis focuses on reducing errors during the processing phases, reducing noise, increasing consistency between different node spacing and frequencies, and, especially, mitigating loss of information. We show that we reduce residuals by targeting the inversion results while increasing resolution on structures without any additional computational cost.

The TOMODEC2005 experiment (Ibáñez et al., 2017) produced data that have been extensively used at Deception Island for velocity (Ben-Zvi et al., 2009; Zandomeneghi et al., 2009), attenuation, and scattering tomography (Prudencio et al., 2013, 2015a; Del Pezzo et al., 2016, 2018). Seismic waves produced at the volcano suffer high scattering, and waveforms recorded at stations can be described by the diffusion equation (Prudencio et al., 2013). In any coda-attenuation study, the waveforms must respect at least one condition: the envelope has to decrease exponentially from its main peak amplitude (Sato et al., 2012). The coda decay is controlled by coda attenuation (or inverse coda quality factor, Q c − 1 ), which, in the diffusive approximation, is a measurement of seismic absorption, i.e., of the energy lost permanently by the interaction of the wave with the medium (Calvet and Margerin, 2013).

The starting dataset for this study is the subset of TOMODEC2005 used to perform attenuation tomography (Ibáñez et al., 2017). It comprises 20283 waveforms of different duration, each detected automatically knowing the time of the shot and corrected for the instrument function. The waveforms were then cut to 10 s from the time of the shot as signals always last less than this duration (Prudencio et al., 2015a). This dataset has been used to image seismic absorption by (Prudencio et al., 2013) with regular node spacing of 3 and 1 km. To verify the dataset quality, we made a first visual check (Valentine and Woodhouse, 2010). We subdivided the dataset into three groups: those with spurious signals (Figure 2A), those primarily characterised by null values (Figure 2B) and those visually identified as “good” (Figure 2C). These visual checks are standard in data selection. Still, the notion of “good” is qualitative, and operators should better follow standard seismological quality criteria, such as signal-to-noise ratio (Valentine and Woodhouse, 2010). However, these criteria do not target the ability of the dataset to perform a specific task, as imaging coda attenuation: such criteria could remove waveforms from datasets unnecessarily.

For each dataset, we obtain four Q c − 1 maps depending on frequencies at constant node spacing of 2 km (Figure 4) and 1 km (Figure 5). The input and output of the chequerboard test for each spatial resolution are shown in Figure 6.

While previous works show a single high-attenuation anomaly in the inner bay of the island (Prudencio et al., 2013; 2015a), our selection procedure highlights anomalies that spread from the centre of the bay to four areas: Telephone Bay, Pendulum Cove, Fumarole Bay and Entrance Point. All located across the caldera margin, some of these are the recognised locations of eruptions of the last century. One corresponds to the area occupied by the small island formed by the Telefon Bay eruption on 4 December 1967 (Smellie et al., 2002). A sub-glacial pyroclastic fissure produced by the 1969 eruption crosses Pendulum Cove (Smellie, 2002). On the contrary, Fumarole Bay is the location of the early volcanic edifice, the oldest exposed formation, and the area characterised by the highest ground temperatures across the island (Smellie et al., 2002). If interpreted with the extensive geochemical and geological records, these maps provide structural constraints to define the complex magmatic evolution of the island (Geyer et al., 2019).

The workflow highlights the need for inversion-specific data-selection strategies for volcanic imaging. If based only on standard attributes like the SNR ratio, the reduction in usable seismic waveforms leads to lower stability and robustness of the solution than with the current workflow (Table 1). A combination of standard and inversion-specific automatic selections followed by a visual analysis on the final dataset is a sequence that can be easily implemented on any active and passive dataset. The only time-consuming step remains the visual analysis. The automatic selection was unable to recognize about 10% of traces that were worsening residuals without adding useful information. While it was largely non-influential for a final interpretation of our results (compare panels d and e) in Figures 5D,E, 4, we cannot ensure that a fully automatic selection would be as successful with another dataset.

The sequence chosen to apply selection parameters is irrelevant if the last automatic selection is performed using the attribute that defines the data quality in the tomographic procedure. In our case, this parameter is the correlation coefficient of the coda decay. It is the quality of a P- or S-wave picks in travel-time tomography. A successful selection will substantially reduce the dataset without deteriorating the resolution in the area of interest, allowing to reconstruct smaller-scale details without redundancy in the data. In the procedure, we used the similarity between maps obtained at different frequencies as an additional marker of solution stability. This qualitative approach can be likely improved using a 2D correlation analysis.

While different frequencies theoretically map anomalies at different scales, these differences are usually minimal when using active data at volcanoes (Prudencio et al., 2013, 2015a,b; Ibáñez et al., 2017). However, a comparison with the local geology could confirm if the residual differences are linked to volcanic structures. The area resolved after the selection is smaller than the one resolved using Dataset 1 (Figure 6B) at 2 km node spacing, and generally interpreted in attenuation tomography studies (Prudencio et al., 2013, 2015a; Del Pezzo et al., 2018). However, our results demonstrate that this illumination likely overestimates the effective resolution potential of the dataset, primarily because the misfit between data and model is measured inadequately. The result can be the failed detection of primary features, as the NNW-SSE trend of the central high-attenuation anomaly. In Figure 3B, we compare the results obtained at the end of our iterations with the contour of the homogeneous bodies obtained by Prudencio et al. (2013), where such a trend is absent, especially at higher frequencies. NNW-SSE is the direction of extension of the Bransfield Strait underlined by the strike of several normal faults and alignment of the eruptive centres between 1967 and 1970 (Martí et al., 1996; Smellie et al., 2002; Geyer et al., 2019). This trend is just visible in velocity tomography maps (Zandomeneghi et al., 2009) and absent from previous attenuation studies (Prudencio et al., 2013; 2015a).

Seismic tomography models remain “just models” if they are not compared with the extensive interdisciplinary knowledge generally available at each volcano. This is particularly true at Deception Island, where researchers have provided some of the most extensive geological, geophysical, tectonic and geochemical reconstructions of magma feeding systems, especially for an ice-capped volcano that has erupted within the last century. Once compared to this extensive interdisciplinary knowledge within a geographic information system (GIS), the new model could provide a decisive piece to the puzzle of how this volcano has evolved, fed, and erupted. Apparently, the maps better fit geological and geochemical knowledge about the volcano (Martí et al., 1996; Geyer et al., 2019). They can thus provide an essential link between what we know from field and modelling studies and what we can image of the volcano interior.

This work is a stepping stone to produce more refined, sounder tomographic models of volcanoes. Active seismic experiments are rapidly becoming the standard in volcano tomography. Here, we demonstrate that efficient data selection strategies can retrieve stable and resolved anomalies that remain blurred without them. Once applied to other available active datasets, the workflow could reveal fine details of the volcanic structures that have remained invisible but could be essential to understand the volcano dynamics. Its application to passive seismicity is certainly more challenging, given the uncertainties associated with volcano seismicity detection and characterisation. At the same time, selection strategies targeting specific tomographic techniques could drastically improve detection of structures at hazardous volcanoes with decades of seismic records, like Mt. Etna (Italy) or Mount St. Helens (US), where data redundancy is most likely.

An inversion-specific semi-automated selection strategy applied to an active seismic dataset recorded at Deception volcano improves stability and robustness of tomographic maps exploiting amplitude information. The combination of standard and inversion-specific selection attributes improves images obtained with different node spacing and at multiple frequencies, showing a consistent picture of the volcano interior. These maps are compatible with existing attenuation maps. The inversion for coda attenuation using the final dataset shows higher resolution on the centre of the caldera at the expense of lateral illumination. The images highlight sites of eruptions and extensional trends, which, if confirmed by comparison with the extensive geological and geochemical records, could better link the present state of the caldera with its geological past.

Volcano tomography using active datasets is now well past its infancy. Our results demonstrate that it cannot simply rely on data-selection strategies taken from global tomography or seismic exploration in oil and gas. Attributes like coda attenuation are not part of standard processing procedures for seismic imaging. Our method focuses on object-oriented data quality, first evidencing the most prominent standard attributed biasing results (the percentage of zeroes across the dataset) and then selecting datasets to optimise the attributes used for the inversion. The results prove that using only standard quality assessments, as the signal-to-noise ratios, leads to trade-offs and unnecessary loss of information. The same procedure can be defined for any attribute used to do imaging, and the method is easily extendable to passive datasets that are more relevant for volcano-monitoring institutions. In the long-range, the adaptation of these workflows to the detection of seismic signals produced by the volcano could help the long-standing effort of the community to recognise pre-eruptive signals. Once combined with the increasing interdisciplinary knowledge available at Deception Island and other volcanoes, their implementation in passive tomography can help better constrain the interpretation of pre-existing structures, current processes and future eruptions by avoiding trade-offs from redundant data.