Imaging the subducted slab under the Calabrian Arc, Italy, from receiver function analysis

Knowledge of the geometry of the subducting plates is critical for constraining tectonic models of the central Mediterranean Sea. Here, plates were reorganized during rapid micro-plate movements driven by slab rollback (Malinverno and Ryan, 1986). This was accompanied by both mountain building and extension as well as the development of a tightly curved orogenic belt, the Calabrian Arc (Fig. 1). Paleogeographic reconstructions show that the Calabrian Arc is the remnant of a once wider subduction zone, running along the entire Italian peninsula and along North Africa (Fig. 1). Neogene rollback of this larger Calabrian Arc opened the southern Tyrrhenian Sea in its wake. Simultaneously, the arc collided with Apulia and Africa to produce the Apennine and Maghrebides orogenic belts (Malinverno and Ryan, 1986; Gueguen et al., 1998; Faccenna et al., 2001; Lucente et al., 2006; Rosenbaum and Lister, 2004). The present Calabrian Arc is the last remnant that is still subducting oceanic lithosphere. The relationship of the subduction zone to the belts of crustal thickening and thinning is needed to evaluate subtle differences in the multiple scenarios proposed for the kinematic history (e.g., Gueguen et al., 1998; Rosenbaum and Lister, 2004).

At present, active subduction in the central Mediterranean is indicated by the occurrence of deep seismicity in the southeastern corner of the Tyrrhenian Sea, beneath the Calabrian Arc, where hypocenters define a 200-km-wide Wadati-Benioff plane continuous down to 500 km depth (Selvaggi and Chiarabba, 1995; Chiarabba et al., 2005). Global Positioning System (GPS) geodesy (D'Agostino and Selvaggi, 2004) shows that the rapid motion of Calabria has halted. It has been suggested that subduction is now waning (Westaway, 1993; Argnani, 2000); however, convergence across the plate boundary must still continue at the slow Africa-Eurasia rate. Alternatively, it has been proposed that subduction and rollback of the lower plate continue, and that only the upper plate is stuck (Gvirtzman and Nur, 1999, 2001). The knowledge of the geometry of the overriding and sub-ducting plates is an important key for interpreting the ongoing tectonics.

The Mesozoic crust of the subducting Ionian Sea plate is imaged by reflection seismics beneath the accretionary prism of the Calabrian Arc (Cernobori et al., 1996). In addition, the steep descent of the slab to depths of >500 km is imaged by seismicity (Selvaggi and Chiarabba, 1995; Chiarabba et al., 2005), tomography (Lucente et al., 1999; Wortel and Spakman, 2000; Piromallo and Morelli, 2003), and reflection seismics (Finetti, 2005a). However, the geometry of the subducting plate beneath Calabria, where it rolls over from the shallow dip of the subduction zone to the steep dip of the downgoing slab is largely unknown.

In this study, we present new data from the CAT/SCAN experiment, a temporary deployment of broadband seismic stations in southern Italy (Steckler et al., 2008; Baccheschi et al., 2007). Teleseisms recorded at 21 stations are used to produce an image of the Ionian Moho across the Sila Plateau and the northern Calabrian region (Fig. 2). We compute a set of teleseismic receiver functions (RFs) for each station and project them along an E-W section. A regional one-dimensional (1-D) S-wave velocity profile, derived from a large data set of RFs recorded at a permanent seismic station, is used to convert our time-domain section to depth. This along-profile representation of the computed RF allows us to identify the location where the overriding and downgoing crusts separate and the shallow subducting crust begins to plunge steeply into the mantle. Using a broader array of RFs, we constructed a map of the slab surface under Calabria in order to highlight the three-dimensional (3-D) geometry of the subducting Ionian plate Moho. These results provide insights into the tectonics of the Calabrian Arc, offering new information on the relationship between the shallow subducting plate geometry and changes of tectonic style observed at the surface and on the northward extent and geometry of the subducting oceanic lithosphere.

The Calabrian Arc rifted off Sardinia ca. 10–12 Ma ago, opening the Tyrrhenian Sea (Malinverno and Ryan, 1986; Faccenna et al., 2001; Lucente et al., 2006). The progressive collision of the Calabrian Arc with the Adria plate produced the Apennines on its NE side, and collision with Africa produced the Maghrebides-Sicilian nappes on its southern side (Fig. 1; Gueguen et al., 1998; Rosenbaum and Lister, 2004). Based on industry seismic data, the Southern Apennines appear to have stopped thrusting in the mid-Pleistocene (Patacca and Scandone, 2001; Ferranti and Oldow, 2005a, 2005b). On the Sicilian side, thrusting has also halted, but a zone of earthquakes north of Sicily attest to continued convergence and possibly incipient subduction reversal (Pondrelli et al., 2004; Billi et al., 2007). The narrow Calabrian Arc thus represents the sole remaining segment still subducting the old Mesozoic oceanic crust of the Ionian Sea.

The incoming Ionian plate has a 5–6-km-thick sediment cover (Makris et al., 1986; de Voogd et al., 1992) including 1.2 km of Messinian evaporites (Montedert et al., 1978; de Voogd et al., 1992). This sediment cover has built a large accretionary prism, whose tip, characterized by active compressional structures (e.g., Finetti, 2005a; Cernobori et al., 1996; Doglioni et al., 1999), is ~250 km offshore from Calabria in the Ionian abyssal plain. The bathymetry of the prism shallows gradually to Calabria, where the imbricate thrust stack is capped by exotic crystalline basement. The Moho of the incoming Ionian plate deepens from 17 km under the Ionian abyssal plain to >30 km as it approaches the eastern coast of Calabria (Ferrucci et al., 1991; Cernobori et al., 1996). Beneath Calabria, data on crustal thickness (Lüschen et al., 1992; Çifçi and Dondurur, 2002; Barberi et al., 2004) are limited. Interpretations show Ionian-crust Moho deepening to 40–45 km, whereas the Tyrrhenian side Moho is 20–25 km deep (e.g., Gvirtzman and Nur, 2001; Finetti, 2005b). The transition between the two Mohos, Ionian and Tyrrhenian, or where the descending crust separates from the overriding-plate crust, is largely unconstrained and varies from model to model (Gvirtzman and Nur, 2001; Dèzes and Ziegler, 2001; Cassinis et al., 2003, 2005; Finetti, 2005b).

The 200-km-wide Calabrian Arc, comprising Calabria and the Peloritani Mountains of NE Sicily, is capped by Hercynian crystalline rocks similar to those appearing in Sardinia (Fig. 1). These units form the elevated (1–2 km) backbone of Calabria, including the squarish Sila Plateau and the linear Catena Costiera (Fig. 2) covered by the CAT/SCAN experiment. These two uplands are separated by the Crati Valley, one of the arc-parallel (longitudinal) extensional basins found along the Tyrrhenian side of Calabria. These longitudinal basins are generally thought to be asymmetric grabens, with the main faults bordering the higher eastern flank (e.g., Ghisetti, 1979). These structural features are present all along the Tyrrhenian side of Calabria and the Apennines and are associated with the postorogenic extensional tectonics on top of a retreating subduction slab (Cifelli et al., 2007, and references therein). In addition, arc-normal (radial) basins or structural inflections (Fig. 2) interrupt longitudinal structures in Calabria (Tortorici et al., 1995).

The origin of these transverse basins, such as the Sibari Plain, is more uncertain and may be related to stretching of the Arc accommodated by transcurrent faulting (Ciaranfi et al., 1983; Van Dijk et al., 2000; Seeber et al., 2008) or to slab tear-faulting episodes during subduction segmentation (Rosenbaum et al., 2008).

Information on the structure of the crust and upper mantle can be derived by the analysis of mode-converted seismic waves. The RF method aims to isolate P- to S-wave conversions associated with crust and mantle seismic discontinuities near the receiver (Langston, 1979). We compute RFs using teleseismic data recorded at 21 broadband seismic stations during the CAT/SCAN experiment from December 2003 to October 2005 (Fig. 2). The stations were deployed across northern Calabria with an average spacing of 10 km. We also include data from the permanent MedNet seismic station TIP, located on the SE flank of the Sila Plateau. We select 586 teleseismic events of M_w ≥ 5.5 from an epicentral distance between 30° and 100° (Fig. 3). Figure 2 shows the station distribution together with the Ionian Moho piercing points of the teleseismic rays used for the RFs. Piercing points were computed using Ionian Moho depths found in this study.

We compute RFs using a frequency-domain method (Di Bona, 1998) based on the deconvolution of the vertical record from both the radial and transverse horizontal records. This technique allows us to compute single RF variance. We select “good” RFs using a threshold for RF variance. Noisy RFs were discarded. We focus on the radial RFs—i.e., on the isotropic bulk subsurface structure. Whereas transverse RFs can be useful for ascertaining the presence of subsurface anisotropic layers or dipping discontinuities; our ability to utilize them is limited by the incomplete back-azimuth coverage at some stations. However, radial RFs have proven to be effective in reconstructing the structure of complex tectonic settings, such as subduction zones (e.g., Bannister et al., 2007).

RFs are band-passed using a Gaussian filter with parameter a = 1, which excludes frequencies higher than ~0.5 Hz. We collect an average of 40 RFs with a high signal-to-noise ratio for each station. To obtain an interpretable image of the crustal and uppermost mantle structure across the subduction, we adopt the method developed in Ferris et al. (2003): we stack together all radial RFs for a single station, producing a “summary” RF for each station (Fig. 4). Summary RFs within 35 km are projected along an E-W profile (i.e., perpendicular to the local slab strike inferred from the deep earthquake distribution, Fig. 2B). Then, the selected summary RFs are binned using a box-shaped filter with 10 km half-width and 50% overlapping scheme. This stratagem emphasizes coherent energy from nearby stations while minimizing effects of scattering and noise, and furnishes a readable time-domain representation of the whole RF data set (Fig. 5). The time-domain image in Figure 5 does not depend on a migration model and gives a clear vision of the geometry of the main subsurface structures.

To convert from time to depth, we apply a 1-D S-wave velocity profile computed at TIP, a permanent seismic station belonging to the MedNet network (Mazza et al., 2008), through a Monte Carlo inversion of its large RF data set. The RF inversion is inherently a nonlinear optimization problem. In order to correctly retrieve the properties of the crust and upper mantle seismic structure, we apply a neighborhood algorithm (NA) inverse method. Such an inversion scheme was developed by Sambridge (1999) in the framework of global directed search algorithms. We refer to Sambridge (1999) for details about NA. This method has been successfully applied to the RF inverse problem (Piana Agostinetti et al., 2002; Sherrington et al., 2004; Bianchi et al., 2008).

We stack RFs coming from the same epicentral distance, and we use these stacks as input for the RF inversion code. Input RF stacks are shown in Figure 6. We use as input for the NA the RFs filtered using different frequencies (low-passed at 0.5, 1.0, and 2.0 Hz) to obtain an S-wave velocity profile that can fit both the gross and fine details of the RF time series. Discrepancy between synthetic and observed RF is measured through a classical chi-squared misfit function. The parameter space searched by the NA comprises five homogeneous layers overlying a half-space. The layers and half-space are separated by planar horizontal seismic interfaces. A priori information about the S-wave velocity at depth is given in the form of minimum and maximum allowed values. S-wave velocity (V_s) boundaries are reported in Figure 7. The initial velocity probability distribution is uniform between the minimum and maximum values. We tune our NA to explore this parameter space in the following way: an ensemble of 900 models that randomly sample the parameter space represents our initial population. Then, the parameter space around the six best-fit models is sampled to produce 36 new samples (i.e., six new samples in the neighborhood of each best-fit model) that are added to the previously sampled population. Next, we select the new six best-fit models inside the entire population and repeat the previous step. This procedure is iterated for 500 resampling steps. As a result, 18,900 models are generated. Finally, we select a best-fit family of sampled models, accepting all models with a misfit equal or <1.25 times the misfit of the best-fit model.

Results from the inversion process are presented in Figure 7, and synthetic RFs computed using the best-fit S-wave velocity model are compared to observed RFs in Figure 6. A first-order value for the error associated with the depth of each of the discontinuities in the bestfit model comes from the standard deviation of the depth, computed from the average depth of that parameter in the best-fit family. Errors computed in this way are <2 km.

The 1-D S-wave velocity model derived from TIP (Fig. 7) is characterized by a strong velocity inversion at ~21 km depth, followed at depth by two positive velocity jumps at ~30 and ~35 km, respectively. From the comparison of the observed and synthetic RFs (Fig. 6), we deduce that the strong negative phase (PS1, hereinafter), observed at ~3 s, is responsible for the S-wave velocity inversion (i.e., an S-wave velocity that decreases with depth) at ~21 km depth, whereas the broader positive phase at ~5 s (PS2, hereinafter) is a combined phase that contains the converted phases from the two deeper interfaces with S-wave velocity increases. The lowermost interface in the 1-D model, at ~35 km (Fig. 7), marks the transition to S-wave velocity values typical for the mantle (>4.0 km/s). We therefore interpret this interface as the Moho of the subducting Ionian plate, and the interface 5 km above it as the top of the subducting Ionian plate crust. The nature of the velocity inversion at ~21 km depth is more obscure. We associate the S-wave velocity inversion (Fig. 7) with the contrast between the Calabrian continental crust and subducting oceanic sediments below. We interpret the thick low-velocity layer at 21–30 km depth as a mixture of both in-place sediments and off-scraped–imbricated sediments from the down-going Ionian plate. Such accreted sediments crop out east of the Sila Plateau (Van Dijk et al., 2000). The 21 km thickness of the overlying crust is consistent with estimates from seismic refraction (Çifçi and Dondurur, 2002; Cassinis et al., 2005; Finetti, 2005b) from the western, Tyrrhenian, side of Calabria. Our interpretation agrees with other RF studies in different sub-duction zones where subducting sediments were imaged as low-velocity materials directly below the overriding crust (i.e., Savage et al., 2007).

The time-depth relationship at TIP enables us to convert local earthquake depths (Chiarabba et al., 2005) to S-wave traveltime by back-migrating the hypocenters. We are then able to plot the earthquakes over the time-domain RF image (Fig. 5) and thus to correlate seismicity to the subsurface structure. Conversely, we use the 1-D S-wave velocity model determined at TIP station as a reference velocity for the other stations. The summary events for each station are migrated to depth with the bulk seismic velocity determined at TIP. A depth-migrated map of interface interpreted as Ionian Moho at all of the stations is shown in Figure 8 (Moho depth and pulse amplitude are presented in Table 1 and Fig. 4), highlighting the complex 3-D geometry of this main mode-converting discontinuity. We observe that Ionian Moho depth values, computed using the 1D S-wave velocity model obtained for TIP, represent a lower bound for the stations that lie along the Tyrrhenian coast. In this area the presence of mantle wedge materials between the stations and the Ionian Moho should increase the Ionian Moho depth values estimated here, increasing the dip angle of the Moho toward the Tyrrhenian Sea. However, our approximation does not influence the overall geometry of the retrieved surface.

The RF data set (Fig. 5) provides a clear and coherent image of the crust and uppermost mantle to depths of ~100 km across the Calabrian Arc. Here we discuss the most striking features in the RF profile, whose continuity can be followed along the entire 90 km length of the profile.

On the eastern, Ionian, side, we observe a strong positive pulse at ~5.5 s (i.e., marked PS2 in TIP data set), which gently dips under the Sila Plateau. Its amplitude decays toward the west. This phase follows a negative pulse ~2.5 s (i.e., PS1 in TIP data set), which becomes broader toward the west. The continuity of these two phases along the profile constrains the presence of two interfaces at depth: The upper one is associated with the S-wave velocity inversion or decrease with depth, whereas the lower one is associated with a strong positive V_s jump. The modeling of TIP indicates that PS2 is associated with a two-step velocity increase. These interfaces dip gently to the west, the direction of subduction.

Beneath the Crati Basin on the western side of Calabria, the arrival time of the PS2 phase increases rapidly, whereas the PS1 phase almost disappears. The increase of the arrival time of the PS2 phase indicates a change in the dip of the lowermost interface. PS1 phase degradation possibly marks a strong change in the nature of the velocity contrast associated with the upper interface owing to the dehydration of the subducting materials (Kawakatsu and Watada, 2007).

Two other pulses are present in our RF image, visible mainly under the Ionian side of the profile: a negative phase, which follows PS2 in time, and a later positive phase at ~13 s. As shown in Figure 6, the negative phase can be modeled as a multiple from our upper interface. The delay of the second positive phase seems to increase under the Calabrian Arc, mirroring the pattern observed for the PS2 pulse, but its arrival time does not coincide with arrival times for PpPs multiple phases from our lower interface. From a simple depth migration, we argue for a depth of >100 km for the converting interface for this phase. We therefore associate this phase with the lithosphere-asthenosphere boundary, but we do not proceed further with interpreting this feature, leaving it to techniques better suited for investigating the lithosphere-asthenosphere boundary (e.g., S-wave RF; Kumar et al., 2007).

Figure 5 also compares the RF image to the local seismicity (Chiarabba et al., 2005). We show events with focal depths >10 km that have depth and horizontal location errors <1 km. In the eastern part of our profile, the deeper local seismicity seems to have been confined to the lowermost part of the crust (3–5 s) between PS1 and PS2 (Fig. 5). Under the Crati Valley, local earthquakes shallow slightly, but other events are deeper (time-delay >5 s) and appear to have followed the dip of the PS2 pulse.

Our interpretation of the RF profile is presented in Figure 9. We interpret PS2 as marking the position of the underthrust lower plate (Figs. 5 and 9). Where it is stronger, beneath the Sila Plateau, it is a combined phase from the velocity jumps at the top and bottom of the Ionian oceanic crust. The steepening and weakening of PS2 to the west correspond to the plunge of the plate into the mantle. Here the signal is primarily from the Ionian Moho, as the upper surface of the crust, and any low-velocity material associated with it lies beneath the mantle wedge. The 10-km-thick, low-velocity zone between PS1 and PS2 is thicker than the incoming sediment thickness, implying that low-velocity, water-rich sediments are accumulating before being dewatered and compacted, and then accreted to the overlying Calabrian crust, much of which is probably composed of such offscraped sediments.

West of the S-point, the gap between the thinner “Tyrrhenian” crust and the top of the slab must be filled by the mantle wedge (Fig. 9). The S-point is the point at which the overriding and downgoing crusts separate and mark a discontinuity in the velocity field that controls the growth of the orogenic wedge (Willett et al., 1993). We also note that west of the S-point the amplitude of PS1 also decreases (Fig. 5). The thick zone of imbricated sediments is probably no longer present. The negative conversion over the slab in this area is probably due to the velocity decrease from the mantle wedge to the top of any remaining sediments, the downgoing plate, or serpentinites associated with this transition.

The steeping of the underthrust Ionian plate into a subducted slab coincides with the transition from the uplifted Sila Plateau to the extensional Crati Valley. The entire Apennine chain exhibits a transition from compression on the Ionian-Adriatic side of the mountain belt to extension on the Tyrrhenian side (Chiarabba et al., 2005). Our RFs reveal that the transition between these two regimes in Calabria overlies the S-point (Willett et al., 1993). Thus in Calabria, the boundary between the thickened “Ionian” crust and the thinned “Tyrrhenian” crust coincides with the change of tectonic style at the surface. Previous estimates of this crustal transition in Calabria (Fig. 5) ranged from near the Tyrrhenian coast to offshore from the Ionian coast (Gvirtzman and Nur, 2001; Dèzes and Ziegler, 2001; Cassinis et al., 2003, 2005; Finetti, 2005b). Based on estimates having the S-point in the Ionian, Gvirtzman and Nur (2001) suggested that flow of hot asthenosphere into a mantle wedge beneath Calabria caused the current rapid uplift (Westaway, 1993; Bordoni and Valensise, 1998). Our data indicate that the asthenospheric wedge beneath Calabria extends only half as far in front of the volcanic arc as Gvirtzman and Nur (2001) estimated and not at all under the Sila Plateau. The thicker crust under Sila removes the large mantle buoyancy anomaly they estimated there, implying that this was not the cause of the Sila Plateau uplift. Instead, the interpreted 2-D image (Fig. 9) also suggests that underplating and/or slab-bending processes may have provided larger contributions to the observed rapid uplift of the Sila Plateau (Westaway, 1993; Bordoni and Valensise, 1998; Ferranti et al., 2006). We also note the similarity of the amplitude and width of the low-velocity zone to the height of the Sila Plateau.

The coincidence of the structural transitions at depth and at the surface suggests that the lower plate has not rolled back relative to the upper plate since Calabria slowed. The slowing of Calabria since the opening of the Marsili oceanic basin probably occurred between 1.6 Ma ago, the oldest estimate (and probably an overestimate) for the end of spreading there (Nicolosi et al., 2006), and 0.66 Ma ago, the end of thrust advance in the Southern Apennines (Patacca and Scandone, 2001). If the lower plate had continued to roll back during the Pleistocene at the 6 cm/a average Pliocene rate of the Calabrian Arc, then the hinge in the subducting plate should have shifted by at least 40–60 km relative to the Calabrian upper plate. In that case, before the tectonics changed, the S-point would have been well offshore in the Ionian Sea. Although not conclusive, the S-point position beneath the surface tectonic transition suggests that the SE motion of the Calabrian Arc and the hinge in the Ionian Sea plate ceased at the same time.

The Moho depth map (Fig. 8) reveals a complex geometry for the subducting plate. The Ionian Moho is shallowest under the Sila Plateau. In the southern part of the study area the Moho dips westward toward the Tyrrhenian coast, almost perpendicular to the Crati Valley (Fig. 2), steepening to become the subducted slab. However, in the northern part of the plateau the dip vector rotates to NNW, normal to the radial Sibari Plain basin. The northern edge of the Sibari Plain corresponds to the surface boundary between the Calabrian lithospheric block and the Southern Apennines. The existence of a lithospheric discontinuity, such as a slab tear or a slab window, was proposed for this surface junction (e.g., Gvirtzman and Nur, 1999; Rosenbaum and Lister, 2004). Older tomographic studies show a slab gap farther north in the Southern Apennines (e.g., Lucente et al., 1999; Piromallo and Morelli, 2003; Montuori et al., 2007), but a recent high-resolution tomographic study (Chiarabba et al., 2008) better locates the slab discontinuity at the Calabria-Apennine boundary. In this case, the NNW deepening of the Ionian Moho is related to this recent tear in the slab, which possibly was caused by the incipient subduction of continental Apulian lithosphere in the Apennines. The tearing of the slab therefore would be associated with SE rollback of the subduction system. The late Pliocene–Pleistocene formation of the Crati and Sibari Basins would be related to extension from the rollback and a clockwise rotation of Calabria (Cifelli et al., 2007). However, we note that the similar depths (Table 1) of the Moho found near the Tyrrhenian coast (Fig. 8) and the 85 km observed at nearby MedNet station CUC (Piana Agostinetti et al., 2008) suggest continuity of the descending slab across the Calabria-Apennines transition at these depths. Indeed, Chiarabba et al. (2008) found the slab to be torn only below 100 km. Alternative explanations for the NNW Moho deepening include the underthrusting of the thicker Apulia plate in the Southern Apennines and thrusting of the Apennines over Calabria.

Receiver functions in Calabria using data from CAT/SCAN provide new constraints on the crust and uppermost mantle structure of the Calabrian Arc, one of the narrowest subduction zones in the world. An E-W profile shows a clear image of the subducting Ionian plate where it bends and starts to descend into the mantle. The depth of the Ionian Moho increases from ~35 km, where the Ionian lower plate is in contact with the Calabrian upper plate crust, to ~80 km where the Ionian crust separates and descends into the mantle. An S-wave velocity profile shows a low-velocity zone, probably made of imbricated sediments, above the oceanic Ionian crust. The change in the subduction zone dip corresponds to the transition from the uplifted Sila Plateau to the extensional Crati Basin at the surface. This suggests that the lower Ionian plate has not continued to roll back relative to the upper Calabian plate, implying that the uplift of the Sila Plateau is more likely due to underplating than to growth of the mantle wedge. However, the 3-D pattern of Moho depths also shows an unexpected deepening northward toward the Southern Apennines. This could be related to a proposed tear in the slab in the Southern Apennines or to thrusting that involves either Apulia, the Apennines, and/or Calabria.

We would like to thank all our colleagues who were involved in the CAT/SCAN experiment. Ignazio Guerra, Anna Gervasi, and the staff of the Dipartimento di Fisica (Universitá della Calabria), in particular, gave fundamental support to the seismic deployment in Calabria. We would also like to thank the members of the Calabrian Arc project for valuable and informative discussions. The manuscript was improved by reviews from Martha Savage and an anonymous reviewer. This work was funded by U.S. National Science Foundation grants EAR99-10554, EAR06-07687, and by INGV. Data collection and archival were facilitated by IRIS. This is LDEO publication 7268.