The Gibraltar slab dynamics and its influence on past and present-day Alboran domain deformation: Insights from thermo-mechanical numerical modelling

The origin and tectonic evolution of the Gibraltar Arc system is the result of a complex geodynamic evolution involving the convergence of the Eurasian and African plates and the dynamic impact of the Gibraltar slab. Although geologic and geophysical data collected in the last few years have increased our knowledge of the Gibraltar Arc region, it is still unclear which are the mechanical links between the Gibraltar slab and the past deformation of the overriding Alboran lithosphere, as well as to which degree this subduction system is presently active. In this study, we use 2D numerical modelling to investigate the impact of the Gibraltar slab dynamics on the deformation of the overriding Alboran lithosphere. Our model simulates a WE generic vertical section at an approximate latitude of 36°N and considers an initial setup at about Burdigalian times (∼20 Ma), when the subduction front position is relatively well constrained by recent tectonic reconstructions. Our modelling shows a switch in the overriding plate (OP) stress state from extensional stresses during the slab rollback to compressional stresses near the trench when the rollback velocity decreases, caused by the change in slab-induced mantle flow. We also find that much of the crustal and lithospheric deformation occur during fast slab rollback and OP extension in the first 10 Myr of evolution, while after that only moderate deformation associated with subduction is predicted. Finally, we find that despite the subduction rollback ceases, the ongoing motion of the deeper portion of the slab induces a mantle flow that causes some amount of west-directed basal drag of the Alboran lithosphere. This basal drag generates interplate compresional stresses compatible with the distribution of intermediate-depth earthquakes in western Alboran.

The origin and tectonic evolution of the Gibraltar Arc system is the result of a complex geodynamic evolution involving the convergence of the Eurasian and African plates and the dynamic impact of the Gibraltar slab. Although geologic and geophysical data collected in the last few years have increased our knowledge of the Gibraltar Arc region, it is still unclear which are the mechanical links between the Gibraltar slab and the past deformation of the overriding Alboran lithosphere, as well as to which degree this subduction system is presently active. In this study, we use 2D numerical modelling to investigate the impact of the Gibraltar slab dynamics on the deformation of the overriding Alboran lithosphere. Our model simulates a WE generic vertical section at an approximate latitude of 36°N and considers an initial setup at about Burdigalian times (∼20 Ma), when the subduction front position is relatively well constrained by recent tectonic reconstructions. Our modelling shows a switch in the overriding plate (OP) stress state from extensional stresses during the slab rollback to compressional stresses near the trench when the rollback velocity decreases, caused by the change in slab-induced mantle flow. We also find that much of the crustal and lithospheric deformation occur during fast slab rollback and OP extension in the first 10 Myr of evolution, while after that only moderate deformation associated with subduction is predicted. Finally, we find that despite the subduction rollback ceases, the ongoing motion of the deeper portion of the slab induces a mantle flow that causes some amount of west-directed basal drag of the Alboran lithosphere. This basal drag generates interplate compresional stresses compatible with the distribution of intermediate-depth earthquakes in western Alboran.

Introduction
The Gibraltar Arc system is located in the convergence zone of the Eurasian and African plates and is part of the tectonic framework of the Western Mediterranean region (Figure 1). It surrounds the Alboran Sea and comprises the Betic Chains in southern Spain and the Rif Mountains in northern Morocco. Several tomographic studies have imaged a subducted slab located under the Gibraltar Arc, Betics and Rif Mountains

Geodynamic context
Although in recent years slab rollback has been postulated as the most likely mechanism to explain the tectonic evolution of the westernmost Mediterranean, there is still debate about its spatialtemporal tectonic evolution. The proposed models can be reduced to three different geodynamic scenarios. Chertova et al. (2014) evaluate them using a 3D numerical modelling approach (see their paper for a further discussion of these three scenarios and their quantitative evaluation). The first scenario proposes that both the Alboran and Algerian basins opened from a single NW-dipping slab initially located south of the Balearic Islands (e.g., Rosenbaum et al., 2002;Spakman and Wortel, 2004;Chertova et al., 2014a;van Hinsbergen et al., 2014;van Hinsbergen et al., 2020). In this setting, subduction is thought to have started at ∼ 35 Ma due to gravitational instabilites and after a 5-10 Myr period of slab rollback to the south, the slab is segmented into two parts. The Betic-Rif part of the slab rotates clockwise about 100°and rolls back towards the Strait of Gibraltar while the other part reaches northeastern Algeria. The second proposed scenario is a subduction process that already started at 35 Ma from a long trench extending from the Gibraltar Arc to the Alps and where the rollback takes place as a continuous subduction front (e.g., Faccenna et al., 2004;Jolivet et al., 2009;Carminati et al., 2012). Finally, the last scenario suggests a model with two opposing subduction zones (e.g., Gelabert et al., 2002;Vergés and Fernàndez, 2012), recently modelled by Peral et al. (2022). In this setting, the Alboran domain is characterized by a slab that dipped to the southeast beneath Africa, and then rotated and retreated westward to the Gibraltar Arc. On the contrary, the Algerian domain is characterized in this scenario by a slab subducting to the northwest and retreating southeastward.
Apart from these scenarios, the recent kinematic reconstruction of the Mediterranean region by Angrand and Mouthereau (2021), which takes the advantage of reconciling previous western Mediterranean models (Vergés and Fernàndez, 2012;Romagny et al., 2020) with recent thermochronological analyses in western Betics (Daudet et al., 2020), also points to sublithospheric processes through slab retreat and delamination having dominated since Neogene (∼ 23 Ma). Regardless of the differences between these models, most reconstructions (e.g., Rosenbaum et al., 2002;Fernàndez et al., 2019;Romagny et al., 2020;Gómez de la Peña et al., 2021;Moragues et al., 2021) agree in that at ca. 20 Ma the trench had already fully rotated to the southwest and the predominantly westward rollback of the Gibraltar slab started taking place (Figure 2). This phase was followed by the cessation of the western migration of the subduction system at the latestmost Miocene (∼ 10-8 Ma), with the NW-SE convergence between Africa and Eurasia starting to dominate the evolution of the region (e.g., Gómez de la Peña et al., 2021). At this point, the deep part of the slab kept moving favoured by slab tearing in the northern edge (e.g., Faccenna et al., 2004;Spakman and Wortel, 2004;Garcia-Castellanos and Villasenor, 2011;Mancilla et al., 2015;Mancilla et al., 2018) followed by possible delamination processes beneath the eastern and central Betics (e.g., Duggen et al., 2005;Mancilla et al., 2013;Negredo et al., 2020), and steepening in the mantle (e.g., Faccenna et al., 2004;Spakman and Wortel, 2004;Garcia-Castellanos and Villasenor, 2011).
The 2D modelling of this work is focused on the last 20 Ma, characterized by the fast westward slab rollback followed by trench retreat cessation, to investigate the influence of subduction dynamics on the deformation of the overriding Alboran lithosphere. The 2D Frontiers in Earth Science 02 frontiersin.org

FIGURE 1
Topography of the Gibraltar Arc region and main geographic names. WAB, West Alboran Basin; EAB, East Alboran Basin; ABB, Algero-Balearic Basin. The shaded area indicates the positive velocity anomaly at 270 km depth found in the tomography study of Bezada et al. (2013) associated with the Gibraltar slab. Orange line denotes the Rif-Betic front (from Gutscher et al. (2012)) and purple line shows the position of the modelled cross-section in this work. The dashed portion of the profile is an extension of the model to prevent boundary effects. White arrows show the convergence direction of the Iberian and African (Nubia) plates. Blue dots indicates the regional intermediate seismicity (40 km < h < 150 km), from Instituto Andaluz de Geofisica Catalogue. The tomographic cross-section has been obtained from Hosseini et al. (2018) using the tomography model UU-P07 (Amaru, 2007).
Frontiers in Earth Science 03 frontiersin.org approach has the advantage of accounting for the evolution of the deformation of the overriding Alboran lithosphere, while extending the modelling to previous evolution stages clearly requires a 3D approach out of the scope of this paper. Moreover, this modelling does not accout for the regions presently affected by post-collisional slab tearing and/or contiental delamination, i.e. the central and eastern Betics.

Governing equations
We use version 2.4.0-pre of the finite-element code ASPECT (Advanced Solver for Problems in Earth's ConvecTion) (Kronbichler et al., 2012;Heister et al., 2017;Gassmöller et al., 2018;Bangerth et al., 2021a; to solve the coupled conservation equations for mass, momentum and energy that govern the behaviour of a viscous fluid and establish how the dependent variables (velocity, pressure and temperature) evolve with time. In this work we use the Boussinesq approximation, which assumes that density variations can be disregarded (incompressible fluid) in all equations except in the buoyancy term of the conservation of momentum. Adiabatic and shear heating are also neglected. With these assumptions, the conservation equations of mass, momentum, energy and compositional fields simplify as follows: where ρ is the density, u is the velocity field, μ is the viscosity, ε = 1 2 (∇u + ∇u T ) is the strain rate tensor, P is the pressure, g is the gravity acceleration, c p is the specific heat, T is the temperature, k is the thermal conductivity, H is the radiogenic heating and the quantities c i are the advected compositional fields that are used to track materials throughout the simulation. We assume that the density satisfies the following relationship: where α is the thermal expansion coefficient and ρ o is the density at a reference temperature T o .

Model setup and boundary conditions
We adopt a 2D box geometry (Figure 3) that simulates a vertical section at a latitude of about 36°N (purple line in Figure 1). The width of the domain is 1320 km, while the height is 660 km to simulate flow in the upper mantle. We use an adaptive mesh refinement (AMR) based on increasing the mesh resolution where strong contrasts on viscosity, density, temperature and/or composition occur. Since the areas that require high resolution change with time, the mesh is adapted every few time steps. The maximum local resolution is 644x644 m. The geometry of the initial setup has been built using the Geodynamic World Builder (GWB) version 0.4.0 (Fraters et al., 2019;Fraters, 2021) and it emulates the situation of the Gibraltar slab at ca. 20 Ma, when the trench strike was roughly N-S and its later migration was mainly westwards (Figure 2). The geometry consists of a 220 km-long subducted slab with a dip angle of 40°. This geometry has been selected after a number of tests to evaluate the influence of the initial geometry on model evolution (Supplementary Figures S5-S7 in the Supplementary Material). The initial configuration and the input thickness values (of initially homogeneous layers) are therefore chosen to best fit the amount of trench migration (∼ 300 km) and the present-day lithospheric structure on average, as well as first order characteristics of the slab geometry beneath Gibraltar. This setting allows for rollback initiation without imposing kinematic conditions or external forces at the lateral sides.
The subducting plate (SP) is assumed to be 120 km thick oceanic lithosphere. The thickness of 120 km is supported by recent observations using S-wave receiver functions (e.g., Molina-Aguilera et al., 2019). We impose a narrow weak zone of 10 km width (Figure 3) characterized by a low viscosity to decouple the subducting plate from the overriding plate. This layer extends to the west of the slab hinge up to x = 350 km and allows subduction and slab rollback to take place. At the top of the westernmost part of the domain, we use oceanic crust that has the same parameters as the weak zone but with a higher viscosity. With this viscous material we impose a slowing down in the subduction rate to account for the cessation of the western migration of the subduction system at the latestmost Miocene. The OP is completely continental and has an initial thickness of 80 km. It consists of 60 km of lithospheric mantle and 20 km of continental crust, where we distinguish between upper and lower crust (10 km thick each). As mentioned above, these initial thickness values have been selected after a number of tests (Supplementary Figure S6) as they best fit the present-day Alboran lithospheric structure and the amount of subduction migration. Finally, the initial trench position is x = 620 km and we use a tracer to track its position and velocity over time.
The initial temperature profile is characterized by a linear increase from 293 K at the surface to 1643 K at the lithosphere-asthenosphere boundary (LAB) and this value is constant in the rest of the mantle. We use Dirichlet temperature boundary conditions at top and bottom boundaries. Velocity boundary conditions are free slip for all boundaries, which means that the velocity perpendicular to the boundaries and shear stresses are imposed to be zero. We made some tests including a true free surface at the top boundary to evaluate the effect of this condition on the evolution of the model (see Section 4 in the Supplementary Material). From these tests we conclude that the use of a free-slip condition in the top boundary does not significantly change the results, while avoiding certain strong mesh deformations and associated numerical instabilities that arise in the free-surface case. Therefore, as we are not focusing on surface processes such as the development of an accretionary wedge or erosion/sedimentation, and after verifying that the free-slip condition leads to results similar to those with a free surface (Supplementary Figure S9), we prefer using the former condition on the top boundary.

Rheology
We use a visco-plastic rheology in which the viscous deformation is given by a composite viscosity that takes into account both diffusion creep (Newtonian rheology) and dislocation creep (powerlaw rheology): where μ diff and μ disl are the diffusion and dislocation viscosities. Both viscosities can be conveniently expressed with one equation but with different parameters (e.g., Hirth and Kohlstedt, 2003).
where the subindex i denotes the deformation mechanism (diffusion or dislocation), A is a prefactor of the equation, d is the grain size, m is the grain size exponent, n is the stress exponent,ε ‖ is the square root of the deviatoric strain tensor second invariant, E and V are an activation energy and activation volume respectively, p is the pressure and R is the gas constant. For the diffusion case (n = 1 and m ≠ 0) there is no dependence of the viscosity on the strain rate tensor, while in the dislocation case (n > 1 and m = 0) there is a power-law dependence. The plastic deformation is modelled by a plastic viscosity μ pl given by: where σ y is the yield stress and is given by the Drucker-Prager criterion: with C being the cohesion and ϕ the friction angle. To combine viscous and plastic deformations, we assume that they are separate processes and the mechanism resulting in the lowest effective viscosity μ eff is favoured: This final effective viscosity μ eff is limited by a predefined minimum viscosity μ min and maximum viscosity μ max to avoid large viscosity jumps and thus ensure the stability and convergence of the simulations. For the upper mantle (asthenosphere and lithosphere) we use a composite rheology from wet olivine (Hirth and Kohlstedt, 2003). These parameters give a composite viscosity of 10 20 Pa s at a depth of 150 km for a strain rate of 10 -15 s −1 . For the continental lower crust we adopt a diffusion rheology from wet anorthite feldspar (Bürgmann and Dresen, 2008). For the continental upper crust, the effective viscosity is given by the lesser between the plastic viscosity and a relatively low viscosity of 10 22 Pa s. In this way we account for a widely distributed deformation, as observed in the Alboran lithosphere. In terms of yielding, we consider a cohesion of 20 MPa and a friction angle of 20°. We do not account for strain softening and therefore cohesion and internal friction angle are kept constant, which inhibits the formation of shear zones. Since we include an initially homogeneous OP crust and lithosphere, including strain softening would require introducing an initial heterogeneity to seed extension (e.g., Huismans et al., 2005). However, in a natural setting, it is more likely that deformation was focused on inherited structures as crustal faults. Since these structures are poorly constrained in our study and locating an initial seed would be arbitrary, we have preferred modelling a distributed pure shear deformation given by a combination of plastic yielding and viscous deformation. Finally, the oceanic crust and weak layer viscosities are chosen to be constant since Frontiers in Earth Science 05 frontiersin.org these are tuning parameters that allow us to reproduce the geologically constrained trench retreat (Table 1).

Results
We have run a large number of experiments to test the influence of the rheological parameters and the initial model geometry (see Supplementary Figures S5-S7). We present here the model that provides the most realistic geodynamic evolution.

Subduction dynamics
The slab initially sinks into the upper mantle due to its negative buoyancy, with the trench retreating westwards and the slab reaching a maximum dip angle of 90°after 3 Myr ( Figure 4A). The slab pull force increases as the slab keeps sinking, increasing the subduction rate and the trench retreat velocity (V T ) up to a maximum of ∼ 7 cm/yr ( Figure 4F). At this time, the high trench retreat velocities result in a slight reduction of the slab dip angle ( Figure 4B). When the slab tip begins to feel the viscosity increase in the upper mantle (at depth of ∼ 400 km, Figure 4C), V T begins to decrease and drops quickly to zero ( Figure 4F). This reduction in V T is also triggered and accelerated by the interaction of the high viscous oceanic crust with the trench. From 10 Myr of evolution onwards, the trench is stagnated and there is no more rollback, although the slab is still sinking in the upper mantle and interacting with the 660-km-depth upper-lower mantle discontinuity (Figures 4D, E and Supplementary Video S1). This viscous interaction leads to the slab flattening horizontally and generates two mantle poloidal flow cells: one clock-wise behind the slab and another counter-clockwise in front of the slab (Figures 4D, E). Additional figures in Supplementary Material (Supplementary Figures S1-S4) illustrate the effect of subduction dynamics on temperature, density, strain rate and horizontal velocity distribution for this model.

Overriding plate deformation
In order to analyse the stress state and deformation of the OP, we plot the horizontal deviatoric normal stresses τ xx (Figure 5). A convention of positive compressive stress is followed. Additionally, we plot the contours of the lithosphere (red) and the crust of the OP (blue) to show their evolution over time.
During the phase of fast slab rollback, extensional stresses spread over the entire OP (negative values in Figure 5A). The highest stresses (> 100 MPa) are found in the OP and in the bending region of the SP at depths > 40 km, where there is horizontal compression. The extension of the OP is caused by the westward trench retreating, and by the basal traction induced by the mantle flow, which increases westwards beneath the OP (Figure 5A and Supplementary Figure S4A) (e.g., Holt et al., 2015;Alsaif et al., 2020). During this period, the OP undergoes roughly uniform thinning and the crust thins up to ∼ 18 km, except in its westernmost part, which thickens to a value of ∼ 26 km due to the downward slab motion (Figure 5A). After 6 Myr of model evolution, compressional stresses develop in a narrow region (∼ 20 km wide) in front of the trench (Figure 5B). At this time, the  crust has thinned to ∼ 15 km to the east and thickened up to ∼ 30 km close to the subduction front. The lithosphere thins uniformly up to ∼ 70 km, whereas near the OP trailing edge it thins to ∼ 55 km. Once the trench retreat starts to decrease, the extension near the trench turns into compressional stresses of ∼ 100 MPa and the extensional stresses are reduced to ∼ 60-80 MPa in the easternmost part of the model ( Figure 5C). There is also a small neutral region (∼ 300 km wide) in the OP at ∼ 180 km from the trench with τ xx ≈ 0. This change in the OP stresses is caused by the basal traction produced by the poloidal flow that drives the OP towards the SP at a faster rate that the trench retreats ( Figure 5C and Supplementary Figure S4C) (e.g., Duarte et al., 2013;Schellart and Moresi, 2013;Holt et al., 2015). From 10 Myr of evolution onwards, only moderate crustal and lithospheric deformation is predicted. The lithosphere slightly thickens due to thermal relaxation while only crustal thickening is predicted near the subduction front (x ∼ 390 km). Additionally, compression near the trench and extension in the eastern part of the OP coexist (Figures 5D, E and Supplementary Video S2). Figure 6A shows the crustal and lithosphere structure resulting from our modelling after 20 Myr of evolution (∼ present-day). We obtain a practically constant Moho depth in the East Alboran Basin (EAB) and West Alboran Basin (WAB) (∼ 15 km), which thickens rapidly towards the Gibraltar Arc (> 30 km, Figure 6A). This abrupt transition occurs over a distance of ∼ 50 km. The maximum crustal thickness of ∼ 35 km is predicted in the westernmost part of the WAB, just above the position of the slab. As the oceanic crust is dragged down with the subducting slab, Moho depths increase from 10 km in the Gulf of Cadiz to > 40 km towards the Gibraltar Strait ( Figure 6A).

Crustal and lithospheric structure
Our modelling predicts large variations in the lithosphere structure. Shallow depths of the LAB (∼ 80-90 km) are obtained for the EAB lithosphere. The LAB gradually deepens to the WAB (∼ 100 km), where the slab subducts almost vertically towards the east (Figure 6A). In the Gulf of Cadiz, we obtain a LAB depth of ∼ 100-110 km which increases sharply towards the Strait of Gibraltar, reaching a depth of ∼ 160 km at 5.5°W due to the influence of the slab. Finally, we obtain an incipient slab necking deformation at ∼ 160 km depth ( Figure 6B).
Heat flow data have been collected from a large number of works (e.g., Polyak et al., 1996;Fernàndez et al., 1998, and references therein) (Figure 6C)

Moho and LAB depths
The Moho depth obtained by our modelling in the EAB is mostly in agreement with the crustal structure given by Fullea et al. (2010) based on 3D modelling of the Atlantic-Mediterranean region, while at longitudes of 3°-4°W we obtain shallower crustal depths compared to those given by Fullea et al. (2010) and Torne et al. (2000) from integrated modelling of gravity, elevation and heat flow data ( Figure 6A). The general pattern is a steep thinning of the crust from the Gibraltar Strait and the westernmost part of the WAB towards the WAB and a constant Moho depth in the rest of the Alboran basin. The compilation of Moho depths by Diaz et al. (2016) also shows large lateral Moho depth gradients, with values ranging from ∼ 35-40 km in the Gibraltar Strait to ∼ 15-20 km in the Alboran domain. This trend is also supported by Palomeras et al. (2017) from the data derived from surface wave tomography. Our results show a maximum crustal thickness of ∼ 35 km in the westernmost part of the WAB, which is slightly less than the receiver function (RF) estimates (> 40 km) (Mancilla and Diaz, 2015) (Figure 6A). To the Gulf of Cadiz, our model predicts a Moho depth of ∼ 10 km, which is consistent with a 7 km thick oceanic crust overlain by a 1-3 km thick layer of Mesozoic to Eocene (roughly 251-33.9 Ma) sediments found by Sallarès et al. (2011) in the central Gulf of Cadiz. The Moho depths of 25 km obtained by Fullea et al. (2010) are much deeper than those found in the present study (note that we are not including the water and sediments column). These differences are likely due to the fact that we are not accounting for the flexural bending of the SP, nor for the development of the accretionary wedge in the Gulf of Cadiz, where sediment thickness reaches 8-10 km (e.g., Gutscher et al., 2012).
The LAB structure obtained in this work is compared with the previous studies of Torne et al. (2000) and Fullea et al. (2010), and with the structure obtained using S-wave receiver functions by Molina-Aguilera et al. (2019) (Figure 6A). The main LAB trends are similar in all these works, with a deeper LAB in the WAB than in the EAB. In the easternmost part of the Alboran domain, we obtain a lithosphere of ∼ 80 km thickness more in agreement with Fullea et al. (2010) Fernàndez et al. (2004) obtained a thin lithosphere of ∼ 90 km near the Gulf of Cadiz, more consistently with our results. Discrepancies in the depth of the LAB result from differences in the methodology and dataset of each work, but the overall trend is consistent. The shape of the slab structure obtained in this work is also in close agreement with the slab geometry imaged by different tomographic studies (e.g.; Amaru, 2007;Bezada et al., 2013;Villasenor et al., 2015;Civiero et al., 2020) ( Figure 6A) Regarding the seismicity, Figure 6B indicates that intermediatedepth earthquakes likely occur in an area of compressional stress accumulation related to the interaction between the subducting slab and the Alboran lithosphere. It is worth noting that intermediatedepth earthquakes occur in an area of relatively cold lithosphere due to the down pull of the slab. This cold geotherm is needed to have brittle deformation instead of ductile aseismic creep, and well agrees with tomographic studies showing that this seismicity concentrates in areas of high seismic velocity (Monna et al., 2013;Villasenor et al., 2015;Sun and Bezada, 2020). Santos-Bueno et al. (2019) obtained focal mechanism of the intermediate seismicity finding that most of them present reverse faulting solutions more abundant in the southern sector, beneath the central Alboran Sea and Northern Morocco. Their P-axis show similar NE-SW orientations, in contrast with the NW-SE convergence direction between Africa and Eurasia plates. They attributed the plausible cause of this orientation to horizontal traction that may originate from a sublithopheric process. The westward basal drag of the Alboran lithosphere found in the present study would be consistent with this interpretation (Figure 5). Our study supports that this seismicity is likely related to interplate compression processes rather than to shallow slab necking proposed by Sun and Bezada (2020).
Our results show that the slab pull and mantle flow interacting with the overriding lithosphere are first-order drivers of crustal and lithospheric deformation. We find simultaneous crustal thickening near Gibraltar and thinning further to the east during the first 10 Myr of evolution. Our model also shows that the main crustal and lithospheric deformation occurs during slab rollback and OP extension in this first 10 Myr, while after that only moderate deformation (some additional crustal thickening close to the subduction front and slightly lithospheric thickening) associated with subduction is predicted (Figure 5). However, it is important to note that our simulations only consider crustal deformation due to trench-perpendicular plate motions and do not consider other features such as the NW-SE relative plate convergence between the African and Eurasian plates that may play a role and are necessary to explain certain properties as the regional trends in the stress field or the present-day velocity fields (e.g., Jiménez-Munt and Negredo, 2003;Pérouse et al., 2012;Neres et al., 2016;Spakman et al., 2018). Additionally, the 2D nature of our numerical modelling does not take into account the resistance of the mantle to the N-NNE slab dragging caused by the absolute surface motion of the African plate (Spakman et al., 2018). This resistance to slab dragging has been proven to affect both crustal and lithospheric deformation and satisfactorily explains some of the main tectonic features in the western Mediterranean in the last 8 Ma, such as shortening in the Rif and extension in the eastern Betics (Spakman et al., 2018). Linking our results with these previous works, slab pull and slab rollback seem to have been the main drivers of lithospheric and crustal deformation prior to the stagnation of the slab in the Strait of Gibraltar, producing most of the current structure, while the relative NW-SE plate convergence and the mantle resistance to slab dragging, have dominated in the last ∼ 8 Ma, without introducing major changes in the global crustal and lithospheric structure. Figure 7 summarises the two stages of the OP deformation based on our results. During the fast sinking phase and westward rollback, there are extensional stresses throughout the entire OP, which contributes to the Alboran basin formation. These stresses are generated by the horizontal basal traction at the base of the OP increasing towards the trench (Figure 5A; Figure 7A and Supplementary Figure S4A). Our results agree with earlier geodynamic studies that found overall OP extension during slab rollback and trench retreat (e.g., Duarte et al., 2013;Chen et al., 2015;Chen et al., 2016;Xue et al., 2022). Additionally, although previous 3D subduction models have shown the importance of toroidal mantle flow during the first phases of subduction (e.g., Piromallo et al., 2006;Schellart, 2008;Stegman et al., 2010) and have suggested that toroidal flow is the primary contributor to OP extension (e.g., Duarte et al., 2013;Meyer and Schellart, 2013;Schellart and Moresi, 2013;Chen et al., 2016), the poloidal mantle flow beneath the OP associated with slab rollback is still large enough to generate extensional stresses. It is worth noting that fixing the lateral boundary of the OP has an impact on its deformation. Previous works have revealed that a fixed OP will undergo higher extensional stresses (e.g., Capitanio et al., 2010a;Chen et al., 2015;Holt et al., 2015). This boundary condition is usually applied in modelling efforts in the Mediterranean region, where "land-locked" basins surrounded by orogens show significant back-arc extension and thinning related to slab rollback (e.g., Royden and Faccenna, 2018).

Alboran domain deformation
The cessation of slab rollback after 10 Myr of evolution results in a change in the style of the mantle flow beneath the OP (Figure 5C; Figure 7B and Supplementary Figure S4C), which drags the OP towards the trench at a faster rate than the trench retreats, causing compressional stresses in the subduction front. The width of the extensional zone in the OP decreases as the mantle flow does not produce the necessary gradients to create sufficient extensional stresses (e.g., Holt et al., 2015).
This pattern of both compression near the trench and extension away from it has already been observed in previous modelling studies, both in 3-D and 2-D (e.g., Capitanio et al., 2010a;Nakakuki et al., 2010;Schellart and Moresi, 2013;Holt et al., 2015;Alsaif et al., 2020). For example, models from Holt et al. (2015) showed that the OP stress state is time-dependent and tends to a more compressive regime when the trench retreat rate is reduced as the slab sinks into the viscous lower mantle, just as obtained in this work.
The coexistence of compressional and extensional stresses at lithospheric levels in the OP is in agreement with different studies of focal mechanisms in the Alboran Sea (e.g., Stich et al., 2006;Martín et al., 2015). The overall observed stress pattern is attributed primarily to two sources: the ongoing NW-SE oblique convergence of the Eurasian and African plates and a secondary stress source that we propose is likely related to the westward basal drag of the Alboran lithosphere. Our simulations suggest that the period of westward slab rollback followed by an abrupt deceleration has had, and still has, a great impact on the deformation pattern of the region.

Is the subduction system presently active?
The question of whether the Western Mediterranean subduction system is still active or not has been a subject of controversy for the last two decades. Based on a marine seismic survey in the Strait of Gibraltar, Gutscher et al. (2002) proposed an active eastward-dipping subduction zone beneath Gibraltar. Later on, Gutscher et al. (2012) clarified that the term "active subduction" is used to specify an ongoing-rollback activity and not in the sense of subduction resulting from horizontal plate motion. They suggested that the subduction does not seem to have completely ceased due to the currently 0.4-0.5 cm/yr WSW motion observed by GPS in the westernmost Rif and the Gibraltar Arc (e.g., Fadil et al., 2006;Vernant et al., 2010;Koulali et al., 2011;Palano et al., 2015). However, several authors claimed that subduction is not presently active on the basis of the absence of a Wadati-Benioff zone (e.g., Platt and Houseman, 2004), on the post-Tortonian deformation of the Gulf of Cadiz Imbricate Wedge being only due to NW-SE Africa-Iberia convergence (Iribarren et al., 2007), on the lack of seismicity in the accretionary wedge, which is sealed by an nearly undeformed package of sediments deposited since the Early Pliocene (Zitellini et al., 2009), and on the lack of significant relative motion between GPS stations located on both sides of the subduction front (Stich et al., 2006). More recent works further point to the inactivity of the subduction system. For example, Santos-Bueno et al. (2019) computed the focal mechanisms for 42 earthquakes at depths from 50 to 100 km beneath the Gibraltar Arc without a relation with downdip stresses, suggesting that there is no signature of ongoing subduction. Based on a new GPS dataset, Civiero et al. (2020) proposed that the subduction system is no longer active as a continuous arc and it has slowed down or stopped. However, they suggest that the slab may be still moving in the mantle and consequently causing material to flow around the slab.
As suggested by Civiero et al. (2020), our simulations show that although the subduction rollback has ceased, the slab is still deforming and sinking in the upper mantle as a result of its negative buoyancy. The deep movement of the slab generates an induced mantle flow that produces a basal traction beneath the OP directed towards the slab (Figure 4E, Figure 7B). The basal traction naturally obtained in our models was required in the neotectonic modelling by Neres et al. (2016) to better reproduce the geodetic data, stress directions and seismicity of the area. According to our understanding, this slabinduced mantle flow is still contributing to the present-day kinematics of the Alboran domain and the Gibraltar Arc region. This idea agrees with previous work suggesting that sub-crustal or sub-lithospheric processes are necessary to drive the surface velocities in the Gibraltar Arc region (Fadil et al., 2006;Pérouse et al., 2010;Baratin et al., 2016;Neres et al., 2016). Therefore, although the subduction process seems to have stopped, the slab is still moving in the upper mantle, which has an impact on the tectonics and kinematics of the region. This terminal phase of subduction is likely related to incipient slab necking. In turn, this necking may prevent stress transfer between the shallow and deep portions of the slab, which is consistent with the absence of seismicity at depths greater than ∼ 150 km. This is in contrast, for example, to the Calabrian subduction system, which comprises a Wadati-Benioff zone about 500 km deep (e.g., Anderson and Jackson, 1987;Chiarabba et al., 2005).

Model limitations
The geometry of the subduction zones has a three-dimensional nature and would therefore require 3D models to perform comprehensive and accurate modelling. Our 2D approach can not handle three-dimensional effects such as the toroidal flow or the propagation of slab pull along the trench direction. Toroidal mantle flow is particularly important during the first stages when the slab sinks into the uppermost mantle and would lead to greater extension of the overriding plate (e.g., Funiciello et al., 2003;Duarte et al., 2013;Meyer and Schellart, 2013;Schellart and Moresi, 2013), while the poloidal flow dominates when the slab tip reaches the more viscous layers of the deep mantle (Funiciello et al., 2003). We hypothesise that predominant toroidal flow during the first stages of slab migration is consistent with crustal extension initiating first in the WAB (Gómez de la Peña et al., 2021). Instead, later dominating poloidal flow could be responsible for the later development of the EAB. This independent evolution of the WAB and EAB through the latemost Oligocene and Miocene (26.5-5.3 Ma) proposed by Gómez de la Peña et al. (2021) is not properly reproduced by our 2D modelling where OP extension is more distributed. Thus, our two-dimensional models are underestimating the extension of the overriding plate. These limitations should be considered when interpreting the results. For example this modelling reproduces a thin continental crust for the EAB, consistently with Medaouri et al. (2014) but the alternative interpretation as magmatic arc crust (e.g.; Booth-Rea et al., 2007;Gómez de la Peña et al., 2018) cannot been discarded due to the mentioned underestimation of extension. This possible underestimation of stretching in the eastern part of the model is also compatible with model-predicted surface heat flow underestimating observations by ∼ 20 mW/m 2 ( Figure 6C).
Other features that cannot be reproduced in our modelling are the effects of NW-SE convergence between Iberia and Africa, which became dominant from 8 Ma, when rollback stalled. Not including this oblique convergence during the phase of fast westards slab rollback should have little influence on the results since the WE component of this velocity (< 0.5 cm/yr) is very small compared to trench retreat velocities up to 7 cm/yr during this phase (Figure 4F). Similarly, the N to NNE dragging of the Gibraltar slab by the African Plate absolute motion (Spakman et al., 2018) cannot be accounted for in this modelling. On the other hand, 3D numerical models are so challenging that are commonly oversimplified and cannot account for a detailed description of the deformation of the OP at crustal levels, which can be done more accurately in 2D models. Future 3D time-evolving models will be able to account for 3D flow and OP deformation in the Betics and Rif since the Miocene.

Conclusion
In this work we have studied, by means of 2D numerical modelling, the impact of the Gibraltar slab on the Alboran domain deformation during the last 20 Myr. Our modelling has also allowed us to gain insights and understanding about the crustal and lithospheric deformation over time, the origin of intermediatedepth seismicity and the ongoing debate on present-day activity of the subduction system. Our results lead to the following conclusions.
• The period of trench retreat is characterized by overriding plate extension, while the cessation of the slab rollback at ca. 10 Ma produces a change in the induced mantle flow that results in compressional stresses near the trench and extension to the east of the overriding Alboran lithosphere. • Much of the crustal and lithospheric structure was shaped during slab rollback and overriding plate extension in the first 10 Myr of evolution, while only minor deformation associated with subduction is predicted after the cessation of slab rollback. Furthermore, this modelling properly reproduces simultaneous crustal thickening near Gibraltar and thinning in the rest of the Alboran basin. • Despite the cessation of subduction rollback, the slab is still sinking in the upper mantle as a result of its negative buoyancy. This sinking is generating an induced mantle flow that produces a west-directed basal shear traction beneath the Alboran domain, as proposed in previous neotectonic studies. • This terminal phase of subduction is likely leading to an incipient slab necking. This necking may prevent stress transfer between the deep and shallow parts of the slab, which is in accordance with the absence of earthquakes deeper than 150 km in the region.

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 in the article/Supplementary Material.

Author contributions
All authors designed the numerical experiments. PG wrote most of the original draft. FM and AN supervised research and reviewed and edited the manuscript. All authors participated in the discussion of the results.

Funding
The present research has been funded by Spanish Ministry of Science and Innovation projects PID 2019-109608GB-I00 (FM and PG), PGC 2018-095154-B-I00 and PID 2020-114854GB-C22 (AN) and FEDERJunta de Andalucia-Conserjeria de Economia y Conocimiento/B-RNM-528-UGR20. PG is supported by the European Social Fund through a youth guarantee contract.