Subduction or delamination beneath the Apennines ? Evidence from regional tomography

In this study we present a new regional tomography model of the upper mantle beneath Italy and the surrounding area derived from the inversion of travel times of P and S waves from the updated International Seismological Centre (ISC) catalogue. Beneath Italy, we identify a highvelocity anomaly which has the appearance of a long, narrow “sausage” with a steeply dipping part down to a depth of 400 km and then expanding horizontally over approximately 400 km. Rather than to interpret it as a remnant of the former Tethyan oceanic slab, we consider that it is made up of the infra continental lithospheric mantle of Adria, which is progressively delaminated, whereas its overlying crust becomes progressively accreted into the Apenninic tectonic wedge.


Introduction
The Mediterranean region is located in the convergence zone between the African and European plates which is characterized by a very complex interaction of different tectonic regimes including subduction, collision, spreading and shear zones (e.g., Faccenna et al., 2004Faccenna et al., , 2014)).In some parts of the region, such as in the Italian Peninsula and its surroundings, all of these regimes are confined to a limited area that results in very complex geological and tectonic structures (Fig. 1).
Today for instance, even the nature of the lithosphere and crust beneath the deep Ionian Basin remains debated, the debate being whether it is made up of a remnant of the for-mer Mesozoic Tethys ocean or by the distal portions of the stretched African continental crust (Dercourt et al., 1985(Dercourt et al., , 1993;;Stampfli and Borel, 2004;Roure et al., 2012).
When going back in the past, the paleogeographic reconstructions of the various platform and basinal domains of Adria, with the Albanides, and Hellenides in the east, and the Apennines and Sicily in the west, are also increasingly ambiguous.Crucial in this respect is whether the Mesozoic Ionian allochtonous units of the Albanides and Hellenides, and coeval Mesozoic basinal series from the Lagonegro and Imerese units of the southern Apennines and Sicily, are considered as lateral equivalents of the modern, still deep water Ionian basin (Roure et al., 2004(Roure et al., , 2012)).
Seismic tomography is one of the key tools which are used to resolve the mechanisms of deep tectonic processes.However, despite the large number of different tomography models, knowledge on the major tectonic processes is often ambiguous and contradictory.In this paper we present a new 3-D seismic model of the upper mantle beneath the Italian region constructed based on similar calculation schemes as those used in Koulakov et al. (2009), but using a considerably larger data set.Based on the distributions of P and S velocity anomalies, we propose a new interpretation for the recent tectonic history of the Italian region.
For a long time, the European part of the Mediterranean has been an attractive region for seismic tomography studies thanks to relatively dense distribution of seismic stations, intensive seismicity and fairly heterogeneous deep structure controlled by complex geodynamic processes.In the early 1990s, the regional mantle structure beneath Europe was studied by the use of travel times of body waves (Spakman, 1990;Spakman et al., 1993;Spakman and Wortel, 2004) and surface waves (Zielhuis and Nolet, 1994).The global models by Bijwaard et al. (1998) and Bijwaard and Spakman (2000) have provided the compatible resolution for the European region with the existing regional models.Later, Piromallo and Morelli (2003) presented another P velocity regional model for the entire European region based on body waves.More recently, Koulakov et al. (2009) have constructed P and S tomography models for the upper mantle beneath Europe, which took into account a newly obtained crustal model by Tesauro et al. (2008).In parallel there were several models of Europe created based on surface-wave data (Boschi et al., 2004(Boschi et al., , 2009;;Marone et al., 2004;Legendre et al., 2012).Adjoint tomography has been used for studying the European region by Zhu et al. (2012).Most of these models are generally consistent with each other, especially for the uppermost mantle.They identify highly contrasting features which detect general lithospheric structures known from geology.However, there are some differences between results obtained based on body and surface waves.Surface-wave tomography usually provides stronger amplitudes of anomalies and smoother lateral and sharper vertical variations.
The Italian Peninsula and surrounding regions have been studied using the regional data of the Italian permanent and temporary networks in many different studies.P wave veloc-ity models of the crust and uppermost mantle based on travel times of body waves from regional events were the focus of many investigations, mostly performed by Italian scientists (e.g.Amato et al., 1993;Alessandrini et al., 1995;Selvaggi and Chiarabba, 1995;Di Stefano et al., 1999;Cimini and De Gori, 2001;Orecchio et al., 2011;Gualtieri et al., 2014).
Both regional and local scale tomography models provide valuable information on the deep processes which explains many observations on surface tectonics.Together with various geological and geophysical data, the tomography results are used for geodynamic interpretation and to reconstruct the scenarios of plate interactions in the European regions.One of the key regions is the Calabrian-Tyrrhenian system where a complex shape of the plate boundary passing through Calabria, the Apennines, the Alps and the Dinarides is identified (Fig. 1).The existence of deep seismicity down to ∼ 400 km beneath Calabria and the Tyrrhenian Sea and of the active volcanism of the Eolian Arc suggests that the subduction processes are active here (e.g., Isacks and Molnar, 1971;Selvaggi and Chiarabba, 1995).It was first proposed by Malinverno and Ryan (1986) that the loop-shaped boundary, leading to back-arc extension in the Tyrrhenian sea (Spadini et al., 1995), was formed due to strongly curved subduction occurring in the narrow zone in front of the Apennines.However, such a strong bending of the subducting plate appears to be mechanically not plausible, as was shown by analogue experiments by Faccenna et al. (1996).Alternatively, Faccenna et al. (2004) have proposed that the complex boundary shape in the Apennines region was formed due to the initial evolution of the western Mediterranean subduction zone (WMSZ) followed by trench retreat and back-arc extension.According to their model, the counterclockwise rotation of the Apennines block started at ∼ 35 Ma due to the opening of the Ligurian-Provençal basin, as also proposed by Mattei et al. (2002), and then continued due to the spreading of the Tyrrhenian basin.However, the final stage, which resulted in the origin of the narrow subduction zone remains not completely clear.Spakman and Wortel (2004) have drawn attention to the segmented nature of the descending slabs in the western Mediterranean region, whereas Govers and Wortel (2005) have pointed out that STEP faults are playing an important role in the dynamics of the underlying subduction dynamics.In this paper we present an alternative scenario based on new tomography models of P and S velocities.We also provide some additional testing results to show that derived seismic structures related to the Apennines region are robust.

Data and tomography algorithm
To construct the upper mantle velocity structure beneath the central Mediterranean, we use P and S travel limes from the International Seismological Centre (ISC) catalogue in a time period from 1964 to 2012.The calculations were performed using generally same computing schemes as in Koulakov et al. (2009) but with the use of a large amount of data accumu-  lated in the ISC catalogue since 2007 (the last year of data in the previous study).Here we use more than 11 million travel times (∼ 8.9 × 10 6 of P and 2.3 × 10 6 of S picks) from more than 150 000 events (72 picks per event).In the previous study by Koulakov et al. (2009), the number of data was several times smaller.
Prior to using the ISC data for tomography, they were reprocessed using the location tools and outlier analysis developed by Koulakov and Sobolev (2006), which resulted in rejection of almost 30 % of the data.All data corresponding to seismic rays traveling inside the selected window, at least partly, are selected for the inversion.They include data from ∼ 66 000 earthquakes located in the study region (red dots in Fig. 2a) recorded by the worldwide stations (blue triangles in Fig. 2b), as well as data from ∼ 90 000 teleseismic events (gray dots in Fig. 2b) recorded by 1915 seismic sta-tions located inside the study area (black triangles in Fig. 2a).The travel times of all seismic rays in the European region were corrected for the crustal model developed by Tesauro et al. (2008).
The inversion is performed in three circular windows covering the study area having a diameter of 1500 km (Fig. 2a), which appears to be the most appropriate size to study the upper mantle, as estimated by Koulakov and Sobolev (2006).In this study we use the windows, which cover a much larger area than the region of interest of this study, in order to avoid boundary effects that may appear at the margin of a circle due to smearing of outside anomalies.Such a window-bywindow approach allows more optimal definition of inversion parameters depending on the data amount in each window, and it provides higher resolution than a global inversion for the entire region.After performing the inversions, the results in all windows are combined in a single model.
The algorithm of tomographic inversion was developed by Koulakov et al. (2002) and then significantly modified by Koulakov and Sobolev (2006).The parameterization of the velocity models is performed using the nodes installed according to the ray density on horizontal levels at depths of 10,25,50,100,150,220,290,360,430,500,570,640,710,780,850,930, and 1000 km.The minimum grid spacing in horizontal levels is set at 30 km.No nodes are installed in areas where there is no data.To avoid an effect of the grid orientation upon the results, the inversions were performed in four grids with different basic orientations (0, 22, 45 and 67 • ) and then combined in a single model.When merging the results computed for different circular windows and differently oriented grids, we calculate a 3-D weight function.It depends on the distance from the nearest parameterization node: if the distance is less than 40 km, the weight is equal to 1; at distances from 40 to 80 km it decreases linearly from 1 to 0. This weight is also scaled depending on the distance to the border of the circular area (from 85 to 100 % of the radius, the scaling factor linearly decreases from 1 to 0).
The inversion was performed simultaneously for the P and S velocity anomalies, source parameters (dx, dy, dz and dt) and station corrections.Although, P and S models are theoretically coupled in inversion through source parameters, this coupling is very weak, and they can be considered independent.There is no other constraint linking the P and S velocity anomalies.The matrix inversion was performed using the LSQR algorithm designed by Paige and Saunders (1982) and Nolet(1987).The quality of the solution was controlled by roughness regularization which was conducted by minimizing the differences of velocity variations between all pairs of neighboring nodes.
The same tomography code was used for studying various regions of the world including several subduction zones, such as the Kuril-Kamchatka and Aleutian (Koulakov et al., 2011), Mariana and Izu-Bonin arcs (Jaxybulatov et al., 2013), Taiwan (Koulakov et al., 2014).In all subduction zones, both P and S velocity models revealed consistent im-ages of the slabs coinciding with the locations predicted by other studies and marked with deep seismicity.

Tomography results
The computed P and S velocity anomalies, which are considered to be the main results of our calculations, are presented in horizontal and vertical sections in Figs. 3 and 4. The anomalies are given in percent in respect to the 1-D reference model AK 135 (Kennett et al., 1995).The model was computed down to 1000 km, while here we show only upper mantle structures down to 700 km.
One of the key inversion parameters which strongly affects the solution is the value of damping which controls the smoothing of the resulting anomalies.Figure 5 presents an example of solution with weaker smoothing parameters equal to 10 and 16 for the P and S models, respectively, as compared to 15 and 25 used for the main model shown in Figs. 3 and 4. It can be seen that this solution adds some smaller patterns and increase the amplitudes of anomalies, especially in the lower part of the model.Although this model appears to be stable, and the size of anomalies remains larger than those resolved in synthetic tests (see next paragraphs), we prefer staying at a conservative side and present a smoothed model as the main result of this study.
This model appears to be generally consistent with the results by Koulakov et al. (2009) and other tomography models.For example, vertical sections in Figure 2.8a in Spakman and Wortel (2004) looks very similar to our result for the P velocity in vertical section A1-B1 having approximately same location.
Before discussing the observed seismic anomalies, we present a series of synthetic tests which were specially performed to assess the robustness of the obtained model in the studied region, which appears to be very important, regarding some inconsistencies with previously published models.First of all, we present the results of the traditional checkerboard test which consists of a reconstruction of periodic rectangular anomalies with the amplitude of ± 4 % based on the existing data configurations.In the example presented in Fig. 6, we present results of two checkerboard models with different sizes of P and S patterns.In the upper two rows, the horizontal sizes of P and S anomalies are 1 • by 1 • and 2 • by 2 • , respectively.In the second case (lower two rows in Fig. 6), the sizes of anomalies are 2 • by 2 • and 4 • by 4 • , for P and S models, respectively.In all cases the sign of anomalies changed with depth every 200 km (200, 400, 600 km etc).The checkerboard anomalies were set in the entire Earth; thus the ray paths corresponding to remote stations or sources were affected by the outside anomalies.The data were perturbed with random noise of 0.3 and 0.6 s of mean magnitude.The results of the checkerboard reconstructions for the P and S models are presented in Figure 6 in horizontal sections at 100, 300 km and 500 km (middle depths of the checkerboard levels).For the P anomalies with 1 • size pat-terns, the correct reconstruction occurs in selected areas with the maximum amount of data.In deeper sections, the resolution is much poorer compared to the result at 100 km depth.For the 2 • size anomalies, the P model is reconstructed almost perfectly in the onshore areas and in the Adriatic Sea.In this test, the P anomalies almost do not lose their amplitude with depth.In the offshore areas, the anomalies are strongly smeared, which is related to lack of available data.The resolution of the S model appears to be poorer mostly because of a smaller amount of data.In addition, most of the S rays correspond to short rays traveling in the uppermost layers; for the deeper parts of the model, the S ray coverage is even poorer.We observe that in the case of 2 • size of patterns, the anomalies below 200 km depth are smeared and lose their amplitudes.Nevertheless, their locations in space are correct, which is important for qualitative interpretation of the results.
The second test shown in Fig. 7 is aimed at studying the capacity of the algorithm to resolve a sausage-shaped anomaly similar to the one we obtain in the observed data inversion.We define the synthetic anomalies along the Profile A1-B1, the same as is used for presenting the main model in Fig. 4, as free-shaped polygonal horizontal prisms with a thickness of 200 km in the direction across the profile.We considered two anomalies corresponding to the Calabrian slab and to the Apennine collision belt.It can be seen that all these anomalies are robustly resolved in P and S models in correct depth intervals, although in the case of the S model, the amplitude of the resolved slab-related anomaly is much weaker compared to the true model.This test is especially important for this study, in view of the inconsistency in depth determination of the Calabrian slab in different tomography models discussed in the next section.

Discussion
The most prominent feature found in our tomography model is a high-velocity anomaly dipping down from the Calabrian arc, which is observed in Profile A1-B1 in Fig. 4 in both P and S anomalies.This anomaly has been detected in most other tomography results for the same region (e.g., Bijwaard and Spakman, 2000;Piromallo and Morelli, 2003;Cimini and Marchetti, 2006).In most of these studies, the positive P anomaly zone extends down to the depth of ∼ 700 km and is interpreted as a slab sinking down to the transition zone.In our results, anomalies with generally higher P velocities along section A1-B1 are also observed down to 700 km depth; however, the anomaly in 300-400 km depth seems to be separated from that at 600-700 km depth.The S velocity model in our result does not show any positive anomaly below 400 km depth.These and other facts make us to propose that these positive anomalies represent two different bodies.One of them is a narrow anomaly which steeply dips from the Calabrian trench down to 400 km depth and then propagates horizontally between 300 and 400 km depth.This anomaly  is clearly seen in both P and S models.The maximum concentration of deep seismicity at ∼ 300 km depth coincides with the bending part of this anomaly where it changes its shape from steeply dipping to horizontal.It is interesting that in Section A2-B2, this high-velocity pattern is seen as an isometrical anomaly.Taking into account possible smearing, the diameter of this body is estimated as 150-200 km.Indications are the same as in Figs. 3 and 4. 2010).Sinking of the narrow body occurs vertically, without strong dependence on the configuration of this body in space, whereas the flat plate slides down along the inclined bottom surface.
The deeper anomaly is seen in Section A1-B1 as a large high-velocity body located between 500 and 700 km depth.It is only observed in the P model; the S anomalies tend to be negative in this depth interval.We propose that this pattern represents the Alpine-Tethys remnant, which subducted in previous stages the paleoocean's closing (as proposed, for example, by Spakman and Wortel, 2004).If the time of stagnation was sufficiently large, no thermal anomaly, which affects the S velocity, is left, but there might be still compositional factors, which increase the P velocity.
The present configuration of anomalies beneath the central Mediterranean originated from a complex series of different geodynamic episodes which occurred in the western Mediterranean in the Cenozoic.According to the reconstruction proposed by Faccenna et al. (2004) shown in a simplified form in Fig. 8, the present configuration of the contact zone between African and European plates evolved from a classical subduction which occurred at about 35 Ma along the present coast of Spain and France (Fig. 8a).The back-arc processes caused the opening of the Ligurian-Provençal basin (Fig. 8b) and subsequently the Tyrrhenian basin (Fig. 8c) which resulted in counterclockwise rotation of the suture zone, which presently forms the Apennines, by approximately 90-100 • .During this rotation, along this suture zone, an ocean-ocean type of subduction took place and resulted in arc volcanism which is still recorded by calcalkaline material.At some moment, this subduction reached the margin of the Adriatic Shelf represented by a transition from oceanic to continental type of the lithosphere; the latter is more buoyant than the oceanic lithosphere.As a result, the oceanic lithosphere in the west was overriding the continental lithosphere of the Adriatic Shelf, and the subduction was transformed into a collision of continental type which caused shortening of the crust on both sides of the suture zone.This has led to active mountain building and strong deformation of the Apennine crust.We assume that the subducted oceanic lithosphere from the Adriatic side was detached and sank.The high P velocity anomaly observed in our tomography model below 500 km depth might be the trace of the remnant lithosphere.In the recent past, the continental crust of the western portion of the Adriatic plate became progressively accreted into the tectonic wedge, leading to the delamination and detachment of its lithospheric mantle which behaved as a subducting slab and remained below the Apennines (Roure et al., 2012), as shown in the scheme in Fig. 9.We propose that the observed high-velocity sausage-shaped anomaly beneath the Apennines represents the sunken part of the Adriatic lithosphere.Actually, its overall lateral and vertical lengths are quite compatible with the restored palinspastic surfaces of the Apenninic upper crustal and sedimentary units.Balanced cross-sections in the southern Apennines account for instance for more than 200 km of shortening among the Apulian and Apenninic platforms and intervening Lagonegro basinal units are currently accreted into the tectonic wedge (Casero et al., 1991;Roure et al., 1991Roure et al., , 2012)).
In Fig. 9 we present our interpretation of the main structures below the Apennines in present time and give a hy- According to the GPS data (e.g., Hollenstein et al., 2003;Nocquet and Calais, 2003), in present days, there is no considerable convergence between Calabria and the African plate (points D and F).The present displacements of Tunisia, Sicily and Calabria are not significantly different.At the same time, most scientists accept the retreat of the Calabrian arc.Thus, the point E tends to be replaced from D to F (Fig. 9).Between the southern and northern Apennines (points D and C), the geodetic observations record considerable shortening.The same conclusion follows from the analysis of stress patterns based on focal mechanisms (Serpelloni et al., 2007) which indicate the extension in SW-NE direction and compression in NW-SE direction.Between the northern Apennines and continental Europe there is an obvious shortening which resulted in the Plio-Quaternary overthrusting of the northern Apennines towards the Po Plain, following pre-Messinian episodes of mountain building in the facing southwestern segment of the Alps (Roure et al., 1989(Roure et al., , 1996;;Turrini et al., 2014).
To explain the detachment of the Adriatic lithospheric mantle, the key process is the NW-SE shortening of the Apennines between points C and D. The crust, which is a weaker part of the lithosphere, has been compressed and reduced its length.For the lithospheric mantle part, shortening was less plausible; instead it has been detached and formed a curved structure sank into the asthenosphere.In the tomography image we observe that the NW part of this body has been completely detached, whereas the SE end of it is still attached to the Calabrian arc.Vertical sinking of the "Calabrian sausage" causes the slab retreat of the contact zone between the newly formed Tyrrhenian Sea and the Ionian-African foreland lithosphere.
An open question relates to a long detachment or tearing of the Adrian mantle lithosphere that would allow the "sausage" to separate and sink into the mantle instead of dangling like a curtain.We hypothesize that the sunk part may correspond to the transitional ocean-to-continent type of the lithosphere with neutral or slightly negative buoyancy, whereas the remnant part of Adria is composed of buoyant continental lithosphere.In the case of regional shortening, they behave differently: the continental Adriatic part remains close to the surface, whereas the transitional part is pressed down and detached as a long sausage.However, this process appears to be too complicated, and its detailed understanding needs accurate three-dimensional thermo-mechanical simulations and more data on the deep structure.

Conclusions
Based on an analysis of the P and S tomographic models of the upper mantle beneath the central Mediterranean region derived from the tomography inversion of the ISC data, we present a new interpretation for the existence of the highvelocity anomaly beneath the Apennines, which was previously interpreted as a subduction zone of the former Tethyan oceanic lithosphere.We found that this anomaly behaves as a long narrow "sausage" with a steeply dipping part down to a depth of 400 km and then expanding horizontally over approximately 400 km.In cross section, this anomaly appears to be 150-200 km thick.We propose that this pattern represents the detached part of the Adriatic mantle lithosphere which was delaminated due to the final episode of the collision along the Apennines.The sinking of this segment was due to the negative buoyancy of the lithosphere material and was additionally triggered by NW-SE shortening of the Apennines.The NW part of this "sausage" was detached, whereas the SE end of it is still connected to the Ionian-African foreland lithosphere.
Author contributions.I. Koulakov, A. Jakovlev and I. Zabelina performed all tomographic calculations and figures preparation.I. Koulakov, S. El Khrepy and F. Roure provided geodynamical interpretation of presented results.I. Koulakov prepared manuscript with contributions of all co-authors.

Figure 1 .
Figure 1.Topography/bathymetry map of the Tyrrhenian Sea and Calabrian regions.Major structural units and sense of transport along main thrust fronts and extensional detachments are shown according to Faccenna et al. (2004).Red lines with arrows indicate thrust and subduction zones.

Figure 2 .
Figure 2. Distribution of data used for tomography in the study region (a) and globally (b).Circles indicate the target areas where the inversions were performed.Black and blue triangles denote stations inside and outside the target circles, respectively.Red and gray dots are the events inside and outside target circles, respectively.

Figure 3 .
Figure 3. P and S velocity anomalies in horizontal sections.

Figure 4 .
Figure 4. P and S velocity anomalies in vertical sections.Locations of the profiles are shown in Fig. 3. Dots indicate projections of the earthquakes located at distances of less than 50 km from the profile.Exaggerated relief along the profiles is presented above each plot.Dotted line indicates the intersection with another profile.
Thus, when viewed in 3-D space, this anomaly looks like a sausage-shaped body of ∼ 800 km long and ∼ 200 km thick penetrating to the upper mantle.It appears to be very different from the images of slabs in other subduction zones which usually behave as flat, conveyor-type plates.Mechanically, a sausage-shaped body in the viscous mantle should behave differently than a rigid flat plate (e.g.Loiselet et al.,

Figure 5 .
Figure 5. Inversion results with the use of smaller damping.Two horizontal and one vertical sections are shown for the P and S models.Indications are the same as in Figs. 3 and 4.

Figure 6 .
Figure 6.Checkerboard tests for two different P and S velocity models.In all cases, the signs of anomalies change at 200, 400, 600 km etc. Dotted lines mark boundaries of the synthetic anomalies.

Figure 7 .
Figure 7. Synthetic test with a synthetic model of realistic configuration.The shape of the synthetic model is highlighted with a contour.The thickness of the anomaly in the direction across the section is 200 km.

Figure 8 .
Figure 8. Simplified plate reconstruction in the western Mediterranean based on Faccenna et al. (2004) with our modifications.Dark blue indicate areas of back arc spreading.Yellow is Adriatic plate with transitional ocean-to-continent structure.Dark brown color in the Apennines and the Alps highlights shortening areas.Indications: Lig -Ligurian Sea, CR -Corsica, SA -Sardinia, Tyr -Tyrrhenian Sea, Io -Ionian Sea, AD -Adriatic Sea.

Figure 9 .
Figure 9. Schematic representation of the origin of the Calabrian sausage background in plots with sections of present configuration corresponds to P velocity anomalies in vertical sections 1 and 2. Black arrows schematically indicate the displacements of blocks; red arrow denotes the direction of the Calabrian trench retreat.Abbreviations: Ap -Apennines; Io -Ionian Sea; Tyr -Tyrrhenian Sea; AD -Adriatic Sea; Sa -Sardinia; Lig -Ligurian Sea.