PARTOS-Passive and Active Ray TOmography Software : description and preliminary analysis using TOMO-ETNA experiment ’ s dataset

In this manuscript we present the new friendly seismic tomography software based on joint inversion of active and passive seismic sources called PARTOS (Passive Active Ray TOmography Software). This code has been developed on the base of two well-known widely used tomographic algorithms (LOTOS and ATOM-3D), providing a robust set of algorithms. The dataset used to set and test the program has been provided by TOMO-ETNA experiment. TOMO-ETNA database is a large, highquality dataset that includes active and passive seismic sources recorded during a period of 4 months in 2014. We performed a series of synthetic tests in order to estimate the resolution and robustness of the solutions. Real data inversion has been carried out using 3 different subsets: (i) active data; (ii) passive data; and (iii) joint dataset. Active database is composed by a total of 16,950 air-gun shots during 1 month and passive database includes 452 local and regional earthquakes recorded during 4 months. This large dataset provides a high ray density within the study region. The combination of active and passive seismic data, together with the high quality of the database, permits to obtain a new tomographic approach of the region under study never done before. An additional user-guide of PARTOS software is provided in order to facilitate the implementation for new users.


Introduction
TOMO-ETNA experiment aimed to enlighten the complex internal structure of Mt.Etna and surrounding areas by developing dense active and passive seismic surveys.Previous chapters of this special volume [e.g.Coltelli et al. 2016, in this volume;Ibáñez et al. 2016aIbáñez et al. , 2016b, in this volume] , in this volume] describe the motivation, performance, data acquired and some different approaches of TOMO-ETNA experiment.The most straightforward procedure to study the internal structure of Mt.Etna volcano using seismic signals is to perform a joint active and passive seismic tomography inversion to derive the 3D velocity structure.As described by Ibáñez et al. [2016b, in this volume] the main success of the experiment was to obtain a large high quality seismic dataset composed by both active and passive seismic sources.
Seismic tomography methods appeared in the early 1970s and from their first application these techniques have been widely used in different scales and environments.Their theoretical and mathematical bases are well known by the scientific community, but still new algorithms and approached are being developed.Implementation procedures on a tomographic code may vary depending on the scope of the study (global, regional or local) and the data available [Koulakov 2012;Koulakov and Shapiro 2015].At present, several tomographic algorithms can be easily found among the international community.Some of these codes were developed for a specific region and later adapted to other regions.For example the code developed by Toomey et al. [1994] was initially applied in the Pacific Rift and later adapted to other regions such as Deception Island in Antarctica [Zandomeneghi et al. 2009].On the other hand, many codes have been developed to be used in a wide range of areas.Some of the most relevant and well known codes are: SIMULPS [Thurber 1983;Eberhart-Phillips 1993] applied in Vesuvius [Lomax et al. 2001] and Mt.Etna [Patanè et al. 2006] among others; FAST [Zelt and Barton 1998] applied in India by Rao et al. [2015]; TO-MODD (Zhang and Thurber [2003]; applied for example by Kato et al. [2010] in Wakayama, Japan; or Alparone et al. [2012] in Mt.Etna, Italy); or LOTOS [Koulakov 2009a].Each of these codes is different in some procedures such as ray tracing, grid parameterization or matrix inversion among others.When only passive seismic signals (earthquakes) are used, its application is mainly limited to regions with high seismicity.
LOTOS software [Koulakov 2009a] has become one of the most widely used codes for tomographic inversion of passive data, reflected by the numerous published works detailed in Table 1.LOTOS code is a friendly widely used tomographic code that permits study from small regions using local events to regional inversions using teleseismic events [Jakovlev et al. 2013].The robustness of the algorithms implemented and the capacity to handle different scales studies (regional or local), together with its friendly use makes of this code one of the most popular and used software among the scientific community.Robustness of the code has been highly proved by the numerous of published works ap- plying this software in a large variety of regions (Table 1).These studies vary from active volcanic areas, dormant volcanoes and non-volcanic areas, and demonstrate the versatility of the code.Passive seismic data comes from local, regional and teleseismic earthquakes.These data are generated within a wide range of depths and therefore give usually a better ray coverage in depth of the studied region.However, number of earthquakes depends on the region.Earthquake location coordinates and origin time are unknown and must be calculated by using the available information: arrival of the seismic waves to the different seismic stations.Local and regional seismicity usually produces clear P-and S-waves that are, then, inverted for obtaining plausible velocity models and source parameters.
The limitation of lack of seismicity is partially solved with the use of active sources, either at sea using airguns or on-land using explosions [Koulakov 2012;Koulakov and Shapiro 2015].At present, active source tomography has been implemented in few codes such as FAST and ATOM-3D [Koulakov 2009b].
ATOM-3D code follows the same structure as LOTOS, but for active datasets.In this case, no relocation of the events is performed, as location and origin time of the active source is a priori known, and it is taken as a fixed parameter.This algorithm has been successfully tested in some active seismic studies such as Tenerife [García-Yeguas et al. 2012] and eastern Greece [Shahrukh et al. 2012].Active data comes from humanmade seismic sources and there is a wide range of methods: chemical reactions, drop of a mass, airgun shots, explosions, etc.The main advantage of the use of active sources is that they permit to register seismic waves in whatever region.Nonetheless, it is very expensive and therefore not always available [Koulakov and Shapiro 2015].Commonly, active sources give very detailed information of the shallower layers of the crust.However, resolution in depth is poor due to the fact that they can only be produced near the surface.The main characteristic of the active sources that make them useful for seismic tomography is that we certainly know the exact time origin and location of both source and receiver, and therefore we only have the ray tracing variable.Clear Pwave arrival can be registered from active sources while S-wave arrivals are not always available.
Recently, some efforts have been focused on the use of both seismic sources to obtain a joint seismic tomographic image.For example, a few algorithms merge both active and passive data [e.g.Rawlinson and Urvoy 2006;Wagner et al. 2007;Battaglia et al. 2008].The major complexity of these new joint procedures is that active and passive seismic data have different characteristics such as its source location, and therefore it is important to take into account these differences in order to better adapt the algorithms for each case.The combination of these two types of sources permits to diminish their lack of information and enhances their valuable information.Active sources can cover areas where no natural seismicity is located.On the other hand, earthquakes may give important information of the deep layers of the crust and of the S-wave properties.By merging active and passive data we obtain better ray coverage in both shallow and deep layers, leading to a higher resolution of seismic tomography models.
In this manuscript, we describe PARTOS, a new tomographic code developed for joint inversion of active and passive datasets.This code was born to solve the necessities above mentioned.

Description of the code
PARTOS code has been developed for simultaneous tomographic inversion of active and passive seismic data.Workflow of the algorithm and general data structures are shown in Figure 1.

Input data
The input data for calculations include the data groups related to active and passive seismic experiments, as well as the starting 1D velocity model and other parameters for calculations.The passive source data include two files: one with geographical coordinates of stations and another with arrival times of P-and S-waves from earthquakes to the stations.The information of the source coordinates and origin times is not strictly required.In case of its absence, the sources can be placed at the center of the network or beneath the station with minimum arrival time.If the locations of events are available a priori, they can be included to the initial data file that may facilitate the calculations.See Appendix 1.For additional information and examples.
For the active source data, there is only one file.Each line in this file corresponds to one ray containing geographical coordinates of the source and the receiver and travel time of seismic wave.The first three lines should correspond to the gathers for which the data were picked.In the cases of offshore surveys with the use of airguns, the gathers are usually the stations, whereas in the cases of using a few explosions and numerous receivers, the gathers are the sources.
The initial data files may contain the two-dimensional map of topography (in case of absence, the flat topography is predefined at zero depth).It is especially important in cases of offshore experiments, when the travel time in the water layer between airguns and sea bottom should be taken into account.
The file with the 1D model of the P-and S-wave velocities should be defined.The velocity values are set at some depth levels; between these levels, the velocity is linearly interpolated.The values of the initial S-velocities are often determined from the P-velocities based on constant Vp/Vs ratio.
In addition, a file with a set of different parameters used for calculations should be defined.This file includes, among others, the grid spacing, the parameters for ray tracing and source locations, weights for the active and passive data, damping parameters for inversion, and visualization settings.Defining these parameters is mostly based on several trials of performing the inversions for real and synthetic datasets.

Source locations and ray tracing
For the passive data, the location of sources is conducted in a same way as in the LOTOS code [Koulakov 2009a].Three main steps conduct the source locations: (i) absolute preliminary locations; (ii) relocation using bending method procedures and; (iii) calculation of source parameters applying station corrections and P and S anomalies.For the preliminary locations, there are two options providing simplified algorithms for computing travel times of seismic rays.The first option presumes computing a table of reference times in the starting 1D model for all combinations of source depths and epicentral distances.In this step, analytical formulas for ray tracing in the 1D model from Nolet [1981] are used.Then the travel times of rays are calculated by linear interpolations of values from this table.For some areas where the depths of events are compatible with the lateral dimensions of the study region, another option with computing travel times along straight ray paths appears to be more efficient.The use of these approximations makes possible conducting calculations of source coordinates for voluminous amounts of data using the grid search method [see Koulakov and Sobolev 2006a].
The second step localizes the sources in the 3D P and S velocity models using the gradient descending method.Unreliable locations are discarded by identifying the outlier time residuals [Alinaghi et al. 2007].In this case, the calculation of the travel times is performed using the bending method of the ray tracing based on the principle of time minimization proposed by Um and Thurber [1987].The localization procedure is performed iteratively.At each iteration, an updated 3D velocity model is used for the localization.
The same algorithm of ray tracing is used for the active-source data.In the case of the offshore shooting using airguns on the sea surface, the ray is computed for the entire paths including the parts in the water layer, where velocity is presumed constant and equal to 1.44 km/s.The calculations are performed using the regular 3D bending ray tracer between the source on the sea surface and the station.

Grid construction
To parameterize the 3D velocity distribution, we construct the mesh with nodes distributed in the study volume according to the ray density.Between the nodes, the velocity anomalies are continuously approximated using a linear interpolation.Because of significant disproportion in the amount of passive and active data, the ray density is computed using different weights for the P and S data.The nodes are placed only in areas with sufficient ray coverage (typically, with the ray density of 0.1 of the average value).In the map view, the nodes are installed regularly with a fixed grid spacing.In the vertical direction, the grid spacing depends on the ray density; this is, the higher the ray density is the finer the grid PARTOS-TOMOGRAPHY SOFTWARE spacing can be.In well covered areas, the spacing is equal to the minimum predefined value; and it becomes larger where the ray density is lower.
Usually, the grid spacing is set smaller than the expected resolution so that every restored anomaly was based on several nodes.This helps to avoid the dependency of the results on the grid configurations.To further reduce this effect, we perform the inversion for several grids with different basic orientations (typically, four grids with orientations of 0, 22, 45 and 67 degrees).After computing results for all grids, they are averaged in one 3D model using a regular mesh.The grid is constructed only one time in the first iteration.In the following iterations, the velocity anomalies are updated in the same nodes.

Matrix calculation and inversion
To perform the tomography inversion, we compute one matrix which includes the pairs of the active and passive data.The matrix elements responsible for the velocity distributions are the first derivatives A ij representing time deviation along the j-th ray due to the unit velocity variations in the i-th node (Equation 1): where Dg j represents the slowness perturbation at the point of the ray caused by slowness anomaly Dv j , at the j-th node [Koulakov and Sobolev 2006b].Their values are computed by numerical integration along the ray paths derived after the steps of locations (for the passive data) and ray tracing (for the active data).For the passive data, each row of the matrix contains also four elements responsible for source parameters (Px, Py, Pz and 1) and one unit element responsible for the station correction (P or S).For the active source data, there is only one additional element responsible for the gather correction (usually, stations, for the airgun surveys and sources, in cases of explosions on land).The data vec- . Checkerboard test calculated for the active dataset (upper frame) and the joint dataset (lower frame) using a 12-km synthetic anomaly. (1) tor contains the time residuals computed for the current velocity model.The inversion is damped by using two additional matrix blocks.The first block, which is responsible for amplitude damping, is a diagonal matrix contains only one non-zero element in each row and zero data vectors.Increasing weight of this block decreases the amplitude of the solution.The second damping block controls the flattening of the solution.Each row of this block contains two nonzero elements with equal values, but opposite signs, corresponding to all combinations of neighboring nodes in the parameterization mesh.Increasing the weight of this block decreases velocity gradients in the resulting models.
The inversion of the large sparse matrix is performed using the LSQR algorithm [Nolet 1981;Paige and Saunders 1982].The most important parameters controlling the inversion with the values optimized for a spe-cific case of the TOMO-ETNA experiment are presented in Table 1.The values of these parameters are determined after several trials with synthetic and observed data.

Combining the velocity models and iterating
The models computed for several differently oriented parameterization meshes are averaged and combined in one model defined on a regular grid.This model is used in the next iteration for performing the ray tracing (for active data) and source re-location (for the passive data).Number of iterations is usually fixed at a predefined value as a compromise between the calculation time and accuracy of the solution.To establish the number of iterations, we use a threshold parameter named dtot that represents the time residuals [Koulakov 2009a].This threshold is usually set at 0.5 and we consider the solution as stable once the dtot value decreases slowly down the threshold (usually after 5-6 iterations).In this case, tuning of the solution properties is conducted using the damping parameters.A regular inversion running in a 'core i7' personal computer takes approximately 2 hours to perform 5 iterations.

Synthetic modeling
Synthetic modeling is an important element of the tomography workflow that allows not only assessing the resolution, but also estimating the optimal values of the parameters for calculations.The synthetic modeling should as adequately as possible represents the realistic procedure of the tomography inversion.The synthetic model is defined as a sum of the 1D reference model (might be different of the starting model used for reconstruction) and synthetic anomalies.The code gives several possibilities of defining seismic anomalies with the use of checkerboard or free-shaped polygonal prisms in horizontal or vertical sections.The calculation of the synthetic data is conducted for the same ray paths as in the case of observed data.The ray tracing is performed using the 3D ray tracer based on the bending algorithm.The synthetic travel times are perturbed with random synthetic noise having similar statistical distribution as the residuals after inversion in the case of observed data.The amplitude of the noise is defined to achieve similar variance reduction during the inversion as during inversion of the experimental data.After computing the synthetic travel times, for the passive data we randomly shift the coordinates and origin times.
The reconstruction of the synthetic model contains the same calculation steps as in the case of ob-served data analysis including the stage of initial source location for the passive data.During the synthetic inversion, we try different values of free parameters to achieve the best correspondence of the reconstructed anomalies with the original model.These parameters are then used to invert the experimental data.

TOMO-ETNA dataset
In order to check PARTOS code, we used the database recorded during TOMO-ETNA experiment [Ibáñez et al. 2016b, in this volume].This dataset consists of both active and passive seismic data recorded during periods of 3 weeks (active data) and 4 months (passive data).The performed calculations are based on the dataset with 184,797 active source rays and 4595 rays from 452 earthquakes [Barberi et al. 2016, in this volume].The active dataset used are those recorded during the refraction experiment of TOMO-ETNA (for more information, see Ibáñez et al. [2016a] in this volume).The area analyzed comprises the whole TOMO-ETNA experiment region: from Aeolian Islands to south Mt.Etna volcano (Figure 2).The P-phase arrivals from both active and passive data used for this seismic tomography were detected and picked using an adapted version an automatic picking algorithm named AMPA [Álvarez et al. 2013;García et al. 2016, in this volume].Only P-phases have been detected and inverted in this study.

Synthetic tests
First, we present the results of synthetic modeling for three data subsets including first only active data, then only passive data and finally the joint dataset.We have performed several checkerboard tests with different configurations of anomalies.To explore the horizontal resolution, we use the model with periodic anomalies in the map view that do not change the shape with depth.Here we present the results for two models with the horizontal spacing of 25 and 12 km.For the larger anomaly size (Figure 3), the shallow structures are similarly well reconstructed for both cases of using only active and combined datasets.Even outside the station area, some anomalies can be reconstructed.For the deeper section of the same model with 25 km anomalies, we can see a considerable difference between the cases with active and joint inversions.Notice the improvement of the reconstruction quality in the central part of the study area where most seismicity occurs.In the second test with 12 km anomalies (Figure 4), the area of the good resolution in the shallow section is smaller compared to the previous case with larger anomalies.In this case, the anomalies can only be resolved beneath the station area.For the deeper section, the effect of adding the passive data ap-   pears to be more prominent: without passive data, no structures can be resolved at 8 km depth.
The second type of the test aims at studying the vertical resolution along one vertical section.The velocity anomalies are defined in the vertical section as alternated positive and negative patterns.In the direction across the section, the length of anomalies is 12 km.The amplitude of anomalies is 15% (Figure 5).Since the rays from active sources are mostly oriented horizontally, they cannot enable good vertical resolution.Passive data also cannot provide good vertical resolution because of the trade-off between velocity and source parameters.Together these two types of data enable considerable improvement in the resolution quality.
In addition to testing the horizontal and vertical resolution, these synthetic analyses were used to estimate the optimal values of the free parameters used for calculations that were then used for calculations of the main models based of experimental data presented in the next section.

Inversion of the experimental data
As in the case of synthetic modeling, the inversions of experimental data were performed using three subsets: with only active sources, only passive sources and the joint dataset.The best resolved area is design by the parameterization grid and the ray density (Figure 6).with this preliminary inversion are the contrast between the volcanic system, with low velocity anomalies up to 40%, and the metamorphic chain located in the northern part of Sicily.Additionally, and even using a large grid, a high velocity anomaly has been detected beneath the central craters of Mt.Etna that corresponds with the 'High Velocity Body' described by many authors in previous tomographic studies on this region [Chiarabba et al. 2004;Patanè et al. 2006;Alparone et al. 2012].It can be seen that the results based on active data have high resolution in a larger area compared to that in the case of using passive data.However, for the deeper section, the passive source data provide better resolution.Correspondingly, the joint dataset benefits the advantages of both datasets and provides better resolution that in each of the individual cases.

Final remarks
PARTOS is a friendly tomography code born from the necessity of combining active and passive seismic data during a tomographic inversion.The robustness of the new code used is widely proved as most of the algorithms applied were adapted from previous versions of LOTOS and ATOM-3D codes, which are widely used among the scientific community during the last decade.New features have been implemented to enhance the accuracy and completeness of the calculations, making the software fast and adaptive to different kind of datasets.
We have tested the code with the data collected during TOMO-ETNA experiment by performing different synthetic tests and real inversions.Results remark that merging two different seismic sources dramatically enhances the resolution and quality of the tomographic results.Indeed, we complement the high resolution in the shallowest layers provided by active source tomography with the deep structure enlightening of the passive tomography.
Acknowledgements.This paper has been partially funded by the following research projects: the European project MED-SUV funded by the European Union's Seventh Framework Program for research, technological development and demonstration under grant agreement No. 308665; the Spanish COCSABO project (COC-DI-2011-08); the European project EUROFLEETS2 (Seventh Framework Programme, grant agreement No. 312762) through transnational access to the research vessels "Sarmiento de Gamboa" operated by CSIC (Spain) and "Aegaeo" by HCMR (Greece); the Geophysical Instrument Pool Potsdam (GIPP) from GFZ (Potsdam) with the project (Seismic TOMOgraphy of ETNA volcano and Eolian Islands, Italy, using active and passive seismic data).We would like to thank the following supporting institutions: Dipartimento Regionale della Protezione Civile, Regione Siciliana; Dipartimento Azienda Regionale Foreste Demaniali, Ufficio Provinciale di Catania; Ente Parco dell'Etna; Unidad de Tecnología Marina -CSIC in Barcelona (Spain); Stato Maggiore Marina (Italian Navy General Staff ), CINC-NAV (Command in Chief of the Fleet) and Marisicilia (Navy Command of Sicily); Coastal Guard of Messina and Riposto; to obtain support and navigation permissions for the oceanographic cruises: Spanish Foreign Office and Italian Foreign Office.This paper has been partially supported by the Spanish projects TEC2015-68752-R (MINECO/FEDER), KNOWAVES and CGL2015-67130-C2-2.We would like to thank all private and public owners of the sites selected to deploy seismic stations for their kind and unselfish disposal to use their properties.This manuscript has been largely improved by the insightful comments of Dr. R. Abreu and Dr. Stich and by the editor José Morales.15.05900 37.39500 -0.00500 15.30630 37.28030 0.01000 5.92000 15.05900 37.39500 -0.00500 15.30990 37.27830 0.01000 5.92000 15.05900 37.39500 -0.00500 15.32050 37.27230 0.01000 6.37000 15.05900 37.39500 -0.00500 15.32230 37.27130 0.01000 5.94000 15.05900 37.39500 -0.00500 15.33990 37.26150 0.01000 6.77000 15.05900 37.39500 -0.00500 15.05200 37.42400 -0.00600 14.94300 37.45600 -0.02600 14.80000 37.43100 -0.04400 14.68700 37. 39900 -0.31700 15.5630 38.4520 131.82

Figure 1 .
Figure 1.Schematic PARTOS workflow for active and passive data.Rectangles represent the data files and rounded blocks depicts the steps of the program.Steps with dotted lines are performed only once in the first iteration.

Figure 2 .
Figure 2. (a).Map of the passive dataset with the stations (blue triangles) that registered the events (red dots).(b) Map of the active source distribution including the seismic network used to register the shots (red dots).

Figure 3 .
Figure 3. Checkerboard test calculated for the active dataset (upper frame) and the joint dataset (lower frame) using a 25-km synthetic anomaly.

Figure 5 .
Figure 5. Vertical Checkerboard test calculated for the active dataset (upper frame), passive dataset (middle frame), and the joint dataset (lower frame) using a 12-km synthetic anomaly.Red triangle shows the location of Mt.Etna.

Figure 6 .
Figure 6.Example of the ray density and node distribution calculation for the TOMO-ETNA dataset.

Figure 7 .
Figure 7. Horizontal layers (2 km b.s.l. in the left column and 8 km b.s.l. in the right column) of the real data inversion of the TOMO-ETNA dataset using a 6-km-spacing grid: active dataset (upper frame); passive dataset (lower frame).Black triangles show the location of Mt.Etna.

Figure 8 .
Figure 8. Real data inversion of the TOMO-ETNA joint dataset using a 6-km-spacing grid.Black triangles show the location of Mt.Etna.
Figures 7 and 8 present the results of inversions for all these cases in two horizontal sections.Figure 9 presents the results of the three different inversions (active, passive and joint datasets) for one vertical profile.It is important that in the central part of the study area, both active and passive datasets, which are absolutely independent, provide consistent structures which are also similar to the results obtained by other authors [e.g., Patanè et al. 2006; Alparone et al. 2012].The main features that are enlighten

Figure 9 .
Figure 9. Vertical section of the real data inversion of the TOMO-ETNA dataset using a 6-km-spacing grid: active dataset (upper frame); passive dataset (middle frame) and joint dataset (lower frame).Red triangle shows the location of Mt.Etna.

Figure A. 1 .
Figure A.1.Root directory structure of PARTOS software.Orange boxes represent folders and white boxes represent files.

Figure A. 2 .
Figure A.2. Directory scheme of the DATA folder.

Table 1 .
Recent most relevant tomographic studies performed using LOTOS code.