# Mantle Dynamics in the Mediterranean and Plate Motion of the Adriatic Microplate: Insights From 3D Thermomechanical Modeling

Christian Schuler, Boris J P Kaus, Eline Le Breton, Nicolas Riel, Anton A Popov

► To cite this version:

Christian Schuler, Boris J P Kaus, Eline Le Breton, Nicolas Riel, Anton A Popov. Mantle Dynamics in the Mediterranean and Plate Motion of the Adriatic Microplate: Insights From 3D Thermomechanical Modeling. Geochemistry, Geophysics, Geosystems, 2025, 26 (3), pp.e2024GC011996. 10.1029/2024GC011996. [https://hal.science/hal-04984023v1](https://hal.science/hal-04984023v1)

HAL Id: hal-04984023

Submitted on 10 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Geochemistry, Geophysics, Geosystems®

## RESEARCH ARTICLE

10.1029/2024GC011996

**Special Collection:**

Alpine mountain belts in 4-dimensions

**Key Points:**

* Systematic 3D thermomechanical simulations of the Mediterranean incorporating Africa-Eurasia convergence and all major subduction zones

* Asthenospheric mantle dynamics in the Mediterranean is dominated by toroidal flow caused by complex subduction dynamics

* Adria's motion depends on Africa-Eurasia convergence, Alpine slab retreat and the distance between the Calabrian and Hellenic trenches

**Supporting Information:**

Supporting Information may be found in the online version of this article.

**Correspondence to:**

C. Schuler, c.schuler@uni-mainz.de

**Citation:**

Schuler, C., Kaus, B. J. P., Le Breton, E., Riel, N., & Popov, A. A. (2025). Mantle dynamics in the mediterranean and plate motion of the Adriatic microplate: Insights from 3D thermomechanical modeling. *Geochemistry, Geophysics, Geosystems*, 26, e2024GC011996. https://doi.org/10.1029/2024GC011996

Received 31 OCT 2024
Accepted 10 FEB 2025

**Author Contributions:**

**Conceptualization:** Christian Schuler, Boris J. P. Kaus, Eline Le Breton, Nicolas Riel

**Formal analysis:** Christian Schuler

**Funding acquisition:** Boris J. P. Kaus, Eline Le Breton

**Investigation:** Christian Schuler

© 2025 The Author(s). Geochemistry, Geophysics, Geosystems published by Wiley Periodicals LLC on behalf of American Geophysical Union. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

<figure><img src="image_2.png" /></figure>

## Mantle Dynamics in the Mediterranean and Plate Motion of the Adriatic Microplate: Insights From 3D Thermomechanical Modeling

Christian Schuler¹,*, Boris J. P. Kaus¹,*, Eline Le Breton², Nicolas Riel¹, and Anton A. Popov¹

¹Institute for Geosciences, Johannes Gutenberg-Universität Mainz, Mainz, Germany, ²Géosciences Rennes - UMR 6118, University Rennes, CNRS, Rennes, France

**Abstract** The motion of the Adriatic microplate is thought to be highly sensitive to the surrounding subduction zones and the convergence of Africa and Eurasia. However, our understanding of the mantle dynamics in the Mediterranean region and its effect on plate motion remains incomplete. Here, we present a large set of 3D thermomechanical models of the entire Mediterranean region over the last 35 Myr to understand what controls the motion of the Adriatic microplate. The simulations take the convergence of the African and Arabian plates with the Eurasian plate into account, along with the dynamics of the subduction systems in the western (Apennines-Calabria), central (Dinarides-Hellenides) Mediterranean and in the Alpine-Carpathian region. Our results demonstrate that the subduction systems around Adria are highly coupled, which gives rise to complex asthenospheric flow in the central Mediterranean. We find that the plate motion of the Adriatic microplate over the last 35 Myr is controlled by the interplay of three main factors: (a) the convergence between the African and Eurasian plates, (b) the retreat of the Alpine subduction zone to the north of Adria, and (c) the distance between the Calabrian and Hellenic subduction zones around Adria. Furthermore, in a system characterized by active convergence between Africa and Eurasia, the slab pull exerted by nearby subduction zones can only notably influence the motion of the Adriatic microplate if these subduction zones are located within a few hundred kilometers of Adria.

**Plain Language Summary** Over the last 35 Myr, different geological processes were active in the Mediterranean. The African and Arabian tectonic plates collided with the Eurasian plate, and as a result of this convergence, different subduction zones were active around the Adriatic microplate, which is situated in the middle of this complex geological system. This study presents 3D thermomechanical simulations of the Mediterranean region over the last 35 Myr. This method lets us simulate how rocks behave and deform over large time and space scales. We look at how different subduction zones affect the motion of the Adriatic microplate. The subduction zones around Adria have opposite polarity and have moved closer together over time, resulting in a complex flow of the subjacent mantle. Our results show that the plate motion of Adria over the last 35 Myr is mainly affected by the African plate pushing it from the south, the retreat of the Alpine subduction zone to the north and of the western Mediterranean and Dinaric-Hellenic subductions to the west and east of Adria. The subduction systems around Adria can notably influence its plate motion only if they are located within a few hundred kilometers of Adria.

## 1. Introduction

The geological evolution of the Mediterranean region is a complex interplay of continental terranes and subduction systems in both time and space. Understanding the interacting processes between different subduction systems, such as in the central Mediterranean with both the Calabrian and Hellenic active subduction zones, remains the subject of an ongoing scientific debate (e.g., Chertova et al., 2014; Faccenna et al., 2001, 2014; Király et al., 2021; Royden & Faccenna, 2018; Wortel & Spakman, 1992, 2000).

The last 35 Myr are of particular interest, as this period approximately marks the onset of collision in the Alps and of the rollback subduction of the Apennines-Maghrebides and the Hellenides, both accompanied by upper plate extension in the western Mediterranean and in the Aegean, respectively (e.g., Faccenna et al., 2014; Le Breton et al., 2021; van Hinsbergen et al., 2020). Consequently, 35 Ma serves as an appropriate starting point to investigate the interactions between the retreating subduction zones in the Mediterranean and the associated mantle dynamics. The Adriatic microplate plays a key role in this setting due to its location between the
**Methodology:** Christian Schuler, Boris

J. P. Kaus, Nicolas Riel

**Software:** Anton A. Popov

**Supervision:** Boris J. P. Kaus

**Validation:** Christian Schuler, Boris

J. P. Kaus, Eline Le Breton

**Visualization:** Christian Schuler

**Writing – original draft:**

Christian Schuler

**Writing – review & editing:** Boris

J. P. Kaus, Eline Le Breton, Nicolas Riel

converging European and African plates to the north and south, respectively (Figure 1). Africa has been converging with Eurasia since the Cretaceous (e.g., Dewey et al., 1989) and today, the Adriatic microplate is almost entirely surrounded by orogenic belts, with the retreating Calabrian and Hellenic subduction zones to the south, the Apennines to the west, the Dinarides to the east and the Alps to the north (Figure 1). As a result, the Adriatic microplate is highly sensitive to the various processes occurring within the Mediterranean (e.g., Király, Holt, et al., 2018; Le Breton et al., 2017; van Hinsbergen et al., 2014). It was long assumed that Adria acted as a promontory of Africa (e.g., McKenzie, 1972; Muttoni et al., 2001). However, geodetic data show that the present-day motion of Adria is clearly different from that of Africa (e.g., D'agostino et al., 2008), with a counterclockwise rotation with respect to Europe and Africa that may have started after approximately 20 Ma (e.g., Le Breton et al., 2017; Rosenbaum et al., 2004; Ustaszewski et al., 2008).

A common approach to explore the geodynamic history of the Mediterranean is to reconstruct the kinematic motion of the plates involved, based on a compilation of geological (structural, metamorphic, magmatic, stratigraphic) and geophysical data, such as seismic tomography models of the mantle (e.g., Handy et al., 2010, 2015; Le Breton et al., 2021; Menant et al., 2016; Rosenbaum et al., 2002b; Schettino & Turco, 2011; Stampfli & Kozur, 2006; van Hinsbergen et al., 2020). These reconstructions provide valuable information about plate motion over geological timescales and how it compares with the structure of the mantle today. However, it can be challenging to make quantitative assumptions about the mantle processes that caused the plate motion and the present-day mantle structure. Thermomechanical modeling attempts to narrow this gap by using physically consistent models that can help understanding the driving forces and their interactions through geological time.

Previous work on the interaction of different mantle processes using 3D numerical and analog modeling has mainly focused on a limited number of pre-selected processes. Király, Holt, et al. (2018) use analog models to investigate the effect of a double-sided subduction zone on plate motion, inspired by the Adriatic microplate. Holt et al. (2017) investigate double-sided subduction zones with numerical modeling. Confal et al. (2018) and Lo Bue et al. (2021) use numerical modeling to study the evolution of the Hellenic slab and the Calabrian slab, respectively. Martinod et al. (2024) use analog models to analyze the influence of slab break-off in the western Mediterranean on Adria's motion and the Alpine dynamics. Duarte et al. (2024) investigate the evolution of the Gibraltar subduction zone and Peral et al. (2022) use numerical modeling to test the hypothesis of opposing subduction systems in the western Mediterranean. Efforts have also been made to investigate the current structure of the upper mantle using instantaneous models (e.g., Faccenna & Becker, 2010; Glerum et al., 2021; Kumar et al., 2022). Other analog modeling studies focus specifically on smaller-scale structures such as slab gabs (Király, Holt, et al., 2018) or detachment processes (Boutelier & Cruden, 2017), as these can be related to the Apennine-Calabrian subduction system and the Alpine subduction zone, respectively.

Nonetheless, a common feature of all these studies is that they focus on rather regional processes while not necessarily taking into account all concurrent processes such as the convergence of Africa and Eurasia or the surrounding subduction zones. Here, we aim to include all the main driving forces in the Mediterranean region, namely the pushing African and Arabian plates (w.r.t. a fixed Eurasian continent), as well as the subduction systems in the western Mediterranean, the eastern Mediterranean and in the Alpine and Carpathian regions. We created 3D models of the Mediterranean region and carried out thermomechanical forward simulations to study the evolution of this area over the last 35 Myr. The dynamics of the subduction zones is predominantly influenced by internal forces, including slab pull forces, the opening of rift systems and rigid plate interactions at the surface. The motion of the African continent is prescribed at the southern boundary and exerts a far-field force on the subduction dynamics.

In this study, we develop a reference model to reproduce the first-order mechanisms governing the Mediterranean region. We perform a parameter analysis to investigate the driving forces behind subduction dynamics and plate motion of the Mediterranean. Special emphasis is placed on the interaction between the various subduction systems, which remains a topic of debate (e.g., Király et al., 2021), as well as on the motion of the Adriatic microplate. Our aim is to identify the primary drivers of Adria's plate motion and establish a connection to its present-day dynamics.

## 2. Tectonic Evolution in the Mediterranean

The Mediterranean region exhibited pronounced tectonic activity throughout its geological history, which was investigated by various kinematic reconstructions (e.g., Handy et al., 2010; Jolivet & Faccenna, 2000; Le Breton
<figure><img src="image_3.png" /><figcaption>**Figure 1.** (a–f) Kinematic reconstruction of the Mediterranean over the last 35 Myr, modified after Le Breton et al. (2021). Beige color shows continental crust, orange color shows thinned continental margin and blue color shows oceanic crust. Back in time, orange also includes crustal domains deformed during Alpine Orogeny. Black lines with triangles show trench position of the various subduction systems. The Adriatic microplate is bound by black lines. LPB stands for Liguro-Provençal Basin, AB for Algerian Basin, Al for Alboran, Ka for Kabylides, Pe for Peloritani, Ca for Calabria, Co for Corsica and Sa for Sardinia.</figcaption></figure>

et al., 2021; Menant et al., 2016; van Hinsbergen et al., 2020). Since approximately 83 Ma, Africa converged with respect to Eurasia and the south-dipping Alpine subduction was initiated, which consumed the Piemont-Ligurian (PL) ocean until about 35 Ma. To the east of the Adriatic microplate, the Hellenic-Dinaric subduction retreated and consumed the Vardar oceanic crust. Subsequently, the retreat of the roll-back slowed down, as thin continental crust was subducted (e.g., Menant et al., 2016).

At approximately 35 Ma, Adria collided with the European continent (e.g., Le Breton et al., 2021; van Hinsbergen et al., 2020). This phase of indentation was accompanied by lateral extrusion toward the east, fast rollback of the Carpathians and extension of the Pannonian Basin (Figure 1; see also e.g., Faccenna et al., 2014; Handy
et al., 2015). From 35 Ma, onset of fast roll-back subduction triggered upper plate extension and the opening of the western Mediterranean basins and extension in the Rhodopes and Aegean Sea (Figures 1a–1c; see also e.g., Faccenna et al., 2014). The subduction system in the western Mediterranean retreated toward Africa. After reaching the African continent at around 15 Ma, the trench retreated both toward Gibraltar and toward Adria (Figures 1d–1f; see also e.g., Rosenbaum et al., 2002a; Spakman & Wortel, 2004). The retreat of this subduction system is reflected in the motion of continental blocks that dock in the Alboran sea and the Kabylies between 15 Ma and 10 Ma (Figures 1d and 1e: see also e.g., Rosenbaum et al., 2002a). East of Adria, the Hellenic trench retreat accelerated since 15 Ma due to the subduction of old Ionian oceanic crust (Figures 1d–1f; see also e.g., Menant et al., 2016).

The motion of the Adriatic microplate before 20 Ma was mainly influenced by the African continent which induced a northwards motion of Adria relative to Europe (e.g., Le Breton et al., 2021; Ustaszewski et al., 2008). After that, the interplay between the advancing African continent, the indentation in the Alpine orogen to the north and the retreating subduction zones both west and east of Adria must have acted on the counterclockwise rotation and independent motion of Adria relative to Africa (e.g., Le Breton et al., 2017; van Hinsbergen et al., 2014). In this study, we aim to test this hypothesis with thermomechanical modeling to learn more on the forces driving the motion of Adria.

# 3. Method

## 3.1. Conservation Equations

To account for lithospheric-scale deformation over long timescales, the incompressible Stokes equations are solved, where the momentum and mass conservation equations are formulated as:

$$ \frac{\partial \tau_{ij}}{\partial x_j} - \frac{\partial P}{\partial x_i} + \rho g_i = 0, \qquad (1) $$

$$ \frac{\partial v_i}{\partial x_i} = 0, \qquad (2) $$

where $x_i$ ($i = 1,2,3$) denote spatial coordinates, $\tau_{ij}$ is the deviatoric stress tensor, $P$ is the pressure, $\rho$ is the density, $g_i$ is the gravitational acceleration vector and $v_i$ is the velocity in each direction.

We employed a simplified density parametrization, in which the density $\rho$ depends on temperature $T$, reference temperature $T_0$, thermal expansivity $\alpha$ and reference density $\rho_0$:

$$ \rho = \rho_0(1 - \alpha(T - T_0)). \qquad (3) $$

Simultaneously with the mass and momentum conservation equations we also solve the energy conservation equation which is given by:

$$ \rho C_p \frac{DT}{Dt} = \frac{\partial}{\partial x_i} \left( k \frac{\partial T}{\partial x_i} \right) + H_{SH} + H_A + \rho A, \qquad (4) $$

where $D/Dt$ denotes the material derivative, $T$ is the temperature. $C_p$ is the specific heat capacity of the material and $k$ the thermal conductivity. $A$ is the radiogenic heat production, the shear heating source term $H_{SH}$ is prescribed in the following way:

$$ H_{SH} = \tau_{ij} (\dot{\epsilon}_{ij} - \dot{\epsilon}_{ij}^{el}), \qquad (5) $$

where $\dot{\epsilon}_{ij}$ is the total deviatoric strain rate and $\dot{\epsilon}_{ij}^{el}$ the elastic strain rate component. The term $H_A$ accounts for adiabatic heating and is given by:

$$ H_A = T\alpha \frac{\partial P}{\partial x_i} v_i, \qquad (6) $$
## 3.2. Rheology

The rheology is defined as visco-elasto-plastic, with the deviatoric strain rate composed of the following components:

$$ \dot{\epsilon}_{ij} = \dot{\epsilon}_{ij}^{\text{el}} + \dot{\epsilon}_{ij}^{\text{dif}} + \dot{\epsilon}_{ij}^{\text{dis}} + \dot{\epsilon}_{ij}^{\text{pl}} = \frac{\hat{\tau}_{ij}}{2G} + A_{\text{dif}} \tau_{ij} + A_{\text{dis}} (\tau_{\text{II}})^{n-1} \tau_{ij} + \dot{\lambda} \frac{\partial Q}{\partial \tau_{ij}}, \quad (7) $$

where $\dot{\epsilon}_{ij}^{\text{el}}$, $\dot{\epsilon}_{ij}^{\text{dif}}$, $\dot{\epsilon}_{ij}^{\text{dis}}$, $\dot{\epsilon}_{ij}^{\text{pl}}$ are the elastic, diffusion, dislocation and plastic strain rate, respectively. G is the elastic shear modulus, $\hat{\tau}_{ij} = \frac{\partial \tau_{ij}}{\partial t} + \tau_{ik}\omega_{kj} - \omega_{ik}\tau_{kj}$ is the Jaumann objective stress rate, $\omega_{ij} = \frac{1}{2}(\frac{\partial v_i}{\partial x_j} - \frac{\partial v_j}{\partial x_i})$ is the spin tensor. $A_{\text{dif}}$, $A_{\text{dis}}$ are diffusion and dislocation constants, respectively, n is the power-law exponent for dislocation creep. $\tau_{\text{II}}$ is the square root of the second invariant of the deviatoric stress tensor (i.e., $\tau_{\text{II}} = (\frac{1}{2}\tau_{ij}\tau_{ij})^{\frac{1}{2}}$). $\dot{\lambda}$ is the magnitude of the plastic strain rate and Q is the plastic flow potential.

The rheological constants of the diffusion and dislocation creep mechanisms, $A_{\text{dif}}$ and $A_{\text{dis}}$, respectively, can be expressed in the following generalized way:

$$ A = B d^{-p} (C_{OH})^r \exp\left(-\frac{E + PV}{RT}\right), \quad (8) $$

Here E stands for the activation energy, R the gas constant, V the activation volume, B the specific prefactor constant, d the grain size, p the exponential factor of grain size, which does only apply for diffusion creep processes, $C_{OH}$ the water content and its exponential factor r. Table 2 summarizes the values of experimentally determined material parameters used to compute both creep constants. Both creep mechanisms are thus temperature- and pressure-dependent.

The magnitude of the plastic strain rate is determined by enforcing the Drucker-Prager yield-criterion:

$$ F = \tau_{\text{II}} - \tau_{\text{Y}}, \quad \tau_{\text{Y}} = \sin(\phi)P + \cos(\phi)c, \quad (9) $$

where $\tau_{\text{Y}}$ denotes the yield stress, $\phi$ is the friction angle and c the cohesion. In incompressible materials, plasticity is dilation-free (e.g., Kaus, 2010), which implies that the plastic flow potential function Q equals $\tau_{\text{II}}$. Plasticity is only active when $F \ge 0$, such that $\dot{\lambda}F = 0$.

## 3.3. Numerical Implementation

The simulations were conducted using the thermo-mechanical code LaMEM (Kaus et al., 2016) where the conservation equations listed above are discretized using a staggered finite difference approach (Harlow & Welch, 1965). The material properties are advected using a Marker and Cell technique (Harlow & Welch, 1965). To solve the nonlinear system of equations LaMEM utilizes the PETSc infrastructure (Balay et al., 2024) employing distributed arrays (DMDA) that enable LaMEM to run on massively parallelized machines.

The effective viscosity of the nonlinear visco-elastic creep rheology is computed in every control volume by numerically solving the following nonlinear scalar equation $r(\eta_{\text{eff}})$ using a bisection algorithm (e.g., Popov & Sobolev, 2008):

$$ r(\eta_{\text{eff}}) = \dot{\epsilon}_{\text{II}}^* - \frac{\tau_{\text{II}}}{2G\Delta t} - A_{\text{dif}}\tau_{\text{II}} - A_{\text{dis}}(\tau_{\text{II}})^n, \quad (10) $$

in which the effective deviatoric stress is evaluated as:

$$ \tau_{\text{II}} = 2\eta_{\text{eff}}\dot{\epsilon}_{\text{II}}^*, \quad (11) $$

and the effective strain rate ($\dot{\epsilon}_{ij}^*$) incorporates a term resulting from advected and rotated deviatoric history stress from the previous time step ($\tau_{ij}^*$) as follows:
$$ \dot{\epsilon}_{ij}^* = \dot{\epsilon}_{ij} + \frac{\tau_{ij}^*}{2G\Delta t}. \qquad (12) $$

History stresses are advected using markers, which is followed by a distance-based averaging of the stress components on the corresponding control volumes of the staggered grid. The stress rotation algorithm is based on the numerical integration of the Jaumann objective stress rate, which gives asymptotically correct results for relatively large finite time steps (Gerya, 2019; Thielmann et al., 2015).

In case that plastic failure occurs, the corresponding effective viscosity is simply evaluated as:

$$ \eta_{\text{eff}} = \frac{\tau_{\text{Y}}}{2\dot{\epsilon}_{\text{II}}^*}. \qquad (13) $$

The discretized system of nonlinear equations is solved with the help of the SNES nonlinear solution infrastructure of PETSc. We utilize a Picard linearization combined with a line search algorithm to facilitate the convergence of nonlinear iterations. The linearized system of algebraic equations is solved using the KSP solver interface of PETSc. We use 5 levels of coupled Galerkin multigrid V-cycle equipped with 5 smoothing sweeps of Richardson iteration with Jacobi preconditioner and damping factor of 0.5. Taking into account the resolution of the simulations (512 × 384 × 128), the processor partitioning (8 × 8 × 4) and the number of multigrid levels (5) the final coarse grid resolution per processor is 4 × 3 × 2. The coarse grid system is solved redundantly on 64 cores using the SuperLU_DIST parallel sparse direct solver (X. S. Li & Demmel, 2003). The nonlinear solver tolerance was set to a relative value of 10⁻². We compared this with simulations that have a relative tolerance of 10⁻⁴, which yielded similar results (see Table 4, Simulation Tol).

## 4. Initial Model Setup

The following is a detailed description of the initial setup for the reference model (REF). The setup was developed by finding a reference model that best fits the first-order mechanisms in the Mediterranean since 35 Ma. It is based on various tectonic reconstructions of the area (e.g., Faccenna et al., 2014; Le Breton et al., 2021; Rosenbaum et al., 2002a) and was continuously updated based on the knowledge gained from running over 500 simulations.

The model setup consists of a rectangular area, with dimensions of 6,000 km in the east-west direction and 3,500 km in the south-north direction. It was shown that for models of this size, the Earth's curvature has a minor effect on the dynamics of the lithosphere and a cartesian model box is therefore sufficient (Macherel et al., 2024). The vertical extension reaches down to a depth of 660 km (Figures 2a and 2b). On top of the crustal units, there is a 50 km thick layer of sticky air to account for a free surface (Crameri et al., 2012), so that the total depth of the model is 710 km. The material parameters of this layer are given in Table 1. The initial crustal configuration of the model features a diverse composition including a 10 km thick continental margin (thinned continental crust), 10 km thick oceanic crust and 35 km thick continental crust and 40 km thick orogenic crust (Figure 2a). Material parameters for these units are given in Table 1.

The initial setup displays several weak zones (WZ) in its crust and lithosphere (Figure 2a). As shown in Table 3 these different weak zones are characterized by different phases and depths. Initially the asthenosphere beneath the crust is homogeneous except for heterogeneities beneath the continental margin and the initial slabs that extend into the mantle. Tilted ellipses are inserted within the mantle lithosphere of the continental margin, creating a weak horizontal layer (Appendix A). These inclusions consist of a weaker mantle phase (Table 1 WZM). The Mediterranean region is characterized by the presence of three distinct subduction systems. The slab surface is composed of 15 km of oceanic crust and is overlaid by a 25 km thick weak zone that imitates the subduction channel. This facilitates the rollback of the slab, enabling it to rapidly reach a realistic geometry. In the western Mediterranean, a slab extends from the northernmost tip of the Adriatic microplate toward the Iberian Peninsula (Figures 2a and 2b). This western Mediterranean slab (WMS), which initially dips to the north, reaches a depth of 150 km. To the north of the trench, continental blocks can be found that are separated by weak material. In the eastern Mediterranean, the northeast-dipping Dinaric-Hellenic slab (DHS) is prescribed (Figures 2a and 2b). The Dinaric slab in the north is 300 km deep, while the Hellenic slab in the south initially reaches a depth of 600 km. The length of the slabs at 35 Ma was determined using the present-day length of subducted lithosphere observed in tomographic images (e.g., Piromallo & Morelli, 2003; Rappisi et al., 2022) and subtracting the amount of plate convergence since 35 Ma
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Initial model setup. (a) Mapview of the setup at 35 Ma and the different phases of the crust. (b) 3D perspective of the model setup. (c) Inflow boundary velocity with the example of the Arabian plate. (d) Vertical temperature profile for continental and oceanic lithosphere.</figcaption></figure>

<table><caption>**Table 1**

Material Parameters of All Model Units of the Reference Simulation</caption><thead><tr><th></th><th>SA</th><th>AS</th><th>WM</th><th>WZ</th><th>WZC</th><th>CC</th><th>OcC</th><th>CM</th><th>OrC</th></tr></thead><tbody><tr><td>ρ [kg/m<sup>3</sup>]</td><td>1,000</td><td>3,300</td><td>3,300</td><td>3,300</td><td>2,800</td><td>2,800</td><td>2,850</td><td>3,300</td><td>3,300</td></tr><tr><td>Dislocation η<sub>eff</sub> = 10<sup>19</sup>Pas</td><td>Dry Olivine</td><td>Wet Olivine</td><td>Wet Olivine</td><td>η<sub>eff</sub> = 10<sup>19</sup>Pas</td><td>Dry Olivine</td><td>Quartzite</td><td>Plagioclase</td><td>Wet Quarzite</td><td>Wet Quarzite</td></tr><tr><td>Diffusion</td><td>—</td><td>Dry Olivine</td><td>Wet Olivine</td><td>Wet Olivine</td><td>Dry Olivine</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>c [MPa]</td><td>—</td><td>20</td><td>20</td><td>—</td><td>10 – 20</td><td>20</td><td>10</td><td>20</td><td>20</td></tr><tr><td>c sof. [MPa]</td><td>—</td><td>2</td><td>2</td><td>—</td><td>1 – 2</td><td>2</td><td>1</td><td>2</td><td>2</td></tr><tr><td>φ [°]</td><td>—</td><td>30</td><td>30</td><td>—</td><td>4.5 – 12.5</td><td>30</td><td>1</td><td>10</td><td>30</td></tr><tr><td>φ sof. [°]</td><td>—</td><td>3</td><td>3</td><td>—</td><td>0.45 – 1.25</td><td>3</td><td>0.1</td><td>1</td><td>3</td></tr><tr><td>A [Wm<sup>-3</sup>]</td><td>—</td><td>2 × 10<sup>-8</sup></td><td>2 × 10<sup>-8</sup></td><td>2 × 10<sup>-8</sup></td><td>2 × 10<sup>-8</sup></td><td>2 × 10<sup>-6</sup></td><td>2 × 10<sup>-7</sup></td><td>2 × 10<sup>-6</sup></td><td>2 × 10<sup>-6</sup></td></tr></tbody></table>

Note. SA is sticky air layer, As is asthenosphere, WM is weak mantle, WZ is weak zone, WZC is weak zone continental margin, CC is continental crust, OcC is oceanic crust, CM is continental margin and OrC is orogenic crust. Dislocation and diffusion creep laws are described in Table 2. * A constant shear modulus G = 5 × 10<sup>10</sup>Pa is used. * A constant conductivity k = 3Wm<sup>-1</sup>K<sup>-1</sup> (SA: k = 100Wm<sup>-1</sup>K<sup>-1</sup>) and heat capacity C<sub>p</sub> = 10<sup>3</sup>Jkg<sup>-1</sup>K<sup>-1</sup> (SA: C<sub>p</sub> = 10<sup>6</sup>Jkg<sup>-1</sup>K<sup>-1</sup>) is used. * A constant thermal expansion α = 3 × 10<sup>-5</sup>K<sup>-1</sup> is used (not used for SA).
<table><caption>**Table 2**

Creep Law Parameters for the Creep Laws Listed in Table 1</caption><thead><tr><th>Dislocation creep</th><th>B [Pa<sup>-n</sup>s<sup>-1</sup>]</th><th>E [Jmol<sup>-1</sup>]</th><th>V [m<sup>3</sup>mol<sup>-1</sup>]</th><th>n []</th><th>p []</th><th>r []</th></tr></thead><tbody><tr><td>Dry Olivine<sup>a</sup></td><td>1.1 &times; 10<sup>5</sup></td><td>530 &times; 10<sup>3</sup></td><td>11.0 &times; 10<sup>-6</sup></td><td>3.5</td><td>0.0</td><td>0.0</td></tr><tr><td>Wet Olivine<sup>a,d</sup></td><td>1.6 &times; 10<sup>3</sup></td><td>520 &times; 10<sup>3</sup></td><td>10.0 &times; 10<sup>-6</sup></td><td>3.5</td><td>0.0</td><td>1.2</td></tr><tr><td>Plagioclase (An<sub>75</sub>)<sup>b</sup></td><td>3.3 &times; 10<sup>-4</sup></td><td>238 &times; 10<sup>3</sup></td><td>0.0 &times; 10<sup>-6</sup></td><td>3.2</td><td>0.0</td><td>0.0</td></tr><tr><td>Quartzite<sup>b</sup></td><td>6.7 &times; 10<sup>-6</sup></td><td>156 &times; 10<sup>3</sup></td><td>0.0 &times; 10<sup>-6</sup></td><td>2.4</td><td>0.0</td><td>0.0</td></tr><tr><td>Wet Quartzite<sup>b</sup></td><td>3.2 &times; 10<sup>-4</sup></td><td>154 &times; 10<sup>3</sup></td><td>0.0 &times; 10<sup>-6</sup></td><td>2.3</td><td>0.0</td><td>0.0</td></tr><tr><td>Diffusion creep</td><td>B [Pa<sup>-n</sup>s<sup>-1</sup>]</td><td>E [Jmol<sup>-1</sup>]</td><td>V [m<sup>3</sup>mol<sup>-1</sup>]</td><td>n []</td><td>p []</td><td>r []</td></tr><tr><td>Dry Olivine<sup>a,c</sup></td><td>1.5 &times; 10<sup>9</sup></td><td>375 &times; 10<sup>3</sup></td><td>10.0 &times; 10<sup>-6</sup></td><td>0.0</td><td>3.0</td><td>0.0</td></tr><tr><td>Wet Olivine<sup>a,c,d</sup></td><td>2.5 &times; 10<sup>7</sup></td><td>375 &times; 10<sup>3</sup></td><td>8.0 &times; 10<sup>-6</sup></td><td>0.0</td><td>3.0</td><td>0.8</td></tr></tbody></table>

<sup>a</sup>Hirth and Kohlstedt (2003). <sup>b</sup>Ranalli (1995). <sup>c</sup>A constant grain size of d = 10<sup>-3</sup>m is used. <sup>d</sup>A water content of C<sub>OH</sub> = 1000H/10<sup>6</sup>Si is used.

(e.g., Faccenna et al., 2001; Le Breton et al., 2021). Present-day tomographic images (e.g., Piromallo & Mor-elli, 2003; Rappisi et al., 2022; Šumanovac et al., 2017) show a slab gap beneath the Dinarides toward the north, which suggests that a slab-breakoff must have occurred. Based on the ages of calc-alkaline, subduction-related, magmatic rocks in the region, the slab breakoff is suggested to have occurred between approximately 37 and 22 Ma (e.g., Handy et al., 2019; Schefer et al., 2011). Therefore the initial setup includes a horizontal weak zone in the mantle at a depth of 200 km in the northern part of the DHS, which initiates a slab break-off beneath the Dinarides. The Alpine-Carpathians slab (AICS) in the north dips steeply southwards, with a curved trench geometry and a depth of 250 km (Figures 2a and 2b). The subduction system in the Alps has been simplified with the introduction of a continuous trench from the Alps in the west to the Carpathians in the east. However, the geological history of this region may have been more complex (e.g., Handy et al., 2010).

The African plate is separated from the Arabian plate by a vertical weak zone (Figure 2 a WZ7). During the simulation time of 35 Myr, both the African and Arabian plates are pushed northwards. The Arabian plate is pushed at a constant velocity of 2.2 cm/yr, the African plate is pushed at a velocity of 0.8 cm/yr in its westernmost part and 1.6 cm/yr in its easternmost part, with a linearly increasing plate velocity in between. For this purpose, the y-velocities at the southern boundary are prescribed as an inflow boundary condition. The inflow velocity is constant in the vertical direction down to a depth of 100 km, and decreases to zero between a depth of 100 and 180 km (Figure 2c). In order to compensate for the inflow of material, an outflow is prescribed at all boundaries between 300 and 660 km depth. The magnitude of the outflow is uniform throughout the domain and precisely compensates for the inflow prescribed at the southern boundary. The top boundary above the sticky air layer is stress-free and all other boundaries are shear stress free and have zero normal velocity components. This type was chosen to mimic the convergence between Africa and Eurasia in a Eurasia-fixed reference frame. Alternative approaches that would prescribe the velocities of the lithospheric plates and incorporate stress-free boundary conditions in the mantle are also conceivable. This could lead to a different timing of subduction dynamics close to the boundaries. However, our research focus is on the Adriatic microplate and once the subduction systems reach the vicinity of Adria they are unlikely to be strongly influenced by the lateral boundary conditions.

The initial temperature structure of the model is calculated based on the half-space cooling model (e.g., Tur-cotte & Schubert, 2002):

<table><caption>**Table 3**

The Parameters for All Weak Zones That Are Part of the Crustal Configuration (Figure 2a)</caption><thead><tr><th></th><th>WZ 1</th><th>WZ 2</th><th>WZ 3</th><th>WZ 4</th><th>WZ 5</th><th>WZ 6</th><th>WZ 7</th><th>WZ 8</th><th>WZ 9</th><th>WZ 10</th><th>WZ 11</th></tr></thead><tbody><tr><th>Phase</th><td>WZ</td><td>WZ</td><td>WZ</td><td>WZC</td><td>WZC</td><td>WZC</td><td>WZ</td><td>WZ</td><td>WZC</td><td>WZC</td><td>WZ</td></tr><tr><th>Depth [km]</th><td>100</td><td>100</td><td>60</td><td>40</td><td>65</td><td>40</td><td>240</td><td>40</td><td>65</td><td>65</td><td>80</td></tr><tr><th>c [MPa]</th><td>—</td><td>—</td><td>—</td><td>20</td><td>10</td><td>20</td><td>—</td><td>—</td><td>20</td><td>20</td><td>—</td></tr><tr><th>φ[°]</th><td>—</td><td>—</td><td>—</td><td>4.5</td><td>5.0</td><td>4.5</td><td>—</td><td>—</td><td>12.5</td><td>12.5</td><td>—</td></tr></tbody></table>
<table><caption>**Table 4**
Different Simulations Are Presented, Along With a List of the Parameters That Were Modified in Each of the Simulations</caption><thead><tr><th>Simulation</th><th>Modified parameters</th></tr></thead><tbody><tr><td>REF</td><td></td></tr><tr><td>Tol</td><td>Relative solver tolerance 10<sup>-4</sup></td></tr><tr><td>Sim1</td><td>WZ 5 fr = 4.5°</td></tr><tr><td>Sim2</td><td>WZ 5 fr = 4.75°</td></tr><tr><td>Sim3</td><td>WZ 5 fr = 4.0°</td></tr><tr><td>Sim4</td><td>WZ 5 fr = 4.85°</td></tr><tr><td>Sim5</td><td>WZ 5 fr = 5.15°</td></tr><tr><td>Sim6</td><td>WZ 5 fr = 5.25°</td></tr><tr><td>Sim7</td><td>No push African/Arabian plate</td></tr><tr><td>Sim8</td><td>African/Arabian plate inflow velocities halved</td></tr><tr><td>Sim9</td><td>African/Arabian plate inflow velocities quartered</td></tr><tr><td>Sim10</td><td>WZ 9 fr = 11.5°</td></tr><tr><td>Sim11</td><td>WZ 9 fr = 13.5°</td></tr><tr><td>Sim12</td><td>WZ 4 fr = 3.5°</td></tr><tr><td>Sim13</td><td>WZ 4 fr = 5.5°</td></tr><tr><td>Sim14</td><td>Oceanic crust t<sub>age</sub> = 80 Myr</td></tr><tr><td>Sim15</td><td>Oceanic crust t<sub>age</sub> = 40 Myr</td></tr><tr><td>Sim16</td><td>Weak zone southwest of Adria</td></tr><tr><td>Sim17</td><td>Weak zone southeast of Adria</td></tr><tr><td>Sim18</td><td>Asthenosphere: V<sub>dis</sub> = 11.5 × 10<sup>-6</sup> m<sup>3</sup>/mol, V<sub>dif</sub> = 10.5 × 10<sup>-6</sup> m<sup>3</sup>/mol</td></tr><tr><td>Sim19</td><td>Asthenosphere: V<sub>dis</sub> = 10.5 × 10<sup>-6</sup> m<sup>3</sup>/mol, V<sub>dif</sub> = 9.5 × 10<sup>-6</sup> m<sup>3</sup>/mol</td></tr><tr><td>Sim20</td><td>Asthenosphere: V<sub>dis</sub> = 8 × 10<sup>-6</sup> m<sup>3</sup>/mol, V<sub>dif</sub> = 7 × 10<sup>-6</sup> m<sup>3</sup>/mol</td></tr><tr><td>Sim21</td><td>Slabs t<sub>age</sub> = 60 Myr</td></tr><tr><td>Sim22</td><td>Slabs t<sub>age</sub> = 100 Myr</td></tr><tr><td>Sim23</td><td>Slabs t<sub>age</sub> = 120 Myr</td></tr><tr><td>Sim24</td><td>Oceanic crust southwest of Adria narrower</td></tr><tr><td>Sim25</td><td>Oceanic crust southwest of Adria wider</td></tr><tr><td>Sim26</td><td>No initial continental crust southwest of Hellenic trench</td></tr><tr><td>Sim27</td><td>No weak mantle phase</td></tr><tr><td>Sim28</td><td>No Alpine-Carpathian slab</td></tr><tr><td>Sim29</td><td>No initial slab tearing Dinaric slab</td></tr><tr><td>Sim30</td><td>Adria's weak zone depth 200 km</td></tr><tr><td>Sim31</td><td>Adria's weak zone depth 160 km</td></tr><tr><td>Sim32</td><td>Adria's weak zone width 50 km</td></tr><tr><td>Sim33</td><td>No Alpine-Carpathian slab and no weak zone north of Adria</td></tr></tbody></table>

$$T = T_{surf} + (T_{Mantle} - T_{surf}) \text{erfc}\left(\frac{d}{2\sqrt{\kappa t_{age}}}\right), \quad (14)$$

where $T_{surf}$ is the surface temperature (0°C) at the surface of the slab which reaches into the asthenosphere. $T_{Mantle}$ = 1350°C is the temperature at the lithosphere-asthenosphere boundary, $d$ is the depth, $\kappa = 10^{-6} m^2 s^{-1}$ is the thermal diffusivity. For oceanic crust $t_{age}$ is 60 Myr and for continental crust and continental margin $t_{age}$ is
120 Myr. A gradient of 0.3 K km⁻¹ is added to account for the adiabatic gradient in the asthenosphere. To ensure a more realistic initial temperature setup and account for radiogenic heating, the energy equation is solved for 0.2 Myr and 20 timesteps before the simulation begins. The resulting temperature structure is displayed in Figure 2b. Figure 2d shows the geothermal gradients for oceanic and continental lithosphere, respectively. A temperature and effective viscosity cross-section is displayed in Appendix B. Phase changes, such as the transition from olivine to wadsleyite, are not considered, as we cannot fully account for them in an incompressible formulation of the Stokes equations. However, this could have an impact on the timing and dynamics of subduction zones.

During the simulation, the effective viscosity is limited to be between 10¹⁹ Pas and 10²⁴ Pas. The bottom temperature is fixed at 1,600°C. The model domain is discretized using a grid of 512 × 384 × 128 grid cells in x, y, z-direction. This results in a spatial resolution of about 12 km × 9 km in the x and y directions, respectively, and a vertical resolution ranging from 6 km in the lower mantle to 4.3 km in the lithosphere. The grid spacing in the z-direction is variable. There are 94 cells between the bottom of the model and –100 km, 24 cells between –100 and 5 km, and 10 cells between 5 and 50 km. Plastic softening is applied for all phases, which means that between 10% and 60% of accumulated plastic strain, the cohesion c and the friction angle φ are reduced linearly by 90%.

# 5. Results

Over 500 3D simulations were conducted to determine a reference model (REF) that best fits the motion of the various subduction-collision fronts, the continental blocks in the Western Mediterranean (Corsica-Sardinia, Al-Ka-Pe-Ca (see Figure 1)) and the motion and rotation of Adria compared to geological reconstructions of the region (Section 2). Indeed, the REF model reproduces well the trench retreat of the different subduction zones (e.g., ca. 1,000 km retreat of the Calabrian and Hellenic subduction zones). We compared the total amount of counterclockwise rotation and northward motion of the Adriatic microplate. The REF model reproduces well the rotation of Adria, while the overall northward motion of Adria is lower in the model (ca. 170 km) compared to the tectonic reconstructions (ca. 280 km). The relation between northward motion of Adria and the subduction dynamics in the Alps is discussed in Section 6.3. More information on the quantitative comparison between the REF model and tectonic reconstructions is available Text S1 and Figure S1 in Supporting Information S1.

The following section provides a detailed description of REF, with particular attention paid to the dynamics around the Adriatic microplate. To visualize the subducting slabs, a 80 km thick layer beneath the slab surface and the crust composed of mantle material is used to track the motion. The results will be displayed with reference to geological time (Ma), with 35 Ma marking the starting point of the simulation. We provide animations of the key simulations in Supporting Information S1.

## 5.1. Reference Simulation

The subduction dynamics of REF is dominated by rollback subduction. Generally, the WMS retreats to the south, the DHS retreats to the southwest and the AICS retreats to the north (Figure 3). The slab retreat is influenced by various parameters, such as the extent and shape of the continental and oceanic crust or the temperature structure of the lithosphere. Below we describe each subduction system in more detail.

At 35 Ma, the WMS extends from the northeast to the southwest and initially reaches a depth of 150 km (Figure 3a). The curved trench retreats at an angle perpendicular to its trench-line, with the central part of the slab retreating at a faster rate than its edges. This results in a curved trench geometry (Figures 3b and 4a). The motion of the trench can be observed by tracking the motion of the continental blocks to the northwest of the trench. The central blocks move in a southerly direction toward Africa, while the western blocks rotate in a clockwise direction toward the west and eventually dock at the Iberian Peninsula at 20 Ma (Figure 3c). The easternmost block, which represents the Corsica-Sardinia continental block, rotates anticlockwise toward the east and subsequently settles into a north-south orientation (Figures 3d–3f). A rift system opens during the retreat of the WMS slab. A north-south oriented rift zone in the center of the Mediterranean is active between 28 Ma and 20 Ma while the WMS is propagating toward Africa (Figures 3b and 3c). Finally around 20 Ma, the subduction zone reaches the African continent, resulting in a slab break-off at approximately 16 Ma (Figure 4b). This break-off propagates to the east and west, subdividing the subduction zone into the Gibraltar slab (GS) to the west and the Apennine-Calabrian slab (ApCS) to the east (Figure 4c). This results in a relatively narrow GS that extends down to
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Evolution of the reference model (REF) over 35 Myr. (a–f) Six different timesteps of a mapview of the 3D simulation, with the continental crust (gray) and oceanic crust (blue) depicted slightly translucent. (a) Shows the abbreviations for the three different slabs. The red color highlights the parts of the crust where significant strain rate (second invariant of the strain rate tensor) occurs. The colors ranging from violet to green indicate the depth of the subducted mantle lithosphere between –120 km and –660 km. See Movie S1 for an animation of the simulation.</figcaption></figure>

660 km depth and retreats to the west (Figures 4h and 4i). As the GS retreats, the gap between the Iberian Peninsula and Africa almost closes. Concurrently, the slab break-off progresses at a faster rate to the east, resulting in the ApCS retreating in that direction (Figure 4c). At 20 Ma, the trench of the ApCS is oriented north-south (Figure 3c). Subsequently, the Ionian oceanic crust begins to be subducted in the south, while the rift zone opens between the Corsica-Sardinia continental block and Adria (Figures 3d–3f), and slab retreat accelerates in this region. The northern part of the ApCS then slows down as it partially consumes the Adriatic continental microplate, while the southern part of the ApCS continues to subduct Ionian oceanic crust resulting in a curved trench geometry of the ApCS (Figures 3e and 3f).

Between 35 Ma and 25 Ma the AlCS retreats northward at a constant velocity across the trench maintaining the initial trench geometry (Figures 3a and 3b). Since 25 Ma the retreat of the eastern part of the slab has accelerated toward the northeast, consuming the oceanic crust in that region (Figure 3c). This leads to a bending of the trench line to the east (Figure 3c). At approximately 21 Ma a slab break-off occurs in the central part of the slab
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Different viewing angles of REF. (a–c) Present the model box from the northwest, (d–f) show a close-up view from the south, and (g–i) display the 3D box from the northeast. The figure shows different parts of the crust, including continental crust (gray), thinned continental crust (orange-brown) and oceanic crust (blue). Subducted lithosphere is shown between a depth of 115 and 660 km with the color scheme identical to that of Figure 3.</figcaption></figure>

(Figure 4d), initiating a slab tearing both to the west and to the east (Figure 4e). Initially the slab tearing propagates at approximately the same velocity in both directions, but then it accelerates to the east around 12 Ma (Figure 4f). It quickly acts along the curved trench resulting in a complete slab break-off in the east at 10 Ma (Figure 4f). In contrast, to the west, slab tearing largely ceases around 10 Ma, resulting in a subduction system that remains stagnant in the northwest of Adria (Figures 4d–4f). As illustrated in Figures 3d–3f, the position of the slab in the west has remained relatively stable over the past 15 Myr.

Simultaneously, the DHS retreats to the southwest (Figures 3a–3c). At 35 Ma, a weak zone vertically separates the DS in the mantle, causing the lower part of the slab to sink rapidly to a depth of 660 km (Figure 4g). North of the Arabian plate, the DHS breaks off. The slab tearing progresses toward the west until approximately 15 Ma (Figures 3b–3d). As the DHS retreats to the southwest, it initially only consumes continental margin (Figure 4g). At 25 Ma, the southern part of the DHS reaches oceanic crust (Figure 3b) and the slab retreat accelerates to the south, resulting in a curved trench (Figure 4h). Additionally, subduction of continental margin to the north continues at a slower rate during the fast slab retreat in the south. Between 21 Ma and 0 Ma the DHS consumes almost the entire oceanic crust southeast of Adria (Figure 4i). At present-day, the trench almost reaches the African continental margin and the southern tip of Adria (Figures 3f and 4i). Although the northern part of the slab reaches Adria along its eastern boundary, the retreat was much slower than in the south. Consequently, the northern slab only extends to a depth of approximately 300–400 km (Figure 4i). South of Adria, the ApCS and the DHS are in close proximity, with less than 100 km of oceanic crust remaining unsubducted (Figure 3f).
### 5.1.1. Mantle Flow

This section provides a detailed description of asthenospheric mantle flow in REF. In general, two types of flow patterns can be observed: poloidal flow, which circulates in a vertical plane, and toroidal flow around the slab edges, which circulates horizontally.

We begin by examining the poloidal flow patterns in the asthenosphere. The direction and magnitude of asthenospheric flow can be observed in Figure 5. The DHS retreats to the southwest, inducing trench perpendicular flow between 35 Ma and 25 Ma. This is particularly pronounced in the southern half of the slab (Figures 5a and 5b). Figure 5g shows a vertical cross section through the model. It can be observed that the retreat of the DHS introduces poloidal flow cells to the southwest and northeast of the slab. The poloidal flow cell to the northeast is located in the lower half of the asthenosphere, above the horizontal slab piece that lies at the bottom of the model. As the slab remains fully connected, the slab pull forces are partly transmitted through the subducting slab to the lower part of the subduction. The greatest velocities are observed in the back-arc of the system, where asthenospheric material ascends. This results in the formation of a poloidal flow cell between the back-arc extensional area and the flat slab at the bottom of the upper mantle. The second poloidal flow cell is situated to the southwest of the slab (Figure 5g). It is located below the lithosphere, which is undergoing a northward motion. This results in the formation of a clockwise rotational flow cell.

The magnitude of mantle flow associated with poloidal flow in REF is typically smaller than horizontal mantle flow associated with toroidal flow (Figure 5). The following section will therefore concentrate on the horizontal mantle flow. The retreat of the different subduction zones introduces horizontal pressure gradients in the asthenosphere, which are the main drivers of mantle flow in the upper mantle. To this end, we introduce a new variable, the dynamic pressure $P_{dyn}$, which in this study is defined as the pressure $P$ minus the average lithostatic pressure per depth level (Figure 6).

Between 25 Ma and 21 Ma the DHS slab breaks off north of the Arabian plate, resulting in a slab tearing from east to west until 14 Ma (Figures 5a–5d). Through the rapid slab tearing, a strong toroidal flow around the southern slab edge evolves which is strongest around 18 Ma (Figure 5c). After 14 Ma a curved trench develops in the south and the magnitude of the toroidal flow gradually weakens until present-day (Figures 5d–5f). Nevertheless, the slab is still retreating, as evidenced by the remaining trench perpendicular flow (Figure 5f). The northern part of the DHS generally exhibits a slower rate of retreat than the southern part, which results in a slab bending between north and south. The WMS first induces a trench perpendicular flow in the back of the subduction zone, caused by the retreat of the slab (Figures 5a–5c). On the western edge of the subduction zone, this is accompanied by a clockwise rotating toroidal flow. This continues until the WMS breaks off at Africa which opens a slab gap (Figures 5d and 6c). The retreat of the slab to the south has resulted in a significant $P_{dyn}$ gradient between the back of the slab and the front (Figures 5d, 5e and 6c, 6d). After the slab gab has opened, a toroidal flow cell develops around the southern edge of the ApCS, resulting in the compensation of horizontal $P_{dyn}$ differences (Figures 6c and 6d). In the western part of the slab gap, in the vicinity of the GS, the asthenospheric flow is nearly perpendicular to the trench (Figure 5e). Following the slab break-off, the toroidal flow accelerates the retreat of the ApCS to the east, ultimately resulting in a curved trench geometry of the ApCS (Figures 5d–5f and 6d–6f).

The ApCS and the DHS retreat toward each other, creating a region of elevated $P_{dyn}$ between the two slabs (Figure 6). From 35 Ma to 25 Ma, the mantle flow between these two slabs is directed northwards. However as the toroidal flow cells south of the ApCS and the DHS are advancing, they cause a southward directed mantle flow underneath the southern part of Adria (Figure 6b), which leads to a switch in asthenospheric flow. After 20 Ma, a transition is observed in the asthenospheric flow beneath Adria. The flow direction undergoes a change from northward in the northern region to southward beneath the southern half of Adria (Figures 6b–6f). This pattern can still be observed at present-day (Figure 6f).

Between 35 Ma and 25 Ma, the AICS retreats to the north inducing a trench perpendicular mantle flow (Figure 5a). After 25 Ma, the slab retreat accelerates to the east, which increases the magnitude of trench perpendicular mantle flow in this region (Figure 5b). This is followed by the slab break-off in the center of the AICS. Consequently, an anticlockwise rotating toroidal flow cell develops in the eastern region of the slab gap (Figures 5d, 5e and 6b–6d). The toroidal flow rapidly accelerates, resulting in the complete break-off of the slab in the east at 9 Ma (Figures 6d and 6e). For approximately 3 Myr, a combination of strong downward-oriented vertical flow and toroidal flow can be observed (Figure 5e). On the western edge of the slab gap no toroidal
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Mantle flow of REF. (a–f) The gray lines outline the edges of the continents, whereas the purple color represents the subducting lithosphere between a depth of 120 and 400 km. Mantle flow directions are displayed by flow lines, colored by flow magnitude. The red arrows indicate the direction of the toroidal flow cells, while the pink line in panel (b) represents the location of the vertical cross-section depicted in panel (g). (g) SW-NE vertical cross-section through the model box. The red arrows indicate the flow direction of the poloidal flow cells, while the black line displays the 1,300°C isotherm.</figcaption></figure>

flow develops. South of the slab, to the west, the $P_{dyn}$ remains relatively high throughout the simulation, due to asthenospheric escape flow from the south. Conversely, to the east, a region with low $P_{dyn}$ develops as the trench retreats at an accelerated rate after 25 Ma (Figure 6). Subsequently, a trench parallel flow establishes south of the AlCS around 15 Ma, which lowers the $P_{dyn}$ difference (Figures 6b–6d). This mantle flow distinctly weakens after approximately 10 Ma when the $P_{dyn}$ is almost entirely equilibrated (Figures 6e and 6f).
<figure><img src="image_8.png" /><figcaption>**Figure 6.** (a–f) Dynamic pressure variations at a depth of 300 km at different timesteps. Arrows indicate the direction of horizontal mantle flow, with colors representing the magnitude of the horizontal flow.</figcaption></figure>

5.1.2. Plate Motion of Adria

The complex subduction dynamics discussed above also leads to an intricate pattern of plate motion. This section
gives a detailed description of the temporal evolution of Adria's plate motion.

Figure 7 displays the plate motion of the Adriatic microplate and the surrounding area like the African continent in the south or the continental blocks within the western Mediterranean. Between 35 Ma and 22 Ma Adria moves in conformity with Africa to the north (Figure 7a). After 22 Ma the motion of Adria and Africa start diverging partially. The southern part of Adria still moves exactly in the same direction as Africa, whereas the northern part of Adria undergoes a change in direction to the northwest (Figure 7b). As illustrated in Figures 7a and 7b, after 22 Ma the northernmost extent of Adria moves approximately in the same northwest-direction as the retreating trench of the AlCS, situated to the north. Overall this leads to an anticlockwise rotation of the Adriatic microplate (Figure 8a).

The southern part of Adria moves in the same direction as Africa until approximately 13 Ma (Figure 8c). Sub-
sequently, the southern part of Adria begins to move differently from Africa, by shifting to the northeast, as the
entire Adriatic plate rotates (Figures 7c–7e). Between 13 Ma and present-day, the change in direction of southern
Adria's plate motion is quite constant (Figures 7c–7e and 8c). Furthermore, both the southern and the northern
parts of Adria rotate to the northeast (Figure 8a). At 13 Ma, the angles between “Adria northeast” (green dot in
Figure 8a) and “Adria south” (blue dot in Figure 8a) were approximately 40° offset from each other. However, the
direction of plate motion of these two parts converged until the present day, when both are moving in approx-
imately the same direction (8c). From 5 Ma, Adria rotates around a rotation pole located to the northwest (purple
dot in Figure 8a), as this part of Adria almost stops moving around this time (8b).

In the initial stages of the simulation, the entire Adriatic microplate moves with the same magnitude as the African
plate (Figure 7a). At around 22 Ma, the northern part of Adria decelerate in comparison to the southern region
<figure><img src="image_9.png" /><figcaption>**Figure 7.** (a–f) Plate motion at a depth of 10 km at different timesteps. Arrows indicate the direction of plate motion, while the magnitude is represented by the color of the arrow. Red lines indicate the trench position of the different subduction systems, while the triangles mark the dipping direction of the slab. (f) Presents the results of the simulation after a simulation time of 38 Myr.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** (a) Visualizes the plate motion of four selected tracers in the crust. Initially, these tracers are located at a depth of 10 km. Different shades of gray depict the initial position of Adria and Adria's position at present-day. (b) Plots the magnitude of plate motion of the four tracers over time. (c) Plots the direction of the tracers over time. The direction is indicated by the angle, where 0° represents a northward direction, -90° indicates a westward motion, and 90° represents an eastward motion.</figcaption></figure>
(Figure 8b). From that point onwards, the velocity constantly decreases until 5 Ma (Figure 8b). The decreasing velocity in the northern part of Adria is in alignment with the observed slowing of the Alpine subduction (Figures 7a–7c and 8b). The southern half of Adria starts to slow down around 18 Ma (Figure 8b). This occurs at a slower rate than in the northern part (Figure 8b). From 18 Ma onwards, the velocity of the southern part of the Adriatic plate motion decreases at a more rapid rate than that of the African plate (Figure 7b).

The REF simulation was run for an additional 5 Myr, and Figure 7f illustrates that the direction of Adria changes abruptly from northeast to southeast at approximately 38 Myr since the beginning of the simulation. This abrupt change happens when the Calabrian subduction zone and the Hellenic subduction zone converge south of Adria, cutting off the connection between the Adriatic microplate and the African continent.

## 5.2. Parameter Study

Over the course of this study, many simulations were conducted with varying configurations in order to identify a model that adequately represents the first-order processes observed in the Mediterranean. While we cannot present all models, we here discuss a representative selection in which specific parameters are changed with respect to REF, to obtain an understanding on how variations in the model parameters affect the outcome. These simulations are summarized in Table 4.

### 5.2.1. Influence of Weak Zones

The weak zones of the model are very sensitive to changes in material properties. Small changes in friction angle quickly lead to changes in subduction dynamics, as observed in Sim12 and Sim13 (Table 4). In Sim12 the friction angle of weak zone 4, which encompasses the Carpathians region, is decreased from 4.5° (REF) to 3.5°. This leads to an accelerated slab tearing. However if the friction angle is increased to 5.5° (Sim13), the slab does not break off entirely, so that part of the slab is still attached after 35 Myr of simulation time. Similar sensitivity tests were conducted for weak zone 9, in the region of the Gibraltar slab, and weak zone 5, at the northern boundary of the African plate. In Sim10 the friction angle of the weak zone 9 is decreased to 11.5° leading to a faster slab retreat of the Gibraltar slab, while an increase of friction angle to 13.5° reduces the retreat velocity in the region and the slab starts to break off in the mantle. The strength of the weak zone 5 is of critical importance, as it determines the retreat velocities of both the Apennine-Calabrian and Dinaric-Hellenic subduction zones. In Sim1, in which the friction angle of weak zone 5 is decreased to 4.5°, the aforementioned subduction zones already converge after 34 Myr, whereas in REF this only happens after 37 Myr. Sim5 has a stronger weak zone 5 ($\phi$ = 5.25°), so that the subduction zones have not yet converged after 40 Myr. Changing the friction angle of the weak zone also affects the interaction of different subduction zones. If the friction angle of weak zone 5 is larger, the propagation of the slab break-off at Africa is slower. This results in a longer period of lateral connection between the GS and the ApCS compared to the REF. In the case of the Gibraltar slab, this results in a generally lower retreat velocity, which leads to a complete slab break-off. This is not the case if the friction angle of weak zone 5 is distinctly lower (Sim1).

### 5.2.2. Mantle Rheology

The simulation results are also highly sensitive to the mantle rheology. A combination of diffusion and dislocation is used and by changing parameters such as the activation volume within the error bars given by Hirth and Kohlstedt (2003), the mantle dynamics changes drastically (see Sim18, Sim19, Sim20). By changing the activation volume from $V_{dis} = 11 \times 10^{-6} \text{ m}^3\text{mol}^{-1}$ to $V_{dis} = 8 \times 10^{-6} \text{ m}^3\text{mol}^{-1}$ and $V_{dif} = 10 \times 10^{-6} \text{ m}^3\text{mol}^{-1}$ to $V_{dif} = 7 \times 10^{-6} \text{ m}^3\text{mol}^{-1}$ the mantle rheology becomes remarkably weaker and thus the slab retreat velocities are significantly higher (Sim20).

### 5.2.3. Subduction of Adriatic Lithosphere North of Adria

A subduction of Adriatic lithosphere is not observed in REF. In Sim28, there is no initial AICS. In this case, the initiation of a subduction process, whereby Adriatic lithosphere is subducted beneath Europe, can be observed to the north of Adria after 2 Myr. A similar behavior can be observed if the weak zone north of Adria (WZ 11) is deeper than in REF. If the depth of the weak zone initially reaches 160 km instead of 80 km (Sim31), a subduction north of Adria is initiated at 25 Myr, after the Alpine subduction has ceased.
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Simulation results of Sim7, without any northward motion of Africa and Arabia. (a–d) Show the subduction dynamics at four different timesteps. The colors ranging from violet to green indicate the depth of the subducted mantle lithosphere (e–h) show the plate motion and the trench locations at the different timesteps. The colors of the arrows indicate the magnitude of plate motion.</figcaption></figure>

### 5.2.4. Extent of Ionian Crust

The extent and shape of the different crustal phases also has a significant impact on the simulation results. Modifying the extent of the oceanic corridor south of Adria influences the subduction dynamics. In Sim24, the oceanic corridor is narrower than in REF, which results in a slower retreat velocity of the ApCS. This, in turn, causes the subduction to cease retreating without reaching the same position as in REF. A wider corridor, however, facilitates the retreat of the Apennine-Calabrian subduction zone (Sim25).

### 5.2.5. No Pushing of the African and Arabian Plates

In Sim7 neither the African plate nor the Arabian plate move northward. The ApCS and the DHS exhibit almost identical dynamic behavior as in REF, however the Adriatic microplate does not move to the north at all. At around 24 Ma, Figures 9a and 9e shows that Adria undergoes an anticlockwise rotation. During this phase plate velocities are generally low. This changes for the first time at approximately 14 Ma, when the motion of the Adriatic plate shifts to a southwest direction (Figures 9b and 9f). This takes place at the same time as the slab break-off of the WMS at the boundary with Africa. The plate velocities of the entire Adriatic microplate increase between 14 Ma and 10 Ma (Figures 9f and 9g). Between 10 Ma and present-day both the ApCS and the DHS retreat toward the south of Adria and the plate motion of Adria changes to the southeast.

## 6. Discussion

### 6.1. Comparison of Reference Model With Tectonic Reconstructions

The reference model is intended to reproduce the first-order processes observed in the Mediterranean, as inferred from geological reconstructions. There are many model parameters that affect the dynamics of the simulations, and particularly the timing of the processes is frequently very sensitive to small changes in the rheology of the mantle. Yet, we can nevertheless learn some useful lessons from the first-order features observed in the simulations. In the following section, the temporal evolution of the subduction in the western Mediterranean, the Alpine-Carpathian subduction and the Dinaric-Hellenic subduction in REF are portrayed closer and placed into context with other work. In order to facilitate a comparison of other reconstructions with REF, Figure 10 illustrates the trench motion of the subduction zones in REF.
<figure><img src="image_12.png" /><figcaption>**Figure 10.** (a–f) Trench Evolution of REF at six different timesteps. The beige color shows continental crust, the orange color shows continental margin, the blue color shows oceanic crust and the brown color shows orogenic crust. The black lines indicate the position of the trench, while the black triangles represent the direction of the slab dip.</figcaption></figure>

In REF at 35 Ma, the WMS extends from Corsica to the Balearic Islands and then retreats to the south and to the west (Figures 10a and 10b). After 20 Ma, the slab breaks off at the northern African margin and the resulting slab gap propagates to the west and the east (Figures 10c and 10d). The subduction rollback of this subduction zone is in general agreement with studies such as Faccenna et al. (2014), Le Breton et al. (2021), Rosenbaum et al. (2002a) or Spakman and Wortel (2004). After the slab break-off, two different subduction systems can be identified: an eastward dipping Gibraltar subduction zone and a westward dipping subduction zone, which retreats toward Adria. The westward retreating Gibraltar subduction zone is supported by numerical studies (Chertova et al., 2014; Duarte et al., 2024) and by analyzing tomography studies of this area (e.g., Gutscher et al., 2002; Spakman & Wortel, 2004). The eastern part of the subduction zone retreats toward Adria, exhibiting a curved trench geometry at the present-day situation (Figure 10f). This is a consequence of accelerated slab retreat in the south, where oceanic crust is subducted. In general, the present-day trench geometry agrees well with that of other reconstructions, such as those presented in kinematic reconstructions like Rosenbaum et al. (2002a),
Faccenna et al. (2014) or the numerical models of Lo Bue et al. (2021). Rosenbaum et al. (2002a) also describe the motion of various continental blocks in the western Mediterranean. The Betic/Rif blocks exhibited a clockwise rotation as they move to the west, a phenomenon also observed in REF (Figure 10). In addition, Rosenbaum et al. (2002a) propose that the Kabylies blocks have moved to the south and accreted at Africa, which is also consistent with the REF model.

The Dinaric and Hellenic subduction system retreated in a southwest direction (e.g., Faccenna et al., 2014; Menant et al., 2016; Ring et al., 2010), which is similar to what can be observed in REF (Figure 10). In the REF model, the Hellenic subduction zone subducts continental crust between 35 Ma and 25 Ma (Figures 10a and 10b), after which the retreat is accelerated while subducting old oceanic crust south of Adria (Figures 10c–10f) (e.g., Menant et al., 2016; Pearce et al., 2012; Ring et al., 2010). The accelerated retreat of the Hellenic slab in the south leads to a curved trench geometry at present-day (Figure 10f). In REF, the Dinaric slab initially is much shorter than the southern Hellenic part of the slab. This is due to the slab being initially separated horizontally in the mantle, similar to what is proposed by Handy et al. (2015). This results in a more rapid retreat in the southern region when compared to the northern region.

Kinematic reconstructions of the Alpine-Carpathian subduction zone indicate that the Alpine-Carpathian slab underwent northward retreat between 35 Ma and 25 Ma, followed by a slab break-off in the center of the slab (e.g., Faccenna et al., 2014; Handy et al., 2015; Le Breton et al., 2021; van Hinsbergen et al., 2020). After 20 Ma, the slab retreated faster in the Carpathian region, which led to a curved subduction front (e.g., Faccenna et al., 2014; Handy et al., 2015). In REF, the Alpine-Carpathian subduction zone initially propagates to the north. Following a slab break-off event in the center, two distinct subduction systems retreat to the west and east, respectively. The western slab then begins to retreat to the northwest, and this process continues until it stops around 10 Ma (Figures 10a–10e). In the east, an acceleration of slab retreat starts around 25 Ma and lasts 10 Myr. The slab break-off propagates in a circular manner to the south, where the slab is finally completely detached at 10 Ma (Figures 10c–10e).

## 6.2. Dynamic Interactions Between Subduction Systems in the Mediterranean

A robust finding in all simulations is that mantle flow in the Mediterranean is predominantly toroidal, consistent with previous work (e.g., Kincaid & Griffiths, 2003; Z.-H. Li & Ribe, 2012; Schellart, 2024; Stegman et al., 2006), which show that toroidal flow is favored by narrow rollback subduction zones. The toroidal flow in the asthenosphere is primarily driven by dynamic pressure variations (e.g., Holt et al., 2017; Pusok & Stegman, 2019). The REF model shows that the different subduction systems active within a few hundred kilometers distance from each other, as is observed in the Mediterranean, result in a highly non-uniform horizontal dynamic pressure pattern. This results in a wide variety of dynamic interactions of the subduction systems, such as trench parallel flow south of the Alpine-Carpathian slab or asthenospheric escape flow underneath Adria. The retreat of the subduction systems is primarily buoyancy-driven. When comparing the retreat of the REF model with Sim7, where there is no convergence between Eurasia and Africa/Arabia, it is evident that the motion of the African and Arabian plates only slightly affects the rollback of the subduction zones in the Mediterranean.

The escape flow under Adria to the north and to the south is caused by the retreat of the subduction systems located to the east and west of Adria. The convergence of these subduction systems results in an increase in dynamic pressure between the slabs, which in turn gives rise to an escape flow in both the north and south directions. This observed trench-parallel flow is consistent with other numerical models of two outward-dipping subduction models (e.g., Di Leo et al., 2014; Király et al., 2021; Király, Holt, et al., 2018). In the north, the escape flow is driven by the increased pressure under Adria. The toroidal flow around the northern edge of the Dinaric slab is weaker than the northward-directed escape flow (Figure 6f). This is consistent with the SKS splitting results of various studies that have analyzed the present-day mantle flow (e.g., Salimbeni et al., 2022; Subašić et al., 2017). In the south, the escape flow has a strong toroidal component around the southern slab edge of the Calabrian slab (Figure 6f). This circular flow pattern around the slab edge has also been observed in other studies (e.g., Civello & Margheriti, 2004; Jolivet et al., 2009; Scarfi et al., 2023).

Kovács et al. (2012), Link and Rümpker (2023) or Petrescu et al. (2020) propose trench parallel asthenospheric flow in the Alpine region, which in our models can be justified by the observed change in flow direction over the past 35 Myr. The flow direction has shifted from a trench-perpendicular orientation prior to 25 Ma to a trench-parallel one in the Alpine region. This phenomenon can be explained by the difference in dynamic pressure between the north of the Adria in the west and the Carpathians in the east. In the western region, the dynamic
pressure is elevated between the Apenninic slab and the Dinaric slab. In contrast, in the eastern region, the dy-
namic pressure is low due to the rapid retreat of the Carpathian slab. This results in the trench-parallel flow in the
Alps, which supports the trench-perpendicular flow in the Carpathians. This, in turn, supports the retreat of the
subduction zone in the Carpathians, which has already been proposed by Kovács et al. (2012).

The slabs in REF are a simplified representation of the actual slab structure observed by seismic studies (e.g., Paffrath et al., 2021; Rappisi et al., 2022; Spakman & Wortel, 2004). As previously stated, the existence of slab windows has been identified east and west of Adria (e.g., Rappisi et al., 2022). With the opening of a slab window, toroidal flow cells develop around the slab edges (e.g., Király et al., 2020; Király, Faccenna, & Funiciello, 2018). This can also be observed in Sim16, where a slab window opens to the southwest of Adria. A consequence of this is an even more complex flow pattern and a change in trench geometry. This change in trench geometry has been investigated by numerical and analog models (e.g., Király et al., 2020; Király, Faccenna, & Funiciello, 2018). In the case of the Apenninic-Calabrian subduction zone this led to the formation of the Apenninic and Calabrian arcs, which explains the difference in trench geometry between REF and the observations. In REF the trench of the ApCS is straight and elongated, while observations show a curved trench of both the Apennine and Calabrian subduction (Faccenna et al., 2004). However, many of these slab windows appear to be a recent feature such as the slab window between the Apenninic and Calabrian slab which likely happened in the last 2 Myr (e.g., Faccenna et al., 2004; Gueguen et al., 1998; Király et al., 2020, 2021). The slab window between the Cyprus slab and the Hellenic slab reaches to about the same depth (Rappisi et al., 2022), with only the slab gap between the Hellenic and the Dinaric slabs reaching a greater depth of approximately 400 km (Rappisi et al., 2022). Given the aforementioned time constraints, it is likely that the overall dynamics of the Mediterranean has not been significantly altered by the presence of slab windows.

**6.3. Plate Motion of Adria and Driving Forces**

The plate motion of Adria over the past 35 Myr is thought to have been significantly influenced by the African plate (e.g., McKenzie, 1972). However, at some point, the direction of Adria's plate motion with respect to Africa changed (e.g., D'agostino et al., 2008). The reference model reproduces fairly well the present-day motion of the Adriatic plate described by geodetic data, which show a counterclockwise rotation of Adria relative to Europe, with a pole located to the northwest of the plate (D'agostino et al., 2008). By examining the plate motion of Adria in REF and the parameter study, we can see that this rotation and independent motion of Adria results from the interplay between the motion of the African plate and the dynamics of the surrounding subduction zones.

It can be observed that while Africa is constantly pushing, Adria begins to rotate with respect to Africa at a certain stage. This proceeds while the Calabrian subduction in the southwest of Adria and the Hellenic subduction in the southeast retreat toward Adria. In order to further understand this, the difference in direction between the southern part of Adria and Africa is plotted against the shortest horizontal distance between the Calabrian trench and Hellenic trench (Figure 11). This allows to compare different simulations and make them independent of time, but dependent on the subduction dynamics around Adria. The rotation of the Adriatic microplate with respect to Africa is initiated only when the trenches undercut a threshold of approximately 1000 km ± 100 km (Figure 11a). Prior to reaching this stage, Adria and Africa move in the same direction. Subsequently, the angle difference between Adria and Africa increases linearly as the trench distance decreases (Figure 11a). This can be explained by an interplay of various factors. The northern part of Adria already starts to slow down prior to the rotation of the southern part (Figure 8b). This can be attributed to the slowing down of the Alpine subduction zone to the north of Adria. During this period, the African plate is moving northwards at a distinctly faster rate than the northern part of the Adriatic plate. As the Apenninic-Calabrian and Hellenic-Dinaric subduction zones retreat toward Adria, the connection between Adria and Africa is becoming increasingly narrow. This results in Africa only pushing the southeastern part of Adria, where the northwestern part of Adria remains fixed (Figures 7b–7d). This geometric constellation forces Adria to rotate anticlockwise, in order to compensate for the northward movement of Africa. The slab pull generated by the surrounding slabs alone is insufficient to pull Adria to the northeast. This can be seen in Sim7 in which Africa is not pushing. In this case, the only factors that determine the plate motion of Adria are the surrounding subduction zones. If these subduction zones are within a few hundred kilometers distance from Adria, Adria's motion will be directed toward the southeast, rather than to the northeast (Figure 9h).
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Comparison of different simulations. The direction of “Adria south” minus the direction of “Africa” (Figure 8a for exact locations) is plotted against the horizontal distance of the Calabrian trench and the Hellenic trench. (a) Simulations with similar behavior are colored in gray, these are Sim1, 2, 6, 10, 11, 13, 14, 21, 22, 25, 27, 28. The plot also shows the REF simulation and Sim3, 4, 12 in color. (b) Shows a comparison of REF with Sim9 and Sim8. (c) Shows a comparison of REF with Sim23 and Sim29.</figcaption></figure>

However, the plate motion of Adria can still be partly influenced by the surrounding subduction zones. In specific cases, the linear trend in angle change of Adria may be interrupted, resulting in a sudden increase in the rate of change of the angle as can be observed in Sim3, Sim4, and Sim12 (Figure 11a). In these simulations, the friction angle of weak zone 4 and 5 is changed. That results in the occurrence of a slab break-off event of the Hellenic slab at the African margin. After this slab breaks-off, the slab pull forces of the Hellenic slab rearrange, which in turn results in an increased plate pull toward the Adriatic microplate. Consequently, the Hellenic slab has a more pronounced influence on the plate motion of the Adriatic microplate. Furthermore, the rate of change of Adria's motion direction is also affected by the velocity of Africa. The slower the motion of Africa, the faster the rotation of Adria with respect to Africa (see Sim8 and Sim9, where the inflow velocities of Africa are halved and quartered, respectively, Figure 11b). If Africa is slower, the slab break-off of the WMS at Africa happens later than in REF. This results in a slower retreat of the Calabrian slab. In that case, the Hellenic slab is always closer to Adria than the Calabrian slab, so that it has a stronger influence on the motion of the Adria and pulls it eastwards.

Another trend can be observed when looking at Figure 11c. Sim29 presents a simulation in which the initial slab tearing in the Dinaric slab is absent, resulting in a significantly larger slab volume. This strengthens the slab pull to the northeast and therefore increases the difference between Adria and Africa, allowing for the independent movement to start earlier. In Sim23 the thermal age of the slabs is increased to 120 Mry, which means that the slabs are cooler and denser than in REF. This accelerates the retreat of the slabs. Nevertheless, the retreat of the Apenninic-Calabrian slab is accelerated to a greater extent than that of the DHS, which in turn increases the slab pull force to the west of Adria and therefore it decreases the difference from Africa compared to REF.

A comparison of the REF simulation with kinematic reconstructions (see Supporting Information S1) reveals a lower absolute northward motion of Adria (e.g., Le Breton et al., 2017; van Hinsbergen et al., 2020). Le Breton et al. (2017) propose that during the rotation of Adria, there was still a northward translation of more than 100 km in the last 20 Ma. However in REF, Adria's motion toward the north ceases almost entirely during the rotational phase after 13 Ma, when the subduction stops retreating in the Alps. This also shows the important influence of the Alpine subduction zone on the northward motion of Adria. We note that a subduction of Adriatic lithosphere beneath Europe (e.g., in Sim31) to the north of Adria would pull Adria further to the north. In our simulations, such a subduction is only initiated when the weak zone to the north of Adria reaches a depth of at least 160 km and if the retreat of the Alpine slab has already ceased (Sim31) or if the Alpine slab in not retreating at all (Sim28). All these scenarios are however not consistent with geological reconstructions of plate convergence and recent interpretation of seismic tomography in the Alps (e.g., Handy et al., 2010, 2021).
## 6.4. Limitations of Our Modeling Approach

The dimensions of the model are finite, unlike nature. The large-scale flow pattern can be influenced by the dimensions of the model and the prescribed outflow at all boundaries of the box. Furthermore, in our models, the lower boundary, at a depth of 660 km, is closed and there are no interactions with a potential lower mantle. However, this could also influence the dynamics, by introducing poloidal flow cells that act across the lower and upper mantle (e.g., Schellart, 2024). Especially in the case of the Hellenic subduction system, the lower mantle may exert an influence on subduction dynamics, given that in REF at 35 Ma the Hellenic slab already reaches down to the 660 discontinuity.

Our model also lacks a westward-moving Anatolian plate, which is likely to change the plate motion of the plate southeast of Adria (e.g., Glerum et al., 2021; Serpelloni et al., 2022). This motion can also contribute to the observed plate motion south of Adria, which in this case is not only related to Africa but also to the Anatolian plate. Nevertheless, it should not significantly affect the dynamics of the model, as this overriding plate is not the primary driving force behind the Hellenic subduction; rather, it is the buoyancy forces of the slab itself. Another limitation is the eastern vertical boundary of the model domain. In our models, the Hellenic slab quickly detaches from the eastern boundary to rapidly roll-back westward. The retreat is enhanced by the slab pull of the detached slab. In reality, the slab remained attached to the east, reducing the slab pull force that facilitates its westward retreat in our models.

It is generally assumed that the convergence between Africa and Europe was active over the last 35 Myr, while the magnitude of this process may have changed over time. In their study Jolivet and Faccenna (2000) propose an increase in the convergence rate between Africa and Eurasia at 30 Ma, followed by a decrease after 20 Ma. Reilinger and McClusky (2011), on the other hand, propose two periods during which the convergence rate slowed down, one occurring approximately after 25 Ma and the other at 11 Ma. In order to account for the convergence of Africa and Eurasia in our simulations, a simplified process was applied by using a constant pushing velocity of Africa over time.

## 7. Conclusions

The present study employed 3D numerical modeling to investigate the subduction dynamics around the Adriatic microplate in the Alpine-Mediterranean region. We were able to reproduce the first-order trench retreats of the major subduction zones in the Alpine-Mediterranean, as well as the motion and rotation of the plates and continental blocks involved. The modeling results highlight the complex interactions between these subduction zones and their influence on mantle dynamics and the kinematics of the Adriatic microplate. It was possible to extensively explore a large parameter space by running more than 500 3D thermomechanical simulations.

The subduction zones around Adria, the Alpine-Carpathian subduction to the north, the Apennines-Calabrian subduction to the west and the Dinaric-Hellenic subduction to the east, exert a strong influence on the mantle flow. The principal findings with regard to mantle flow are as follows:

1. Toroidal flow is the dominant mode of mantle dynamics in the Mediterranean.

2. The retreat of the Dinaric-Hellenic subduction and the Apenninic-Calabrian subduction results in the formation of an escape flow under Adria to the north and to the south.

3. After 20 Ma, a trench parallel flow is observed south of the Alpine slab, driven by the dynamic pressure difference between the high dynamic pressure under Adria and the low dynamic pressure south of the fast retreating Carpathian slab.

4. The two toroidal flow cells, one situated around the southern Calabrian slab edge and the other around the southern Hellenic slab edge, induce a strong southward-directed flow under southern Adria.

Furthermore, the influence of surrounding subduction zones and the convergence between Africa and Eurasia on the plate motion of Adria was investigated. It was demonstrated that the northern part of Adria is primarily influenced by the interaction between the African plate motion and the subduction system in the Alps. Once the Alpine slab has ceased its northward retreat, the motion of Adria toward the north also ceases despite the continued northward movement of Africa.

An examination of the entire Adriatic microplate reveals that it begins to rotate with respect to Africa at a specific point. This rotation can be attributed to a number of factors. First, the northern part of the Adriatic plate decelerates
as a result of the ceasing of the Alpine slab. Conversely, Africa continues to move northwards, while the Calabrian subduction to the west and the Hellenic subduction to the east are approaching Adria, effectively narrowing the connection between Adria and the African continent. As the distance between the two trenches narrows (at about less than 1,000 km), the influence of Africa's push reduces to the southeastern portion of Adria. Since Adria is fixed in the northeast at this point, this results in a counterclockwise rotation of the microplate. The slab pull of the subduction zones surrounding Adria alone appears to not be sufficient to be the primary driver of Adria's plate rotation. Nevertheless, the closer the trenches are to the Adriatic microplate, the more profound the influence of the surrounding subduction zones on Adria's motion becomes. Abrupt changes in subduction dynamics, such as slab break-offs, can significantly increase the impact of the subduction zone on Adria's plate motion.

In conclusion, in a system where there has been a continuous convergence between Africa and Eurasia, the plate motion of Adria is primarily influenced by the interplay of three factors: (a) the convergence between the African and Eurasian plates, (b) the retreat of the Alpine subduction zone to the north, and (c) the distance between the retreating Calabrian and Hellenic subduction zones around Adria.

## Appendix A: Heterogeneities in the Lithosphere

The heterogeneities in the mantle lithosphere are created by the insertion of horizontal ellipses, with horizontal minor axes measuring 100 km in length and vertical minor axes measuring 20 km in length. In order to place them the model is subdivided into horizontal control areas of 100 km edge length each. In each of these control areas the x and y center points are randomly positioned within a 10 km inner square. The depth z is randomly assigned a value between 48 and 78 km.

## Appendix B: Temperature and Viscosity Structure

Here, we show the temperature and effective viscosity structure of a representative cross-section through the REF simulation at 10 Ma Figure B1.

<figure><img src="image_14.png" /><figcaption>**Figure B1.** Cross-sections through the REF model. (a) Mapview of REF simulation at 10 Ma (see Figure 3b). The purple line indicates the location of the cross-section. (b) Temperature cross-section. (c) Effective viscosity cross-section.</figcaption></figure>

## Acronyms

REF Reference model

WMS Western Mediterranean slab

AICS Alpine-Carpathian slab

DHS Dinaric-Hellenic slab

ApCS Apennine-Calabrian slab

GS Gibraltar slab
## Acknowledgments

We thank Claudio Faccenna for helpful discussions on Mediterranean tectonics and Andrea Piccolo who helped with the visualization. This paper greatly benefited from the constructive reviews of Manuele Faccenda and an anonymous reviewer. Paraview was used for the visualization of the modeling results (https://www.paraview.org/). The authors gratefully acknowledge the computing time granted on the supercomputer MOGON 2 at Johannes Gutenberg University Mainz (hpc.uni-mainz.de). CS acknowledges financial support provided by the German Research Foundation (DFG), SPP project Mountain Building in 4-Dimensions (4D-MB) through Grant KA3367/10-1. ELB acknowledges financial support from the German Research Foundation (DFG) through Grant BR4900/4-1. Open Access funding enabled and organized by Projekt DEAL.

## Data Availability Statement

The software to run the simulations is available on zenodo via Schuler et al. (2024). This repository contains the version of LaMEM (Kaus et al., 2016) and the geomIO package (Bauville & Baumann, 2019) that were used in the reference simulation, along with the input scripts that were used to run the simulation.

## References

Balay, S., Abhyankar, S., Adams, M. F., Benson, S., Brown, J., Brune, P., et al. (2024). *PETSC/TAO users manual revision* 3.21 (Tech. Rep.). Argonne National Laboratory (ANL).

Bauville, A., & Baumann, T. S. (2019). *geomio: An open-source matlab toolbox to create the initial configuration of 2-D/3-D thermo-mechanical simulations from 2-D vector drawings*. *Geochemistry, Geophysics, Geosystems*, *20*(3), 1665–1675. https://doi.org/10.1029/2018gc008057

Boutelier, D., & Cruden, A. (2017). Slab breakoff: Insights from 3D thermo-mechanical analogue modelling experiments. *Tectonophysics*, *694*, 197–213. https://doi.org/10.1016/j.tecto.2016.10.020

Chertova, M., Spakman, W., Geenen, T., Van Den Berg, A., & Van Hinsbergen, D. (2014). Underpinning tectonic reconstructions of the western mediterranean region with dynamic slab evolution from 3-D numerical modeling. *Journal of Geophysical Research: Solid Earth*, *119*(7), 5876–5902. https://doi.org/10.1002/2014jb011150

Civello, S., & Margheriti, L. (2004). Toroidal mantle flow around the Calabrian slab (Italy) from SKS splitting. *Geophysical Research Letters*, *31*(10). https://doi.org/10.1029/2004gl019607

Confal, J. M., Faccenda, M., Eken, T., & Taymaz, T. (2018). Numerical simulation of 3-D mantle flow evolution in subduction zone environments in relation to seismic anisotropy beneath the eastern mediterranean region. *Earth and Planetary Science Letters*, *497*, 50–61. https://doi.org/10.1016/j.epsl.2018.06.005

Crameri, F., Schmeling, H., Golabek, G., Duretz, T., Orendt, R., Buiter, S., et al. (2012). A comparison of numerical surface topography calculations in geodynamic modelling: An evaluation of the ‘sticky air’ method. *Geophysical Journal International*, *189*(1), 38–54. https://doi.org/10.1111/j.1365-246x.2012.05388.x

D'agostino, N., Avallone, A., Cheloni, D., D'anastasio, E., Mantenuto, S., & Selvaggi, G. (2008). Active tectonics of the Adriatic region from GPS and earthquake slip vectors. *Journal of Geophysical Research*, *113*(B12). https://doi.org/10.1029/2008jb005860

Dewey, J., Helman, M., Knott, S., Turco, E., & Hutton, D. (1989). Kinematics of the western mediterranean. *Geological Society, London, Special Publications*, *45*(1), 265–283. https://doi.org/10.1144/gsl.sp.1989.045.01.15

Di Leo, J., Walker, A., Li, Z.-H., Wookey, J., Ribe, N. M., Kendall, J.-M., & Tommasi, A. (2014). Development of texture and seismic anisotropy during the onset of subduction. *Geochemistry, Geophysics, Geosystems*, *15*(1), 192–212. https://doi.org/10.1002/2013gc005032

Duarte, J. C., Riel, N., Rosas, F. M., Popov, A., Schuler, C., & Kaus, B. J. (2024). Gibraltar subduction zone is invading the atlantic. *Geology*, *52*(5), 331–335. https://doi.org/10.1130/g51654.1

Faccenna, C., & Becker, T. W. (2010). Shaping mobile belts by small-scale convection. *Nature*, *465*(7298), 602–605. https://doi.org/10.1038/nature09064

Faccenna, C., Becker, T. W., Auer, L., Billi, A., Boschi, L., Brun, J. P., et al. (2014). Mantle dynamics in the mediterranean. *Reviews of Geophysics*, *52*(3), 283–332. https://doi.org/10.1002/2013rg000444

Faccenna, C., Becker, T. W., Lucente, F. P., Jolivet, L., & Rossetti, F. (2001). History of subduction and back arc extension in the central mediterranean. *Geophysical Journal International*, *145*(3), 809–820. https://doi.org/10.1046/j.0956-540x.2001.01435.x

Faccenna, C., Piromallo, C., Crespo-Blanc, A., Jolivet, L., & Rossetti, F. (2004). Lateral slab deformation and the origin of the western mediterranean arcs. *Tectonics*, *23*(1). https://doi.org/10.1029/2002tc001488

Gerya, T. (2019). *Introduction to numerical geodynamic modelling*. Cambridge University Press.

Glerum, A. C., Spakman, W., van Hinsbergen, D. J., Thieulot, C., & Pranger, C. (2021). Sensitivity of horizontal surface deformation to mantle dynamics from 3D instantaneous dynamics modeling of the eastern mediterranean region. *EarthArXiv eprints*, X5FW59.

Gueguen, E., Doglioni, C., & Fernandez, M. (1998). On the post-25 ma geodynamic evolution of the western mediterranean. *Tectonophysics*, *298*(1–3), 259–269. https://doi.org/10.1016/s0040-1951(98)00189-9

Gutscher, M.-A., Malod, J., Rehault, J.-P., Contrucci, I., Klingelhoefer, F., Mendes-Victor, L., & Spakman, W. (2002). Evidence for active subduction beneath Gibraltar. *Geology*, *30*(12), 1071–1074. https://doi.org/10.1130/0091-7613(2002)030<1071:efasbg>2.0.co;2

Handy, M. R., Giese, J., Schmid, S. M., Pleuger, J., Spakman, W., Onuzi, K., & Ustaszewski, K. (2019). Coupled crust-mantle response to slab tearing, bending, and rollback along the Dinaride-Hellenide orogen. *Tectonics*, *38*(8), 2803–2828. https://doi.org/10.1029/2019tc005524

Handy, M. R., Schmid, S. M., Bousquet, R., Kissling, E., & Bernoulli, D. (2010). Reconciling plate-tectonic reconstructions of alpine Tethys with the geological–geophysical record of spreading and subduction in the alps. *Earth-Science Reviews*, *102*(3–4), 121–158. https://doi.org/10.1016/j.earscirev.2010.06.002

Handy, M. R., Schmid, S. M., Paffrath, M., Friederich, W., & the AlpArray Working Group. (2021). Orogenic lithosphere and slabs in the greater alpine area-interpretations based on teleseismic p-wave tomography. *Solid Earth*, *12*(11), 2633–2669. https://doi.org/10.5194/se-12-2633-2021

Handy, M. R., Ustaszewski, K., & Kissling, E. (2015). Reconstructing the Alps–Carpathians–Dinarides as a key to understanding switches in subduction polarity, slab gaps and surface motion. *International Journal of Earth Sciences*, *104*, 1–26. https://doi.org/10.1007/s00531-014-1060-3

Harlow, F. H., & Welch, J. E. (1965). Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface. *Physics of Fluids*, *8*(12), 2182–2189. https://doi.org/10.1063/1.1761178

Hirth, G., & Kohlstedt, D. (2003). Rheology of the upper mantle and the mantle wedge: A view from the experimentalists. *Geophysical monograph-american geophysical union*, *138*, 83–106.

Holt, A., Royden, L., & Becker, T. (2017). The dynamics of double slab subduction. *Geophysical Journal International*, *209*(1), 250–265.

Jolivet, L., & Faccenna, C. (2000). Mediterranean extension and the Africa-Eurasia collision. *Tectonics*, *19*(6), 1095–1106. https://doi.org/10.1029/2000tc900018

Jolivet, L., Faccenna, C., & Piromallo, C. (2009). From mantle to crust: Stretching the mediterranean. *Earth and Planetary Science Letters*, *285*(1–2), 198–209. https://doi.org/10.1016/j.epsl.2009.06.017

Kaus, B. J. (2010). Factors that control the angle of shear bands in geodynamic numerical models of brittle deformation. *Tectonophysics*, *484*(1–4), 36–47. https://doi.org/10.1016/j.tecto.2009.08.042
Kaus, B. J., Popov, A. A., Baumann, T., Pusok, A., Bauville, A., Fernandez, N., & Collignon, M. (2016). Forward and inverse modelling of lithospheric deformation on geological timescales. *Proceedings of nic symposium, 48*, 978–983.

Kincaid, C., & Griffiths, R. (2003). Laboratory models of the thermal evolution of the mantle during rollback subduction. *Nature*, **425**(6953), 58–62. https://doi.org/10.1038/nature01923

Király, Á., Faccenna, C., & Funiciello, F. (2018). Subduction zones interaction around the adria microplate and the origin of the apenninic arc. *Tectonics*, **37**(10), 3941–3953. https://doi.org/10.1029/2018tc005211

Király, Á., Funiciello, F., Capitanio, F. A., & Faccenna, C. (2021). Dynamic interactions between subduction zones. *Global and Planetary Change*, **202**, 103501. https://doi.org/10.1016/j.gloplacha.2021.103501

Király, Á., Holt, A. F., Funiciello, F., Faccenna, C., & Capitanio, F. A. (2018). Modeling slab-slab interactions: Dynamics of outward dipping double-sided subduction systems. *Geochemistry, Geophysics, Geosystems*, **19**(3), 693–714. https://doi.org/10.1002/2017gc007199

Király, Á., Portner, D. E., Haynie, K. L., Chilson-Parks, B. H., Ghosh, T., Jadamec, M., et al. (2020). The effect of slab gaps on subduction dynamics and mantle upwelling. *Tectonophysics*, **785**, 228458. https://doi.org/10.1016/j.tecto.2020.228458

Kovács, I., Falus, G., Stuart, G., Hidas, K., Szabó, C., Flower, M., et al. (2012). Seismic anisotropy and deformation patterns in upper mantle xenoliths from the central Carpathian–Pannonian region: Asthenospheric flow as a driving force for cenozoic extension and extrusion? *Tectonophysics*, **514**, 168–179. https://doi.org/10.1016/j.tecto.2011.10.022

Kumar, A., Cacace, M., Scheck-Wenderoth, M., Götzé, H.-J., & Kaus, B. J. (2022). Present-day upper-mantle architecture of the alps: Insights from data-driven dynamic modeling. *Geophysical Research Letters*, **49**(18), e2022GL099476. https://doi.org/10.1029/2022GL099476

Le Breton, E., Brune, S., Ustaszewski, K., Zahirovic, S., Seton, M., & Müller, R. D. (2021). Kinematics and extent of the Piemont–Liguria basin–implications for subduction processes in the alps. *Solid Earth*, **12**(4), 885–913. https://doi.org/10.5194/se-12-885-2021

Le Breton, E., Handy, M. R., Molli, G., & Ustaszewski, K. (2017). Post-20 ma motion of the adriatic plate: New constraints from surrounding orogens and implications for crust-mantle decoupling. *Tectonics*, **36**(12), 3135–3154. https://doi.org/10.1002/2016tc004443

Li, X. S., & Demmel, J. W. (2003). Superlu Dist: A scalable distributed-memory sparse direct solver for unsymmetric linear systems. *ACM Transactions on Mathematical Software*, **29**(2), 110–140. https://doi.org/10.1145/779359.779361

Li, Z.-H., & Ribe, N. M. (2012). Dynamics of free subduction from 3-D boundary element modeling. *Journal of Geophysical Research*, **117**(B6). https://doi.org/10.1029/2012jb009165

Link, F., & Rümpker, G. (2023). Shear-wave splitting reveals layered-anisotropy beneath the european alps in response to mediterranean subduction. *Journal of Geophysical Research: Solid Earth*, **128**(9), e2023JB027192. https://doi.org/10.1029/2023jb027192

Lo Bue, R., Faccenda, M., & Yang, J. (2021). The role of adria plate lithospheric structures on the recent dynamics of the central mediterranean region. *Journal of Geophysical Research: Solid Earth*, **126**(10), e2021JB022377. https://doi.org/10.1029/2021jb022377

Macherel, E., Räss, L., & Schmalholz, S. M. (2024). 3D stresses and velocities caused by continental plateaus: Scaling analysis and numerical calculations with application to the Tibetan plateau. *Geochemistry, Geophysics, Geosystems*, **25**(3), e2023GC011356. https://doi.org/10.1029/2023gc011356

Martinod, J., Daou, A.-C., Métral, L., & Sue, C. (2024). Did subduction in the western mediterranean drive Neogene alpine dynamics? Insights from analogue modeling. *BSGF-Earth Sciences Bulletin*, **195**(1), 5. https://doi.org/10.1051/bsgf/2024001

McKenzie, D. (1972). Active tectonics of the mediterranean region. *Geophysical Journal International*, **30**(2), 109–185. https://doi.org/10.1111/j.1365-246x.1972.tb02351.x

Menant, A., Jolivet, L., & Vrielynck, B. (2016). Kinematic reconstructions and magmatic evolution illuminating crustal and mantle dynamics of the eastern mediterranean region since the late cretaceous. *Tectonophysics*, **675**, 103–140. https://doi.org/10.1016/j.tecto.2016.03.007

Muttoni, G., Garzanti, E., Alfonsi, L., Cirilli, S., Germani, D., & Lowrie, W. (2001). Motion of Africa and Adria since the permian: Paleomagnetic and paleoclimatic constraints from northern Libya. *Earth and Planetary Science Letters*, **192**(2), 159–174. https://doi.org/10.1016/s0012-821x(01)00439-3

Paffrath, M., Friederich, W., Schmid, S. M., Handy, M. R., & the AlpArray and AlpArray-Swath D Working Group. (2021). Imaging structure and geometry of slabs in the greater Alpine area—A P-wave travel-time tomography using AlpArray seismic network data. *Solid Earth*, **12**(11), 2671–2702. https://doi.org/10.5194/se-12-2671-2021

Pearce, F. D., Rondenay, S., Sachpazi, M., Charalampakis, M., & Royden, L. H. (2012). Seismic investigation of the transition from continental to oceanic subduction along the western hellenic subduction zone. *Journal of Geophysical Research*, **117**(B7). https://doi.org/10.1029/2011jb009023

Peral, M., Fernandez, M., Vergés, J., Zlotnik, S., & Jiménez-Munt, I. (2022). Numerical modelling of opposing subduction in the western mediterranean. *Tectonophysics*, **830**, 229309. https://doi.org/10.1016/j.tecto.2022.229309

Petrescu, L., Pondrelli, S., Salimbeni, S., & Faccenda, M. (2020). Mantle flow below the central and greater alpine region: Insights from SKS anisotropy analysis at Alparray and permanent stations. *Solid Earth Discussions*, **2020**(4), 1–26. https://doi.org/10.5194/se-11-1275-2020

Piromallo, C., & Morelli, A. (2003). P wave tomography of the mantle under the alpine-mediterranean area. *Journal of Geophysical Research*, **108**(B2). https://doi.org/10.1029/2002jb001757

Popov, A., & Sobolev, S. V. (2008). Slim3d: A tool for three-dimensional thermomechanical modeling of lithospheric deformation with elasto-visco-plastic rheology. *Physics of the Earth and Planetary Interiors*, **171**(1–4), 55–75. https://doi.org/10.1016/j.pepi.2008.03.007

Pusok, A. E., & Stegman, D. R. (2019). Formation and stability of same-dip double subduction systems. *Journal of Geophysical Research: Solid Earth*, **124**(7), 7387–7412. https://doi.org/10.1029/2018jb017027

Ranalli, G. (1995). *Rheology of the earth*. Springer Science & Business Media.

Rappisi, F., VanderBeek, B., Faccenda, M., Morelli, A., & Molinari, I. (2022). Slab geometry and upper mantle flow patterns in the central mediterranean from 3D anisotropic p-wave tomography. *Journal of Geophysical Research: Solid Earth*, **127**(5), e2021JB023488. https://doi.org/10.1029/2021jb023488

Reilinger, R., & McClusky, S. (2011). Nubia–Arabia–Eurasia plate motions and the dynamics of mediterranean and middle east tectonics. *Geophysical Journal International*, **186**(3), 971–979. https://doi.org/10.1111/j.1365-246x.2011.05133.x

Ring, U., Glodny, J., Will, T., & Thomson, S. (2010). The hellenic subduction system: High-pressure metamorphism, exhumation, normal faulting, and large-scale extension. *Annual Review of Earth and Planetary Sciences*, **38**(1), 45–76. https://doi.org/10.1146/annurev.earth.050708.170910

Rosenbaum, G., Lister, G. S., & Duboz, C. (2002a). Reconstruction of the tectonic evolution of the western mediterranean since the oligocene. *Journal of the Virtual Explorer*, 8(January), 107–130. https://doi.org/10.3809/jvirtex.2002.00053

Rosenbaum, G., Lister, G. S., & Duboz, C. (2002b). Relative motions of Africa, Iberia and Europe during alpine orography. *Tectonophysics*, **359**(1–2), 117–129. https://doi.org/10.1016/s0040-1951(02)00442-0

Rosenbaum, G., Lister, G. S., & Duboz, C. (2004). The mesozoic and cenozoic motion of Adria (central mediterranean): A review of constraints and limitations. *Geodinamica Acta*, **17**(2), 125–139. https://doi.org/10.3166/ga.ga.17.125-139