STEFANO SOLARINO

ISTITUTO NAZIONALE DI GEOFISICA E VULCANOLOGIA, INGV, ITALY, STEFANO.SOLARINO@INGV.IT

Abstract

Seismic tomography serves as a potent means to investigate the Earth’s internal structure. This method closely resembles the one used in medicine and relies on variations (in path, amplitude or velocity) that seismic rays experience when passing through an anomalous body, such as a rock with higher density or an area subjected to high temperature. In the case of seismic tomography, earthquakes act as emitters while seismic stations function as receivers. Unlike in medicine, where a specialized machine ensures a uniform distribution of emitters and receivers, seismic rays do not originate from all directions, sometimes necessitating techniques to achieve an acceptable distribution. In such instances, installing temporary seismic stations can be beneficial to increase the density of ray intersections and enhance the method’s resolution power. While it has been frequently claimed in literature that this operation enhances seismic tomography, this is the first time where an experiment conducted specifically quantifies this enhancement, allowing a critical analysis of the associated advantages. The present study takes advantage of the large availability of seismic stations resulting from three projects recently conducted in the Alps in northern Italy.

Introduction

Seismic tomography stands out as a highly effective method for studying the Earth’s interior. It offers the unique ability to delve into depths beyond the reach of other investigative techniques, providing detailed insights into the structure of the crust and mantle in three dimensions. This knowledge is essential for a wide range of applications, from understanding seismic phenomena to mapping deep faults, from prospecting for valuable resources like minerals and oil to reconstructing the geodynamic evolution of a region. By harnessing seismic waves, whether naturally occurring or artificially induced, this approach mirrors the widely employed Computer Tomography (hereinafter CT) scanning technology in medical diagnostics in which a narrow beam of X-rays is aimed at a patient and quickly rotated around the body, producing signals that are processed by a series of x-ray detectors, located directly opposite the X-ray source. As the X-rays pass through the patient, they are transmitted to a computer to generate cross-sectional images, or “slices.” These slices are called tomographic images and can give a clinician more detailed information than conventional X-rays. Once a number of successive slices are collected by the machine’s computer, they can be digitally “stacked” together to form a three-dimensional (3D) image of the patient that allows for easier identification of basic structures as well as possible tumors or abnormalities (https://www.nibib.nih.gov/science-education/science-topics/computed-tomography-ct).

Seismic tomography employs a similar principle, but instead of X-rays, it utilizes seismic waves emitted by an earthquake and detected by seismic stations positioned on the Earth’s surface. The ray theory, extensively used in seismology and seismic exploration for modeling of high-frequency seismic body waves, allows the use of rays instead of waves. This simplifies the tracking of the expanding wavefront; in particular it helps us to track the paths in a medium along which energy propagates. The rays’ properties, such as their arrival times at the stations, are influenced by anomalies encountered along their path in comparison to a standard Earth model. By employing suitable algorithms, it becomes feasible to speculate on the source of these propagation variations and subsequently reconstruct the three-dimensional composition of the surveyed volume.

There are, however, important differences between medical and seismic tomography. In medical tomography, the human body is illuminated by rotating the ray source and using a series of detectors. This ensures a uniformly resolved image, where any observed anomalies are solely attributed to the presence of medical abnormalities. Factors related to the tomographic investigation process do not influence these abnormalities. In contrast, seismic tomography does not provide uniform coverage of the volume being investigated. It depends on the location of seismic events and stations. This lack of uniform illumination results in either a complete absence of information for a specific portion of the volume or in an uncertainty in recognizing the three-dimensional structure. The latter is especially cunning because the algorithm may yield an image that seems trustworthy until one delves into the specifics of its computation.

In order to address these challenges, two strategies have been adopted: either by increasing the amount of seismic ray sources or by enhancing the number of receiving stations.

The addition of rays from induced explosions in specific locations selected to cover areas devoid of seismic activity was a common practice in the 1980s and 1990s (Ansorge et al., 1992; Aichroth et al., 1992). However, this controlled source experimentation faced numerous drawbacks associated with the use of explosives, such as exorbitant costs, security concerns, and challenges in conducting explosions in densely populated regions. The application of such techniques is confined to the shallow crustal regions, limiting their utility to studies of the shallow subsurface. Furthermore, active seismic methods have predominantly been employed along profiles, resulting in a primarily 2D representation. To address this constraint and achieve a quasi-3D image, researchers have occasionally conducted multiple perpendicular profiles in small-scale areas (see, for example, Gasparini and Tomoves Working Group, 1998). Improving data accuracy is also achievable by integrating seismic instruments for capturing “lost or not detected” rays.

In the present paper, after a short explanation of the principles of seismic tomography, we discuss how the installation of additional seismic stations in the area to be investigated can enhance ray coverage. This discussion draws insights from recent seismic projects conducted in the Western Alps, Northern Italy.

Seismic tomography

Since the pioneering seismic tomography experiment conducted by Aki and Lee in 1976 and Aki et al. in 1977, the technique for investigating the Earth’s interior has evolved significantly. This evolution has been closely tied to the increasing availability of seismograms and advancements in computer technology, allowing for more sophisticated analyses. Curtis and Snieder (2002) provide a comprehensive overview of the history of seismic tomography.

Over time, seismic tomography has undergone various changes, including modifications in the size of the analyzed area, the depth explored by the method, and the volume of data processed during tomographic procedures. Methodological improvements have enabled the incorporation of additional information beyond the arrival times of seismic waves, such as, for example, wave amplitudes. The advent of powerful computers has facilitated the utilization of larger datasets in experiments and the implementation of complex algorithms for precise seismic ray path reconstructions. This has replaced the simplistic geometric models used in the early stages of tomography, offering valuable insights into the Earth’s internal structure.

As a result, a multitude of tomographic approaches have emerged, differing in the scope of the investigated volume, the type of data utilized (ranging from basic P and S wave arrival times to complete seismograms), and the computational strategies employed to solve the problem. Nonetheless, all these methods share a common foundation: the ability to refine an initial Earth model based on general knowledge and imposed a priori from the user through iterative adjustments that consider the constraints imposed by the input data.

Generally speaking, the volume under investigation is divided into layers, with each layer consisting of blocks or nodes. During each iteration, the algorithm updates the values of the unknowns, such as, for example, the velocity values for P and S phases for each node or block within the model, and then compares the current solution with the previous one in order to seek enhancements relative to a predefined parameter that is specific to the type of tomography being utilized. If an enhancement is achieved, the algorithm retains this solution and continues to adjust it in subsequent iterations with the aim of further improving the reference parameter. Conversely, if the solution deteriorates, the algorithm reverts back to the previous solution and explores alternative options. The calculation process is terminated based on an end-iteration criterion, which is triggered either when a predetermined threshold is met in the reference parameter or when further improvement becomes unattainable due to constraints inherent in the initial data.

Local Earthquake Tomography (LET)

The technique known as LET (Local Earthquake Tomography) is employed to examine regional-scale crustal rock volumes and it is a widely used optimization technique for the joint reconstruction of P- and S-wave velocities and hypocentre parameters (e.g. Thurber 1983; Eberhart-Phillips 1986). This method relies on the seismic rays generated by earthquakes within the area under investigation. During each step of the process, the earthquakes’ locations are identified, the paths of the rays originating from each hypocenter are traced, and the velocity model is established according to the new earthquake locations. This involves one direct and two inverse problems. LET is a non-linear problem and several inversion methods have been developed to solve the problem in a linearized way or by consecutive iterative linearized steps. The comprehensive explanation of the technique and related issues fall beyond the scope of this paper. Readers who are interested in the tomographic process can explore a wealth of literature that discusses various aspects such as optimizing the direct problem, handling large amounts of data effectively, and utilizing mathematical algorithms to solve problems. One good starting point is the paper by Fichtner et al., 2024.

The primary focus of the present paper is to delve into the elements associated with the enhancement due to data availability. This is achieved by examining case studies of recent tomographic experiments conducted specifically in the northwestern Alps (Italy) area. The abundance of rays, particularly when it involves a more uniform sampling of the regions under study, plays a crucial role in constraining the unknown parameters, therefore, having a greater number of rays is essential for accurate results in this specific type of tomography. The comparison of the two panels shown in Figure 1 highlights the impact of deploying additional seismic stations on enhancing the number of rays available for addressing the tomographic procedure.

Figure 1. With an increased number of seismic stations (red triangles in panel b), there are more rays available for each earthquake (yellow stars) in the volume being examined. These rays, in turn, sample areas that would otherwise not be crossed by rays (panel a) and increase cross-firing, essential for reconstructing anomalies in a tomographic inversion.

The Cifalps and Alparray projects in the Western Alps

The Alps arose as a result of the collision of the African and Eurasian tectonic plates, in which the Tethys Ocean, which was formerly in between these continents, disappeared. The western Alps formed during Oligocene-Miocene times, as a consequence of the convergence and indentation of the Adriatic microplate toward Europe. The tectonic reconstruction of this part of the chain has been, and still is, under debate. It is essential to have a thorough understanding of the composition and distribution of crustal and sub-crustal rocks in order to obtain a reliable geological history. Over the past twelve years, studies have been carried out to enhance knowledge of the Western Alps. We here briefly describe those that involved the seismic stations used in this paper (Figure 2).

The projects Cifalps (China-Italy-France transect along the Alps) (Zhao et al., 2015) and Cifalps2 (Paul et al., 2022) took place in 2012-2013 and 2018-2020 respectively. They involved the installation of 56 seismic stations along two alignments, one running in a W-E direction and the other with a NW-SE orientation. The geometry of the installation was aiming at receiver function studies. However, in order to exploit the potential of seismic tomography, a few stations were deployed around the transect to improve the density of seismic stations and decrease the interdistance between sites (Figure 2).

The Alparray project (Hetényi et al., 2018) involved 24 institutions that collaborated to install the largest seismic temporary network ever in Europe. Between 2016 and 2019, a total of 312 stations were installed to monitor the seismic activity across the entire Alpine chain. Even after this period, several stations continued to record passive seismicity until 2022, and eventually, a few of them were converted into permanent installations.

Between 2018 and 2020, approximately 80 temporary stations operated in the Western Alps region as a result of the concurrent timing of two projects Cifalps and Alparray (Figure 2). The subsequent paragraph will focus on the analysis of the extent to which these instruments contribute to ray coverage and the enhancements they bring to seismic tomography.

Figure 2. Projects recently conducted in the Western Alps. Blue triangles: seismic station installed in the Cifalps project; red triangles Cifalps 2; green triangles Alparray.

A discussion on improved ray coverage

The impact of adding seismic stations is examined by comparing the tomographic reconstruction of a synthetic model. Two sets of reconstructions are compared: in one scenario data from existing seismic stations in the area was used and in the second scenario data from additional temporary stations were added. The histogram for each station is shown in Figure 3, which illustrates the enhanced ray coverage in scenario 2.

Figure 3. Histograms of ray counts. Each bar shows the number of rays per station. In the left panel only permanent seismic stations are shown, while in the right panel the “improved” distribution is shown. The bars are color-coded to indicate varying quantities: green signifies less than 50 rays, purple represents a range between 50 and 100, while red denotes more than 100.

The figure shows that the increase in seismic receiver density mainly affected the central area of the study region. This sector is that of the Po Plain, where the thick sedimentary layer of the terrain, intense human activity, and technical characteristics of seismometers have discouraged, and in some cases, prevented the installation of seismic stations until the end of the last century. The Cifalps projects successfully covered the area through strategic site selections, providing unprecedented data.

In seismic tomography experiments, the resolution power of the association between data and method, and therefore the method’s ability to provide a reliable image, is often evaluated using synthetic tests. Typically, a form of anomalies is imposed (for example, a checkerboard pattern alternating positive and negative anomalies with respect to the velocity of the layer) to choose the parameters for tomographic inversion. The parameter set that yields the image most similar to the synthetic one, often as a result of visual analysis, is then utilized for the carrying on the experiment. Generally, given a specific data distribution, the parameters governing node spacing, the number of iterations for the inversion process, and the mathematical parameters guiding the algorithm are adjusted according to the result of the tests. In our case, we will take the opposite approach, keeping all parameters constant but varying the amount of data (in this case, by adding seismic stations) to assess and potentially quantify the impact of adding seismic stations and, consequently, rays.

Figure 4. Model with synthetic anomalies later used to evaluate reconstruction capability of tomography. The model features two rectangles, one with a P-wave velocity value about 18% less than the layer velocity (red) and a larger one with velocity some 18% greater than the layer’s (yellow).

The experiment of reconstructing anomalies, which were imposed a priori, was carried out on the synthetic velocity distribution illustrated in Figure 4. Two rectangles were placed in the shallower layers (down to 5 km depth), one with P-wave velocity of 5 km/s (red vertical rectangle at the top of the figure, about 18% less than the layer velocity) and another larger one with a velocity of 7 km/s (yellow rectangle at the bottom of the figure, about 18% more than the average velocity). By employing the Simulps code (Thurber, 1983), a total of 1000 earthquakes were inverted on a custom dataset. Initially, only the permanent stations were utilized, yielding the results shown in Figure 5.

Figure 5. Reconstruction of the imposed anomalies using only data from permanent stations.

The tomographic inversion was carried out on a 36 by 36 by 13 node grid, with two iterations performed. Subsequently, the stations associated with the aforementioned projects were added, resulting in the image presented in Figure 6. By comparing Figures 5 and 6, several differences become apparent, which are exclusively due to the increased number of rays in the second case, as per the conducted test. The first significant difference is the reduced presence of smearing in the layers at depths ranging from 10 to 25 km in the tomographic image obtained with more rays (Figure 6). Smearing, as is well-known, poses a challenging problem in interpreting a tomographic image. The term smearing refers to a contamination by anomalies in a certain layer of the upper and lower layers. This effect can only be identified if the three-dimensional extent of a real anomaly is known, which is of course not. Thus it is impossible to distinguish the artifact from the actual distribution of anomalies and a great effort must be made to avoid smearing or at least minimizing it.

Figure 6. Reconstruction obtained using all available stations, including thus those installed in the temporary projects Cifalps and Cifalps2.

Our test proves that an increased number of stations does not avoid smearing, but it contributes to diminish it. This is evident, for example, in the layer at 15 km depth, where the greenish “stain” due to the leaking of the big rectangle turns to yellow color and becomes less regular when all stations are used.

Figure 7. The three upper layers of the model at 0, 3 and 5 km depth are here compared. In top panel, the imposed anomalies; in the middle panel, the reconstruction with permanent stations; in the bottom panel, anomalies reconstructed using all available data. Anomalies are in percentage with respect to layer velocities.

Figure 8. The three layers of the model at 10, 15 and 25 km depth are here compared. In top panel, the imposed anomalies; in the middle panel, the reconstruction with permanent stations; in the bottom panel, anomalies reconstructed using all available data. Anomalies are in percentage with respect to layer velocities.

The improvement introduced by the presence of more rays is even more evident when observing the velocities expressed as a percentage. The anomalies in this instance are shown as blue rectangles (positive values, greater velocity than the layer) and red rectangles (negative values). Figures 7 and 8 present a comparison between the a priori imposed anomaly (upper panel), the one reconstructed using only the data from permanent stations (middle panel), and the one obtained by adding the stations installed during the projects shown for different layers (bottom panel) for two sets of layers (0,3 and 5 km depth in Figure 7 and 10, 15, and 20 km depth in Figure 8). It is clearly evident that the blue rectangle anomalies in layers at 0 and 3 are much better reconstructed using all data, as the white spots in the rectangle of the middle panel confirm. The red anomaly of the layer at 5 km depth is better reconstructed as well. In Figure 8, the smearing is less pronounced; not only is the blue rectangle smaller in size in the bottom panel compared to the middle one, but its color is much lighter (5-7% in the bottom panel compared to more than 11% in the middle one). It should be noted that extreme values were used in this test. In fact, an anomaly of 18% is quite rare or at least less common in nature. This means that in a real case and for more moderate values, the tomography obtained with all stations would not show smearing, unlike that obtained with only temporary stations. This comparison confirms the enhancement achieved through the availability of more seismic rays and that a higher distribution of rays leads to a decrease in artifacts, sometimes even significantly. It should be noted that the presented results are obtained after two iterations and could therefore improve if the number of iterations were higher. This means that in a real case and for more moderate values, the tomography obtained with all stations would not likely show smearing, unlike that obtained with only temporary stations.

One of the main limitations of LET is that the maximum depth that can be investigated depends on the depth of the earthquakes used for inversion, which is the earthquakes occurring within the volume under study. The inclusion of additional seismic stations does not substantially impact this limitation, as evidenced by Figures 9 and 10, which display the theoretical path of hypocenter-station rays. It is important to note that this is a simplification, as the ray tracing process in inversion is not purely geometric (a straight line), but rather involves the utilization of specific algorithms that consider various variables and much more complicated ray paths. Nevertheless, the reconstructions in Figures 9 and 10 provide a rough estimation of the potential improvement in sampling the deeper layers. In particular, Figure 9 shows the comparison of ray coverage for P-waves in the case where only permanent stations are used (left panel) and where all available stations are used instead (on the right side). Figure 10 shows the same comparison but for S-waves.

Figure 9. Comparison of ray coverage for P rays with depth. Left: only permanent stations are used, right all stations.

Figure 10. Same as Figure 9 but for S rays.

Conclusions

While it is intuitive that the utilization of a greater number of seismic stations would enhance the level of detail in the tomographic image (Fichtner et al., 2024), this is the first instance in which this improvement has been quantified.

In this research, we have indeed taken advantage of the availability of seismic stations installed during projects for studying crustal and sub-crustal structure to observe the effect of a larger number of rays on seismic tomography LET. Comparing the reconstruction of a synthetic model obtained with the two sets of stations reveals that a greater number of rays allow for a more precise definition of the anomalies imposed a priori and reduces the smearing effect for deeper layers, albeit at the cost of increased computational burden. In fact, when incorporating more data, one must also consider aspects related to the expansion of data matrices and issues concerning the mathematical solution of the problem. These aspects may involve additional mathematical techniques, which are not taken into account in this paper (Virieux et al., 2024).

It is necessary to observe that the number of rays alone does not justify an improvement, on the contrary, as previously mentioned, it could even represent a burden in calculations. This happens because both the number of samples and the direction from which the rays that intersect the node originate are significant in any tomography. The more evenly distributed they are, the higher the resolution of that particular node will be. Hence, the improvement observed in this experiment is applicable solely to the data distribution resulting from the positioning of the stations used in this test, and it may vary considerably in other scenarios.

Finally, it is established that the addition of data does not enhance, or at least not significantly, the coverage of deep rays, which remains the main constraint of an LET tomographic experiment.

Aichroth, B., C. Prodehl and H. Thybo,1992. Crustal structure along the Central Segment of the EGT from seismic-refraction studies. Tectonophysics, 207 (1–2), 43-64. https://doi.org/10.1016/0040-1951(92)90471-H

Aki, K. and W. Lee, 1976. Determination of three-dimensional velocity anomalies under a seismic array using first P arrival times from local earthquakes: 1. A homogeneous initial model. Journal of Geophysical Research, Solid earth and planets, 81 (23), 4381-4399. https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/jb081i023p04381

Aki, K., A. Christoffersson and E. S. Husebye, 1977. Determination of the three-dimensional seismic structure of the lithosphere. Journal of Geophysical Research, Solid earth and planets, 82 (2), 277-296. https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/JB082i002p00277

Ansorge, J., D. Blundell and St. Mueller, 1992, Europe’s lithosphere – seismic structure. In Blundell D., Freeman R. and Mueller St. eds. A continent revealed. The European Geotraverse. Press Syndacate of the University of Cambridge. https://epic.awi.de/id/eprint/35095/2/the_european_geotraverse.pdf.

Curtis, A. and R. Snieder, 2002. Probing the Earth’s interior with seismic tomography. Chapter 52. In Lee, W., H. Kanamori, P. C. Jennings and C. Kisslinger, eds. International Handbook of Earthquake and Engineering Seismology, Part A. Academic Press 81A, Int. Geophys. Series, 861-874. https://www.geos.ed.ac.uk/~acurtis/assets/Papers/Curtis_Snieder2002.pdf

Eberhart-Phillips, D., 1986. Three-dimensional velocity structure in northern California Coast Ranges from inversion of local earthquake arrival times, Bulletin of the Seismological Society of America, 76 (4), 1025-1052. https://doi.org/10.1785/BSSA0760041025

Fichtner, A., B.L.N. Kennett, V. C. Tsai, C. H. Thurber, A. J. Rodgers, C. Tape, N. Rawlinson, R. D. Borcherdt, S. Lebedev, K. Priestley, C. Morency, E. Bozdağ, J. Tromp, J. Ritsema, B. Romanowicz, Q. Liu, E. Golos, F. Lin, 2024. Seismic Tomography 2024. Bulletin of the Seismological Society of America. https://doi.org/10.1785/0120230229

Gasparini, P. and TomoVes Working Group, 1998. Looking inside Mt. Vesuvius. Eos, Transactions, American Geophysica Union. 79 (19), 229-232. https://agupubs.onlinelibrary.wiley.com/doi/epdf/10.1029/98EO00165

Hetényi, G., I. Molinari, J. Clinton et al., 2018. The AlpArray Seismic Network: A Large-Scale European Experiment to Image the Alpine Orogen. Surveys in Geophysics. 39. https://link.springer.com/article/10.1007/s10712-018-9472-4.

Paul, A., M.G.Malusà, S. Solarino, S. Salimbeni, E. Eva, A. Nouibat, S. Pondrelli, C. Aubert, T. Dumont, S. Guillot, S. Schwartz and L. Zhao, 2022. Along-strike variations in the fossil subduction zone of the Western Alps revealed by the CIFALPS seismic experiments and their implications for exhumation of (ultra-) high-pressure rocks. Earth and Planetary Science Letters, 598. https://doi.org/10.1016/j.epsl.2022.117843.

Thurber, C., 1983. Earthquake Locations and Three-Dimensional Crustal Structure in the Coyote Lake Area, Central California. Journal of Geophysical Research, 88 (B10), 8226-8236. http://dx.doi.org/10.1029/JB088iB10p08226

Virieux, J., A. Paul, M. Langlais, G. Janex, P. Guéguen, A. Helmstetter and L. Stehly, 2024. Assessing the reliability of local earthquake tomography for crustal imaging: 30 yr of records in the Western Alps as a case study, Geophysical Journal International, 236 (1), 99–118. https://academic.oup.com/gji/article/236/1/99/7285801

Zhao, L., A. Paul, S. Guillot, S. Solarino, M.G. Malusà, T. Zheng, C. Aubert, S. Salimbeni, T. Dumont, S. Schwartz, R. Zhu and Q. Wang, 2015. First seismic evidence for continental subduction beneath the Western Alps. Geology 43, 815–818. https://doi.org/10.1130/G36833.1

About the Author

Stefano Solarino, PhD, is Senior Researcher at the Osservatorio Nazionale Terremoti (National Earthquake Observatory) at the INGV, Italy. Stefano’s expertise lies in the fields of seismotectonics and seismic tomography. He holds a position as a member of the Scientific Council for the Genoa Science Festival. Furthermore, he is the founder and head of the Working Group “Increasing the impact of scientific communication in seismology” within the European Seismological Commission. He has taken part in many field projects as the person in charge of data acquisition operations. Stefano is engaged in scientific outreach for numerous years in schools and public settings. He has coordinated and joined educational events. He has written, acted, and sung in a performance about Earth’s history, merging his singing passion with outreach endeavors.