Abstract
Lithospheric slab tearing, the process by which a subducted lithospheric plate is torn apart and sinks into the Earth’s mantle, has been proposed as a cause for surface vertical motions in excess of 100 s of meters. However, little is known about the mechanisms that help initiate and control the propagation of slab tearing and the associated uplift. This study aims to explore these processes by means of 3D thermo-mechanical geodynamic modelling of a slab retreat oblique to a continental margin, using the Gibraltar Arc region (Betic Cordillera) as a scenario for inspiration. Our results suggest that the obliquity of the continental passive margin relative to the subduction trench leads to an asymmetric distribution of subduction forces and strength, facilitating the initiation of slab tearing. The model results predict a lateral migration of the tearing point at a velocity ranging between 38 and 68 cm/yr for a sublithospheric-mantle viscosity of up to 1e+22 Pa s. This fast slab tearing propagation yields uplift rates of 0.23–2.16 mm/yr above the areas where the subducted slab is torn apart, depending on mantle viscosity. Although a more detailed parametric exploration is needed, this range of uplift rates is compatible with the uplift rates required to overcome seaway erosion along the Atlantic-Mediterranean marine corridors during the Late Miocene, as proposed for the onset of the Messinian Salinity Crisis.
Highlights
• A subduction trench arriving diachronously to a continental passive margin facilitates slab tearing of the retreating slab.
• The viscosity of the sublithospheric mantle and the amount of shortening prior to tearing are the key controls on the slab-tearing dynamics.
• Surface uplift rates of 0.23–2.16 mm/yr for the overriding plate are obtained, compatible with those required for blocking the seaways of the Mediterranean Sea as during the onset of Messinian Salinity Crisis.
1 Introduction
The perception that some areas of the continental crust have subsided/uplifted at rates that cannot be explained by crustal thickening or fault activity alone, has pointed the need to other epeirogenic surface uplift mechanisms (England and Molnar, 1990) including dynamic topography or isostatic responses to mass motions. Slab breakoff is among the deep-seated mechanisms invoked to justify the long-wavelength, high rates of surface uplift observed in some continental collision settings (Spakman et al., 1988; Wortel and Spakman, 1992; 2000; Davies and von Blanckenburg, 1995; Van Der Meulen et al., 1998). Uplift in these scenarios is purportedly driven by the cancellation of the negative buoyancy force, the same force that drives slab pull and subduction (e.g., Garcia-Castellanos et al., 2000; Conrad and Lithgow-Bertelloni, 2002; Billen, 2008; Boonma et al., 2019; Jiménez-Munt et al., 2019). Slab breakoff is a process happening at depth within the mantle consisting of the detachment of a subducted oceanic lithospheric slab from the more buoyant continental lithosphere during continental collision. The concept of slab breakoff, as inferred from seismic tomography, was first proposed for the geodynamical evolution of the Mediterranean by Spakman et al. (1988) and Wortel and Spakman (1992). Slab breakoff was then used to explain post-collisional magmatism and exhumation of high-pressure rocks in the European Alps by Davies and von Blanckenburg (1995). Garzanti et al. (2018), and references therein, gave a comprehensive global overview of where slab breakoff has been invoked to explain changes in plate kinematics and tectonic deformation in regions including the Alps (Spakman et al., 1988; Wortel and Spakman, 1992; Davies and von Blanckenburg, 1995; Sinclair, 1997; Fox et al., 2015), the Mediterranean region (Carminati et al., 1998; Wortel and Spakman, 2000; Spakman and Wortel, 2004; Rosenbaum et al., 2008; Chertova et al., 2014; van Hinsbergen et al., 2014; Spakman et al., 2018), the Anatolia-Zagros orogen (Şengör et al., 2003; Faccenna et al., 2006), and Himalaya and Tibet (Jiménez-Munt and Platt, 2006; Van Hinsbergen et al., 2014; Wu et al., 2014; Liang et al., 2016; Qayyum et al., 2022). Many of these studies confirm the hypothesis that a complete slab break off is preceded by the lateral propagation of a slab tear or detachment within the subducted plate (Yoshioka and Wortel, 1995), ascribing short-lived, long-wavelength exhumation events or sudden pulses in sediment supply to this process. This lateral propagation of tearing makes the process of breakoff more gradual and prone to be influenced by the 3D tectonic configuration. The lateral tearing of lithospheric subducted slabs has been invoked to explain a diversity of observations, such as the seismicity pattern in the Vrancea slab (Mitrofan et al., 2016) or the seismic tomography in the South Iberian margin, Apennines and Hellenic arc (Wortel and Spakman, 2000). Previous 3D numerical models of subduction that exhibited slab tearing, such as those by van Hunen and Allen (2011); Duretz et al. (2014); Magni et al. (2017), all investigating a passive margin parallel to the subduction trench, which ended up producing slab tearing in the slab interior rather than on the slab edge. How does this initial tectonic setting condition the slab tearing in the first place? How can other configurations affect the progress of the tearing? How much does slab tearing contribute to the buoyancy-driven isostatic surface uplift? Our study addresses these questions inspired by the geological region of the Gibraltar Arc in the westernmost Mediterranean.
Based on seismic tomographic imaging, Spakman and Wortel (2004) suggested that slab tearing might have occurred in the Gibraltar Arc region as a consequence of the continental convergence between Africa and Iberia and the subsequent slab rollback (Elsasser, 1971; Karig, 1971) declined during early Miocene. The uplift observed in the internal Betic basins after late Tortonian are best constrained from the present elevation of tectonically undeformed Miocene marine sediment in that region (marine to non-marine transition), often above 600 m elevation (Garcés et al., 1998; Iribarren et al., 2009). The biostratigraphic studies by Krijgsman et al. (2018) and van der Schee et al. (2018) revised the age of the uplift and provided a new age constraint on the western Betic intramountain basins to be older than 7.51 Ma (late Tortonian). This has been interpreted as the result of a westward migration of a lateral tear of the hanging slab seen in tomography underneath the Gibraltar Arc (Spakman and Wortel, 2004; Garcia-Castellanos and Villaseñor, 2011; Bezada et al., 2013). Slab tearing occurred in previous thermomechanical models by Chertova et al. (2014) and Spakman et al. (2018), but the cause and timing of slab tearing was not investigated by them. The timing and the topographic response to this tearing process remains poorly constrained.
The uplift of the intramountain basins within the Betics and Rif (Figure 1), possibly caused by the slab tearing, has been linked to the closure of the marine gateway across the Gibraltar Arc during the Late Miocene which led to the partial desiccation of the Mediterranean Sea, (Messinian Salinity Crisis; MSC) (Garcia-Castellanos and Villaseñor, 2011; Coulson et al., 2019; Capella et al., 2020). According to this hypothesis, once the marine gateways closed off due to the regional epeirogenic uplift, it triggered the massive salt accumulation of the MSC (5.96–5.33 Ma), possibly the most abrupt environmental change on Earth since the beginning of the Tertiary.
FIGURE 1
This study aims to understand the lithospheric slab tearing process and the resulting surface uplift using 3D thermomechanical modelling. For this purpose, we design a model setup of convergence between two tectonic plates, one of them being a continental margin attached to an oceanic plate which subducts under a second plate (Figure 2, we try both continent and oceanic structure for this one). This model setup is inspired by the tectonic setting of the westernmost Mediterranean. The tectonic kinematics of this region for the last 35 Myr are controversial (Faccenna et al., 2004; Spakman and Wortel, 2004; Vergés and Fernàndez, 2012) in terms of subduction polarity, the direction of slab retreat, and the exact timing. To study the tearing process, our model setup (Figure 2C) follows those tectonic models implying a subducted slab oblique to the Iberian margin and laterally ending towards the east (Spakman and Wortel, 2004; Vergés and Fernàndez, 2012; Van Hinsbergen et al., 2014; Angrand and Mouthereau, 2021). It is important to note that our study will focus on the tearing process itself during Middle to Late Miocene rather than on the previous evolution of the slab retreat. We address this by testing how different model parameters act on the initiation of the slab tearing and its propagation along the trench and the resulting surface elevation changes (due to viscous flow, temperature evolution, and dynamic topography). In addition, we will provide insights into how the slab tearing dynamics fit within the realm of the western Mediterranean, and how these surface vertical motions may have helped restricting the connections with the Atlantic Ocean during the Messinian Salinity Crisis.
FIGURE 2
2 Methods
The modelling in this work is carried out using a 3D thermo-mechanical coupled numerical code, “I3ELVIS” (Gerya, 2013). The code is based on finite-differences and marker-in-cell numerical schemes (Harlow and Welch, 1965; Gerya and Yuen, 2003).
2.1 Governing equations
The governing physical laws such as conservation of mass, conservation of momentum, and heat equation are discretised on a staggered Eulerian grid.
(1) The conservation of mass is described by the continuity equation of an incompressible viscous medium:where is the flow velocity and where i (as well as j shown hereafter) is index, which denotes spatial directions i, j = (x,y, z) in 3D using Einstein notation.
(2) The conservation of momentum is described by the Stokes equation for creeping flow, which states:where is the deviatoric stress tensor, denotes the total pressure (Pa), and is the gravitational acceleration (m/s2). In this form, we consider the extended Boussinesq approximation, where the density in the buoyancy term varies locally as a function of temperature (T), pressure (P), and composition (c).
(3) The heat balance in a convective medium is described by the heat conservation equation, which states:
where is the heat capacity (J/K), is temperature (K), denotes thermal conductivity (W/(m·K), which is a function of temperature, pressure, and composition, i.e., k(T, P, c), is radioactive heat production; is adiabatic heat production/consumption; is latent heat production/consumption (due to phase transformation), and is the shear heat production (a product of deviatoric stress and strain rate). More details are provided in Gerya and Yuen (2003), Gerya and Yuen (2007).
For each time step, a fourth order Runge-Kutta scheme is used to spatially advect the markers. The multi-grid method is used to speed up the convergence of the Gauss-Seidel iterative solver.
2.2 Density model and phase changes
The rocks’ densities vary with temperature T (K) and pressure P (Pa) according to the equation of state:where is the reference density at = 1 MPa and = 298.15 K, the coefficient of thermal expansion = 2×10−5 1/K, and the coefficient of compressibility = 6×10−12 1/Pa. Our models take phase transition of olivine in the mantle into account. As the dry olivine is subjected to greater pressure at depths, it first undergoes exothermic phase transition (∼410 km) and transforms into wadsleyite (Katsura and Ito, 1989). At a greater depth and pressure, the wadsleyite exothermically transforms into ringwoodite (∼520 km), which decompose (endothermically) into bridgmanite (silicate perovskite) at an even greater depth (∼660 km) (Ito et al., 1990). The eclogitization of the subducted oceanic crust (basaltic and gabbroic) is taken into account by linearly increasing the density of the crust with pressure from 0% to 16% in the P-T region between the experimentally determined garnet-in and plagioclase-out phase transitions in basalt (Ito and Kennedy, 1971).
2.3 Visco-plastic rheology
A composite visco-plastic rheology is used that captures the general dynamics of continental convergence and the underlying mantle flow (
Gerya and Yuen, 2003;
2007). The ductile rheology is charactesized by a ductile viscosity
that is approximated by a combination of effective viscosities for diffusion
and dislocation
creep. We assume that the total strain rate is the sum of the strain rates due to diffusion creep and dislocation creep, and that they act in parallel and independently from each other.
(1) In the crust, we assume constant grain size and and is computed as:
where
Ris gas constant (8.314 J/(K·mol)),
Pis pressure (Pa),
Tis temperature (K),
is the second invariant deviatoric strain-rate tensor,
σcris the critical stess,
Ais the pre-exponential factor (Pa
n·s),
Edenotes activation energy (J/mol),
Vis activation volume (J/Pa), and
nis the stress exponent of the viscous creep (
Table 1).
(2) In the mantle, the ductile creep is implemented with grain size growth and reduction processes. In the case of the mantle, the composite rheology in Eq. 5 still stands,
where,
is the experimentally determined pre-exponential factor for diffusion creep (Pa·s) and
is pre-exponential factor for dislocation creep (Pa
n·s),
his grain size (m),
mis the grain size exponent (
Table 1).
TABLE 1
| Material | Density ρ0 (kg/m3) | Thermal conductivity (W/m·K) at TK, PMPa | Flow Law | Flow law parameters |
|---|---|---|---|---|
| Upper continental crust (Felsic) | 2750 | 0.64 + 807/(T + 77) | Wet quartzite Ranalli (1995) | A = 1.97 × 1017 Pans, n = 2.3, E = 1.54 × 105 J/mol, V = 0 J/(mol·MPa), σcr = 3 × 104 Pa, C = 3 MPa, μ0 = 0.3, μ1 = 0 |
| Lower continental crust (Gabbro) | 3000 | 1.18 + 474/(T + 77) | Plagioclase An75 Ranalli (1995) | A = 4.80 × 1022 Pans, n = 3.2, E = 2.38 × 105 J/mol, V = 0 J/(mol·MPa), σcr = 3 × 104 Pa, C = 3 MPa, μ0 = 0.3, μ1 = 0 |
| Upper oceanic crust (Basalt) | Wet quartzite Ranalli (1995) | A = 4.80 × 1022 Pans, n = 3.2, E = 2.38 × 105 J/mol, V = 0 J/(mol·MPa), σcr = 3 × 104 Pa, C = 3 MPa, μ0 = 0, μ1 = 0 | ||
| Lower oceanic crust (Gabbro) | Plagioclase An75 Ranalli (1995) | A = 4.80 × 1022 Pans, n = 3.2, E = 2.38 × 105 J/mol, V = 0 J/(mol·MPa), σcr = 3 × 104 Pa, C = 3 MPa, μ0 = 0.6, μ1 = 0.0 Mod4: μ1 = 0.3 | ||
| Mantle | 3300 | 0.73 + 1293/(T+77) × exp (0.000004P) | Dry olivine Hirth and Kohlstedt (2003) | m = 3, Adiff = 1.50 × 1015 Pa s, Ediff = 3.75 × 105 J/mol, Vdiff = 0.7 J/(mol·MPa), Adisl = 1.10 × 1016 Pans, n = 3.5, Edisl = 5.30 × 105 J/mol, Vdisl = 2.6 J/(mol·MPa), C = 3 MPa, μ0 = 0.6, μ1 = 0.0 Mod3: Vdiff = 0.8 J/(mol·MPa), Vdisl = 3.0 J/(mol·MPa) Mod4: μ1 = 0.3 |
| Mantle weak zone | m = 3, Adiff = 1.50 × 1015 Pa s, Ediff = 3.75 × 105 J/mol, Vdiff = 0.7 J/(mol·MPa), Adisl = 1.10 × 1016 Pans, n = 3.5, Edisl = 5.30 × 105 J/mol, Vdisl = 2.6 J/(mol·MPa), C = 3 MPa, μ0 = 0, μ1 = 0 |
Material properties used in the numerical experiments. The flow law include: A is the pre-exponential factor; E denotes activation energy; V is activation volume; n is the stress exponent; m is grain size exponent; σcr is critical stress or the assumed diffusion-dislocation transition stress; C is the rock compressive strength at p=0 MPa; μ0 and μ1 are the initial and final internal coefficients, respectively parameters (Karato and Wu, 1993; Ranalli, 1995; Hirth and Kohlstedt, 2003; Turcotte and Schubert, 2014). The subscripts “diff” and “disl” indicate that those parameters are associated with diffusion and dislocation creep processes, respectively. Mod3 has higher values for mantle activation volume than other models. Mod4 has higher final internal friction coefficient for lower oceanic crust and the mantle. Other properties for all rock types include: heat capacity Cp = 1000 J/(kg·K), thermal expansion α = 2 × 10−5 1/K; and compressibility β = 6 × 10−12 1/Pa.
The viscous rheology is combined with a brittle rheology to compute an effective visco-plastic behavior adopting a Drucker-Prager yielding criterion (e.g., Ranalli, 1995) that effectively accounts for a strain weakening:where is the internal friction coefficient (; are the initial and final internal friction coefficient, respectively), is the rate of faults weakening with integrated plastic strain ( is the upper strain limit for the fracture-related weakening), C is the rock compressive strength at p = 0, t is time (s), is the plastic strain rate tensor (Table 1).
2.4 Model setup
The continental convergence is modelled with an incoming continental block, overriding a subducting oceanic plate, which is connected to a stationary continental block through a passive margin (model box located in Figure 1B). The model setup is based on geodynamic settings where an oceanic subducting plate retreats towards an oblique continental passive margin, inspired on the Neogene-to-present southern Iberian margin (e.g., Spakman and Wortel, 2004; Vergés and Fernàndez, 2012; van Hinsbergen et al., 2014). The 3D model domain (Figure 2) measures 1,500 km × 780 km × 1,200 km, with a resolution of 4.6 km× 3.0 km × 4.6 km, in the x, y (vertical), and z directions, respectively (22.1 million elements, 6 markers per element). The 40-km thick continental crust consists of the upper (20 km) and lower (20 km) continental crust, and thinning toward the ocean. The 8-km thick oceanic crust also consists of the upper (basaltic, 3 km) and lower (gabbroic, 5 km) oceanic crust. Partial melting and melt extraction processes are neglected in order to keep our models simplified and the interaction between the slab tearing and the crust isolated.
We use the “wet quartzite” viscous flow law for both the upper continental crust and the upper oceanic crust, while the “Plagioclase An75” describes the behaviour of the lower continental crust and lower oceanic crust (Table 1; Ranalli, 1995). The lithospheric mantle, asthenospheric mantle, and mantle “weak zone” (less viscous mantle material) all have a “dry olivine” rheology (Table 1; Hirth and Kohlstedt, 2003). The weak zone is prescribed in front of the incoming continental plate (Figures 2A,C) to facilitate the subduction initiation (Table 1). Two weak zones are also prescribed on both the eastern and western side of the incoming continental plate as a way to decouple the moving plate from the surrounding oceanic domain (Figures 2A,C).
In the first stage of the experiment, the convergence velocity is prescribed as constant at a vertical plane at x=1,389 km between depths of y=21–147 km (upper lithospheric mantle). This plane is labelled “ridge” in Figure 2C. The first stage of the experiment involves a forced convergence until the slab reach 200 km depth. After the first stage, the obtained thermo-mechanical state is then used as an initial setup for the second stage. In this second stage, the prescribed convergence rate is either removed, so that the slab sinks due to its own weight (Mod1-reference, Mod2, Mod3, and Mod4), or reduced to a lower value of 4 mm/yr (still pushing from the “ridge” at x=1,389 km) resembling the relative convergence rate between Africa and Iberia in the western Mediterranean (Mod5) (Macchiavelli et al., 2017).
2.5 Boundary conditions
The velocity boundary conditions are free slip for the left (x=0 km), right (x= 1,500 km), front (z=0 km), back (z=1,200 km), and top boundaries (y=0 km). The bottom boundary (y=780 km) is open/permeable, allowing rock materials to flow in or out of the model domain. For the bottom boundary, an external outflux boundary implies zero shear stress conditions and constant normal velocity to be satisfied at ∼300 km below the base of the model domain, which ensures the mass conservation within the computational domain (e.g., Burg and Gerya, 2005; Gerya et al., 2008; Li et al., 2013).
The elevation of the lithosphere is calculated dynamically as an internal free surface through a buffer layer of “sticky air” (ηair = 1018 Pa s, ρair = 1 kg/m3) (Gerya and Yuen, 2003; Schmeling et al., 2008; Crameri et al., 2012). It is 22 km thick on top of the continental plate and 25 km on top of the oceanic plate. We implemented a simplified erosion condition in our model, where instantaneous sedimentation limits a trench depth to 8 km below the water level and the instantaneous erosion is prescribed at 8 km above the initial continental crustal surface where rock markers change into sticky-air markers.
The continental geotherm is prescribed as a linear variation from 0°C at the model surface (y≤22 km, air) to 1,344°C the lithosphere-asthenosphere boundary (y=110 km) (Jiménez-Munt et al., 2019; Kumar et al., 2021). The initial thermal structure of the oceanic lithosphere is calculated using the half-space cooling model (e.g., Turcotte and Schubert, 2014) based on a slab age of 110 Ma. The initial adiabatic temperature gradient of 0.5°C/km is prescribed in the asthenospheric mantle (Katsura et al., 2010). For the initial thermal boundary conditions, the upper boundary has a fixed value of 0°C, and zero horizontal heat flux across the vertical boundaries. Similar to the velocity boundary, the bottom thermal boundary is permeable such that the temperature and vertical heat fluxes can vary along the lower boundary. This implies that the constant temperature condition is satisfied at ∼300 km below the bottom of the model box (e.g., Burg and Gerya, 2005; Gerya et al., 2008; Li et al., 2013).
3 Results
All of the experiments incorporate two stages. The first stage consists of forced convergence until the slab reaches 200 km depth. The time when this happens is defined as t=0 in our models. The second stage (t>0) uses the output from the first stage as an initial setup and continues on as the slab undergoes either free-rollback (Mod1-reference, Mod2, Mod3, and Mod4) or fixed convergence velocity (Mod5). Note that time “t” (in Myr) is relative to the beginning of stage 2 for each model.
3.1 Reference model (Mod1-reference)
After the slab has reached the depth of 200 km, the prescribed convergence rate is stopped. As the dense slab continues to sink due to its own negative buoyancy relative to the surrounding mantle (Table 1), the continental margin (Iberia) starts to bend downward and the incoming block overrides the passive margin. At t=3.84–4.10 Myr, the lithospheric thinning/necking started on the slab’s easternmost side (z=800 km) at 120 km depth (Figure 3). Immediately after this necking develops into tearing at 4.24 Myr, the incoming continental block stops due to the collision with the lower-plate continent. The tearing point propagates westward reflecting in the tilted angle of the slab’s top edge as shown in Figure 3E and Figure 4.
FIGURE 3
FIGURE 4
While the slab is fully attached, the down-dip motion of the slab induces corner flows, and the large slab body induces a large flow around the slab’s edges (Figure 5). The initiation of slab necking and tearing on the eastern side is inherent to our model setup because the slab’s easternmost part is the only region in which the continental-continental convergence occurs (Figure 4a1). As the oceanic plate is fully subducted in the eastern side, slab retreat is inhibited and subduction ends in that area due to the lower density of the upper continental plate continent. This initial slab tearing increases the slab pull on the still attached slab and the slab tear starts propagating westward. Towards the west, the subsequent lithosphere tearing is purely a result of the tearing process that has been set in motion from the east and not a tearing following continental collision. The different amounts of exhumed oceanic crust in the forearc wedge (Figure 4b2, c3) appear to be depending on how large the remaining oceanic domain is in between the southern incoming plate and the northern continental margin. We observe a larger amount of exhumed oceanic crust in the westernmost side (Figure 4c3).
FIGURE 5
The tearing localises where high stress and strain rate are caused between the buoyant continental part of the plate (upward force) and the negatively buoyant hanging slab (downward force). This initiation of slab tearing is reflected by a sharp surface uplift along the collisional belt, with uplift rate ranges from 0.23 mm/yr to 2.16 mm/yr throughout the tear propagation. The rise in elevation caused by the detachment of the slab affects an ever-increasing area (Figures 4A–C), migrating westward and reflecting the tear propagation. The slab is completely detached after t=5.75 Myr (tearing duration of ∼1.65 Myr).
3.2 Influence of model parameters
3.2.1 Effect of no incoming continental block (Mod2)
The reference model (Mod1-reference) incorporates an incoming buoyant continental block implemented to create a continental-continental collision, which the diachronous collision then led to a one-sided slab tearing. We now move on to investigate how the absence of this incoming buoyant continental block affects the subduction dynamics. Rheologically, Mod2 mimics Mod1-reference but the absence of an incoming continental block creates a continental-oceanic arc. At t=3.48 Myr, the retreating intra-oceanic subduction trench reaches the continental passive margin, after which the trench continues to retreat. After 0.5 Myr, high topography developed over the trench (Supplementary Figure S1A). Here, on the eastern side of the slab, the accumulation of crustal materials above the trench prevented the trench from retreating any further and led to the initiation of slab tearing at t=4.66 Myr. The slab is completely detached by t=5.70 Myr (tearing duration of ∼1.04 Myr). The uplift rate during the tear propagation ranges from 0.71 mm/yr to 1.35 mm/yr. The lack of incoming continental block thus does not prevent the initiation of slab tearing.
In Mod2, where the overriding plate does not incorporate a buoyant continental block, the slab-tearing dynamics are similar to the Mod1 (reference model), likely because both models have the same mantle rheological setup. Mod2’s lack of a buoyant continental block on the overriding plate does not affect the rate of trench retreat that is thus mainly controlled by the negative buoyancy of the oceanic slab and the asthenospheric mantle viscosity. In Mod1-reference, the presence of an incoming buoyant continental block does limit the extent of the forearc region, as illustrated in Figure 6. A less dense body (relative to the surrounding mantle) rises up the subduction channel and thrusts under the overlying crustal materials (Figures 6C,E,F). Uplift of the lighter mantle material in this region is partly driven by the extension related to roll back. The absence of a continental block in Mod2 allows the crustal material to spread farther compared to Mod1, where the spreading is limited by the buoyant continental block.
FIGURE 6
3.2.2 Role of mantle viscous rheology
The subduction in the reference model is spontaneous, in the sense that the slab sinks by its own weight causing the subsequent trench retreat. However, the velocity of this process is highly controlled by the viscosity of the mantle and the role of this parameter must be quantively explored. To this purpose we design two model setups (Mod3 and Mod4) in which viscosity is increased relative to the reference setup Mod1.
Mod3 adopts a more viscous mantle, which can be achieved by increasing the ductile viscosity through increasing the activation volume of the mantle (both Vdiff and Vdisl). This results in slowing down the down-going slab due to the increased resistance of the higher-viscosity sublithospheric mantle. In the model setup Mod1-reference, the prescribed activation volume for the dislocation creep is Vdisl=2.6 J/(mol·MPa) and for diffusion creep Vdiff=0.7 J/(mol·MPa), and in model Mod3 we use Vdisl=3.0 J/(mol·MPa) and Vdiff=0.8 J/(mol·MPa) (Table 1). The increased activation volume implies a stronger mantle viscosity increases with pressure (and therefore with depth, Supplementary Figure S2). The evolution of the subduction is similar to the reference model but with much slower rate. For example, when the slab has reached 450 km depth, the slab in Mod1-reference has a maximum downward velocity of 20 cm/yr (t=3.05 Myr) where as Mod3’s maximum downward velocity is 8 cm/yr (t=6.87 Myr). The slab tearing in Mod3 initiated at around t=11.08 Myr as oppose to t=4.34 Myr in the reference model. The surface topography above the initiation of tearing exhibits an elevation of ∼1.5 km (Supplementary Figure S1D), which is similar to Mod1-reference. The uplift rate during the tear propagation ranges from 0.75 mm/yr to 1.68 mm/yr. In Mod3, the slab tear initiated at t=11 Myr and the slab completely detached by t=12.95 Myr (tearing duration of ∼ 2 Myr).
The fast down-going slab, together with the fast trench retreat velocity in the reference model, causes segments of high stress (4–5 MPa) and high strain-rate (10−14–10−12 1/s) to develop at the depth of greater than 120 km which led to a deeper tearing depth. In Mod3, with a more viscous mantle, the ambient mantle offers less resistance to the incoming of the down-going slab leading to a more gradual and shallow stress build-up focussing within the bending zone of the slab. This shallow stress focussing, at a depth of less than 100 km, led to a shallower tearing compared to the reference model.
Another way to increase the viscosity of the mantle is to increase the brittle viscosity, i.e., the upper visco-plastic limit of viscosity. In Mod4, we increase the final internal friction coefficient (μ1 in Eq. 11) for the lower oceanic crust and the mantle from 0 in Mod1 to 0.3 in Mod4 (Table 1). By increasing this coefficient, we decrease the strain weakening by a factor of two and we significantly increase the effective visco-plastic viscosity of deformed cold lithospheric mantle at high pressure/depth. As a result, the slab fails to sink down into the ambient mantle on its own due to a high resistance to local brittle/plastic deformation. This lack of slab’s downward velocity also led to the termination of trench retreat all together (Supplementary Figure S3). The slab only reached a depth of 300 km. This rheological setup shows that a lithosphere-averaged effective viscosity higher than 1024 Pa s should prevent slab roll-back and tearing.
3.2.3 Effect of fixed the convergence velocity (Mod5)
In model Mod5, we take the first stage of the reference model (initial push, slab reaches the depth of 200 km) and enters the model into the second stage with an imposed constant convergence. The imposed plate convergence velocity is set to 4 mm/yr to mimic the average convergent velocity between the Iberian and African plates (Macchiavelli et al., 2017). This velocity is much slower than the velocity resulting from the free hanging slab in previous models, consequently it is reducing the slab retreat. Such slow velocity makes thermal diffusion more relevant relative to advection (Boonma et al., 2019) warming the slab up and lowering its density and viscosity, all of which lead to thermal erosion and lithospheric dripping (Supplementary Figure S4). No slab tearing occurs in this model, but instead a lithospheric dripping takes place (Supplementary Figure S4).
The enhanced thermal diffusion that the slab experiences and the long time that the slab spends hanging among the sublithospheric mantle both allow an arcuate (in plan-view) deformed lower-viscosity slab to develop. In the models dominated by slab roll-back (Mod1-reference, Mod2, Mod3, and Mod4), subduction and trench migration stop once the slab reached the passive margin and the tear has started. In Mod5, however, the continuous pushing of the incoming continental block creates a band of high elevation over the arcuate trench (Supplementary Figure S1D).
4 Discussion
4.1 Geometry of the passive margin and slab-tearing dynamics
In this study, we have explored the role in tearing dynamics of an oblique continental passive margin adding asymmetry to the collision process. Magni et al. (2014) used also 3D geodynamic numerical modelling to explain the formation of backarc basins and later a slab window with an along-trench variation in slab buoyancy, localizing extension within the overriding plate. The results from the present study show that the obliquity of the continental passive margin can also introduce an asymmetry in the dynamics of the subducted slab, initiating slab tearing because it promotes a laterally diachronous consumption of the subducted oceanic plate, and facilitating slab tearing in one end of the plate and its lateral propagation. The obliquity of the passive margin relative to the trench axis causes an initial localized contact of the upper plate continent with the edge of the subducted slab in the trench. The resistance of the buoyant continent to subduction, together with the oppositely-directed slab pull, causes the initial break-off of the oceanic part of the plate. This is different from models implementing a continental margin parallel to the subduction trench direction, where the lateral propagation of tearing is initiated differently. Averaging over 500 km distance (from z=300–800 km), the tearing velocities are 42.6 cm/yr in Mod1-reference, 67.6 cm/yr in Mod2, and 37.6 cm/yr in Mod3 (Table 2). Our tear propagation rates fall well within the range of previous estimations: 7–45 cm/yr from the Carpathians’ depocenter migration by Meulenkamp et al. (1996) from the evolution of the Carpathian-Pannonian system whose geological model is constructed using regional chronostratigraphic sequences; and 10–80 cm/yr from 3D numerical modelling of continental collision by van Hunen and Allen (2011). A 3D stress model accounting for rheology, tear length, and force distribution, by Yoshioka and Wortel (1995) obtained a wide range of tear propagation rates between 2 and 94 cm/yr.
TABLE 2
| Model | Description | Incoming continent | Slab detachment | Slab tear propagation (cm/yr) | Tearing duration (Myr) | Uplift rate from tearing (mm/yr) |
|---|---|---|---|---|---|---|
| Mod1 | Reference model | Yes | Yes | 42.6 | 1.65 | 0.23—2.16 |
| Mod2 | No incoming continental block | No | Yes | 67.6 | 1.04 | 0.71—1.35 |
| Mod3 | Higher ductile viscosity of the mantle | Yes | Yes | 37.6 | 2 | 0.75—1.68 |
| Mod4 | Higher brittle strength of the mantle | Yes | No | - | - | - |
| Mod5 | Fixed the convergence velocity | Yes | No | - | - | - |
Comparison between model results.
The 500 km wide slab takes about 2 Myr to completely detach, which is fast compared to the timescales needed for subduction. This velocity is in a similar order of magnitude than the slab tear propagation in the Betics, where the ages of emersion of the intraorogenic basins (the transition from marine to continental conditions) suggests a shift of 300–400 km migration of the uplifting region in a period lasting between 3 and 6 Myr (the ranges indicate uncertainty; Garcia-Castellanos and Villaseñor, 2011). In our results, the main parameter controlling the timing of the tearing is the mantle rheology. The viscosity of the sublithospheric mantle in Mod3 is ≤ 1022 Pa s whereas it is ≤ 1021 Pa s in Mod1-reference. The more viscous sublithospheric mantle in Mod3 slowed down the sinking slab, hence the slowest tear-propagating velocity. This tear propagation rate is faster than that obtained in Chertova et al. (2014) and Spakman et al. (2018) based on numerical models and seismic tomography. Because the upper mantle viscosity in their model is similar to ours, the difference might rise from the lack of a sharp jump to higher viscosity below the 640 km discontinuity (lower mantle) in our model, which facilitates slabs to sink down to the bottom of our model box, at −780 km.
The ranges of tearing depth from our models are 80–150 km on the eastern side and 170–200 km on the western side. These ranges are similar to previous numerical modelling studies such as 80–240 km from Freeburn et al. (2017), 95–140 km from Schellart (2017), 100–400 km from Gerya et al. (2004), and 120–145 km from Duretz et al. (2014). The difference in tearing depth could plausibly come from the different rheological setting in each numerical model, as well as different tectonic setup. In the models presented here, the slab’s eastern side has a shallow tearing depth, determined by the weakness in the transition zone between the continental and the oceanic lithosphere. Although tearing is triggered by the first subduction of the continental-oceanic boundary of the upper plate in the East, the tearing thereafter propagates through the subducted oceanic lithosphere, and the transition zone stops having a role. As the tearing propagates westwards, the tearing depth becomes deeper. Two mechanisms favouring tear propagation and possibly making it grow deeper are: (i) the negative buoyancy of the already detached portion of the slab hangs from an ever-decreasing length of the rest of the slab; and (ii) the mantle flow through the slab tear window (Figure 5). A similar pattern is reflected in the tearing depth. On the easternmost side, the tearing tends to occur within the subducted continental lithosphere portion, such that the detached slab pinches out some continental crust. Because the tearing depth is deeper in the west, breakoff tends to localise within the subducted oceanic lithosphere.
4.2 Surface uplift: Isostasy vs. dynamic topography
Two known contributions to surface topography are the isostatic equilibrium of the lithosphere floating atop the mantle and the dynamic topography (Forte et al., 1993) caused by the buoyancy-driven mantle convection exerting vertical stress onto the lithosphere. Dynamic subsidence is generally caused by downward mantle flow (downwelling), while dynamic uplift is caused by upward mantle flow (upwelling).
We aim at differentiating topographic changes relate to tectonic deformation from vertical epeirogenic motions related to mantle/slab tearing by separating these two components of topography in our models. We calculate the isostatic effect adopting a compensation depth of 150 km (128 km below crustal surface, where the lowest values of viscosity are predicted, roughly matching the depth of initial tearing). Figure 7 shows model Mod3’s modelled density distribution and the evolution of the modelled elevation, and the two contributions, the isostatic and the dynamic components. This isostatic elevation is due to the density changes at crustal and lithosphere scales, without accounting the dynamics of the flow associated to slab subduction.
FIGURE 7
The dynamic topography is then obtained from taking the isostatic effect away from the modelled elevation. The dynamic topography shown in Figure 7 appears to be reflecting properly the mantle dynamics. From the dynamic topography results, we can identify the slab pull effect with a subsidence and the corner flow and mantle upwelling with an uplift. Prior to slab tearing, the mantle flowing upwards in the subduction channel corresponds with the high dynamic topography (Figures 7A,B). As tearing begins, the tearing gap allows the mantle flow to go through and this channel upward flow is reduced (Figures 7C,D; Figure 5). While the slab is still attached, the slab-pull force is transmitted up to the crust, producing the subsidence on the passive margin (Figure 7A, x=250–300 km). As the slab starts necking and tearing, this transmission of slab-pull force vanishes, reducing the aforementioned subsidence (Figures 7C,D). The mantle convection decreases when the detached slab is at a depth of 450–660 km (Figure 7D).
We also set out to look at the time-response of surface topography to tearing in the mantle and the possible temporal delay involved. The one-to-one (instantaneous) interpretation has been widely utilised by previous studies (Lithgow-Bertelloni and Silver, 1998; Boschi et al., 2010; Faccenna and Becker, 2010; Faccenna et al., 2014; Gvirtzman et al., 2016; Heller and Liu, 2016; Austermann and Forte, 2019; Ávila and Dávila, 2020). However, the tearing process in our models occurs at a relatively fast velocity, which may make it difficult to capture and quantify this delay. Figure 8 displays the modelled evolution of the topographic response as the slab tearing laterally propagates westward. The incoming continental block collides with the passive margin and subsequently stops. The deformation associated to this continental-continental collision prior to tearing increases elevation by >1 km on the eastern side (z=800 km) (Figures 8A,E). As the tearing process starts at the eastern end, the topography increases by another 300 m. As tearing propagates westwards, this additional elevation propagates in the tearing direction (Figures 8B–E), attaining 500–1,000 m. The increase in surface elevation does not only occur above the tearing point but also in a nearby area, as shown in Figures 8E,F, reflecting the flexural lateral transmission of stresses in response to slab unloading. This is also due to the mantle flow: as the tear gap opens, it triggers higher poloidal flow, inducing trenchward mantle movement. This poloidal flow then induces a basal drag that drives trenchward motions under the two converging plates. The trenchward motion exerts compressional force to the relatively immobile subduction zone hinge, in addition to the opposing related to convergence, leading to an uplift of 0.3–0.8 km even before the arrival of the tear (Figure 8E). As the tearing propagates further westward, the high topography on the eastern side starts to subside by as much as 0.2 km (Figure 8E).
FIGURE 8
Figure 7 shows how the flow within the mantle changes during the subduction and tearing stages. As the slab subducts further, its volume in the mantle increases, obstructing the mantle flow and giving rise to corner flow in the mantle wedge. This corner flow increases the velocity of mantle convection by 3–10 cm/yr (Figure 7A), supporting the overlaying crust (Figures 7A,B). As the tearing of the slab initiates, it opens up a new pathway for the mantle to quickly flow through to replace the volume previously taken up by the slab (Figure 7C). Tearing causes the unload of the overlying crust (Figure 7C) leading to a rapid isostatic surface uplift as a prominent signature of slab detachment. Meanwhile, the velocity of the corner flow in the mantle wedge is reduced during the slab tearing (Figures 7B,C), decreasing the dynamic contribution to topography at the trench. As the detached slab sinks further down, the bottom of the slab hits the depth of 660–700 km discontinuity, where it becomes flatter (e.g., Figure 3F) while the mantle convection velocity slows down to pre-rollback values (4 cm/yr). This slower velocity reduces the elevation while the crust and the lithospheric mantle begin to thermally equilibrate. Overall, the surface uplift rates observed in our models, as a response to the slab tearing, range from 0.23 to 2.16 mm/yr. These values fall within surface uplift rates obtained from previous numerical modelling studies ranging from as low as 0.10 mm/yr to as high as 2.65 mm/yr (Andrews and Billen, 2009; Duretz et al., 2011).
4.3 Comparison to the tearing process in the Gibraltar Arc
The Neogene tectonic evolution of the western Mediterranean is a topic of controversy and numerous geodynamic mechanisms have been proposed. The currently prevalent models invoke a slab rollback that ends by reaching the south Iberian margin. Some authors consider an originally NW-dipping subduction zone located along the SE borders of Corsica-Sardinia and Balearic promontories (e.g., Rosenbaum et al., 2002; Spakman and Wortel, 2004; van Hinsbergen et al., 2014) where the slab experiences a clockwise rotation of about 180° dragged by roll-back until reaching the southern Iberian margin and creating the arcuate Betic-Rif orogen. An alternative model considers an originally SE-dipping subduction under Africa, then retreating to the northwest and originating the arcuate Betic-Rif fold and thrust belts on Iberian and NW African plates (Vergés and Fernàndez, 2012). Regardless of the precise geodynamic mechanism and kinematics responsible for the formation of the slab, what is relevant to the present study is that these alternative tectonic evolution models agree on the present SE dipping of the slab that rolled back towards the South Iberian margin prior to 6 Ma. Any tectonic model implying 1) a subducting plate laterally ending towards the East, and 2) its earlier consumption by subduction in that edge, would equally justify our initial setup. Therefore, according to our results, any such model would be consistent with slab tearing starting from the eastern end of the slab and leading to the present slab geometry observed in seismic tomography (Figure 9; Spakman and Wortel, 2004; Bezada et al., 2013; Palomeras et al., 2014).
FIGURE 9
Several 3D numerical modelling studies have addressed the geodynamic evolution of the last 35 My of the western Mediterranean (Chertova et al., 2014; Spakman et al., 2018; Capella et al., 2020; Peral et al., 2022). A thermo-mechanical model by Negredo et al. (2020) suggests that sharp lateral thermal and rheological variations across STEP faults following slab-tearing, favours the occurrence of continental delamination. In the present study, the Gibraltar Arc is used as a reference because the relationship between slab tearing and surface uplift has previously been described in this region (Garcia-Castellanos and Villaseñor, 2011; Spakman et al., 2018). In our model, slab tearing initiates at the eastern slab edge, which according to the tectonic evolution envisaged by Vergés and Fernàndez (2012) may have resulted from a subduction polarity change under the transition between the Betic Cordillera and the Balearic Promontory. The uplift rates obtained (0.23 mm/yr to 2.16 mm/yr) are consistent with the rates considered necessary to compensate erosion of a sustained inflow during the first stage of the Messinian Salinity Crisis (Garcia-Castellanos and Villaseñor, 2011). The prevalent gypsum precipitation during this period suggests a reduced connectivity between the Atlantic Ocean and the Mediterranean Sea that has been linked to a competition between tectonic uplift of the seaway and the erosion of the seaway by the inflowing waters. This inflow from the Atlantic Ocean is needed to compensate for evaporation in the Mediterranean and to explain the gypsum accumulation from 5.97 Ma (e.g., Andreetto et al., 2021). This tectonic-erosion competition would explain the persistent but restricted water inflow from the Atlantic into the Mediterranean Sea required to explain the precipitation of gypsum during the first stage of the MSC. The model by Garcia-Castellanos and Villaseñor (2011) also leads to a critical uplift rate in the range of a few mm/yr needed to close the seaways across the Gibraltar Arc. Coulson et al. (2019) incorporated to that model the effect of glacioisostasy in response to sea level changes, lowering the required critical uplift rate to < 1.5 mm/yr.
Regarding the time span of this uplift, previous stratigraphic studies in the intrabetic basins suggested that the age of marine to continental transition (basin emersion) in the region is in the range of 11–8 Ma, younger towards the West (Garcés et al., 1998; Krijgsman et al., 2000; Iribarren et al., 2009). The last of these basins’ emersion was previously dated at 5.3 Ma, which indicates that the duration of tear propagation in the Betics (across 400 km distance) could last as long as 3–5 Myr, yielding a tearing rate of 80–133 km/Myr. Nevertheless, more recent biostratigraphic studies by Krijgsman et al. (2018); van der Schee et al. (2018) dated the age of the latest preserved marine sediment in the Guadalhorce Corridor (Ronda, Antequera and Arcos sedimentary basins) to 7.51 Ma (late Tortonian; westernmost age in Figure 1B). This would narrow down the plausible tear propagation duration to 1–3 Myr, implying a tearing rate of (133–400 km/Myr), a priori more compatible with the outcomes of our reference model (370–670 km/Myr). However, the available constraints on the timing of uplift show more complexity. The Alborán Sea coastal uplift reported by Guerra-Merchán et al. (2014) may have propagated further to the W until 3.2 Ma, leaving marine sediment at present elevations locally above 150 m. Furthermore, in the East, the Lorca Basin may not have become continental as soon as 11 Ma (Figure 1B) but during Messinian times (Carpentier et al., 2020). Moreover, Duggen et al. (2003) link the topographic uplift to the change in volcanic geochemistry between 6.3 and 4.8 Ma, from a subduction-related to a sublithospheric-mantle source. Seismic focal mechanisms suggest that tearing is today ongoing underneath the city of Málaga (Ruiz-Constán et al., 2011; Mancilla et al., 2015), although probably at much slower rates than those during the Messinian, since that location is close to the Guadalhorce Corridor: if the tearing propagation rates were similar to those mentioned above, it would have already affected the entire slab, leading to a break off that is not visible in the available seismic tomography. In summary, the uncertainties in the measured chronology of the intrabetic basin uplift and in the modelling methodology (both in mantle’s viscosity and the initial setup) impede a more detailed comparison between modelled and observed tearing propagation rates.
The uplift of the intramountain basins within the Betics in southern Iberia is higher on the eastern side (Iribarren et al., 2009), where the slab is interpreted to be detached earlier based on seismic tomography (Figure 9B), relative to the western Betics, where the slab still remains attached to Iberia (Figure 9A). Our models predict a similar trend and geometry, with earlier and higher uplift on the eastern side due to a combination of tectonic shortening deformation and to slab tearing, and a smaller and later uplift in the western end of the slab, where it remains attached (at t=4.24 My; Figures 9C,D).
5 Conclusion
Our results support the idea that the lateral ending of the subducting plate and the obliquity of its subduction relative to the continental passive margin favours the initiation of lateral slab tearing. The obliquity favours the consumption of the subducting plate earlier at the slab edge, determining the first contact between slab and continental margin in one side and leading to incipient slab necking. In the setups explored in this paper, this obliquity leads to slab tearing velocities of ∼37.6–67.6 cm/yr (for a lower-mantle viscosity of up to 1022 Pas) and a surface uplift of 0.5–1.5 km across the forearc region throughout the tearing process.
Within the limited range of model runs presented, lithospheric tearing along the subducting plate does not relate specifically to the process of continental collision or to the obliquity between plate motions, but to the obliquity between the continental margin of the subducting plate and the subduction trench. This obliquity determines in our model which edge of the subducting plate starts tearing first.
The slab tearing depth increases as it propagates along the slab, shallower on the side where the tear initiated (80–150 km in our setups) and deeper tear on the other edge (170–200 km). The speed of slab tearing is geologically fast (370–670 km/Myr, with a duration og 3 Myr in our model scenarios). The key controls on the speed of detachment process are the viscosity of the sublithospheric mantle (low viscosity implying a faster tearing) and the amount of shortening/oceanic subduction prior to tearing (determining the slab pull of the subducted oceanic slab).
We obtain uplift rates ranging from 0.23 to 2.16 mm/yr, compatible with uplift rates needed to achieve an equilibrium between seaway uplift and seaway erosion proposed as responsible for the closure of marine gateways that reduced the water-flow from the Atlantic Ocean into the Mediterranean Sea during the first stage of the Messinian Salinity Crisis.
Statements
Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: http://doi.org/10.5281/zenodo.4637879.
Author contributions
KB: Methodology, Software, Investigation, Formal analysis, Validation, Writing—Original Draft, Visualisation. DG: Conceptualisation, Resources, Funding acquisition, Supervision, Validation, and Writing—Review and Editing. IJ: Conceptualisation, Resources, Funding acquisition, Supervision, Validation, and Writing—Review and Editing. TG: Methodology, Software, and Writing—Review and Editing.
Funding
This work has been supported by EU Marie Curie Initial Training Network “SUBITOP” (674899-SUBITOP-H2020-MSCA-ITN-2015), the Spanish Government national research program (GeoCAM, PGC 2018–095154-B-I00) and the Generalitat de Catalunya grant (AGAUR 2021 SGR 00410). We thank the Laboratorio de Geodinámica at GEO3BCN-CSIC as well as the Euler Cluster at the Scientific Computing centre at ETH Zürich for providing the computing facilities.
Acknowledgments
We gratefully thank Wim Spakman and Nicolas Riel, and the editors Frederic Mouthereau and Guillermo Booth-Rea for the careful and critical reviews that contributed to clarify and enrich this manuscript. Rob Govers also provided useful comments on a previous version of this manuscript.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart.2023.1095229/full#supplementary-material
References
1
AndrewsE. R.BillenM. I. (2009). Rheologic controls on the dynamics of slab detachment. Tectonophysics464 (1–4), 60–69. 10.1016/j.tecto.2007.09.004
2
AndreettoF.AloisiG.RaadF.HeidaH.FleckerR.AgiadiK.et al (2021). Freshening of the Mediterranean Salt Giant: controversies and certainties around the terminal (Upper Gypsum and Lago-Mare) phases of the Messinian Salinity Crisis. Earth-Sci. Reviews216. 10.1016/j.earscirev.2021.103577
3
AngrandP.MouthereauF. (2021). Evolution of the Alpine orogenic belts in the Western Mediterranean region as resolved by the kinematics of the Europe-Africa diffuse plate boundary. Bsgf - Earth Sci. Bull.192, 42. 10.1051/bsgf/2021031
4
AustermannJ.ForteA. M. (2019). The importance of dynamic topography for understanding past sea-level changes. Past. Glob. Chang. Mag.27 (1). 10.22498/pages.27.1.18
5
ÁvilaP.DávilaF. M. (2020). Lithospheric thinning and dynamic uplift effects during slab window formation, southern Patagonia (45˚-55˚ S). J. Geodyn.133, 101689. 10.1016/j.jog.2019.101689
6
BezadaM. J.HumphreysE.ToomeyD.HarnafiM.DavilaJ.GallartJ. (2013). Evidence for slab rollback in westernmost Mediterranean from improved upper mantle imaging. Earth Planet. Sci. Lett.368, 51–60. 10.1016/j.epsl.2013.02.024
7
BillenM. I. (2008). Modeling the dynamics of subducting slabs. Annu. Rev. Earth Planet. Sci.36, 325–356. 10.1146/annurev.earth.36.031207.124129
8
BoonmaK.KumarA.Garcia-CastellanosD.Jiménez-MuntI.FernándezM. (2019). Lithospheric mantle buoyancy: The role of tectonic convergence and mantle composition. Sci. Rep.9 (1), 17953. 10.1038/s41598-019-54374-w
9
BoschiL.FaccennaC.BeckerT. W. (2010). Mantle structure and dynamic topography in the Mediterranean Basin. Geophys. Res. Lett.37 (20), 20303. 10.1029/2010GL045001
10
BurgJ. P.GeryaT. V. (2005). The role of viscous heating in Barrovian metamorphism of collisional orogens: Thermomechanical models and application to the Lepontine Dome in the Central Alps. J. Metamorph. Geol.23 (2), 75–95. 10.1111/j.1525-1314.2005.00563.x
11
CapellaW.SpakmanW.HinsbergenD. J. J.ChertovaM. V.KrijgsmanW. (2020). Mantle resistance against Gibraltar slab dragging as a key cause of the Messinian Salinity Crisis. Terra nova.32 (2), 141–150. 10.1111/ter.12442
12
CarminatiE.WortelR.SpakmanW.SabadiniR. (1998). The role of slab detachment processes in the opening of the Western–central mediterranean basins: Some geological and geophysical evidence. Earth Planet. Sci. Lett.160 (3–4), 651–665. 10.1016/S0012-821X(98)00118-6
13
CarpentierC.VenninE.RouchyJ. M.CornéeJ. J.Melinte-DobrinescuM.HibschC.et al (2020). Ages and stratigraphical architecture of late Miocene deposits in the Lorca Basin (Betics, SE Spain): New insights for the salinity crisis in marginal basins. Sediment. Geol.405, 105700. 10.1016/j.sedgeo.2020.105700
14
ChertovaM. V.SpakmanW.GeenenT.van den BergA. P.van HinsbergenD. J. J. (2014). Underpinning tectonic reconstructions of the Western Mediterranean region with dynamic slab evolution from 3-D numerical modeling. J. Geophys. Res. Solid Earth119 (7), 5876–5902. 10.1002/2014JB011150
15
ConradC. P.Lithgow-BertelloniC. (2002). How mantle slabs drive plate tectonics. Science298 (5591), 207–209. 10.1126/science.1074161
16
CoulsonS.PicoT.AustermannJ.PowellE.MouchaR.MitrovicaJ. X. (2019). The role of isostatic adjustment and gravitational effects on the dynamics of the Messinian salinity crisis. Earth Planet. Sci. Lett.525, 115760. 10.1016/j.epsl.2019.115760
17
CrameriF.SchmelingH.GolabekG. J.DuretzT.OrendtR.BuiterS. J. H. H.et al (2012). A comparison of numerical surface topography calculations in geodynamic modelling: An evaluation of the ‘sticky air’ method. Geophys. J. Int.189 (1), 38–54. 10.1111/j.1365-246X.2012.05388.x
18
DaviesJ. H.von BlanckenburgF. (1995). Slab breakoff: A model of lithosphere detachment and its test in the magmatism and deformation of collisional orogens. Earth Planet. Sci. Lett.129 (1–4), 85–102. 10.1016/0012-821X(94)00237-S
19
DuggenS.HoernleK.van den BogaardP.R€upkeL.Phipps MorganJ. (2003). Deep roots of the Messinian salinity crisis. Nature422, 602–606. 10.1038/nature01553
20
DuretzT.GeryaT. V.MayD. A. (2011). Numerical modelling of spontaneous slab breakoff and subsequent topographic response. Tectonophysics502 (1–2), 244–256. 10.1016/j.tecto.2010.05.024
21
DuretzT.GeryaT. V.SpakmanW. (2014). Slab detachment in laterally varying subduction zones: 3-D numerical modeling. Geophys. Res. Lett.41 (6), 1951–1956. 10.1002/2014GL059472
22
ElsasserW. M. (1971). Sea-floor spreading as thermal convection. J. Geophys. Res.76, 1101–1112. 10.1029/jb076i005p01101
23
EnglandP.MolnarP. (1990). Surface uplift, uplift of rocks, and exhumation of rocks. Geology18 (12), 1173–1177. 10.1130/0091-7613(1990)018<1173:SUUORA>2.3.CO;2
24
FaccennaC.BeckerT. W. (2010). Shaping mobile belts by small-scale convection. Nature465 (7298), 602–605. 10.1038/nature09064
25
FaccennaC.PiromalloC.Crespo-BlancA.JolivetL.RossettiF. (2004). Lateral slab deformation and the origin of the Western Mediterranean arcs. Tectonics23 (1). 10.1029/2002tc001488
26
FaccennaC.BellierO.MartinodJ.PiromalloC.RegardV. (2006). Slab detachment beneath eastern Anatolia: A possible cause for the formation of the north anatolian fault. Earth Planet. Sci. Lett.242 (1–2), 85–97. 10.1016/J.EPSL.2005.11.046
27
FaccennaC.BeckerT. W.MillerM. S.SerpelloniE.WillettS. D. (2014). Isostasy, dynamic topography, and the elevation of the Apennines of Italy. Earth Planet. Sci. Lett.407, 163–174. 10.1016/j.epsl.2014.09.027
28
ForteA. M.PeltierW. R.DziewonskiA. M.WoodwardR. L. (1993). Dynamic surface topography: A new interpretation based upon mantle flow models derived from seismic tomography. Geophys. Res. Lett.20 (3), 225–228. 10.1029/93GL00249
29
FoxM.HermanF.KisslingE.WillettS. D. (2015). Rapid exhumation in the Western Alps driven by slab detachment and glacial erosion. Geology43 (5), 379–382. 10.1130/G36411.1
30
FreeburnR.BouilholP.MaunderB.MagniV.van HunenJ.HunenJ. V.et al (2017). Numerical models of the magmatic processes induced by slab breakoff. Earth Planet. Sci. Lett.478, 203–213. 10.1016/j.epsl.2017.09.008
31
GarcésM.KrijgsmanW.AgustíJ. (1998). Chronology of the late turolian deposits of the fortuna basin (SE Spain): Implications for the messinian evolution of the eastern Betics. Earth Planet. Sci. Lett.163 (1–4), 69–81. 10.1016/S0012-821X(98)00176-9
32
Garcia-CastellanosD.VillaseñorA. (2011). Messinian salinity crisis regulated by competing tectonics and erosion at the Gibraltar arc. Nature480 (7377), 359–363. 10.1038/nature10651
33
Garcia-CastellanosD.TorneM.FernàndezM. (2000). Slab pull effects from a flexural analysis of the Tonga and Kermadec trenches (Pacific Plate). Geophys. J. Int.141 (2), 479–484. 10.1046/j.1365-246X.2000.00096.x
34
GarzantiE.RadeffG.MalusàM. G. (2018). Slab breakoff: A critical appraisal of a geological theory as applied in space and time. Earth-Science Rev.177, 303–319. Elsevier. 10.1016/j.earscirev.2017.11.012
35
GeryaT. V.YuenD. (2003). Characteristics-based marker-in-cell method with conservative finite-differences schemes for modeling geological flows with strongly variable transport properties. Phys. Earth Planet. Inter.140 (4), 293–318. 10.1016/j.pepi.2003.09.006
36
GeryaT. V.YuenD. (2007). Robust characteristics method for modelling multiphase visco-elasto-plastic thermo-mechanical problems. Phys. Earth Planet. Interiors163 (1–4), 83–105. 10.1016/j.pepi.2007.04.015
37
GeryaT. V.YuenD. A.MareschW. V. (2004). Thermomechanical modelling of slab detachment. Earth Planet. Sci. Lett.226 (1–2), 101–116. 10.1016/j.epsl.2004.07.022
38
GeryaT. V.PerchukL. L.BurgJ. P. (2008). Transient hot channels: Perpetrating and regurgitating ultrahigh-pressure, high-temperature crust-mantle associations in collision belts. Lithos103 (1–2), 236–256. 10.1016/j.lithos.2007.09.017
39
GeryaT. V. (2013). Three-dimensional thermomechanical modeling of oceanic spreading initiation and evolution. Phys. Earth Planet. Inter.214, 35–52. 10.1016/j.pepi.2012.10.007
40
Guerra-MerchánA.SerranoF.HlilaR.El KadiriK.Sanz de GaldeanoC.GarcésM. (2014). Tectono-sedimentary evolution of the peripheral basins of the Alboran Sea in the arc of Gibraltar during the latest Messinian-Pliocene. J. Geodyn.77, 158–170. ISSN 0264-3707. 10.1016/j.jog.2013.12.003
41
GvirtzmanZ.FaccennaC.BeckerT. W. (2016). Isostasy, flexure, and dynamic topography. Tectonophysics683, 255–271. 10.1016/j.tecto.2016.05.041
42
HarlowF. H.WelchJ. E. (1965). Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface. Phys. Fluids8 (12), 2182. 10.1063/1.1761178
43
HellerP. L.LiuL. (2016). Dynamic topography and vertical motion of the U.S. Rocky Mountain region prior to and during the Laramide orogeny. Bull. Geol. Soc. Am.128 (5–6), 973–988. 10.1130/B31431.1
44
HirthG.KohlstedtD. (2003). “Rheology of the upper mantle and the mantle wedge: A view from the experimentalists,” in Inside the subduction factory. 83–105. 10.1029/138GM06
45
IribarrenL.VergésJ.FernàndezM. (2009). Sediment supply from the betic–rif orogen to basins through Neogene. Tectonophysics475 (1), 68–84. 10.1016/j.tecto.2008.11.029
46
ItoK.KennedyG. C. (1971). “An experimental study of the basalt-garnet granulite-eclogite transition,” in The structure and physical properties of the Earth’s crust. Editor HeacockJ. G. (American Geophysical Union), 14, 303–314. 10.1029/GM014p0303
47
ItoE.AkaogiM.ToporL.NavrotskyA. (1990). Negative pressure-temperature slopes for reactions formign MgSiO 3 perovskite from calorimetry. Science249 (4974), 1275–1278. 10.1126/science.249.4974.1275
48
Jiménez-MuntI.PlattJ. P. (2006). Influence of mantle dynamics on the topographic evolution of the Tibetan Plateau: Results from numerical modelling. Tectonics25, TC6002. 10.1029/2006TC001963
49
Jiménez-MuntI.TorneM.FernàndezM.VergésJ.KumarA.CarballoA.et al (2019). Deep seated density anomalies across the iberia-africa plate boundary and its topographic response. J. Geophys. Res. Solid Earth124 (12), 13310–13332. 10.1029/2019JB018445
50
KaratoS. I.WuP. (1993). Rheology of the upper mantle: A synthesis. Science260 (5109), 771–778. 10.1126/science.260.5109.771
51
KarigD. E. (1971). Origin and development of marginal basins in the Western Pacific. J. Geophys. Res.76, 2542–2561. 10.1029/jb076i011p02542
52
KatsuraT.ItoE. (1989). The system Mg 2 SiO 4 -Fe 2 SiO 4 at high pressures and temperatures: Precise determination of stabilities of olivine, modified spinel, and spinel. J. Geophys. Res. Solid Earth94 (B11), 15663–15670. 10.1029/JB094iB11p15663
53
KatsuraT.YonedaA.YamazakiD.YoshinoT.ItoE.SuetsuguD.et al (2010). Adiabatic temperature profile in the mantle. Phys. Earth Planet. Interiors183 (1–2), 212–218. 10.1016/j.pepi.2010.07.001
54
KrijgsmanW.GarcésM.AgustıJ.RaffiI.TabernerC.ZachariasseW. J. (2000). The ‘Tortonian salinity crisis’ of the eastern Betics (Spain). Earth Planet. Sci. Lett.181 (4), 497–511. 10.1016/s0012-821x(00)00224-7
55
KrijgsmanW.CapellaW.SimonD.HilgenF. J.KouwenhovenT. J.MeijerP. T.et al (2018). The Gibraltar corridor: Watergate of the messinian salinity crisis. Mar. Geol.403, 238–246. 10.1016/j.margeo.2018.06.008
56
KumarA.FernàndezM.VergésJ.TorneM.Jiménez‐MuntI. (2021). Opposite symmetry in the lithospheric structure of the Alboran and Algerian basins and their margins (Western Mediterranean): Geodynamic implications. J. Geophys. Res. Solid Earth126 (7). 10.1029/2020JB021388
57
LiZ. H.XuZ.GeryaT.BurgJ. P. (2013). Collision of continental corner from 3-D numerical modeling. Earth Planet. Sci. Lett.380, 98–111. 10.1016/j.epsl.2013.08.034
58
LiangX.ChenY.TianX.ChenY. J.NiJ.GallegosA.et al (2016). 3D imaging of subducting and fragmenting Indian continental lithosphere beneath southern and central Tibet using body-wave finite-frequency tomography. Earth Planet. Sci. Lett.443, 162–175. 10.1016/J.EPSL.2016.03.029
59
Lithgow-BertelloniC.SilverP. G. (1998). Dynamic topography, plate driving forces and the African superswell. Nature395 (6699), 269–272. 10.1038/26212
60
MacchiavelliC.VergésJ.SchettinoA.FernàndezM.TurcoE.CascielloE.et al (2017). A new southern north atlantic isochron map: Insights into the drift of the iberian plate since the late cretaceous. J. Geophys. Res. Solid Earth122 (12), 9603–9626. 10.1002/2017JB014769
61
MagniV.FaccennaC.van HunenJ.FunicielloF. (2014). How collision triggers backarc extension: Insight into Mediterranean style of extension from 3-D numerical models. Geology42, 511–514. 10.1130/G35446.1
62
MagniV.AllenM. B.van HunenJ.BouilholP.HunenJ. V.BouilholP.et al (2017). Continental underplating after slab break-off. Earth Planet. Sci. Lett.474, 59–67. 10.1016/j.epsl.2017.06.017
63
MancillaF. d. L.Booth-ReaG.StichD.Pérez-PeñaJ. V.MoralesJ.AzañónJ. M.et al (2015). Slab rupture and delamination under the Betics and Rif constrained from receiver functions. Tectonophysics663, 225–237. 10.1016/j.tecto.2015.06.028
64
MeulenkampJ. E.KováčM.CichaI. (1996). On Late Oligocene to Pliocene depocentre migrations and the evolution of the Carpathian-Pannonian system. Tectonophysics266 (1–4), 301–317. 10.1016/S0040-1951(96)00195-3
65
MitrofanH.AnghelacheM. A.ChiteaF.DamianA.CadicheanuN.VişanM. (2016). Lateral detachment in progress within the Vrancea slab (Romania): Inferences from intermediate-depth seismicity patterns. Geophys. J. Int.205 (2), 864–875. 10.1093/gji/ggv533
66
NegredoA. M.MancillaF. d. L.ClementeC.MoralesJ.FulleaJ. (2020). Geodynamic modeling of edge-delamination driven by subduction-transform edge propagator faults: The westernmost mediterranean margin (central betic orogen) case study. Front. Earth Sci.8, 6. 10.3389/feart.2020.533392
67
PalomerasI.ThurnerS.LevanderA.LiuK.VillasenorA.CarbonellR.et al (2014). Finite‐frequency Rayleigh wave tomography of the Western Mediterranean: Mapping its lithospheric structure. Geochem. Geophys. Geosystems15, 140–160. 10.1002/2013gc004861
68
PeralM. F.VergésJ.ZlotnikS.Jiménez-MuntI.Jimenez-MuntI. (2022). Numerical modelling of opposing subduction in the Western Mediterranean. Tectonophysics830, 229309. 10.1016/j.tecto.2022.229309
69
QayyumA.LomN.AdvokaatE. L.SpakmanW.van der MeerD. G.van HinsbergenD. J. J. (2022). Subduction and slab detachment under moving trenches during ongoing India-Asia convergence. Geochem. Geophys. Geosystems23, e2022GC010336. 10.1029/2022GC010336
70
RanalliG. (1995). Rheology of the Earth. 2nd ed. Springer. ISBN: 978-0-412-54670-9.
71
RosenbaumG.ListerG. S.DubozC. (2002). Reconstruction of the tectonic evolution of the Western mediterranean since the oligocene. J. Virtual Explor.8. 10.3809/jvirtex.2002.00053
72
RosenbaumG.GasparonM.LucenteF. P.PeccerilloA.MillerM. S. M. S. (2008). Kinematics of slab tear faults during subduction segmentation and implications for Italian magmatism. Tectonics27 (2). 10.1029/2007TC002143
73
Ruiz-ConstánA.Galindo-ZaldívarJ.PedreraA.CélérierB.Marín-LechadoC. (2011). Stress distribution at the transition from subduction to continental collision (northwestern and central Betic Cordillera). Geochem. Geophys. Geosyst.12, Q12002. 10.1029/2011GC003824
74
SchellartW. P. (2017). A geodynamic model of subduction evolution and slab detachment to explain Australian plate acceleration and deceleration during the latest Cretaceous-early Cenozoic. Lithosphere9 (6), 976–986. 10.1130/L675.1
75
SchmelingH.BabeykoA. Y.EnnsA.FaccennaC.FunicielloF.GeryaT.et al (2008). A benchmark comparison of spontaneous subduction models-Towards a free surface. Phys. Earth Planet. Interiors171 (1–4), 198–223. 10.1016/j.pepi.2008.06.028
76
ŞengörA. M. C.ÖzerenS.GençT.ZorE. (2003). East Anatolian high plateau as a mantle-supported, north-south shortened domal structure. Geophys. Res. Lett.30 (24). 10.1029/2003GL017858
77
SinclairH. D. (1997). Flysch to molasse transition in peripheral foreland basins: The role of the passive margin versus slab breakoff. Geology25 (12), 1123–1126. 10.1130/0091-7613(1997)025<1123:FTMTIP>2.3.CO;2
78
SpakmanW.WortelR. (2004). “A tomographic view on western mediterranean geodynamics,” in The TRANSMED atlas. The mediterranean region from crust to mantle (Berlin, Heidelberg: Springer Berlin Heidelberg), 31–52. 10.1007/978-3-642-18919-7_2
79
SpakmanW.WortelM. J. R.VlaarN. J. (1988). The hellenic subduction zone: A tomographic image and its geodynamic implications. Geophys. Res. Lett.15 (1), 60–63. 10.1029/GL015i001p00060
80
SpakmanW.ChertovaM. V.Van Den BergA.van HinsbergenD. J. J. J. (2018). Puzzling features of Western Mediterranean tectonics explained by slab dragging. Nat. Geosci.11 (3), 211–216. 10.1038/s41561-018-0066-z
81
TurcotteD.SchubertG. (2014). Geodynamics. Cambridge University Press. 10.1017/CBO9780511843877
82
Van Der MeulenM. J.MeulenkampJ. E.WortelM. J. R. (1998). Lateral shifts of Apenninic foredeep depocentres reflecting detachment of subducted lithosphere. Earth Planet. Sci. Lett.154 (1–4), 203–219. 10.1016/s0012-821x(97)00166-0
83
van der ScheeM.van den BergB. C. J.CapellaW.SimonD.SierroF. J.KrijgsmanW. (2018). New age constraints on the Western betic intramontane basins: A late tortonian closure of the Guadalhorce corridor?Terra nova.30 (5), 325–332. 10.1111/ter.12347
84
Van HinsbergenD. J. J. J.VissersR. L. M. M.SpakmanW. (2014). Origin and consequences of Western Mediterranean subduction, rollback, and slab segmentation. Tectonics33 (4), 393–419. 10.1002/2013TC003349
85
van HunenJ.AllenM. B. (2011). Continental collision and slab break-off: A comparison of 3-D numerical models with observations. Earth Planet. Sci. Lett.302 (1–2), 27–37. 10.1016/j.epsl.2010.11.035
86
VergésJ.FernàndezM. (2012). Tethys–atlantic interaction along the iberia–africa plate boundary: The betic–rif orogenic system. Tectonophysics579, 144–172. 10.1016/j.tecto.2012.08.032
87
WortelR.SpakmanW. (1992). “Structure and dynamics of subducted lithosphere in the Mediterranean region,” in Proceedings of the koninklijke nederlandse akademie van Wetenschappen/C, 95, 325–347.
88
WortelM. J. R.SpakmanW. (2000). Subduction and slab detachment in the mediterranean-carpathian region. Science290 (5498), 1910–1917. 10.1126/science.290.5498.1910
89
WuF.-Y.JiW.-Q.WangJ.-G.LiuC.-Z.ChungS.-L.CliftP. D. (2014). Zircon U-Pb and Hf isotopic constraints on the onset time of India-Asia collision. Am. J. Sci.314 (2), 548–579. 10.2475/02.2014.04
90
YoshiokaS.WortelR. (1995). Three-dimensional numerical modeling of detachment of subducted lithosphere. J. Geophys. Res.100 (B10), 20223–20244. 10.1029/94jb01258
Summary
Keywords
geodynamical modeling, dynamic topography, mantle dynamics, subduction, Gibraltar Arc
Citation
Boonma K, García-Castellanos D, Jiménez-Munt I and Gerya T (2023) Thermomechanical modelling of lithospheric slab tearing and its topographic response. Front. Earth Sci. 11:1095229. doi: 10.3389/feart.2023.1095229
Received
10 November 2022
Accepted
31 March 2023
Published
17 April 2023
Volume
11 - 2023
Edited by
Guillermo Booth-Rea, University of Granada, Spain
Reviewed by
Wim Spakman, Utrecht University, Netherlands
Frederic Mouthereau, Université Toulouse III Paul Sabatier, France
Updates
Copyright
© 2023 Boonma, García-Castellanos, Jiménez-Munt and Gerya.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Ivone Jiménez-Munt, ivone@geo3bcn.csic.es
This article was submitted to Solid Earth Geophysics, a section of the journal Frontiers in Earth Science
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.