Subduction or delamination beneath Apennines ? Evidences from regional tomography

Subduction or delamination beneath Apennines? Evidences from regional tomography I. Koulakov, A. Jakovlev, I. Zabelina, F. Roure, S. Cloetingh, S. El Khrepy, and N. Al-Arifi Trofimuk Institute of Petroleum Geology and Geophysics SB RAS, Prospekt Koptyuga, 3, 630090, Novosibirsk, Russian Federation Novosibirsk State University, Novosibirsk, Russia, Pirogova 2, 630090, Novosibirsk, Russia IFP-Energies Nouvelles, Rueil-Malmaison Tectonics Group, Utrecht University, the Netherlands King Saud University, Riyadh, Saudi Arabia, P.O. Box 2455, Riyadh 11451, Saudi Arabia National Research Institute of Astronomy and Geophysics, NRIAG, 11421, Helwan, Egypt


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 (Figure 1).Today for instance, even the nature of the lithosphere and crust beneath the deep Ionian Basin remains debated, being either made up of a remnant of the former 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 3D seismic model of the upper mantle beneath the Italian region constructed based on generally same calculation schemes as in Koulakov et al. [2009], but using a considerably larger dataset.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 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 nineties, 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 Pvelocity 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. e.g. [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 velocity 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, Apennines, Alps and Dinarides is identified (Figure 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 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 down going 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 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 accumulated 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 Figure 2A) recorded by the worldwide stations (blue triangles in Figure 2B), as well as data from ~90,000 teleseismic events (grey dots in Figure 2B) recorded by 1915 seismic stations located inside the study area (black triangles in Figure 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 (Figure 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 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-by-window 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.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 degrees) and then combined in a single model.When merging the results computed for different circular windows and differently oriented grids, we calculate a 3D weight function.It depends on the distance from the nearest parameterization node: if the distance is less than 40 km, the weight equal 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 algorithm of tomographic inversion was developed by
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 Kurile-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 images 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 as the main result of our calculations, are presented in horizontal and vertical sections in Figures 3 and 4. The anomalies are given in percent in respect to the 1D 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 Figures 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 Figure 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 Figure 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 km, 400 km 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 km, 300 km and 500 km (middle depths of the checkerboard levels).For the Panomalies with 1º size patterns, 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 Figure 7 is aimed at studying the capacity of the algorithm to resolve a sausage-shaped anomaly similar to that we obtain in the observed data inversion.We define the synthetic anomalies along the Profile A1-B1, same as used for presenting the main model in Figure 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 Figure 4 in both Pand 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.Thus, when viewed in 3D space, this anomaly looks as 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., 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 of closing the paleoocean [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 central Mediterranean originated from a complex series of different geodynamic episodes which occurred in Western Mediterranean in Cenozoic.According to the reconstruction proposed by Faccenna et al. [2004] shown in a simplified form in Figure 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 (Figure 8a).The back-arc processes caused the opening of the Ligurian-Provençal basin (Figure 8b) and subsequently the Tyrrhenian basin (Figure 8c) which resulted in counterclockwise rotation of the suture zone, which presently forms the Apennines, by approximately 90-100 degrees.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 calcalcaline 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 in 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 Figure 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 Lago Negro basinal units are currently accreted into the tectonic wedge [Casero et al., 1991;Roure et al., 1991Roure et al., , 2012]].
In Figure 9 we present our interpretation of the main structures below the Apennines in present time and give a hypothetical reconstruction to the recent past.Sections correspond to the location of the Profiles A1-B1 and A2-B2.In Section A1-B1 parallel to the Apennines, reference points are used for the interpretation: A is located in the continental part of Europe; B marks the suture zone in the Alps; C is a point in Po plain in Northern Italy; D marks the southernmost part of the Apennines, E marks the Calabrian arc, F is a point in the African Plate.
According to the GPS data [e.g., Hollenstein et al., 2003, Nocquet andCalais, 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 (Figure 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 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 daggling there like a curtain.We hypothesize that the sunk part may correspond to the transitional ocean-tocontinent 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 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 thermomechanical simulations and mode data on the deep structure.

Conclusions
Based on an analysis of the P-and S tomographic models of the upper mantle beneath the Central Mediterranian region derived from the tomography inversion of the ISC data, we present a new interpretation for the existence of the high-velocity 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 Apennines.The NW part of this «sausage» was detached, whereas the SE end of it is still connected to the Ionian-African foreland lithosphere.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 grey tots are the events inside and outside target circles, respectively.

Figure 1 .
Figure 1.Topography/bathymetry map of the Tyrrenian Sea and Calabrian regions.Major

Figure 2 .
Figure 2. Distribution of data used for tomography in the study region (A) and globally (B).

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

Figure 5 .
Figure 5. Inversion results with the use of smaller damping.Two horizontal and one vertical

Figure 6 .
Figure 6.Checkerboard tests for two different P-and S-velocity models.In all cases, the signs

Figure 7 .
Figure 7. Synthetic test with a synthetic model of realistic configuration.The shape of the

Figure 8 .
Figure 8. Simplified plate reconstruction in western Mediterranean based on Faccenna et al.,

Figure 9 .
Figure 9.Schematic representation of the origin of the Calabrian Sausage Background in plots

Figure 1 .
Figure 1.Topography/bathymetry map of the Tyrrenian 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 grey tots 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 Figure 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.

Figure 5 .
Figure 5. Inversion results with the use of smaller damping compared to the main model shown in Figures 3 and 4. Two horizontal and one vertical sections are shown for the P and S models.Indications are the same as in Figures 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 km, 400 km, 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 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 Apennines and Alps highlights shortening areas.Indications: Lig -Ligurian Sea, CR -Corsica, SA -Sardinia, Tyr -Tyrrenian 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 arrow schematically indicate the displacements of blocks; red arrow denotes the direction of the Calabrian trench retreat.Abbreviations: Ap -Apennines; Io -Ionian Sea; Tyr -Tyrrenian Sea; AD -Adriatic Sea; Sa -Sardinia; Lig -Ligurian Sea.