It has been suggested that post-orogenic lithospheric removal in continental back arcs promotes extension and surface subsidence. However, the surface response of this process and its primary difference from “classical” back-arc opening have remained uncertain. Here, the back-arc extension process with varying continental mantle lithosphere thickness and thermal heterogeneities is studied by using thermomechanical subduction experiments. The experiments illustrate that models with only slab retreat result in minor surface subsidence and extension in the back-arc region. Alternatively, there is notable extension due to the slab retreat and a localized high-temperature zone in the back arc with uniform lithospheric thickness. Models with advecting mantle (after lithospheric removal) in the extending back arc predict rifting (stretching factor β > 2) and surface subsidence (>1.5 km) in the center of the basin. The results of this work suggest that lithospheric removal may be an important trigger for continental back-arc development rather than slab retreat alone causing lithospheric extension and subsidence. The findings help explain rift formation and subsidence in the Aegean Sea–west Anatolia, and possibly other Mediterranean back arcs, such as the Alboran Sea and the Pannonian Basin.


Retreating ocean-continent subduction systems are associated with thin and hot continental back-arc lithosphere of the overriding plate. It has been recognized that the thinning of the back arc occurs due to the dynamic slab pull exerted in the crust by the oceanic plate that causes trench retreat, surface extension, sedimentary basin formation, and high heat flow (Le Pichon et al., 1981; Faccenna et al., 1996). However, the geodynamics of this extension and the primary factor for the uplift and/or subsidence in the back arcs are not well understood—for example, the dominant role of induced mantle convection preceding rifting (Toksöz and Hsui, 1978) versus vertical stresses due to the slab-pull forcing that produces dynamic subsidence (Mitrovica et al., 1989; Husson, 2006). In addition to this, there is still controversy about why some back-arc regions are anomalously topographically high, such as the Canadian Cordillera (Hyndman and Currie, 2011), or tectonic subsidence (1.5–2.5 km) has occurred in the Mediterranean back arcs during the past 15–20 m.y. (e.g., Aegean Sea, Pannonian Basin, Alboran Basin) (Le Pichon et al., 1981; Royden et al., 1983; Watts et al., 1993). In this work, a new class of numerical models is presented to account for the back-arc basin subsidence.

Many geological observations suggest that back-arc extension may have been initiated or promoted by lithospheric thinning after post-orogenic thickening. Two major geodynamic hypothesis have been suggested as lithospheric removal mechanisms (Göğüs and Pysklywec, 2008): (1) mantle lithosphere delamination (Bird, 1979); and (2) Rayleigh-Taylor type convective removal (Houseman et al., 1981). Conceptual removal models have been put forward to explain rapid subsidence, anomalous heating, exhumation of high-pressure rocks, tectonic mode switching (from compression to extension), and high stretching factors at the Pannonian Basin (Horvarth, 1993; Houseman and Gemmer, 2007), the Tyrrhenian Sea (Channell and Mareschal, 1989; Faccenna et al., 1996), the Alboran Basin (Docherty and Banda, 1995; Platt, 2007), and Aegean Sea–west Anatolia (Dewey, 1988; Seyitoğlu and Scott, 1996). Accordingly, all of these back arcs contain a Neotethyan suture that marks the Alpine continental collision and the shortened, thickened lithosphere.

In spite of the potential geodynamic feedback between the lithospheric removal and slab retreat–roll-back driven back-arc extension, proposed models (references given above) have addressed each of these mechanisms in isolation. In this work, thermomechanical numerical models are used to investigate the surface response of hot mantle upwelling (after lithospheric removal) in the retreating subduction zone as a coupled geodynamic process to explain the subsidence and rifting in continental back arcs. A series of computational geodynamic experiments is conducted to test the variation of surface topography and crustal thickness for different back-arc lithospheric thicknesses and thermal fields. The experiments are carried out in two-dimensional procedures for high-resolution purposes, while model predictions are carefully interpreted in the context of three-dimensional geodynamic processes.

As a case study, the numerical results are interpreted in the context of the last 20 m.y. of geodynamic evolution of the Aegean Sea–west Anatolia back arc, where the cause of extension is still not well understood (Figs. 1A and 1B). Although the extension is mainly controlled by the southward retreat of the Aegean slab (e.g., Meulenkamp et al., 1988; Okay and Satır, 2000), there is also strong evidence for lithospheric removal (Gessner et al., 2013) (Figs. 1C and 1D). For instance, seismological studies that suggest extensive low-speed anomalies beneath the Aegean and western Anatolia also support the hypothesis for lithospheric removal and subsequent asthenospheric upwelling (Wortel and Spakman, 2000; Piromallo and Morelli, 2003). Corroborating the geological and seismological interpretations, synthesized petrological work by Pe-Piper and Piper (2006) indicates that the widespread distribution of igneous rocks in the Aegean originated through decompression melting of upwelling asthenospheric mantle. While these asthenosphere-derived rocks are commonly well correlated with areas of low-speed seismic waves (warm), their ascendant commonly postdates the subduction (arc-related) magmatism. Such change in the magma source (subduction to asthenospheric mantle) has also been interpreted through geochemical analysis of volcanic rocks in western Anatolia where mantle lithosphere delamination has been suggested (Aldanmaz et al., 2000).


The configuration of the numerical models is designed as a general representation of the retreating subduction of the oceanic lithosphere (pro-plate) under the continental plate (retro-plate) at the plate boundary region (see the GSA Data Repository1 for numerical method, reference model setup, and material properties, including Figure DR1 and Table DR1). This configuration provides a reasonable physical approximation for the inception of slab retreat–roll-back (e.g., Hellenic slab retreat) where the subducting slab has reached appropriate density and length scale tens of millions of years after the initiation of subduction (Meulenkamp et al., 1988). It is worthwhile to note that the purpose of this work is not to explore the origin and/or cause of instabilities and model exactly the geodynamic evolution of continental back arcs, such as western Anatolia–Aegean, but rather the focus is testing the surface and crustal response of the previously foundered lithosphere (during orogenic collapse) in the retreating subduction systems. Three selected numerical experiments are presented—from simple back-arc extension to a more complex setup—with varying thermal heterogenities and thicknesses of the back-arc lithosphere.

Figures 2A and 2B show the evolution of oceanic slab retreat and resulting back-arc extension experiment EXP-1 with uniform temperature field and thickness along the continental mantle lithosphere (i.e., without any thermal or structural perturbation introduced). By time t = 7.6 m.y., the oceanic lithosphere (pro-plate) is subducting under the continental lithosphere at the trench boundary with a relatively steep angle, and the trench has retreated 60 km (with respect to t = 0). At this stage, the amount of subduction is significant; however, there is only ∼4 km of crustal thinning and extension, distributed evenly along the subduction back arc. The surface topography remains almost uniform (2 km) along the entire back arc except for a sudden drop near the subduction trench.

Further into the experiment, by t = 19.3 m.y., the oceanic slab hits the bottom of the solution box while dragging down portions of continental mantle lithosphere. This type of double-sided subduction, “ablative subduction,” is akin to that proposed by Tao and O’Connell (1992) as an alternative model to one-sided subduction. The total retreat of the trench is ∼130 km, and the crust is thinned by ∼10 km uniformly to the trench. In addition to this, thermal subsidence involved with this model is small (change in Moho temperature is negligible), and topographic subsidence of 500 m occurs with respect to the model edge due to the response of a few kilometers of thinned crust closer to the arc. A noteworthy prediction of the model is that although the subduction retreat process alone can stretch and thin the crust, the effect of this process is considerably small for typical back-arc extension. The first-order results indicate that deformation (thinning and subsidence) does not localize if the back-arc lithosphere initially has a uniform temperature field, thickness, and rheological properties.

Alternatively, a series of experiments is conducted with a pre-existing localized high-temperature (1375 °C) zone in the center of the back arc at 100 km depth (200 km wide) and uniform lithospheric thickness (EXP-2). The purpose of running these experiments is to evaluate the surface response of localized extension in the back arc through an injection of a heat and/or magma source behind the volcanic arc, for instance, upwelling of partially molten mantle. The geodynamic evolution of the retreating subduction is comparable with that of the previous experiment; however, by t = 7.6 m.y., overall extension is more amplified with ∼10 km of crustal thinning (Figs. 3A and 3B). This is because localized heating in the back arc would effectively weaken the lithosphere and induce dynamic stretching. At this time also, the surface has subsided nearly 800 m with respect to early stages of the model (e.g., t = 1.2 m.y.), but a broader region of surface high (∼2.3 km) still persists with respect to the left edge of the model, again owing to the preexisting thermal uplift.

As the trench retreat progresses with 160 km of total retreat, by t = 19.3 m.y., the extension and lithospheric thinning become more significant in the central back-arc lithosphere. At the same locus of the total crustal thinning of ∼20 km, the surface basin is depressed ∼700 m with respect to the left edge of the solution space. There is also an ∼200 m negative deflection at the left edge of the box as well as ∼5 km of crustal thinning. Clearly, the extension is higher in magnitude and less broadly distributed in this experiment compared with the previous experiment.

The reference model (EXP-3) shown in Figure 4A presents three stages during the model evolution. The initial model setup in the back arc includes a high-temperature anomaly (i.e., 1375 °C is computed at a depth of 100 km), and 75-km-thick, 200-km-wide mantle lithosphere is substituted by sub-lithospheric mantle. This high-temperature zone is an approximation of the replacement of cold mantle lithosphere with hot sub-lithospheric mantle that advects heat at shallow lithospheric levels after the inferred lithospheric removal and during collapse of the orogen. Such a removal event likely occurred prior to or contemporaneous with the slab roll-back–retreat along the Alpine collision zone under the continental back arc in the Mediterranean (Dewey, 1988).

At time t = 1.2 m.y., the oceanic lithosphere is in the early stage of slab retreat under the continental lithosphere at the trench boundary. The hot back-arc area (foundered mantle lithosphere) (x = 400–800 km) is represented by a peak surface topography of ∼3.2 km elevation (uplift due to the heating) and a maximum crustal thickness of 40 km (Fig. 4B). The initial Moho temperature is 1100 °C.

By t = 7.6 m.y., the amount of subduction has not changed much with respect to the previous experiments for the same time, but the subduction trench has retreated ∼100 km compared to t = 0. Deformation in the Lagrangian grids suggests that the crust is thinned significantly at the hot back-arc zone. At this location, the crust has thinned by ∼16 km and the surface has subsided 1.4 km since t = 0. This extension is induced by mechanical and dynamic stretching working in tandem at the center of the lithospheric gap: (1) the retreating trench, and (2) the flow of hot and buoyant asthenospheric mantle softening the crust.

By t = 19.3 m.y., the trench has retreated ∼220 km, and the subducting slab is fully coupled with the overriding plate while both are consumed into the sub-lithospheric mantle by way of ablative subduction. At this stage, dynamic stretching is associated with high extension factors (e.g., β > 2, rift basins) where model predictions suggest ∼600 km of crustal extension and 30 km of crustal thinning. The surface topography is characterized by 1.5 km of negative surface deflection in the center of the “rift” with respect to the flanking high topography that persists (∼2 km) with thicker crust of ∼44 km to the edge of the solution box, where further deformation is negligible (Figs. 4A and 4B).


The results of the reference experiment (EXP-3) are consistent with interpretations of geological and geophysical observations in the Aegean–western Anatolia and may generally be applicable to other Mediterranean-type back-arc basins:

  • (1) Seismic tomography interpretations that suggest a zone of low-speed seismic waves at shallow depths under the Aegean Sea–west Anatolia (Wortel and Spakman, 2000; Piromallo and Morelli, 2003) and petrological explanations that suggest transition from subduction-related magmas (calc-alkaline) to the igneous rocks of asthenospheric mantle origin (alkaline chemistry) (Aldanmaz et al., 2000; Pe-Piper and Piper, 2006) are in good agreement with the modeled asthenospheric upwelling under the continental back arc.

  • (2) Estimated 600 km of extension in the Aegean based on the geological reconstructions (Jolivet and Brun, 2010) and geodetic estimates of 2 cm/yr extension (Aktug et al., 2009) for western Anatolia are consistent with the model’s (EXP 3) stretching predictions.

  • (3) At first order, the 2 km of subsidence in the Sea of Crete (Le Pichon et al., 1981) and the pattern of increasing surface elevation to western Anatolia (Menderes massif) along a NNE-SSW cross section (Y-Y′ in Fig. 1B) spatially correlate well with the modeled (t = 19.3 m.y.) surface topography from maximum subsidence to flanking topography, respectively.

The results of this work may provide insight into the geodynamic evolution of other Mediterranean back arcs where lithospheric “dripping” and/or delamination have been proposed and the causes and/or timing of the extension have still been disputed. For example, seismological studies in the Pannonian Basin interpret the presence of warm mantle beneath this basin to be the result of asthenospheric rise (Wortel and Spakman, 2000). According to Harangi et al. (2006), geochemical analysis of igneous rocks also suggests upwelling mantle in the various parts of the basin (e.g., Carpathians) where mantle lithosphere most likely has been removed. Moreover, the basin has been extending with crustal stretching factors of ∼2 and has subsided 2.5 km in the last 17 m.y. (Royden et al., 1983). The reference model predictions (EXP-3) with β > 2 can account for the degree of crustal thinning and subsidence inferred by Royden et al. (1983).

The Alboran Basin is another example of a very thin (60 km) and warm lithosphere where a young ocean basin starts to develop in conjunction with 1.5 km of predicted subsidence (Watts et al., 1993; Platt, 2007). The amount of extension and subduction trench retreat predicted by the reference model (∼220 km total retreat) by t = 19.3 m.y. is comparable with the estimated ∼200 km western retreat of the Betic-Rif arc of the Alboran Basin in the last 20 m.y. (Faccenna et al., 2014).

The investigation of appropriate rheological and thermal properties for the initiation of mantle lithosphere removal in continental back arcs is not the focus of this contribution. Such model parameterizations have been made in other geodynamic experiments. For instance, Currie et al. (2008) suggested that subduction-induced corner flow can also drive small-scale gravitational instabilities under the continental back arc. However, the surface response of these lithospheric instabilities may not be comparable to the results given in this work, because the observed size of instabilities under the back arc were not sufficiently large to produce a wide region of thinned mantle lithosphere in the time scales during which such a process has been proposed for a large number of areas given in this work (e.g., Mediterranean back arcs).


This work’s numerical experiments on subduction retreat demonstrate that a uniform temperature field and mantle lithosphere thickness in the continental back arc (without any perturbation) are associated with minor crustal extension and subsidence (experiment EXP-1). In alternate slab retreat scenarios (EXP-2), with a localized high-temperature anomaly in the back arc and uniform lithospheric thickness, the extension and subsidence are amplified and are localized in the center of the back arc.

Post-orogenic lithospheric removal in continental back arcs may promote rifting (β > 2) and subsidence (>1.5 km) in the extensional back-arc area (EXP-3). The models predict that rifting induced by advecting mantle flow on both sides of the back arc can stretch the lithosphere significantly. The rift center is associated with maximum subsidence and thinning; however, on the rift flanks, the surface depression and extension are less pronounced because of thicker crust. It is reasonable to interpret that, depending on the thickness and pre-rifting history of the lithosphere as well as its unique geological evolution, back-arc tectonics are more complex than the simplified idea that the slab retreat alone can cause a high degree of extension and subsidence. The primary finding of this work is that any type of lithospheric removal can act as a trigger and facilitator for rift development in continental back arcs.

I acknowledge Russell Pysklywec for assisting with the SOPALE modeling code, which was originally developed by Phillip Fullsack at Dalhousie University (Nova Scotia, Canada) under the direction of Chris Beaumont. The paper was improved by constructive reviews from Klauss Gessner, Laurent Husson, and an anonymous reviewer. I thank Hasan Sözbilir, Tadashi Yamasaki, Celâl Şengör, Greg Houseman, and Claudio Faccenna for discussions and the motivation to develop the ideas presented in this work.

1GSA Data Repository item 2015014, numerical method, material properties, and Figure DR1 (model setup), is available online at www.geosociety.org/pubs/ft2015.htm, or on request from editing@geosociety.org or Documents Secretary, GSA, P.O. Box 9140, Boulder, CO 80301, USA.