The formation and evolution of a backarc basin are linked to the dynamics of the subduction system. The opening of the central Mediterranean basins is a well-documented example of backarc extension characterized by short-lived episodes of spreading. The underlying reasons for this episodicity are obscured by the complexity of this subduction system, in which multiple continental blocks enter the subduction zone. We present results from three-dimensional numerical models of laterally varying subduction to explain the mechanism of backarc basin opening and the episodic style of spreading. Our results show that efficient backarc extension can be obtained with an along-trench variation in slab buoyancy that produces localized deformation within the overriding plate. We observe peaks in the trench retreating velocity corresponding first to the opening of the backarc basin, and later to the formation of slab windows. We suggest that the observed episodic trench retreat behavior in the central Mediterranean is caused by the formation of slab windows.

Backarc basins (BABs) form in convergent margins where the overriding plate is in an extensional regime (e.g., Uyeda and Kanamori, 1979; Carlson and Melia, 1984). Their formation is kinematically associated either with trench retreat (e.g., Japan, Lau-Kermadec, and Philippine basins) or with motion of the upper plate away from the trench (e.g., the Mariana Islands). What triggers backarc extension and why it is episodic (lasting only for tens of millions of years; e.g., Karig, 1971) are debated. In several cases, such as the Woodlark Basin (Papua New Guinea) or Havre Trough BAB (New Zealand), it has been noted that backarc extension takes place by an adjacent collisional zone inducing sharp rotation and curvature of the adjacent region (Wallace et al., 2009). Such a setting is typical for the Mediterranean. In the eastern Mediterranean, for example, backarc extension in the Aegean mainly occurred during Arabia-Eurasia collision and accelerated after the associated slab tearing along western Turkey in the east, and near Corinth in the west (Faccenna et al., 2006; Jolivet et al., 2013). Moreover, in the central Mediterranean, the retreat of the Calabrian slab led to the formation of two main BABs. Geological records show how this process occurred episodically, with pulses of greater extension and oceanic spreading (e.g., Malinverno and Ryan, 1986; Faccenna et al., 2001a) (Fig. 1). The first episode of extension occurred between 30 and 16 Ma and led to the formation of the Liguro-Provençal basin. After a pause of ∼5 m.y., extension shifted to the east of Sardinia to form the Tyrrhenian basin (where rifting started ca. 12–10 Ma). Here, pulses of oceanic spreading migrated eastward and formed the Vavilov Basin (ca. 4.3–2.6 Ma; Kastens and Mascle, 1990), followed by the Marsili Basin (ca. 2 Ma; Patacca et al., 1990; Nicolosi et al., 2006). At depth, the slab broke off and tore where collision occurred; first in northern Africa, ca. 10–12 Ma, as evidenced by the change in the nature of volcanism of the Mogods volcanic complex (Faccenna et al., 2005), and later in the Apennines. The lateral tearing of the slab along the continental margins adjacent to the oceanic subduction accommodates the Calabrian slab rollback (Carminati et al., 1998; Wortel and Spakman, 2000). The present-day configuration of the central Mediterranean subduction zone is supported by tomographic images (Piromallo and Morelli, 2003; Fig. 1) that show two slab windows next to the narrow Calabrian slab. The nature of the episodic behavior of BAB opening is debated. It has been postulated that slowing of extension can be driven by buoyant continental crust entering the trench (Malinverno and Ryan, 1986), by the interaction between the slab tip and the 660 km seismic discontinuity (Faccenna et al., 2001b), or by the collision-related lateral slab tear and subsequent development of vigorous toroidal flow around the slab (Faccenna et al., 2005, 2007; Guillaume et al., 2010).

Inspired by the geological cases described here, and in particular by the Mediterranean example, we developed three-dimensional (3-D) numerical models to explore how laterally varying collision may influence the formation and style of backarc extension. Our results show that BABs form next to collisional systems and that they develop episodically if the boundary between oceanic and continental lithosphere ruptures.

Previous studies helped our understanding of the mechanisms of the BAB formation through 2-D (Arcay et al., 2008; Gerya and Meilick, 2011) and laterally heterogeneous 3-D models (Capitanio and Replumaz, 2013; Li et al., 2013), and the episodicity in backarc tectonic regimes (Clark et al., 2008; Baitsch-Ghirardello et al., 2014). However, until now, all these aspects have been studied separately.

We present a set of 3-D numerical models of subduction using the finite element code Citcom (Moresi and Gurnis, 1996; Zhong et al., 2000; van Hunen et al., 2005) to explore the role of collision on BAB opening. Conservation of mass, momentum, thermal energy, and composition describe material flow, deformation, and thermal state. Continental collision is simulated using a viscoplastic rheology that combines diffusion creep, dislocation creep, and a pseudo-brittle deformation simulating Byerlee’s law (van Hunen and Allen, 2011), while the converging plates are decoupled through a mobile weak zone (Magni et al., 2012; see the GSA Data Repository1 for details on the numerical method).

To investigate the role of lateral continental colliding blocks on the side of an oceanic setting, we vary the width of the subducting oceanic plate from 500 to 1200 km (Fig. 2) within a 3300 × 3960 × 660 km regional 3-D model, in which the overriding plate is fully continental, and the subducting plate is oceanic in the middle and continental at the sides (Fig. 2). For every model we compute the potential energy release rate (E′pot) and the energy dissipation rate (E′diss) to investigate how energy is distributed within the system (for details on how energy is computed, see the Data Repository). The models are not designed to simulate specifically any particular case, but rather to capture the dynamics of the system. All of these models show a common behavior during the first phase of evolution, when only oceanic lithosphere subducts and the entire slab slowly rolls back, causing trench retreat. When the continental blocks arrive at the trench, subduction locally slows and the trench starts to advance, while the oceanic part of the plate continues to subduct and sustains trench retreat. These opposite trench motions cause large strains in the overriding plate and a strong curvature of the trench and slab. The largest values of E′diss in the overriding plate are mainly localized where ongoing oceanic subduction causes significant overriding plate deformation (Fig. 3A). For oceanic plates larger than 500 km, the overriding plate behind the oceanic subduction stretches and thins to accommodate the retreat until a BAB is formed ∼10 m.y. after initial collision (Fig. 3B). During this process, the retreating velocity increases rapidly, from 1–2 cm/yr to 4–5.5 cm/yr in just a few million years (Fig. 4A). After the opening of the BAB, E′pot integrated over the entire lithosphere slightly increases in response to this episode of faster trench retreat and associated slab sinking (Fig. 4B). This is because the sinking of cold, dense slabs is responsible for the bulk of the potential energy release that powers the subduction dynamics (and, more generally, plate tectonics). Meanwhile, the E′diss in the overriding plate drops (Fig. 4B), because most energy is now dissipated in the rapidly deforming subducting plate. At depth, the slab is highly deformed with a curved concave shape. The middle, oceanic part of the plate is still subducting, and slab pull is still acting there, whereas near the continental blocks, subduction almost stops and the slab starts breaking off in a necking process. At this point, rollback slows. If the oceanic plate is wider than 800 km, the oceanic slab pull is large enough to stretch and rupture the slab at the ocean-continent boundaries at a depth of ∼230 km (Fig. 3D). These ruptures then rapidly propagate laterally, allowing, in only a few million years, the formation of two wide slab windows, where subslab mantle material can flow through (Fig. 3E). This process corresponds to a second peak of acceleration in the trench retreat velocity (Fig. 4A) and peak of increased E′pot (Fig. 4B) that are not observed in narrower oceanic plates. The remaining oceanic slab becomes narrower (∼400 km wide in the final stages) and the retreat slows.

A very different behavior is observed for the model with a 500-km-wide or narrower oceanic plate, where the retreat is very slow (<1 cm/yr), the slab deforms slightly, and no BAB forms.

We put the behavior of subducting ocean-continent plates in a wider context by investigating simpler models without along-strike variation in the subducting plate. Subduction of a purely oceanic plate (i.e., no continents at the sides) results in a very slow trench retreat and no BAB formation (Fig. DR1 in the Data Repository). In a model with one large continental block colliding along the entire trench after the subduction of oceanic lithosphere, the initial phase is very similar to the previous models. However, once collision occurs, subduction slows and the slab starts to advance slowly (Fig. 4A). At this point, the break-off starts in the center of the slab ∼19 m.y. after collision and quickly propagates laterally (van Hunen and Allen, 2011) until the entire slab detaches from the shallow part of the plate. In both models without lateral variation (purely oceanic and laterally uniform continental lithosphere), no BAB is formed.

These models highlight a number of essential ingredients for the formation of BABs. Localized deformation arising from rheological heterogeneities seems to be a key factor in the formation of a BAB. The model with only oceanic subduction shows that although the viscous dissipation in the overriding plate is high, it is not localized enough to open a BAB. In such a setting, BAB formation is still possible through other physical factors that are not taken into account in our models, such as the weakening of the overriding plate due to the presence of fluids and melt (Arcay et al., 2008; Gerya and Meilick, 2011).

When a continental block collides, the trench motion changes from retreating to advancing with respect to the upper plate and, where the trench remains parallel to the margin, the slab keeps advancing for few million years after collision (Magni et al., 2012). In the case of a continent of a finite size, a torque will be created between the continental advancing sides and the oceanic retreating one. The resulting localized tensile stresses in the overriding plate provide a means for localized stretching, which can lead to the opening of a BAB. The larger the oceanic plate, the greater the potential energy available to the system (Fig. 4B). Therefore, larger oceanic plates produce a larger and faster rollback, rapidly inducing a BAB, whereas smaller oceanic plates may result in only a small amount of retreat, which is not vigorous enough to produce a BAB.

Slab pull also plays a crucial role in the formation of the slab windows at the ocean-continent boundaries, which are observed in the models only for a wide oceanic plate. Similar to BAB formation, slab tears at ocean-continent boundaries form under high tensile stresses, and can be generated only if the amount of slab pull is high enough. Once the oceanic slab has started tearing away from the continental block, it reinvigorates its retreating course, leading to a second trench retreat phase.

The numerical model setup is, by definition, a simplified and more generic version of reality. As such, the shape and geometry of the continental block are not designed to specifically reproduce the natural system as in the central Mediterranean. However, these models may serve to help us understand the dynamics of this complex subduction zone. There, the first phase of backarc extension ends when continental lithosphere enters the trench on both sides of the oceanic slab while the width of the oceanic domain narrows (Fig. 1). The onset of the second phase of extension occurs only a few million years after the first pulse. The first slab window forms ca. 12–10 Ma in northern Africa, and propagates laterally westward beneath Sicily until the middle Pleistocene, opening a 500-km-wide slab window. In the middle Pleistocene, a subduction window forms beneath the Central Apennines. Our model indicates that the opening of those slab windows is a necessary condition for the second phase of extension in the Tyrrhenian that lasted for ∼5–7 m.y. The opening of the central Mediterranean is somewhat different from the simplified, perfectly symmetric numerical model presented here, but the rates, time scales, and dimensions of the model dynamics are quite similar to the natural example, indicating that the model represents the primary features of subduction forces, trench retreat, and slab deformation and tearing in the central Mediterranean. We furthermore propose that the presented model can, more generally, fit several other natural cases of backarc spreading associated with lateral variations along trench (such as the eastern Mediterranean and western Pacific Ocean). As in the case presented here, in these regions fast backarc extension is accompanied by large-scale rotation, caused by lateral collision of positively buoyant indentors, such as continental plates or plateaus (e.g., Tonga, Aegean, Papua New Guinea, New Zealand, the Mariana Islands, Vanuatu, and Fiji), ridges (e.g., the Ryukyu Islands), and seamount chains (e.g., Tonga) (Wallace et al., 2009, and references therein).

We thank the editor J.B. Murphy and reviewers T. Gerya, W. Spakman, and S. Ellis for their constructive comments. This research was supported by the European Young Investigators (EURYI) Awards Scheme (EUROHORCS/ESF, including funds from the National Research Council of Italy) and by the European Research Council (ERC Starting Grant 279828). Models presented in this paper have been realized thanks to CASPUR HPC Standard Grant 2012.

1GSA Data Repository item 2014181, details on the numerical method, energy computation, and distribution, is available online at, or on request from or Documents Secretary, GSA, P.O. Box 9140, Boulder, CO 80301, USA.