SENSITIVITY OF HORIZONTAL SURFACE DEFORMATION TO MANTLE DYNAMICS FROM 3D INSTANTANEOUS DYNAMICS MODELING OF THE EASTERN MEDITERRANEAN REGION ANNE GLERUM* Helmholtz Centre Potsdam-GFZ German Research Centre for Geosciences, Germany

Geodetically estimated surface motions contain contributions to crustal deformation from coupled geodynamic processes active at all spatial scales and constitute key data for lithosphere dynamics research. Data interpretation methods should therefore account for the full range of possible processes, otherwise risking misinterpretation of data signal and incorrect estimation of lithosphere rheology, stress, or deformation fields. Here we explore the sensitivity of surface deformation to sub-lithospheric processes such as viscous plate-mantle and slab-mantle coupling, variations in slab pull, and buoyancy-driven mantle flow. To this end, we perform 3D instantaneous-dynamics numerical modelling of an elaborately structured compressible crust-mantle system designed for the Eastern Mediterranean Aegean-Anatolian region. We first determine a reference model driven by the absolute motions of the major plates, regional slab pull, a 3D mantle buoyancy field, and modulated by plate boundary coupling and mantle viscosity. The RMS motion data fit of ∼5.9 mm/yr of predicted and observed Aegean-Anatolian horizontal surface motions demonstrates that the bulk amplitude of surface motion can be explained by these combined mantle processes. Next, by systematically perturbing reference model features, we assess the crustal sensitivity to each geodynamic driver and to mantle rheology. We find significant changes in crustal velocity gradient amplitudes, often between 10% and 40% of the reference model, with slab morphology effects of up to 93%. This demonstrates the key importance of carefully accounting for each process in modelling lithosphere dynamics. For the Aegean-Anatolia region, we present geodynamic evidence that the Aegean slab pull is the primary driver of the crustal motion field, as was previously suggested from kinematic analysis. Key points • Regional data-driven instantaneous-dynamics modeling of compressible flow predicts surface deformation of the Aegean-Anatolian region. • Sensitivity of predicted crustal deformation to subducting plate morphology, mantle buoyancy structure, mantle rheology and absolute plate motions is investigated. • Modeled Aegean-Anatolian crustal deformation is most sensitive to slab morphology, substantially modulated by sublithospheric mantle rheology and buoyancy, and carries least signal of absolute plate motions.


Introduction
Deformation of Earth's crust is not only an expression of lithospheric drivers but also of geodynamic processes acting in the underlying mantle. These drivers of deformation include plate tectonic forcing (ridge push, plate coupling, slab pull), plate-mantle coupling (Forsyth and Uyeda, n.d.), slab-mantle coupling (Uyeda and Kanamori, 1979), and stress in the continental lithosphere from internal density variation and topography (Fleitout andFroideveaux, 1982, 1983). Geodetic observations are an important source of information on surface deformation and, implicitly, on the various geodynamic drivers presently deforming the crust and mantle. Many investigations therefore employ numerical modelling to identify the driving processes of surface deformation by optimizing the fit between observed and predicted surface motions (Flesch et al., 2001, Hu et al., 2004, Ghosh et al., 2008, Jadamec and Billen, 2010, 2012, Nijholt et al., 2018, Osei Tutu et al., 2018a. When comparing different modeling studies of the same region, the best-performing forward models resulting from such data fitting procedures can however show disagreement about the geodynamic processes that cause deformation. For example, studies of the Aegean subduction system of the Eastern Mediterranean region have put forward different dominant drivers of the upper plate deformation in Aegean-Anatolian region, such as slab rollback and slab-induced mantle flow (e.g. Sternai et al., 2014), slab pull and mantle drag from upwellings (e.g. Boschi et al., 2010, Becker and Faccenna, 2011, Faccenna et al., 2014, or lateral gradients in gravitational potential energy (Özeren and Holt, 2010, Carafa et al., 2015, England et al., 2016. These approaches differ foremost in the representation or treatment of the geodynamic processes that may contribute to surface motions, which may affect their modelled contribution to the overall data fit and consequently may bias the determination of preferred models (e.g. Oreskes et al., 1994).
The significance of lithospheric contributions to surface motion and deformation is most clear, particularly, that of gravitational potential energy gradients due to lateral density variations (e.g Fleitout andFroidevaux, 1982, England et al., 2016), and of tectonic plate coupling across plate boundary faults (e.g. Warners-Ruckstuhl et al., 2012). Much less clear is the particular regional significance of the contributions of larger-scale mantle processes comprising the viscous coupling of the lithosphere to buoyancy-driven mantle flow, slab pull, and the viscous slab-mantle coupling as driven by absolute plate motion. The magnitude of these geodynamic processes, and thus the effect these may have on crustal deformation, is determined by the properties of the mantle system, such as the mantle rheology, 3D slab morphology, and the 3D makeup of the mantle density and temperature fields, which are all still rather uncertain. Here, we aim to clarify the significance of sublithospheric processes and mantle properties for driving crustal deformation. To this end, we take the extensively-studied Aegean-Anatolian region as our natural laboratory and investigate the central question of how sensitive observed surface motions and associated crustal deformation are to each of the deep geodynamic processes.
To estimate the sensitivity of crustal motion predictions to mantle drivers, we focus on surface motions as observed by the Global Navigation Satellite System (GNNS) in the Aegean-Anatolian region of the Eastern Mediterranean. We design regional numerical experiments in a 3D domain of 40 x 40 degrees, centred on the Aegean plate and including the entire mantle below the region. These experiments include detailed representations of the laterally varying crust-lithosphere system (e.g. Ganas and Parsons, 2009) and vary the far-field plate motions (e.g. Petricca et al., 2013), the morphology of the subducting plate, and the seismic-tomography-based mantle buoyancy field (e.g. Boschi et al., 2010). We model compressible mantle flow employing varying nonlinear viscoplastic rheologies for which material properties including phase changes are precomputed with thermodynamics software. Our experiments focus on computing the sensitivity of surface motions and the associated velocity gradient field to changes in plate motions, slab-mantle coupling, slab rollback, slab tearing, mantle flow, and rheology, all with respect to a reference model constrained by GNNS observations. This allows us to rank the investigated mantle drivers in order of importance for explaining crustal deformation. Apart from the generic importance of our findings for the deep causes of plate boundary deformation, we can also infer the main drivers of Aegean-Anatolian crustal deformation in the overall plate-convergent setting of the eastern Mediterranean region.
In the following, we first describe the geodynamic setting of the Eastern Mediterranean that informed our model design. A description of the numerical modelling methods and model setup follows, together with an explanation of the postprocessing of the model data and its comparison to GNNS observations. We then present four sets of modeling experiments and individual discussions of their results. Lastly, we discuss the specific geodynamic drivers of deformation in Aegean-Anatolian platelet of the Eastern Mediterranean region.

Geodynamic setting of the eastern Mediterranean
In this section, we first summarise the present-day geodynamic complexity and structural heterogeneity of the eastern Mediterranean region, as inherited from a long geological evolution caught between the converging African/Arabian and Eurasian tectonic plates. Next, we determine constraints on the morphology of the involved subducting slabs that are relevant for the setup of our numerical modeling of the present-day geodynamics. Finally, we present the eastern Mediterranean crustal flow field as observed by GNSS studies. Determining the sensitivity of this flow field for various geodynamic drivers is this study's main goal.
Africa-Europe convergence was accommodated by (at least) two subduction zones in Anatolia (Şengör and Yilmaz, 1981, Menant et al., 2016. The Central Anatolian back-arc basin was governed by E-W extension  and subsequently N-S shortened and thickened during oroclinal bending (Lefebvre et al., 2013, Isik et al., 2014, Advokaat et al., 2014. In the last ∼10 Myr, central-eastern Anatolia did not experience much internal deformation while moving westwards relative to Eurasia along the North Anatolian Fault (NAF) zone (e.g. Şengör et al., 2005).
Due to the crustal nappe orogen, lithospheric thicknesses in the Aegean-Anatolian region are small (60-80 km, e.g. Endrun et al., 2011). The lithosphere has been further thinned by extension, which in the Aegean region has been widely attributed to solely rollback of the Aegean slab relative to the upper Eurasian plate (Le Pichon and Angelier, 1979, Wortel and Spakman, 2000. This however assumes that the Eurasian plate is mantle-stationary. The Global Moving Hotspot Reference Frame (GMHRF) suggests that Eurasia moved NE-ward during the Cenozoic at ∼1 cmyr −1 (Fig. 2; Doubrovine et al., 2012), such that part of the Aegean extension could also be related to upper plate escape (e.g. Lallemand et al., 2005, Schellart, 2008, van Hinsbergen et al., 2015. Gögüş Figure 1. Plate tectonic setting of the Nubia-Arabia and Eurasia convergence zone. Red lines denote the plate boundaries of Bird (2003), black arrows indicate the relative plate velocities with respect to Eurasia (Kreemer et al., 2014). While our 3D model domain spans the area shown, the dashed black box outlines the area of the eastern Mediterranean of which we present model results. Abbreviations: BS -Balearic Sea, Ls -Ligurian Sea, TS -Tyrrhenian Sea, IS -Ionian Sea, Hs -Hellenides, AS -Aegean Sea, HT -Hellenic Trench, NAT -North Aegean Trough, NAF -North Anatolian Fault, EAF -East Anatolian Fault, DSF -Dead Sea Fault. et al. (2017) suggested that the newly grown mantle lithosphere below Central Anatolia was removed by dripping following the post-Eocene oroclinal bending.
In short, the laterally variable subduction-collision evolution of the eastern Mediterranean has culminated into a rheologically and compositionally strongly heterogeneous crust-mantle system. Our numerical experiments focus on the mantle and include the first-order thickness and composition variation of the heterogeneous lithospheric system, while we test for the effect of sublithospheric mantle structure (e.g. slabs) on crustal motions.
Velocity scale: 10 mm/year  Kreemer et al. (2014) shown here in the mantle frame defined by the absolute Global Moving Hotspot Reference Frame (GMHRF) motion of Nubia of the last 10 My (Doubrovine et al., 2012). Specifically, we use as our mantle anchor the absolute GMHRF Euler pole of the Nubian plate as averaged over the past 10 My, which lies at 33.3 • W, 36.7 • N, with a counterclockwise rotation of 0.161 • My −1 . The GNSS field is obtained by applying this Euler rotation to the GNSS velocity vectors in the Africa-fixed frame as determined by Kreemer et al. (2014).

2.2.
Morphology of the subducting Nubian lithosphere. Several morphological features of the Aegean slab that have been interpreted from tomographic models may affect the impact of slab dynamics on crustal deformation. Spakman et al. (1988) inferred partial horizontal detachment of the Aegean slab below SW Greece and suggested a tear propagating southeastward. Later tomographic models did not clearly corroborate such a detachment (Wortel and Spakman, 2000, van der Meer et al., 2018 and point to one continuous slab under the Aegean region reaching a depth of about 1500 km (Spakman et al., 1993, Bijwaard et al., 1998, van Hinsbergen et al., 2005.
Vertical segmentation of the Aegean slab was proposed in the region of the Kefalonia Fault and Gulf of Corinth (Govers and Wortel, 2005, Suckale et al., 2009, Bocchini et al., 2018, thought to originate from a lateral change in slab density upon transitioning from oceanic to continental crust (Suckale et al., 2009, Royden andPapanikolaou, 2011), although Halpaap et al. (2018) suggest a laterally smooth transition between the two types of lithosphere. We note that directly north of the Kefalonia zone, a slab is still imaged across the entire upper mantle (e.g. Spakman et al., 1993, Suckale et al., 2009, Handy et al., 2019.
In contrast, at the eastern end of the Aegean slab an actual slab edge is imaged that has been interpreted as a vertical slab rupture (De Boorder et al., 1998, van Hinsbergen et al., 2010, Biryol et al., 2011, Govers and Fichtner, 2016 associated with a Subduction-Transform Edge Propagator (STEP) fault (Govers and Wortel, 2005) of which the Pliny-Strabo trenches were considered the surface expression (Zachariasse et al., 2008, van Hinsbergen et al., 2010, Biryol et al., 2011,Özbakir et al., 2013. However, Bocchini et al. (2018) showed that the slab window instead opened as a triangular gap, suggesting a relationship with the kink in the trench when going from the Aegean to the Antalya-Cyprus slabs. This tear-region below western Anatolia would have separated the Aegean slab from the Antalya and Cyprus slabs (Biryol et al., 2011, van der Meer et al., 2018.
The upper mantle section of the Antalya and Cyprus slabs likely consists of > 200 My old (Granot, 2016) oceanic lithosphere of the eastern Mediterranean basin and may be detaching N to NW of Cyprus (Portner et al., 2018). To the west, a connection with the Antalya slab is debated (Biryol et al., 2011, Koç et al., 2016, McPhee et al., 2018. In our experiments, we will test the sensitivity of surface motions to horizontal and vertical tearing of the Aegean slab, and even absence of this slab, while keeping the geometry of the Antalya-Cyprus slabs the same. 2.3. The eastern Mediterranean crustal flow field. To properly account for the dynamic coupling between slab, mantle and plates, we use a modern deep-mantle reference frame, specifically the GMHRF of Doubrovine et al. (2012), to represent observed as well as predicted plate motions and crustal velocities relative to the mantle (Fig. 2). In this frame, the absolute motion of Nubia just south of the Aegean region is ∼1.2 cmyr −1 to the NNE, Arabia is moving ∼N at ∼2 cmyr −1 at the plate boundary with eastern Anatolia, and Eurasia is moving to the NE at ∼1 cmyr −1 in the southwest Black Sea region, i.e. away from the Africa-Eurasia trench at about half the Aegean roll back rate as indicated by the GNSS vectors (Fig. 2). The NW-SE relative motion between Africa and Eurasia is ∼5 mmyr −1 (see also Fig. 1).
Aegean and Anatolian motions differ greatly from the motions of the larger plates (Mc-Clusky et al., 2000, Reilinger et al., 2006, Nocquet, 2012, Kreemer et al., 2014). An overall counterclockwise (ccw) rotation (Fig. 2) is facilitated by strike-slip along both the NAF and the East Anatolian Fault (EAF) (Fig. 1, Reilinger et al. 2006). The rotational motion field accelerates towards the Hellenic trench (Reilinger et al., 2006, Le Pichon and Kreemer, 2010, Nocquet, 2012 from about 13 mmyr −1 to 20 mmyr −1 (Fig. 2). The idea of the present-day westward Anatolian motion being caused by a 'push' from the Arabian plate (e.g. McKenzie, 1972) is invalidated by this acceleration towards the trench and the almost pure strike-slip motion in the absence of fault-normal compression along the EAF (Reilinger et al., 2006). The ccw motion of the Aegean-Anatolian region and the upper plate extension in the Aegean region have otherwise been associated with the following processes, which may work in parallel, and which may represent an incomplete list of drivers: asthenospheric drag (e.g. McKenzie, 1978, Le Pichon and Kreemer, 2010, gravitational spreading due to elevation differences between the Aegean Sea and the Mediterranean Sea and within the Aegean and Anatolian regions (e.g. Le Pichon and Angelier, 1979, Gautier et al., 1999, England et al., 2016 and Hellenic slab rollback (Le Pichon and Angelier, 1979, Spakman et al., 1988, Wortel and Spakman, 2000, Reilinger et al., 2006, Ring et al., 2010, Royden and Papanikolaou, 2011.
We note that adopting the GMHRF Africa motion of the past 10 My as our mantle anchor is not crucial for our experiments. Important is that we define a mantle frame such that the motions of plates and slabs and the implied viscous coupling between slab, plates and mantle, are physically consistent. Alternatively, one can adopt a reference plate, e.g., Eurasia, to be fixed and use plate motions relative to Eurasia as well as impose a uniform mantle wind on all sides with the amplitude of, and direction opposite to, the Eurasian absolute plate motion. We prefer to define a mantle frame with corresponding absolute plate motions as it is physically more intuitive.

Methods
Our workflow for estimating the sensitivity of surface observables to mantle dynamics is centered around the modeling of instantaneous compressible flow. We first describe our numerical modeling tool ASPECT; then the model setup concerning the preparation of model input from observational data to include the complex dynamic processes discussed in Section 2; and last the analysis of the resulting predictions of mantle and crustal deformation through comparison to observations and the reference model.

Instantaneous flow modeling with ASPECT.
ASPECT is an open source, massively parallel, finite element code (Kronbichler et al., 2012 based on the libraries deal.II (Bangerth et al., 2007) (Finite elements), Trilinos (Heroux et al., 2005, Heroux, 2017) (linear algebra) and p4est (Burstedde et al., 2011) (Mesh generation). We use ASPECT (version 1.5.0-pre) to solve the "isothermal compression" formulation (see Bangerth et al., 2017 of the conservation equations of momentum, mass and energy and the advection equations for each Eulerian compositional field c i used to track different materials: Here˙ is the strain rate tensor 1 2 (∇v + (∇v) T ), C represents the compressibility 1 ρ ∂ρ ∂P and ρ is the full density, dependent on the full pressure P and temperature T through look-up tables generated by the thermodynamics software Perple X as outlined in Section 3.2.4. Other symbols are defined in Section 3.2.4 and below. We neglect the contributions of latent heat and radioactive decay to temperature, as they can be considered second order effects on the short timescales we are interested in. The translation of observational data to initial and boundary conditions required to solve equations (1)-(4) is described in Section 3.2.
3.2. Model setup. The model setup defines the driving forces of the instantaneous crustmantle flow. From the geodynamic setting of the eastern Mediterranean (Section 2), we identified five potentially-key influences on the instantaneous dynamics of the region: the geometry and nature of the plates and their boundaries; the geometry of the subducting slabs; the mantle structure in terms of temperature, here determining the density distribution; the rheology of the plates and mantle; and the motions of the plates.
3.2.1. Model domain. Our model domain constitutes a part of a 3D spherical shell, for which it is convenient to rotate the conventional geocentric coordinate system such that the equator of the new coordinate system runs through the center of our target region, the Aegean plate. Therefore, we have transformed all data used for creating the model setup in the region between (1 • E, 21 • N), (41 • E, 16 • N), (62 • E, 52 • N) and (-6 • E, 60 • N) (Fig. 1) to a new coordinate system in which the point (25 • E, 40 • N) above the Aegean slab represents the new origin through which the new equator runs eastward with an azimuth of N100 • E. The computational domain then runs from -20 • to 20 • in both longitude and latitude in the new coordinate system (Fig. 3a) and reaches a depth of 2900 km, as shown in Fig. 3d.
The model domain is discretized using a radially refined hexahedral mesh, going from a mesh size of 1.25 • x 1.25 • x 86.56 km in the lower mantle to 0.16 • x 0.16 • x 8.09 km in the crust.
3.2.2. Synthetic plate and slab geometry: initial composition conditions. We construct a synthetic 3D geometry of the plates, plate boundaries and slabs as shown in Figs. 3a-3d in the stand-alone program Regional Tectonic Geometry Builder (the successor of the tool developed by Pranger, 2014). Our choices are based on the overall tectonic setting (e.g. Faccenna et al., 2014, Fig. 1; Section 2), earthquake hypocenters (NCEDC, 2014), topography (Amante and Eakins, 2009) and the plate boundaries of Bird (2003) and Le Pichon and Kreemer (2010).
We mimic the plate boundary fault zones by constructing discrete, 30 km wide weak zones to the depth of the thickest lithospheric plate in the model. This is a practical choice for simulating major lithosphere cutting faults in the computational domain and may provide a fair approximation for some major faults, e.g. the North Anatolian Fault (Frederiksen et al., 2015, Taylor et al., 2016, Papaleo et al., 2018, but may overestimate the width of subduction channels (e.g. Vannucchi et al., 2012). Subduction channels are implemented for the Dinarides and the Aegean and Cyprus slabs to dip parallel to the subducting plate surface as interpreted from earthquake hypocenters. We note that the rheology of these model plate boundary faults is tunable such that a relatively wide, higher viscosity fault zone may reproduce the regional dynamic coupling that would occur across a lower viscosity, thin fault zone. The traces of the plate boundaries are shown in blue in Fig. 3a, their 3D structure is presented in Fig. 3b.
Crust and mantle lithosphere thicknesses (35 and 120 km, resp.) and regions of different thicknesses within plates and boundaries (Tab. 1) were abstracted from Moho and LAB maps (Carafa et al., 2015, Motavalli-Anbaran et al., 2016 and the ocean floor age map of Müller et al. (2008). We include a 200 My old oceanic region in the Ionian and Levantine Seas (e.g. Meier et al., 2007, Suckale et al., 2009, Speranza et al., 2012, D'Alessandro et al., 2016, a very young oceanic region in the Liguro-Provencal basin (Bayer et al., 1973) and the Tyrrhenian Sea (Nicolosi et al., 2006) and thinned continental lithosphere in the Aegean-Anatolian region and the Black Sea (green and red areas in Fig. 3a). Transitions in crustal and lithospheric thickness from the plates to these regions and vice versa are smoothed with an error function with a halfwidth of 60 km.
Instead of working with the tomographically imaged slab structure, we prefer to work with synthetic slab geometries which gives us control on the slab morphology and its hypothesized tearing. Based on the seismic tomography model UU-P07 of Amaru (2007, available at http://www.atlas-of-the-underworld.org, van der Meer et al., 2018) and earthquake hypocenters (NCEDC, 2014), we constructed 3D volumes of the Aegean slab and the Antalya-Cyprus slab in the upper mantle by picking the slab surface on trench-(semi)perpendicular slices (Fig. 3b, similar to Ganas and Parsons, 2009, Jadamec and Billen, 2010, Alisic et al., 2012. The volume of the slab is then determined by assuming a uniform thickness of 120 km for the at least 200 My old lithosphere involved in subduction (Fig. 3c). Our synthetic slabs are therefore a smooth, condensed form of the tomographic slabs seen in the UU-P07 model. In several experiments, we superpose on this volume lateral and vertical tears to test the sensitivity of surface motions to slab detachment and slab window hypotheses.
The Adria slab underneath the Dinarides is confined to the uppermost ∼150 km in the UU-P07 tomography model and is approximated by a tilted plate boundary up to a depth of 120 km only. Seismic tomography anomalies otherwise recording this slab, as well as the lower-mantle signature of the Aegean and Cyprus slabs, are incorporated in the thermal structure of the mantle (see Section 3.2.3).
The resulting crustal and lithospheric volumes are then represented on a grid of 0.13 • x 0.13 • x 6.7 km resolution. ASPECT reads in these gridded volumes as compositional fields c i using a trilinear interpolation between grid points.
Note that we do not include pressure gradients that may result from surface topography and detailed internal lateral density contrasts. These have been investigated to some extent before in thin-sheet approximations of lithosphere deformation (e.g. Carafa and Bird, 2016, England et al., 2016,Özbakir et al., 2017. Instead we keep the lithosphere laterally homogeneous, apart from variations in crust and lithosphere thickness and in reference density between oceanic and continental domains. Hence, deep lateral gravitational potential energy contrasts resulting from undulations in Moho and lithosphere-asthenosphere boundary depth are intrinsical to our model. 3.2.3. Plate, slab and mantle structure: initial temperature conditions. The prescribed initial temperature conditions depend on the radial depth z of a point r and the type of the compositional field c i at r (i.e. continental vs. oceanic plate, slab or mantle). The lithospheric temperature profile T L (r, c i ) consists of a linear temperature gradient in a continental plate, the plate cooling model inside an oceanic region (Schubert et al., 2001)    a plate cooling profile in the subducting plates according to McKenzie (1970), see Eq. (14) in App. A. For the sublithospheric mantle, we first construct a 1D background temperature profile T bg (z), consisting of a uniform temperature T a that is prescribed to the depth of the deepest lithospheric plate z Lmax and is the adiabatic temperature at this depth, and an adiabatic reference profile T adiabat (z) that is perturbed by a thermal boundary layer at the coremantle boundary (CMB) δT tbl (z) (see App. A for detailed equations and Fig. 26a for a typical temperature profile).
We further perturb the background temperature T bg (z) with spatially varying temperature anomalies by converting the seismic anomalies of different tomographic models to temperature anomalies through depth-dependent seismic velocity derivatives profiles, under the assumption that all seismic velocity anomalies derive from lateral temperature anomalies (as done by for example Steinberger and Calderwood, 2006, Alisic et al., 2012, Zhu et al., 2013, and Adam et al., 2014. Before conversion, the velocity anomalies of each tomography model are interpolated on a grid of 0.5 • x 0.5 • x 10 km resolution. The conversion of velocity to temperature anomaly for each data point is performed with the depth-dependent anharmonic and anelastic profiles of Karato (2008, Fig. 26g&h) according to: where V is the seismic velocity and B represents either the shear or the compressional wave type (Karato, 2008). We do not add δT (r) to the upper 170 km, where we defined our synthetic lithosphere temperature distribution, nor to the lower 200 km of the domain (transitions are smoothed with a hyperbolic tangent with a 25 km halfwidth), as especially in the crust and the D" layer above the CMB, variations in seismic velocities may also incorporate a significant compositional contribution (lithospheric compositional differences are incorporated through our synthetic plate system; Section 3.2.4). We also refrain from adding temperature anomalies to the background temperature profile in the slab area to cut out contributions from the tomographic slab. See Fig. 3e for a representative slice through the initial temperature field.
In order to assess the effects on surface velocities from different mantle buoyancy fields, we experimented with temperature fields derived from 9 tomography models, which differ in wave type and amplitude and smoothness of the imaged mantle structure (see Tab. 7). The models were obtained from their respective online sources and can be inspected at the https://www.earth.ox.ac.uk/~smachine/cgi/index.php (Hosseini et al., 2018) and http://www.atlas-of-the-underworld.org (van der Meer et al., 2018) websites.

Material parameters.
Temperature T -and pressure P -dependent material parameters thermal expansivity α, density ρ and specific heat c P are read in by from (P ,T ) tables. Compressibility is calculated from the tabulated density as 1 ρ ∂ρ ∂P . We computed tables for oceanic and continental crust, pyrolite lithospheric mantle and pyrolite mantle. The latter two (P ,T ) tables were constructed with the Gibbs energy gridded minimalization option of the thermodynamics software Perple X (Connolly, 2009), using the databases of Holland and Powell (1998) and Lithgow-Bertelloni (2005, 2011), respectively. Material constituents were chosen according to the pyrolite mantle composition of Workman and Hart (2005).
The pressure and temperature conditions for which our computations would require tabled material properties of the crust (especially <500 K and <3 kbar) are usually not in thermodynamic equilibrium and we found that Perple X calculations produced results not applicable to our modeling. Therefore we constructed crustal (P ,T ) tables ourselves (see App. B). Typical profiles of α, c P and ρ are given in Fig. 26, while a representative crosssection of the density can be found in Fig. 3f. Thermal conductivity is assumed constant at 4.7 Wm −1 K −1 (e.g. Jacobs and Van den Berg, 2011).
The density contrast with the surrounding unperturbed mantle in the synthetic Aegean slab resulting from the prescribed temperature distribution and the used (P ,T ) table lies predominantly between 40 and 80 kgm −3 , conform Cloos (1993), with a peak of 120 kgm −3 at 400 km depth. The synthetic slab density is higher than the density anomalies generated by the other mantle temperature heterogeneities, which is warranted by the findings of Lu and Grand (2016), who showed that maximally 88% of the subducting plate mass is recovered in tomographic inversion.
3.2.5. Rheology of the crust and mantle. We implement a viscoplastic rheology, including diffusion and dislocation creep, Peierls creep and Drucker-Prager plasticity for the plates and upper mantle (Garel et al., 2014, Glerum et al., 2018. Lower mantle deformation is restricted to diffusion creep only (Čížková et al., 2012), and we smooth the transition from a composite rheology to diffusion creep only with a hyperbolic tangent with a half-width of 40 km around the depth of the 660 km phase transition.
Diffusion, dislocation and Peierls creep can be conveniently formulated with one equation (Karato and Wu, 1993, Karato, 2008, Garel et al., 2014: where in case of diffusion creep, n = 1, while for dislocation creep n > 1 and for Peierls creep n = 20. The effective deviatoric strain rate is defined as˙ e = 1 2˙ ij˙ ij . See Tab. 2 for the definition of other symbols. The effective plastic viscosity is given by (Glerum et al., 2018): with σ y is the full pressure-dependent yield stress. All four deformation regimes act simultaneously (Karato, 2008) under the same deviatoric stress, so their contributions to the effective viscosity are combined as (van den Berg et al., 1993, e.g.): A typical viscosity profile resulting from the parameter values in Tab. 2 is given in Fig. 26f. Plate boundaries, which we mimic with 30 km wide low viscosity zones, are assigned a fixed viscosity (default value = 1·10 21 Pas) that can be varied laterally along the fault trace to tune the coupling between plates. The old oceanic crust in front of the Aegean subduction zone also has a fixed viscosity of 1·10 21 Pas.
We note that our implementation of crust/lithosphere rheology as laterally homogeneous (apart from thinned regions and plate boundaries) and rather stiff (with a strong core of ∼3·10 24 Pas, see Fig. 26f), i.e. of realistic strength (e.g. Burov, 2011), and importantly without internal weaknesses such as lithosphere-scale faults, avoids intra-plate forcing of surface deformation but at the expense of creating a low-pass filter effect on the surface motion response to sublithospheric processes, as was also observed in the analogue models of Sembroni et al. (2017). We will refer to this as the low-pass filter effect of the lithosphere.
3.2.6. Plate motions: boundary conditions. Plate motions in our new coordinate system are prescribed on the lateral domain boundaries down to a depth of 130 km (i.e. along the whole depth of the reference lithosphere) through bilinear interpolation between tabulated values. Plate motions were precomputed at the tabulated points with a resolution of 0.03-0.1 • x 5 km from relative or absolute Euler poles Ω as v = Ω × r. Instantaneous plate motion poles relative to Nubia (Africa) computed from GNSS measurements are taken from Kreemer et al. (2014). Absolute plate motions are constructed by adding to these poles the instantaneous Euler pole of Nubia describing plate rotation relative to the deep lower mantle, for which we again adopt the average Euler pole of Africa during the last 10 My in the GMHRF of Doubrovine et al. (2012), see Tab. 3.
Along the model domain edges, transitions in plate rotation over a plate boundary are smoothed with a hyperbolic tangent with a halfwidth of 0.27 • latitude. Below the prescribed lithosphere velocity, the lateral boundaries are free-slip. While the CMB velocity boundary condition is no-slip, the surface is free-slip to allow for lateral deformation of the crust. Composition and temperature are fixed along the top and bottom boundaries.
3.3. Model analysis. Through the instantaneous flow modeling with ASPECT, we obtain predictions of the flow field and the strain rate. This output is first compared to GNSS observational data (Kreemer et al., 2014) of the whole eastern Mediterranean (Fig. 2) and the Aegean region specifically (Fig. 4) in order to establish a reference model with an overall good fit to the GNSS data. This fit is determined with Eqs. 9&10, where the observed velocities are first interpolated on a regular grid of 0.8 • resolution spanning the Aegean-Anatolian region by solving a linear inverse problem as described in Spakman and Nyst (2002, light blue vectors in Fig. 4). These surface motions on the regular grid are compared with our model predictions computed by evaluating the finite element solution at the same   Kreemer et al. (2014) lateral locations at 5 km depth.
where n = 96, the number of points of comparison. Only points on the Aegean-Anatolian plate are included in this analysis.  Kreemer et al. (2014). This data together with synthetic data for the Eurasian and Nubian plate constructed with the Euler poles of Kreemer et al. (2014) are inverted (Spakman and Nyst, 2002) to obtain the model GNSS predictions in light blue on a regular grid of 0.8 • x 0.8 • . The error ellipses of these predictions are computed from the model covariance.
Our aim is not to produce an instantaneous model that fits the surface motions in some optimal sense. For our purposes, it is sufficient to obtain a model that provides a reasonable data fit such that the model can be accepted as a physical representation of the crust-mantle system under investigation. We then vary the model setup, thus affecting the force balance in the model w.r.t. this reference model, to estimate the change in surface motion, i.e. the sensitivity of surface motion to the particular change in model feature. To quantify this change, we compare the predicted crustal velocity field with that of the reference model rather than with the observed surface motions. This cuts out any observed data misfits from the analysis. We plot the difference in velocity ∆v = v p − v o at the regular grid points in the Aegean-Anatolian regions on top of the velocity magnitude difference field after rotating the reference model velocities to same frame of reference through an experimentdependent, deformation-free uniform Euler rotation. A more quantitative comparison is made by evaluating the velocity gradient field, which is spatially more sensitive than ∆v and defines the crustal deformation in terms of the sum of the strain-and rotation rate. Using the inversion method of Spakman and Nyst (2002), we estimate the velocity gradient field associated with v o and with ∆v from each experiment in the nodes of a surface mesh composed of spherical triangles that covers the A-A region including the diffuse plate boundary regions. The velocity gradient model determined from v o (the Reference model) is shown in Fig. 5 and the corresponding model quality information is provided in Supplementary information Fig. 27. The effective velocity gradient 1 2 ∇ ij v∇ ij v provides the scalar magnitude of the tensor field, where ∇ ij = ∂ i ∂ j is the spatial derivative with summation implied by repeated indices. We adopt as sensitivity measure of deformation change the relative change in this scalar measure of field strength: The term α S is used to suppress S in places where the effective velocity gradient of v o (Fig. 5) is near zero and to suppress effects of amplitude errors in the effective velocity gradient field (Fig. 27). We use α S = 2.5·10 −9 yr −1 ≈ 0.8·10 −16 s −1 , which is generally much larger than the model amplitude errors. We indicate the range of sensitivities obtained for the points included (see Fig. 7) as < S min − S max >.
We also report the velocity of the Nubian plate just in front of the trench (SP vel) and that of the Aegean plate closely behind the trench (OP vel) based on the averages of two point values on each plate.

Experiments
We first build a Reference model that provides a first-order fit to the GNSS data; in effect already establishing the largest sensitivity of surface motions to mantle processes. This reference model then constitutes the basis for comparison with crustal flow predictions from four classes of subsequent experiments (Tab. 4), in which we vary: (1) the prescribed absolute plate motions; (2) the slab geometry in terms of horizontal and vertical tears; (3) the mantle temperature (i.e. density) structure; (4) the sublithospheric mantle rheology.  Fig. 3: Plate velocities are prescribed in the previously defined GMHRF motion frame of Nubia motion of the past 10 My, hereafter called Nubia mantle frame; the synthetic slab is a continuous feature with edges under western Anatolia and northwest Greece; the sublithospheric mantle temperature structure is converted from tomographic model P06 CSloc with the seismic derivative profile of Karato (2008) in Fig. 26g. P-wave model P06 CSloc (Amaru, 2007) is the most detailed seismic model of the eastern Mediterranean out of the set listed in Tab. 7. The model has been determined by Amaru (2007) by inverting the same data set as used for model UU-P07 but relative to a 3D reference model called CUB-S20RTS, in itself a combination of model CUB Ritzwoller, 2002, Ritzwoller et al., 2002) for the upper few hundred kilometers and model S20RTS (Ritsema et al., 2004). This leads to a model comparable to UU-P07 in well-resolved regions, while in poorly-resolved regions the 3D reference model dominates. The rheology is defined by the parameters listed in Tab. 2, including the A-family parameters ofČížková et al. (2012) for the lower mantle.
In summary, the included dynamic drivers of flow are: absolute plate velocities, implicit lateral Gravitational Potential Energy (GPE) gradients due to smooth variations in crust and lithosphere thickness, the negative buoyancy of the synthetic slab and topography of the phase change boundaries, and the mantle buoyancy field. The resisting forces constitute shear stress in the slab related to deformation, the mantle shear stress governed by mantle rheology and the shear stress related to the imposed viscosity of the plate boundary faults.
Starting from this model setup, we tune the plate boundary viscosity around the Aegean-Anatolian plate. We obtain a good fit with the GNSS velocities (ARMS = 17.6 • , MRMS = 5.9 myr −1 ) by lowering the viscosity of the western section of the NAF w.r.t the general plate boundary viscosity of 10 21 Pas to 10 20 Pas. The viscosity of the subduction interface of the Cyprus slab and that of the EAF is increased to 5·10 22 Pas. The Aegean subduction interface remains at 10 21 Pas, a value thatČížková and Bina (2013) found to allow for the penetration of slabs into the lower mantle (like the Aegean slab).
Although a better fit might be found by tuning other parameters, adding more detail to the crustal structure in terms of predefined faults (conform e.g. Nyst and Thatcher, 2004, Floyd et al., 2010, Nijholt et al., 2018, or including more lithospheric heterogeneities, our goal is to find a geodynamically acceptable benchmark that we can use to test the effect of deviations from the reference state on the crustal motion field. 4.1.1. Reference model predictions. Figure 6 presents the computed instantaneous flow field for the reference model setup. The data fit is illustrated by the overall crustal velocity pattern (Fig. 6a) providing a good match with the observed GNSS velocities. Direct comparison of model predictions and observations in the zoom-in of crustal velocities (Fig. 6b) highlights the first-order fit between our model predictions and the GNSS data, particularly the overall rotation and acceleration of the Aegean-Anatolian velocity field towards the trench. The predicted acceleration is however less strong than that observed, as exemplified by the increased misfits around the western continuation of the NAF through the North Aegean Trough (NAT) onto mainland Greece. These misfits could result from the modeled plate boundary fault connecting the continuations of the NAF across the northern Aegean (Le Pichon and Kreemer, 2010) with the lithosphere extension at the Gulf of Corinth, while the actual distribution of lithosphere weaknesses is more complex and accommodated by faults that are not included in our model. The smaller increase in our velocity predictions towards the trench is perhaps also due to the smooth thickness variations we introduced in the Aegean Sea and Anatolia and the oceanic crust south of the Aegean Sea not providing enough weakness and GPE to generate the observed amount of extension. Our misfit in terms of the magnitude of the overriding plate velocity (MRMS = 5.9 mmyr −1 ) is however comparable to the misfit obtained by the thin viscous sheet models of England et al. (2016), where motion is driven only by internal GPE variations and by tuning of boundary forces (5.1-4.7 mmyr −1 ) but not by slab rollback, and those of Ozbakir et al. (2017) of 4.1 mmyr −1 , who focused on detailed fault distributions.
Below the lithosphere at 150 km depth (Fig. 6c), we obtain high velocities underneath the thinned continental lithosphere of Aegea-Anatolia (directed westward and then turning towards the trench) and the Tyrrhenian Sea. As for the strong upwelling underneath the Arabian plate, they are associated with low-velocity bodies in the P06 CSloc tomography model. The decrease in velocity in central Anatolia is in part due to the smoothly reducing the tomography-derived temperature anomalies to zero around the synthetic slab (Section 3.2.3). Decreasing the interval over which the anomalies are reduced does not significantly affect model predictions. Generally, the high variability in the sublithospheric flow field is not recorded by the crustal flow field, which shows only smooth velocity variations (apart from in the plate boundaries). This demonstrates the low-pass filter effect of the (laterally homogeneous) strong lithosphere rheology we employ.
Down to a depth of around 450 km, a small region of upward flow underneath east Anatolia exists, together with down-and westward motion underneath central Anatolia. Underneath western Anatolia, flow is directed SW towards the Aegean slab. In the transition zone (Fig. 6c, right), larger-scale, slower velocity patterns are seen, with material downwelling underneath the whole Aegea-Anatolia region, markedly faster in and around the Aegean slab.
The vertical sections perpendicular to the Aegean trench ( Fig. 6d&6e) show rollback of the slab together with W-SW inflow above the slab. Rollback velocities as recorded by the overriding plate lie around 15 mmyr −1 (OP vel in Tab. 4), lower than what the GNSS data indicates (20 mmyr −1 ; Fig. 4). The slab itself shows deep sinking velocities of up to 60 mmyr −1 . The viscosity of the slab is ∼1·10 22 -5·10 23 Pas, whereas the depth average upper mantle viscosity lies between 6·10 20 and 2·10 21 Pas. Our slab viscosity agrees with estimates from time-dependent and instantaneous numerical modeling as well as fits to observational data ranging from 5·10 21 to 1·10 24 Pas or 100-1000 times upper mantle viscosity (as summarized by Alisic et al., 2012, and references therein). Alisic et al. (2012) also found that slab strengths at the higher end of this range are required to fit plateness and plate motion in global instantaneous flow models.
Sublithospheric mantle velocities above the Aegean slab reach up to 10 cmyr −1 and velocities in some mantle upwellings reach 20 cmyr −1 in the asthenospheric mantle, even though plate velocities are an order of magnitude smaller. Billen (2010, 2012) demonstrated that mantle velocities around the slab could reach values of 80 cmyr −1 , while plate motions did not exceed 10 cmyr −1 in their instantaneous numerical models of the eastern Alaska plate boundary system. In the transition zone, the flow velocity reduces as a result of increasing viscosity, although strong velocities around the Aegean slab are associated with slab rollback. In the lower mantle, flow speeds are generally below 2 cmyr −1 (Fig. 6d).
In summary, we have established a reference model that in absence of detailed lateral density heterogeneity in the crust-lithosphere system and surface topography self-consistently incorporates the dynamic forcing of the lithosphere-mantle system including the 'remote' plate tectonic forcing imposed by absolute plate motions and as such captures the overall crustal flow field of the system. 4.2. Mantle reference frames. In our first set of experiments relative to the reference model, we focus on the sensitivity of the crustal flow field to changing the mantle reference frame while conserving the relative plate motions. We vary the mantle reference frame (MRF) between Eurasia-fixed (Eurasian plate velocity is constrained to be zero relative to the mantle on the model sides), Nubia-fixed (prescribed Nubia velocity is zero on the model sides) and an MRF in which Nubia rotates twice or four times as fast, experiments similar to those of Chertova et al. (2014b). As these various plate motion models are obtained by a uniform Euler rotation of the plate motions prescribed for the MRF, the relative plate motions and slab pull stay the same. The mantle drag below the plates however changes and the slab is dragged differently through the mantle, leading to different viscous interactions with the ambient mantle. Slab dragging is the lateral transport of the slab through the mantle that is forced in the direction of the absolute motion of the subducting plate at the trench , here Nubia trench motion in the MRF. If the impact of viscous mantle resistance is negligible, we would not expect a surface motion field different from that of the reference model rotated to the new MRF.
We compare the resulting crustal flow fields of the Aegean-Anatolian region with that of the Reference model. Slices of the velocity field at 5 km depth in the models' own MRF as well as the velocity difference with the reference model in this frame are given in Figs. 8-11, together with cross-sections perpendicular to the Aegean trench in the individual MRFs. Figure 7 shows the sensitivity of each model relative to the Reference model motions (after being rotated into the new MRF), while Table 4 lists the Subducting Plate (SP) and Overriding Aegean-Anatolian Plate (OP) motions in the respective model frames.
4.2.1. Eurasia fixed. Crustal motions predicted for the Eurasia-fixed MRF are shown in Fig. 8a. Velocity differences with the rotated Reference model (Fig. 8b) show a systematic pattern of W to SW motion with amplitudes of up to 0.8 mmyr −1 in Arabia and the Aegean-Anatolian plate. Although the motion sensitivities w.r.t. the reference model are accordingly small (Fig. 7), this is ∼16% of the amplitude of the NW directed absolute motion of Nubia in this MRF. This sensitivity to the Eurasia-fixed MRF has several causes. First, in the Eurasia-fixed MRF there is no NE directed stress transfer -a pull by Eurasiaacross the NAF as opposed to in the Nubia MRF of the Reference model, allowing Aegea-Anatolia to flow more parallel to this zone of weakness. Second, in the current MRF, the absolute Nubian plate motion is perpendicular to the SSW-directed rollback and more or less parallel to much of the trench pushing the slab to the NNW ( Fig. 8a and 1). Consequently, the SSW rollback due to slab pull is less counteracted by Nubian plate advance to the trench, leading to stronger rollback expressed by an increased southward motion (c) Predicted velocity at 150 km and 600 km depth. As the velocity field is semi-transparent, the 3D velocity vectors pointing into the view plane are lighter colored than the upward vectors.
(d) Vertical slice through Aegean slab of the velocity field with surrounding mantle velocity vectors and temperature contours. Location of the cross-section is indicated in Fig. 6a. The velocity vectors are interpolated on a regular grid.
(e) Vertical slice through Aegean slab of the viscosity and slab velocity vectors. Location of the cross-section is indicated in Fig. 6a. The velocity vectors are interpolated on a regular grid.   (Fig. 5). This measure is based on determining the effective velocity gradient field of the the velocity difference field w.r.t the reference model obtained from each model experiment as is displayed in the upper right panels of the experiment result figures that follow, e.g. Fig. 8. See Section 3.3 for a description of the procedure. component of the overriding plate, while Aegean motion also obtains a larger westward component due to NNW directed trench-parallel slab dragging. Increased rollback is also exemplified by faster flow toward the slab in the mantle wedge compared to the Reference model (Fig. 8c vs. Fig 6d). The differential motion of the Arabian plate w.r.t. the reference model stems from the strong plate contact between Anatolia and Arabia, which transfers the lateral crustal stress due to increased slab pull acting on the Aegean plate, while the weaker Dead Sea fault decouples the Nubian and Arabian plate. 4.2.2. Nubia fixed. Figure 9 shows the motion results obtained for the Nubia fixed model. Differences with the reference velocity field occur in the Nubian and Aegean-Anatolian plate on the order of 1 mmyr −1 (Fig. 9b), which is ∼20% of the SE-directed motion of Eurasia in this MRF. The sensitivity measures are nevertheless small (max. 8 %, Fig. 7).
In the Nubia-fixed MRF, the subducting plate velocity near the trench is practically zero (see top part of the slab in Fig. 9d and the SP velocity in Tab. 4) and slab motion is only driven by the pull of the Aegean slab. Slab rollback is maximum in this setting and results in some internal deformation of the Nubian plate close to the trench (Fig. 9a). Compared to the Eurasia fixed model, slab rollback as recorded by the overriding plate has turned southward. Evidently, the larger W-component of crustal motion in the Aegean predicted in the Eurasian-fixed MRF is caused there by the NW directed absolute motion of Nubia. In the Nubia-fixed MRF there is no slab dragging, allowing for slab rollback to dominate all southern Aegean motion. In addition, in this MRF there is a SE-directed push across the northern plate boundary zone resulting from the SE motion of Eurasia (Fig. 9a). This adds a SE component to the crustal flow field south of the plate boundary compared to that of the Eurasia-fixed MRF. Also, below the Aegean crust, mantle flow to the SSW has increased as a result of rollback.
Comparison between all three models indicates that velocity patterns in the sublithospheric mantle away from the slab are similar in the different MRFs. In particular, we note a broad zone of downward flow of ≤20 mmyr −1 in the lower mantle below the Aegean slab that results from the lower mantle continuation of the slab (e.g. Spakman et al., 1993, van der Meer et al., 2018. 4.2.3. Nubia X 2. Even though the relative motions between Nubia and Eurasia and Nubia and Arabia remain as in the reference model, increasing the absolute motion of Nubia by a factor 2 ( Fig. 10) affects the coupling of the plates across the plate boundaries and the interaction with the underlying mantle. To retain the relative motion between e.g. Eurasia and Nubia, Eurasia motion also changes amplitude and assumes a more northerly direction. Figure 10a shows that in this MRF the Aegean region is almost mantle-stationary, showing little to no motion (compare OP motion in Tab. 4). Differences w.r.t. the reference model after rigid Euler rotation into the Nubia X 2 MRF are also most notable in this region and in the plate boundary zone to the south, reaching up to 0.8 mmyr −1 .
In the Nubia X 2 model, rollback of the slab is effectively cancelled by the advancing motion of Nubia dragging the slab northward, as is visible from the subvertical slab vectors in Fig. 10d. Slab pull now leads to a near vertical sinking of the slab compensated by smaller mantle material inflow from the north (Fig. 10c). The smaller shear tractions of (a) Velocity field at 5 km depth in model MRF of motion. In green the current model predictions, in brown the Reference model predictions.  this flow acting on the thin lithosphere atop may also contribute to the southern Aegean region being mantle-stationary. Now that slab rollback is absent, the surface motions in Anatolia are mostly driven by Nubian plate push, while the Aegean region is in the 'shadow' of the laterally stagnant slab. Northward advance of the slab due to the Nubia motion is apparently resisted by the mantle, as Nubia crustal motion south of the slab is smaller than in the Reference model (i.e. the difference flow field in Fig. 10b is directed southward).
4.2.4. Nubia X 4. In the Nubia X 4 MRF, Nubia motion is ∼50 mmyr −1 (Tab. 4). The overall flow field is associated with a slab rolling forward and with northward trench advance, as can be seen in Fig. 11d, where the upper part of the slab is advancing northward through the mantle, while the lower part of the slab is sinking vertically. The slab tip advances along the 660 km viscosity jump. Slab sinking still excites a trenchward flow below the south-Aegean lithosphere, although less than in the previous models and more westward (Fig 11c). Where plate motions are similar in direction to the sublithospheric mantle flow, mantle velocities are higher compared to the reference model due to shear Figure 10. Nubia X 2 model results. See caption of Fig. 8 for further explanation.
coupling with the lithosphere. Mantle resistance against NNE slab dragging causes deformation in the Nubian plate as exemplified by the outward radiating differential flow field to the south of the trench. This pattern is probably related to lateral deformation of the slab, evidenced by lower slab viscosity, due to slab-mantle coupling during its northward transport. Differences in crustal velocity compared to the Reference model are now considerably larger (Fig. 11): differential motions up to 3.0 mmyr −1 occur in the Aegean region, while differences of about 1 mmyr −1 and 2 mmyr −1 are seen in the Arabian and Nubian plates, respectively. The sensitivity measures are the largest so far, <∼3-22>%, with one outlier of 35% (Fig. 7). 4.2.5. Discussion -Absolute plate motions and the coupling between lithosphere, slab, and mantle. Figures 8-11 demonstrate that the same relative plate motions occurring in distinctly different MRFs produce a crustal response that deviates from the Reference model Figure 11. Nubia X 4 model results. See caption of Fig. 8 for further explanation.
in both horizontal motions and the velocity gradient field due to changes in plate, platemantle and slab-mantle coupling. In the slow-plate motion MRFs of Eurasia fixed and Nubia fixed, changes in plate motions lead to differences in internal plate deformation of both the overriding and the subducting plate of up to 1 mmyr −1 relative to the motions in the Reference model due to a.o. varying interactions of slab rollback and subducting plate advance (S =<∼0-8>%). In the mantle MRFs with faster plate motions, the differences increase to several millimeters per year and S =<0-35>%. The differential motions (Figs. 8b-11b) prove to be regionally systematic and significant, but are intriguingly small in amplitude when compared to cmyr −1 -scale OP motions driven by the regional coupling between the Aegean slab and the mantle (see for example the dramatic kinematic difference between the two MRFs of Fig. 10a&10c and Fig. 6d versus the difference field of crustal motions in Fig. 10b). Generally, much deformation is taken up by the plate boundary zones surrounding the A-A plate while the internal deformation is smooth and of low amplitude although increasing with absolute plate motion. The smoothness of differential motion patterns we attribute to the low-pass filter effect caused by the strong unbroken lithosphere we implemented (Section 3.2). The first conclusion we draw from this set of experiments is that observations of crustal motions, when interpreted in the proper MRF, can lead to the identification of the style of subduction. This is shown in our experiments by the large sensitivity of the absolute crustal flow field (Figs. 8a-11a) for deep processes ranging from pure slab rollback, when slab pull is the only forcing, to various interactions between slab pull and slab dragging. In the models Reference, Nubia X 2 and Nubia X 4, the Nubia plate was given increasing subducting-plate advance velocities. The predicted trench motion (retreat, stationary, advancing, respectively) agrees with the findings of Schellart (2005), who investigated the relation between subducting plate velocity and trench migration with analogue modeling.
The second conclusion from these experiments is that both mantle and surface deformation are sensitive to changes in the dynamic balance between slab dragging by the subducting plate advance velocity, slab pull, and the viscous coupling with the ambient mantle. Interaction of the plates with the mantle, even away from the slab, leads to changes in mantle flow as indicated by the absolute strain rate difference field at 200 km depth in Fig. 12. As strain rate is invariant for uniform Euler rotations of the MRF, a direct comparison between all our models is warranted and shows strong, detailed differences, particularly in and around the slab. These differences may also be expressed in the mantle anisotropy. Within the plate boundary that surrounds the A-A plate, geologically significant strain rate differences on the order of 10 −15 s −1 occur. Strain rate within the A-A plate is generally at least one order of magnitude lower while most of the differential motion pattern shows regional rotation of the relatively stiff A-A plate being accommodated by the enclosing weak plate boundary. The sensitivity estimates of the regional-scale velocity gradient vary between <∼0-34>% of the Reference model regional velocity gradient amplitude (Fig. 7) and indicate modest but significant impact of absolute plate motion on crustal deformation as a result of predominantly slab-mantle interaction.
Our third conclusion concerns the more common case of crustal motions being analyzed in a relative plate motion frame for which our experiments show that the crustal motion field and its internal deformation cannot unambiguously distinguish between slab rollback, a stationary slab, or an advancing slab despite the significantly different surface responses. Our results underline that interpretations of crustal motions in relative plate motion frames are only useful for identifying the nature of crustal deformation, e.g. the sense of fault motion and style of distributed deformation, but cannot identify its deep causes such as slab rollback or advance. For instance, rotating the crustal velocity field of Nubia X 2 into the reference motion frame would lead to a motion pattern comparable in style to that of the reference data set, suggesting slab rollback instead of stationary subduction. Equally, a MRF should be used for interpretation of mantle anisotropy caused by the accumulation of strain, i.e. through the integration of strain rate, in time-dependent mantle flow models (e.g. Faccenda, 2014, Confal et al., 2018. 4.3. Slab geometry. Our reference model considers a continuous slab. In the next set of experiments, we test the crustal motion response to partial and complete detachment, Figure 12. Strain rate difference field at 5 and 200 km depth for experiment set 1. The strain rate difference is computed as the absolute difference between the considered experiment and model Reference. Crustal strain rate differences are predominantly located in the plate boundaries, the overriding plate and the areas in front of the trenches. For model Nubia X 4 the Nubian and Arabian plates also display internal differences. In the mantle, differences are largest in and around the Aegean and Cyprus slabs. complete absence of the slab, and vertical tearing of the slab (see Tab. 4). Although these experiments are not intended to precisely simulate any existing slab geometry for the Aegean region, they are motivated by proposed slab tearing processes in the Mediterranean (e.g. Spakman et al., 1988, Wortel and Spakman, 2000, Spakman and Wortel, 2004, and Section 2). We focus on deep detachment (below 200 km), which was predicted to occur in old lithosphere like the Mediterranean basin by Duretz et al. (2011) and van Hunen and Allen (2011).

Slab detached.
For the Slab detached model, all synthetic Aegean slab material between 200 and 300 km depth is removed. Figure 13a shows that this complete detachment leads to a decrease of velocity of the overriding plate towards the trench. The same effect can be seen in the velocity difference plot (Fig. 13b; the difference in the Aegean-Anatolian region is at most 3 mmyr −1 ; see also Tab. 4). Differences in the velocity do not extend into the Eurasian plate, but do show a stronger response (∼6 mmyr −1 ) in the western Hellenic subduction zone. The plate coupling is thus affected strongest. Figure 13c shows that trenchward inflow above the slab is reduced compared to the Reference model. Sinking of the slab remnant is similar in direction and magnitude however; even though sinking is no longer resisted by the connection of the slab to the Nubian plate, it is hindered by the viscosity increase around 660 km depth. Downward motion of the still attached part of the slab is reduced due to the missing slab pull of the detached slab. Although a strong south-directed motion is still induced directly under the Aegean plate (Fig. 13c&17), the reduced slab pull probably causes the reduced rotational motion of the overriding plate.

Slab partially detached.
When the slab is not completely horizontally detached, but only in the west along the trench up to the island of Crete, the difference in surface response with the Reference model is smaller (Fig. 13f). Model Slab partially detached predicts a smaller change in velocity magnitude than model Slab detached and a distinct change in direction. Aegean motion is more southward than in the reference model, directed towards the part of the slab that is still attached to the east of Crete.
Mantle flow into the subduction wedge is still smaller than in the Reference model (Fig. 17). The motion of the upper 200 km of the slab (not shown) has a slightly larger vertical component than in the completely detached model, but it is still smaller than for the reference setup.

4.3.3.
Slab detached and removed. The relatively small differences between the reference model and the previous detachment models (up to ∼50%, but for the bulk of the evaluated points smaller than 21%, Fig. 7) prompt the question whether the flow induced by the detached slab remnant is an important contributor to the overriding plate motion or that the slab pull generated by the upper 200 km of slab connected to the Nubian plate is dominant, or both.
Removal of the detached synthetic slab segment reduces the trench-normal velocity of the overriding plate by 1-2 mmyr −1 (Fig. 14b and Tab. 4). The vertical model slices show that the inflow above the slab is greatly reduced and the attached part of the slab sinks vertically into the mantle (as in model Slab detached, Fig. 14). Nevertheless, the sensitivity measures are comparable to that of model Slab detached with the slab remnant left in: sensitivity S = <∼1-51>%. Hence, in our models Aegean motion is primarily controlled by the top 200 km of the subducting plate, whereas the sinking of the deeper part of the slab excites an asthenospheric flow toward the slab that has only a secondary effect on crustal motion.

No slab.
When the Aegean slab is completely removed below the depth of the lithosphere (120 km), overriding plate motion towards the trench is reduced by about a factor three w.r.t. the Reference model (Fig. 15). Differences with the reference model are up to 10 mmyr −1 , two-third of the reference OP motion. In the western and eastern margins of the subduction zone, differences are largest, even exceeding 10 mmyr −1 . The difference field in effect outlines the region that is subject to slab pull, mainly the OP and the trench. The substantial response to the absence of the slab is also reflected in the sensitivity measures: S = <∼5-91>% (Fig. 7). The vertical velocity slices show no fast flow underneath Aegea due to the lack of a slab and its rollback.
Absence of the slab leaves the overall crustal flow field of the Aegean-Anatolia region subject to the relative plate motions between the three major plates (by stress transmitted across the plate boundary fault zones) and to tractions at the base of the lithosphere. The relative plate motion between Nubia and Eurasia is NW-SE convergent at ∼5 mmyr −1 and Arabia provides NNW motion. Under eastern Anatolia, WSW directed tractions again occur (Fig. 17), but underneath the Aegean and western Anatolia, mantle motion is now NNW, with even northward flow underneath mainland Greece. Together these lead to a strongly decelerating westward motion of the Aegean-Anatolian plate. 4.3.5. Vertical tear Kefalonia. As a last experiment in this set, we include a vertical slab tear at the Kefalonia fault zone. As slab pull is now more concentrated on the southeastern part of the slab (a similar, but smaller effect was seen for model Slab partially detached ), this results in increased (OP velocity = 17 mmyr −1 ) and more southward motion of Aegea (Fig. 16). Largest differences with the reference model are found right above the slab tear.
The southwestward inflow of mantle material above the slab and under the Aegean Sea is of smaller magnitude, as can be seen in the comparison of mantle flow of this set of experiments in Fig. 17.
(a) (b) Figure 16. Vertical tear Kefalonia model results. See caption of Fig. 8 for further explanation.
4.3.6. Discussion. Negative slab buoyancy in combination with subducting plate motion determines the rollback of a slab, which in turn induces toroidal flow around and poloidal flow from underneath the slab (e.g. Funiciello et al., 2003, Piromallo et al., 2006, Schellart and Moresi, 2013, Chen et al., 2016, and references therein). Although of reduced magnitude (Fig. 17), a slab remnant from complete horizontal tearing is still able to generate toroidal flow (Fig. 13c). Removal of the slab remnant greatly diminishes slab-induced mantle flow, but has a mmyr −1 -scale effect on the crustal deformation. In fact, in our models the slab pull and rollback of the top 200 km of the slab exert the most control on the surface deformation of the overriding plate (about half of Aegean motion). The area most strongly affected by the pull of the slab is the trench and the Aegean part of the overriding plate (up to 60% of the crustal motion, but Anatolia and Adria experience deformation due to the slab as well. Despite the simplicity of our geometric implementation of slab tears (see the studies of Duretz et al. (2012Duretz et al. ( , 2014 for 2D and 3D time-dependent modeling of the process involving viscous necking and sometimes simple shearing), we can conclude that the changes in horizontal motions associated with complete deep slab detachment (>200 km depth) are around 1.5 mmyr −1 relative to 15 mmyr −1 rollback in the Reference model (maximal differences are ∼3 mmyr −1 ). From comparison of models Slab detached and removed and No slab, we infer that shallower detachment will have a much greater signal in the crustal flow field. The 200 km-depth slices in Fig. 17 show that differences in slab pull are sensed by the sublithospheric mantle in a wide radius around the Aegean slab, although most strongly in the region directly above the slab. We cannot easily distinguish between the relative controls of traction from this slab-induced asthenospheric flow and of trench suction on the overriding plate. Numerical models of oceanic subduction indicated an active role of the former in thermally modulating the lithosphere and shearing its base Moresi, 2013, Sternai et al., 2014). Analogue models of Chen et al. (2016) even pointed to an order-of-magnitude larger control of mantle flow than trench suction, as backarc extension coincided with a locally large trench-normal gradient in sublithospheric mantle velocity. In our models, the asthenospheric flow is governed by the sinking of the deeper (below ∼200 km depth) part of the slab. This is revealed first by the Slab detached model, where the slab is entirely detached, but a large part of the sublithospheric flow towards the slab under the Aegean is still generated, and second by model Slab detached and removed that with only a 200-km slab shows no strong sublithospheric flow, but where the rollback of this short slab explains a large component of the slab rollback motion. Hence, we find that even though the SW-directed sublithospheric flow in the Reference model is faster than the SW motion of Aegea, the direct effect of slab pull by far dominates over mantle tractions. We note that we observe this in a laterally highly curved subduction system, which differs geometrically from the models in which toroidal flow appeared to be a dominant driver of backarc basin motion. Moreover, different coupling of the plates over the subduction interface as defined by the subduction zone viscosity could also explain differences in the transmission of slab motion to the overriding plate.
When no slab is present (model No slab), rollback cannot contribute to the deformation of the overriding plate. The crustal flow field of the Aegean-Anatolian region in this case is completely determined by lateral plate interactions and plate-mantle coupling. This results in near zero motion relative to the mantle for the Aegean region. In contrast, extension of the Aegean region is increased as compared to the reference model when the slab is torn vertically, which places the Aegean slab edge underneath Kefalonia. This allows for more southward rollback, as trench retreat is no longer hindered by the connection to the Eurasian plate.
We conclude from these experiments that changes in slab geometry due to tearing processes have a distinct horizontal surface motion signature on the scale of several mmyr −1 , compared to the reference model. This 'detachment signal', limited to the overriding Aegean-Anatolian plate and its subduction zones, is significantly larger than the precision of GNSS-derived surface motions (typically 0.1-0.5 mmyr −1 in the motion field of Kreemer et al. (2014)). In case of deep slab detachment, this signature is is however small compared to the continuing rollback. In terms of strain-and rotation rates, as measured by the effective velocity gradient, the impact of detachment is very high, with changes of up to 92% of reference model values. We also showed that if no slab is present, the southward motion of the Aegean region is strongly diminished and the overall motion of Aegean-Anatolian plate is subject to the relative plate motions and the sublithospheric flow under eastern Anatolia. This indicates that the overall geodetic motion field of the region is strongly sensitive to rollback of the Aegean slab.

Mantle heterogeneity.
It is well known that seismic tomography models of the mantle differ (e.g. Boschi, 2002, Schaeffer andLebedev, 2013), which has led to studies trying to quantify what tomographic models have in common (e.g. Lekic et al., 2012, Shepard et al., 2017. Here we are not concerned with why these models differ, but we rather exploit the fact that they do to generate various data-based temperature and density models of the mantle. This strategy was followed earlier by e.g. Boschi et al. (2010), who noted a dependence of surface deformation on the tomography model used in simplified models of mantle flow. We test this dependence in a third set of experiments by 1) removing any tomographic sublithospheric density variations from the whole mantle (model No Tomo) or from the lower mantle (model No LM tomo); 2) creating a spread in the mantle buoyancy field by converting 8 different tomographic models to temperature anomalies (Tab. 4); and 3) converting P06 CSloc with different end member scaling profiles. Slices of the predicted velocity fields at 200 and 5 km depth are presented together in Fig. 18 and 19.

No tomo and No LM tomo.
Removing the density perturbations derived from mantle tomography models (model No tomo) isolates the contribution of slab pull to the regional sublithospheric (Fig. 18) and crustal (Fig. 19) flow pattern, apart from that induced by lithospheric thickness variations. This slab pull results in inflow above the slab of up to 70 mmyr −1 , increasing towards the slab (similar to Billen 2012, Alisic et al. 2012). The inflow is fed by toroidal flow around the slab's edges, strongest around the eastern edge through the window between the Aegean and Cyprus slab. Additional flow derives via the thinner Anatolian region from below the Arabian plate, although due to the absence of the high temperature body underneath eastern Anatolia, this flow is reduced and more westward as compared to the Reference model. Despite significant slab-pull induced flow, the rollback of the slab and resultant inflow above it are about 25% smaller than in the Reference model. This could perhaps be attributed to the absence of the positive density anomaly in the lower mantle associated with the deep extension of the Aegean slab (Spakman et al., 1993, van Hinsbergen et al., 2005Fig. 20), which may exert a downward suction on the upper mantle slab, although of small amplitude, as has been proposed by Faccenna et al. (2013) (compare bottom row of Fig. 20 and Fig. 6). However, model No LM tomo (Fig. 19) demonstrates that the lower mantle contribution to the crustal flow field is small (S = <∼0-4>%, compared to <∼5-23>% for No tomo; Fig. 7) when upper mantle heterogeneity is included.
The crustal flow field of No tomo shows a change in the direction of motion in Anatolia, while in the Aegean region, the magnitude of velocity is most affected (about 4 mmyr −1 smaller, Fig. 20). The more northward motion of eastern Anatolia could correspond to the smaller and mostly westward flow of the mantle underneath. The slower motion of Aegea most likely results from the reduced slab rollback now that the density-induced mantle flow is not enhancing SW slab motion. Alternatively, the reduced mantle drag of Aegea impedes trench retreat through coupling of the OP and the slab over the plate contact. Because of the discrepancy between crustal and sublithospheric mantle motion directions, the former is more likely.  Interestingly, differences in plate velocity are predominantly seen in the Aegean-Anatolian region. Mantle traction on the bottom of the major plates from the P06 CSloc model does therefore not seem to alter the deformation of these strong entities.

4.4.2.
Mantle structure from different tomographic models. While individual models are discussed in App. D, we can draw several general conclusions from the mantle and crustal flow fields predicted by the 8 tomographic models in Figs. 18 and 19, respectively. Except for model SL2013 S40RTS, the experiments lead to crustal motion difference fields that show a smooth clockwise rotation of the Aegean-Anatolian region that reduces with distance from the trench. Maximum differences are ∼8-9 mmyr −1 , ∼60% of reference OP motion (Tab. 4). In most models, this differential pattern is characterized by reduced Aegean motion and a more northward motion of Anatolia. From inspection of the mantle flow fields at 200 km depth, these changes result from interaction of the slab with the surrounding mantle flow in two ways: 1) the eastern Anatolian upwelling provides a backward push to the Aegean slab; and 2) NNE flow from underneath Nubia interacts with the slab and reduces the lateral backward motion of the slab. Due to the relative absence of the Anatolian upwelling and the NNE-directed flow below Nubian in the reference model, the Aegean slab sinks more vertically changing the style of rollback as compared to experiments based on several other tomographic models. Coupling over the subduction interface facilitates the slab acting as a stress guide transmitting forcing from the surrounding mantle to the overriding plate. Distinctly increased flow around the slab's eastern edge through the Aegean-Cyprus slab window is seen in a few models (LLNL G3D and GAP P4 ) working against trenchward motion of the OP and model SL2013 S40RTS features strong flow of material from underneath the Black Sea contributing to slab rollback.

4.4.3.
Seismic velocity anomaly conversion -K08 min and K08 max. The conversion of the P06 CSloc tomography model to temperature anomalies for the Reference model is done using the ∂lnV P ∂T (z) profile of Karato (2008). When we take into account the uncertainties given by Karato (Fig. 26g&h) to include the range of profiles found in the literature (models K08 min and K08 max ), the Aegean-Anatolian motions differ by at most 1 mmyr −1 (Fig. 19). Sensitivity measures are therefore very small, with maximum values an order of magnitude smaller than for the different tomographic models (S = <0-4>%; Fig. 7). When using the minimum seismic parameter profile, seismic anomalies are translated into larger temperature anomalies and mantle velocities double. However, flow patterns are not changed, explaining the minor effect on the surface deformation. The difference in surface response is even smaller for the maximum seismic velocity-temperature derivative profile, which reduces mantle velocities, but again maintains the flow pattern (in agreement with Warners-Ruckstuhl et al., 2012).

4.4.4.
Discussion. The response of crustal motion to modifications of the buoyancy-driven mantle flow field can be significant, up to 9 mmyr −1 , with a sensitivity of maximally 36% (Fig. 7). Although this response is most pronounced in the Aegean-Anatolian plate that is surrounded by plate boundary faults, the motion of the Arabian plate is sensitive to mantle flow directions as well through stress transfer over the Aegean-Anatolian and Arabian plate boundary and through coupling with the underlying mantle. The sensitivity of Arabia motion to mantle structure was also demonstrated by Faccenna et al. (2013), who found that Arabia's speed increased when lower mantle temperature anomalies and/or positive upper mantle temperature anomalies were added to their global instantaneous models.
The control of mantle flow on OP motion in our experiments has 2 sources: 1) Basal drag exerted by the mantle on the lithosphere directly translating into changes in overlying crustal motion and 2) interaction of the mantle with the subducting Nubian plate transmitted by the slab to the OP via coupling over the subduction interface. Our models allow us to tentatively distinguish between these sources. In the absence of mantle structure, slab-pull induced flow underneath Aegean-Anatolia reaches 70 mmyr −1 , about half that of when P06 CSloc heterogeneity is included, yet OP motion is reduced much less (≤27%), indicating there is at least not a one-to-one translation of mantle to plate motion (as already established in experiment sets 1 and 2). Moreover, like model No tomo, most experiments with a different mantle structure lack the upwelling material beneath Anatolia that is channeled towards the trench by the thinner OP lithosphere, leaving the supraslab mantle flow to be controlled by slab-pull induced motion. However, these experiments predict larger sensitivity measures, pointing to other flow sources interacting with the slab. The approximately northward flow from underneath Nubia hindering slab rollback seems responsible for reduced Aegean motion. Similar responses to remotely-induced trench-oblique mantle flow were inferred for the Gibraltar slab in the western Mediterranean in time-dependent models by Chertova et al. (2018). This is corroborated by the differently scaled models K08 max and K08 min that although greatly changing the magnitude of the flow, do not alter the flow pattern. The balance of mantle flow tractions interacting with the slab therefore remains similar and leads to much smaller crustal motion differences than the other experiments in this set.
The fit between surface observations of plate motion and predictions from geodynamic models based on different tomography-derived mantle heterogeneities has been investigated with viscous sheet models with prescribed bottom mantle tractions from global flow models (e.g. O'Connell, 2001, Warners-Ruckstuhl et al., 2012), coupled 3D upper/lower mantle flow codes (e.g. Ghosh and Holt, 2012, Steinberger et al., 2010, Osei Tutu et al., 2018a, Wang et al., 2019 and 3D upper and whole mantle models (Lithgow-Bertelloni and Guynn, 2004, Boschi et al., 2010, Becker and Faccenna, 2011, Faccenna et al., 2014, Adam et al., 2014, Steinberger, 2016. In agreement with our results, plate motions were found to be quite sensitive to mantle buoyancy forces, and tomography models with more detailed structure or seismically defined slabs provide a better fit with the observational data (Becker and O'Connell, 2001, Boschi et al., 2010, Warners-Ruckstuhl et al., 2012, Adam et al., 2014. Torque balance studies suggest a 20-50% contribution of horizontal tractions deriving from deep mantle flow to plate motions (e.g. Becker and O'Connell, 2001, Ghosh et al., 2008, Warners-Ruckstuhl et al., 2012. Global 3D computations indicate a first-order control of mantle drag from active upwellings on Arabia/Eurasia convergence Faccenna, 2011, Faccenna et al., 2013). Becker and O'Connell (2001) found that including Lower Mantle (LM) density heterogeneities in their global flow models always improved the fit between predicted and observed plate motions, with 40% of plate driving torques coming from the lower mantle. This seems at odds with the low sensitivity (S = <0-4>%) of the No LM tomo model. However, said study employed radially symmetric viscosity profiles that did not reach the low asthenospheric values that our pressure-, temperature-and strain rate dependent rheology including nonlinear dislocation creep exhibits (1·10 20 vs. 1·10 19 Pas; see Fig. 26f). Low asthenospheric values can decouple LM and lithosphere (e.g. Ghosh et al., 2008), as can a larger viscosity gradient at 660 km depth (Conrad and Lithgow-Bertelloni, 2002). Alisic et al. (2012) moreover demonstrated that depth-and temperature-dependent rheology combined with yielding complicates the relation between plate motions and LM heterogeneity in global instantaneous models. While the inclusion of LM structure speeds up most plates connected to low-temperature LM bodies, some plates slow down when connected to large-scale high-viscosity structures (like the Nazca plate in Fig. 9 of Alisic et al.). 4.5. Sublithospheric mantle rheology. The predictions of mantle flow and its coupling to surface deformation are governed by the nonlinear flow laws for mantle rheology (e.g. Coltice et al., 2019). Because upper and lower mantle rheology is still highly uncertain (Karato, 2010, King, 2016a and nonlinear rheologies introduce numerical complexities, many studies are restricted to radially symmetric viscosity profiles (e.g. Steinberger and Calderwood, 2006, Boschi et al., 2010, Warners-Ruckstuhl et al., 2012. Here we test the influence of sublithospheric mantle viscosity on surface deformation (Tab. 4) without changing the rheological parameters of the subducting plate. Slices of the mantle velocity field at 200 km depth are presented collectively in Fig. 21  4.5.1. Two-layer viscosity. For a laterally constant, two-layer radial mantle viscosity profile of 1·10 21 Pas in the upper mantle and 5·10 22 Pas in the lower mantle (conform Faccenna andBecker, 2010, Boschi et al., 2010), model predictions indicate that the motion of the Aegean and west Anatolian crust is reduced by up to 7 mmyr −1 (Fig. 22). It is also directed more trenchnormal, i.e. less southward, than in the Reference case. The vertical slices show that the upper mantle viscosity of 1·10 21 Pas greatly reduces the magnitude of sublithospheric mantle flow and prohibits backwards motion of the slab (only the deep slab sinks vertically and there is no retreating horizontal velocity component). Flow patterns are strongly altered as well (Fig. 21a). Interestingly, underneath the Mediterranean basin, flow remains northward until it hits the Hellenic and Cyprus subduction zones, while flow in the Reference model is SSE under the Nubian shore and deflected eastward in the Mediterranean Sea by the slab (Fig. 17). The poloidal and toroidal inflow of material above the slab is also reduced w.r.t. model No tomo; the dynamic support of the slab by the surrounding mantle shear stress counteracts the slab's negative buoyancy, reducing rollback and rollback-induced mantle flow, which in effect considerably slows the motion of the overriding plate. These results show a significant sensitivity of crustal flow to mantle rheology.

4.5.2.
Mantle viscosity X 0.5. When the nonlinear sublithospheric mantle viscosity is reduced by a factor 2, mantle flow velocities markedly increase, but mantle flow patterns are not affected (Fig. 21b and 23c). The reduced viscosity of the surrounding mantle provides less support for the Aegean slab and this affects the strength of the subducting plate, which exhibits a lower viscosity despite unchanged rheological parameters. The slab sinking velocities also double, inducing an increase in overriding plate motion of up to 5 mmyr −1 and a stronger change in direction of motion towards the trench (Fig. 23). Increased WSW flow underneath eastern Anatolia possibly contributes to the change in direction of the crustal flow field above. The crustal velocity differences are limited to the OP and the Aegean trench. In the latter, motion towards the OP contact is increased. 4.5.3. Mantle viscosity X 2. Doubling the sublithospheric mantle viscosity has, not surprisingly, an effect opposite to model Mantle viscosity X 0.5. Similar to model Two-layer viscosity, it leads to vertical slab sinking due to a higher resistance of the surrounding mantle to lateral slab motion (Fig. 24c&24d), and inflow above the slab is decreased. As a result, motion of Aegea decreases by up to 4 mmyr −1 . However, contrary to model Twolayer viscosity temperature anomalies do feed back into the flow field through the viscosity formulation, and thus mantle patterns are not changed as much.

Mantle viscosity B.
A smaller viscosity jump around the 660 km phase transition does not affect surface deformation as much as the previous rheological changes (Fig. 25 and Tab. 4), but the magnitude of Aegean-west Anatolian motion does increase by up to 1.5 mmyr −1 . This is likely due to increased slab sinking (as evidenced by the faster motions of the slab and overlying mantle in the vertical slices of Fig. 25) now that the lower mantle provides less resistance. 4.5.5. Discussion. The viscosity of the mantle plays a large role in determining flow speeds in the mantle (Fig. 21). The large effect mantle viscosity has on the flow magnitudes is however not one-to-one translated to the OP velocities; the crustal sensitivity in our set of experiments is generally <3 mmyr −1 , with a maximum absolute response of 7 mmyr −1 . This sensitivity is significant compared to the motions of the major plates. The largest differences in velocity occur in the plate boundary zone between the NAT and the Gulf of Corinth for changes in upper mantle rheology. The velocity in this complex zone above the slab scales with the velocity of the OP. Again a large role is played by the Aegean slab in transferring mantle stresses to the crust. Whether directly through coupling over the subduction interface or indirectly through to changes in rollback-induced asthenospheric flow, the lateral motion of the slab controls the deformation of the OP: when the slab does not rollback, Aegean motion decreases and when rollback increases, so does Aegean motion. Moreover, models Mantle viscosity X 0.5 and Mantle viscosity X 2 demonstrate that the region of the OP that is affected is predominantly the western part, which we can link to the increase, respectively, decrease Figure 23. Mantle viscosity X 0.5 model results. Note the different velocity field and vector scaling in the vertical slices as compared to previous models. See caption of Fig. 8 for further explanation.
of slab rollback. This simple scaling of the composite mantle rheology with a factor 0.5 or 2 changes the subduction regime from fast rollback to near-stationary subduction. This again indicates that for OP deformation, the resistance of the mantle to slab motion is more important than the horizontal tractions generated by mantle flow. Otherwise eastern Anatolia would be more responsive to the strong changes in the velocity of the upwelling underneath (Fig. 21). Interestingly, we noticed an significant difference in crustal response between adopting linear upper mantle rheology or composite, nonlinear, rheology. When the mantle viscosity follows a radially symmetric profile, a change in surface motion outside the Aegean-Anatolian plate is seen as well. Both Adria and the Arabian plate show differential motion w.r.t. the Reference model. Here, vertical and horizontal mantle tractions could be the cause, as mantle flow is reduced when the temperature in the hot upwellings underneath the heel of Italy and Arabia no longer affects the viscosity of the mantle. Jadamec and Billen (2012) compared Newtonian and composite mantle (and slab) rheologies for their 3D instantaneous models of flow around the Alaskan slab. As we do, they found that strain rate-dependent mantle rheology leads to much larger mantle-wedge velocities. Similar to Fig. 20d, low-viscosity ellipses are generated in Jadamec and Billen's models directly above the subducting plate and underneath the overriding lithosphere, as well as a lower-viscosity region around the entire slab. These zones decrease the coupling between mantle and overlying lithosphere. In our models, increased coupling due to a constant upper mantle viscosity does however not counteract the effect of reduced slab pull on the overriding plate.
Lowering the viscosity of the top part of the lower mantle (model B-family) increases instantaneous slab rollback. This agrees withČížková and Bina (2013), who showed with time-dependent models that although lower mantle viscosity is not the primary control on trench retreat, higher values do impede slab rollback. An increased viscosity jump also reduces slab sinking velocities and trench retreat in the free-plate subduction models of Garel et al. (2014). Figure 25. B-family model results. See caption of Fig. 8 for further explanation. Note the different scaling of the bottom left velocity field as compared to previous models.
We note that the crustal response obtained in the experiments may underestimate the actual surface response to deep processes, because we used a two-layer approximation of the lithosphere. This approximation on average retains a strong core of lithosphere strength on the order of 1·10 23 -1·10 24 Pas (see for example the radial viscosity profiles in Fig.  26). The vertical heterogeneity of lithosphere viscosity determines the relation between surface deformation and deformation of the lithosphere at depth (Flesch and Bendick, 2012, Bendick and Flesch, 2013, Sembroni et al., 2017. In our case, they will be more strongly coupled, as we did not implement a lower crustal weak layer between the shallow crust and lithospheric mantle (Burov, 2011), as its existence is unverified for the region. Also, our A-A lithosphere does not contain internal weaknesses on which deformation can localize other than smooth thickness variations. Such localization may affect the amplitude and sense of fault motion, crustal block rotations or regional crustal shear, shortening or extension, even on the scale of mmyr −1 (e.g. Spakman et al., 2018). A simple experiment where we uniformly decrease the nonlinear crustal viscosity by a factor 10 shows that the velocity of the A-A region increases to 3 cmyr −1 in the Aegean, exceeding GNSS velocities in the absolute frame by ∼1 cmyr −1 , and that the Aegean region moves more westward than observed.
Similarly, a different rheological parameterization of the subducting plate could also affect model predictions as, on the one hand, weaker plates facilitate the bending required for subduction, while on the other hand, stronger slabs better transmit the slab pull that drives plate motion (Elsasser, 1969, Conrad and Lithgow-Bertelloni, 2002, Stegman et al., 2006, Wu et al., 2008, Garel et al., 2014, Chertova et al., 2014.
Lastly, we approximated plate boundaries as ∼30 km wide zones of reduced viscosity relative to the surrounding lithosphere. We exploit a natural trade-off between fault width and strength assuming their product is some measure of integrated plate boundary strength. As the plate boundary zones accommodate stress and motion transfer across the faults and our implementation suppresses strong local changes, this also contributes to the overall spatially smooth crustal motion responses we observe. We have kept the lateral plate boundary viscosity distribution that we obtained for the reference model constant throughout subsequent model runs, thereby fixing the coupling between the interacting plates. As our experimental focus is on slab-mantle coupling, plate-mantle coupling and mantle flow forcing of surface motion, this ensures that the response to changes in the dynamics of the sublithospheric system can be adequately investigated.

Discussion
We have presented four sets of regional full-mantle instantaneous-dynamics models centered around the Aegean-Anatolian (A-A) plate and the Aegean subduction zone that indicate that the observed crustal deformation, i.e., the velocity field and its associated velocity gradient field, carries information on both shallow and deep geodynamic processes. Discussions of each set of experiments separately are presented in Section 4 and the relative sensitivity obtained from each single experiment is summarized in Fig. 7; here we first discuss other modeling approaches towards explaining crustal deformation of the A-A plate and next what we can infer from our experiments about the drivers of A-A deformation.
5.1. Other modeling approaches towards Aegean and Anatolian crustal deformation. Many studies have focused on simplifications and/or subsets of driving forces in the Eastern Mediterranean and the Middle East. For example, time-dependent 3D models of asymmetric continental collision and oceanic subduction (Sternai et al., 2014) suggested that slab rollback and slab-induced mantle flow may drive Aegean extension and Anatolia motion, but the generic design of such models hampers their application to the natural system. A whole-system approach was followed by Boschi et al. (2010), Faccenna and Becker (2010), Becker and Faccenna (2011), Faccenna et al. (2013, who computed global instantaneous mantle flow and surface motions based on different tomographic models using a radial viscosity profile and a uniform-thickness lithosphere without a crust. These global investigations draw attention to the effects of mantle flow on surface motion, as is synthesized in the review by Faccenna et al. (2014). In their mantle flow approach, Arabia-Eurasia convergence accounts for approximately 40% of Anatolia motion, while mantle drag and slab pull account for the remainder. Further, in several models the WSW to SW direction of Aegean crustal flow was observed but with strongly underestimated magnitude, which the authors attribute to the absence of GPE forcing. Our model setup differs predominantly in its regional character, compressibility, lithospheric thickness variations and the use of nonlinear visco-plastic rheology. From our mantle rheology experiments, we suggest that the linear mantle rheology on which the above studies are based, and used more often in global mantle studies of lithosphere deformation (e.g. Lithgow-Bertelloni and Guynn, 2004, Ghosh et al., 2008, Ghosh and Holt, 2012, Wang et al., 2015, may provide a possible explanation for the crustal motion misfits observed by Faccenna and Becker (2010) and Faccenna et al. (2014) for the Aegean region (see Discussion of Section 4.5). Overall, the impact of mantle upwellings on lithosphere motion as modeled in the latter studies is lower in our experiments, which we attribute in part to our implementation of a rheologically strong, but realistic, lithosphere and to the use of nonlinear mantle rheology, reducing the viscous lithosphere-mantle coupling. Another class of models focuses on modeling of the lithosphere separately. Based on such models, GPE gradients associated with surface topography and crustal structure have been put forward as an important driver of surface motion in the eastern Mediterranean region (Özeren and Holt, 2010, Carafa et al., 2015, England et al., 2016. In contrast, the lithosphere sheet models of Jiménez-Munt and Sabadini (2002), Jiménez-Munt et al. (2003),Özbakir et al. (2017 suggest a minor influence of lateral GPE variation. Using a regional model domain for the lithosphere,Özeren and Holt (2010) combined the GPE body forcing derived from heterogeneous lithospheric structure with tuning of the stress imposed on the domain boundary to optimize the fit to crustal deformation derived from GPS observations. A good fit required significant SW pull along the Hellenic trench from the Kefalonia fault to about 33 • E and an increasingly large NNE push further eastward and an on average weak lithosphere viscosity of 10 21 -10 22 Pas, while they found that boundary stress and GPE forcing contributed comparably to crustal deformation. Basal tractions were not considered byÖzeren and Holt (2010) because earlier global studies (Ghosh et al., 2008(Ghosh et al., , 2009 showed that basal tractions directly below the Aegean do not significantly affect the lithosphere above. In contrast, Carafa et al. (2015) proposed that boundary and basal tractions drive Aegean motion, while GPE gradients control Anatolia motion. The basal tractions suggest forcing by plate-mantle interaction, but do not identify the responsible mantle processes. The lithosphere sheet model of England et al. (2016) is forced by topography-based GPE gradients, does not include mantle tractions, assumes uniform crust and lithosphere thicknesses and requires compressional tractions along the eastern A-A region boundary and tensional tractions on the Hellenic boundary for fitting the regional crustal motions. GPE variations and boundary tractions combined reduce the GPS motion RMS of ∼20 mmyr −1 (in a Eurasia-fixed frame) to ∼4.7 mmyr −1 . This is ∼3.2 mmyr −1 smaller compared to a model in which GPE forcing is not included, suggesting that the largest part of the observed motion is explained by the combination of plate boundary forces and thin-sheet viscosity. The average lithosphere viscosity obtained from the data fitting procedure is comparable in magnitude to what is obtained byÖzeren and Holt (2010). This is 1-2 two orders of magnitude smaller than the strong cold core of the lithosphere in this study, whereas our stronger lithosphere still allows for a comparable GPS motion misfit of ∼5.9 mmyr −1 achieved for our Reference model. In sum, the various approaches to modeling Eastern Mediterranean surface flow led to disparate propositions of dominant drivers of surface motion. Our sensitivity study has demonstrated the potential of modeled mantle processes and properties to modulate the A-A surface motions and associated crustal deformation. Taken together, these stress the importance of properly representing the mantle drivers in numerical models of crustal deformation, in particular slab pull and the interaction of the slab with the surrounding mantle in subduction systems. Through systematic testing of sensitivities, either manually like in this study or, preferably, automatically through for example adjoint methods (Colli et al., 2018, Reuber et al., 2018, Crawford et al., 2018, proposed mantle drivers for each specific region can then be ranked or excluded without a priori assumptions about their importance.

Geodynamic drivers of Aegean-Anatolian deformation.
Although we did not strive to find an optimal model that minimizes the misfit of the Reference model motions with the GNSS observations, we identify the relative importance of modeled mantle drivers of surface deformation from our sensitivity analysis. The Reference model misfit of ∼5.9 mmyr −1 is obtained with an elaborate crust-mantle dynamics model with a synthetic slab morphology derived from tomography and seismicity and by tuning the stress transfer across the plate boundary zone surrounding the A-A region. The largest modeled contribution to A-A plate motion and deformation are the gravitational pull of the Aegean slab, plate interactions and the nonlinear viscous interaction of the slab and other plates with the mantle. The dominance of slab-related processes is most notably demonstrated by removing the slab in model No slab. The additional sensitivity of the crustal deformation to mantle drivers summarized in Fig. 7 shows that any sudden changes in slab morphology will have the second largest impact on surface deformation. The A-A crust is further sensitive to the (far-field) mantle buoyancy field, which controls mantle flow around the slab and viscous drag at the base of the lithosphere in concert with the locally slab-induced flow. The deformation sensitivity to various absolute motions (while keeping the same relative plate velocities) is significant but smallest when compared to most other sensitivities (Fig. 7) because of the low velocity of the Nubia plate (∼1 cmyr −1 ), although the mantle response around the slab can be large (Fig. 12). Viscous slab-mantle coupling and crustal deformation sensitivity expectedly increase with the slab dragging speed. Our experiments thus indicate that the current crustal deformation of the A-A region is for the larger part controlled by slab rollback, and, of secondary importance, slab-mantle coupling, plate interaction and mantle flow. All investigated mantle processes and properties however cause a significant imprint on crustal deformation.

Conclusions
Through numerical instantaneous-dynamics modelling of an elaborate compressible crustmantle system, we systematically investigated the sensitivity of horizontal surface deformation to absolute plate motions, slab morphology, buoyancy-driving mantle flow, and sublithospheric mantle rheology for the natural laboratory of the Eastern Mediterranean region. All experiments were conducted with the same crust-lithosphere system without surface topography, thus excluding topography-related GPE contributions to crustal deformation, but including GPE contributions from laterally varying crust and lithosphere thickness. In addition, we adopted a stiff lithosphere rheology that smoothens the model's crustal response. The four sets of experiments demonstrate the large, but variable, impact of each of the modelled mantle processes and the uncertainty in mantle properties on crustal deformation of the Aegean-Anatolian region. We quantified this impact with the sensitivity shown in Fig. 7, defined as the change in amplitude of the effective velocity gradient field of the shallow crust relative to the reference model in Fig. 5.
We conclude the following: • The reference model, which is in terms of 3D mantle buoyancy field, structural makeup, and material properties the most elaborate used so far, achieves an RMS data misfit of ∼5.9 mmyr −1 with respect to GNNS observations (RMS 17.4 mmyr −1 ) of horizontal surface motions. This misfit is comparable to those obtained in previous lithosphere modeling studies, and it demonstrates a large integral sensitivity of the modeled Aegean-Anatolian crustal flow field to mantle processes. This misfit can possibly be reduced by including more lateral rheological heterogeneity in the crust and lithosphere and forces internal to the lithosphere such as lateral gradients in GPE from surface topography. • The modeled horizontal crustal motions of the Aegean-Anatolian plate result for the major part from Aegean slab pull, with important contributions from viscous coupling of the slab to regional mantle flow and plate boundary stress-transfer. The sensitivity of the crustal flow field for slab pull and for remotely-excited mantle flow that impacts the Aegean slab proves large (0.5-1.0 cmyr −1 scale and determining up to 93% of the reference velocity gradient field). • Predominantly the top 200 km of the modeled slab causes the horizontal crustal motions associated with slab rollback. Changes in slab morphology strongly affect surface motions (sensitivity of up to 68%). • The crustal flow field is also sensitive to the regional patterns of buoyancy-driven mantle flow, which drives the surrounding major plates as well as provides basal drag to the Aegean-Anatolian microplate. The exact scaling of seismic velocity anomalies to temperature and thus density has a much smaller impact however (Fig. 7). • The sensitivity of crustal deformation to absolute plate motion is ranked lowest for the Eastern Mediterranean region (∼ mmyr −1 scale, sensitivity less than 10%), but our experiments suggest it could prove important in subduction systems with faster absolute motion of the subducting plate.
• All modeled mantle drivers of crustal deformation are strongly modulated by the mantle rheology. A significant difference in crustal deformation (sensitivities up to ∼45%) was found when using a depth-dependent viscosity profile as compared to the nonlinear visco-plastic rheology used for the reference model. • Horizontal crustal motions are generally sensitive to mantle processes, even when deformation of the lithosphere is suppressed by absence of local weaknesses. The actual crustal response of a region of interest is expected to be larger in amplitude and less smooth in pattern than modeled here with a laterally strong lithosphere. • The crustal flow field is strongly sensitive to the style of subduction, which we varied from slab rollback to slab advance. This leads to diagnostic flow patterns, but only when viewed in the mantle reference frame. Deformation in the crust and in the mantle surrounding the slab may reveal the impact of the viscous coupling between slab, plates and mantle due to slab dragging. • The inclusion of mantle processes in any modeling that aims at interpreting surface motions requires using a mantle frame of reference for both model and data to ensure that the viscous resistance between slab and mantle and between plates and mantle is properly incorporated, as this viscous coupling depends on their respective relative motion. This is particularly important when only kinematic boundaries are used. A mantle frame of reference ensures that for example prescribed, or inverted, basal tractions relate directly to the mantle. • Because crustal motions and deformation may contain signals of mantle processes, any modeling that does not include mantle processes for predicting observations can at best provide a data-sensitivity test for the geodynamic processes under consideration. • Lastly, because of the intrinsic connections between processes in the crust-mantle continuum, trade-offs may exist between the model implementations of geodynamic drivers as a result of uncertainties in absolute plate motions, crust-mantle structure, rheology, temperature and density. Such trade-offs may also lead to non-linear responses when striving to include as many drivers as possible, which calls for systematic sensitivity studies as a basis for data fitting using forward model searches.
Other, independent data sets with different process sensitivity, such as seismological observations of strain-induced mantle anisotropy and surface topography, may bring independent constraints.
Appendix A. 3D initial temperature conditions The construction of our 3D initial temperature distribution starts with a background temperature profile T bg that varies with depth z (Eq. 13). The 1D profile consists of a fixed surface temperature T top , a uniform temperature T a that is prescribed to the depth of the deepest lithospheric plate z Lmax and is the adiabatic temperature at this depth, and an adiabatic reference profile T adiabat (z) that is perturbed by a thermal boundary layer at the CMB δT ltbl (z) (see Tab. 5). T adiabat (z) is computed below the Nubian continent.
where we follow Steinberger and Calderwood (2006) for the calculation of the perturbed adiabatic temperature and the values T top = 285 K and T adiabattop = 1613 K. However, we neglect the latent heat jump in the adiabatic profile that Steinberger and Calderwood (2006) do include. Inside the compositional fields representing the crust and lithospheric mantle, we overwrite T bg (z) with a lithospheric temperature T L (r, c i ) based on position r and plate type. We prescribe a linear temperature gradient in a Continental Plate (CP), the plate cooling model inside an Oceanic Plate (OP) (Schubert et al., 2001) and a plate cooling profile in the Subducting Plates (SP) according to McKenzie (1970): Tab. 5 for the definition of used symbols. In the sublithospheric mantle, we perturb the background temperature with spatially varying temperature anomalies to capture the structure of the mantle as described in Section 3.2.3.

Appendix B. Crustal (P,T) lookup tables
We created material properties lookup tables for crustal (P,T) conditions by extending a small python program supplied with ASPECT, using the following expressions for density ρ(T, P ) = ρ 0 (1 − α (T − T adiabat )) exp (βP ) , thermal expansivity (Schmeling et al., 2003) and specific heat For parameter definitions and values, see Tab. 6. V (T,P ) V 0 (T ) is approximately a third-order Birch-Murnaghan equation of state (Schmeling et al., 2003). The specific heat formulation is a third-order fit to the line MDG in Fig. 2 of Jacobs and Van den Berg (2011), which is taken from Stixrude and Lithgow-Bertelloni (2005). Figure 26. a-f) In red, depth profiles of temperature, pressure, thermal expansivity, specific heat, density and viscosity for an experiment (model No tomo) without sublithospheric mantle heterogeneities. The profiles are taken in stable Eurasia, far from the plate boundary zones. In black, profiles from the literature for reference; the largest differences arise from the presence c.q. absence of crustal and thermal boundary layers. g&h) Compressional and shear wave velocity temperature derivatives with depth. Solid lines indicate the total of anharmonic and anelastic contributions to the derivatives. The solid red line of Karato (2008) represents our reference profile, with black dashed lines indicating a 20% (10%) uncertainty in the lower (upper) mantle. Dashed and dotted colored lines delineate one standard deviation for the other profiles. Table 6. Parameter definitions and values used in the computation of the oceanic and continental crust (P,T) lookup tables.
Appendix D. Supplement 2 -Flow predictions of experiment set 3 Here we present the crust and mantle flow predictions for each of the models with mantle structure from tomographic models other than UU-P06 CSloc as discussed in Section 4.4.2. Details of the tomographic models can be found in Tab. 7. D.1. CUB S20RTS. The CUB S20RTS model shows reduced overriding plate motion towards the trench (about 2 mmyr −1 , Fig. 28). This coincides with smaller upper mantle velocities in the eastern Anatolian region (Fig. 18). Compared to the model without sublithospheric mantle heterogeneity, toroidal flow around the slab's eastern edge is also reduced. Rollback of the slab seen in the vertical slice is less than in the reference model, as is the downward flow underneath the slab tip. This could be due to the smaller and more laterally spread out cold lower mantle anomaly. The former probably results from the smaller, but still present, SW push against the slab from upwelling underneath eastern Anatolia.

D.2.
TX2015. The TX2015 model shows little upper mantle flow towards the trench underneath Anatolia and as such resembles the model without mantle heterogeneities (Fig. 18). Indeed even w.r.t to the latter model, mantle flow beneath Anatolia seems somehow hindered and inflow above the slab derives mainly from toroidal flow. This leads (a) Velocity field at 5 km depth in model MRF of motion. In green the current model predictions, in brown the Reference model predictions.
(b) Magnitude and vectors of the velocity difference field at 5 km depth (current model minus the Reference model). For a proper comparison between the predicted crustal flow field and that of the reference model, we first rotate the crustal motions of the Reference model into the new plate motion MRF with the same rigid-plate Euler rotation as used in creating the new MRF, before computing the difference between the two flow fields. Note the different velocity scales.
(c) Vertical slice through the Aegean slab showing the velocity field together with black velocity contours every cmyr −1 exceeding the color bar. Temperature isocontours at 1600, 1700, 1800, 1900 and 2000 K are shown in white. The velocity vectors, interpolated onto a regular grid, are not plotted inside the slab. Black lines at 410 and 660 km depth are included for scale.   Table 7. Summary of the seismic tomography models used in this study. All models are sampled at a resolution of 0.5 • x 0.5 • x 10 km before conversion to temperature anomalies. to up to 7 mmyr −1 less trenchward motion of the overriding plate (Fig. 29, upper panels). As can be seen from the bottom panels in Fig. 29, slab motion is subvertical and flow patterns in the upper mantle are smooth. Notably, flow underneath Nubia is directed towards the slab (∼NNE). D.3. S40RTS. In model S40RTS, upper mantle flow from below the Arabian peninsula is conveyed underneath Anatolia to the above-slab region (rotational, westward motion as seen in Fig. 18). It is diverted in the uppermost mantle by the Cyprus slab, while at deeper depths flow is allowed to be more northwest until close to the slab. There is no additional mantle upwelling underneath eastern Anatolia however. Also, the Aegean slab rolls back less than in the reference model. Therefore trenchward motion of the overriding plate is reduced by up to 4 mmyr −1 (Fig. 30) and has a larger northward component. D.4. SL2013 S40RTS. In this model, the upper mantle heterogeneities derive from the SL2013 tomographic model, while the lower mantle anomalies are still derived from the S40RTS model. Overriding plate velocity increases compared to S40RTS to Reference model values (see Tab. 4). Sensitivity measures show both lower and higher values for SL2013 S40RTS (Fig. 7). Southward motion originating from the Black Sea combined with a westward-turning flow coming from underneath Arabia results in a counterclockwise rotational difference of up to 2 mmyr −1 in the whole Aegea-Anatolia-Arabia region (Fig. 31). The flow field underneath the lithosphere shows a very fast upwelling underneath the northwestern edge of the Black Sea that directly feeds the region above the Aegean slab (see Figs. 18 and 31) and probably causes the more southward motion of Aegea. Southwestward flow away from the back of the slab is comparable to the Reference model, while toroidal flow components are reduced.
D.5. SAVANI. For the SAVANI model, the crustal flow pattern in the Anatolian region shows more northward directions, while in the Aegean Sea directions correspond well with the reference model, but the magnitude of velocity is lower (Fig. 32). Aegean motion is reduced by up to 6 mmyr −1 (Tab. 4). The mantle flow field at 200 km depth shows stronger toroidal flow than model No tomo and strong northward flow under the northern part of Arabia. There is little southwestward flow directly behind the slab around the 410 discontinuity, but NNE flow underneath the Mediterranean basin and Nubia parallel to slab rollback (bottom images of Fig. 32). This push against the slab and the absence of flow from eastern Anatolia reduce Aegean motion, while flow from the Arabia upwelling deflects Anatolian crustal motion northward. D.6. LLNL G3D. The LLNL G3D model shows a very strong northward inflow of mantle material through the gap between Aegean and Cyprus slab (Fig. 18). Rotational flow around the western slab edge is also stronger than in the reference model. There is no contribution to inflow above the slab from the Anatolian region however. Slab motion is reduced and mostly downward, perhaps due to the NNE mantle flow from underneath Nubia that deflects around the slab. This drop in rollback compared to the reference model and the absence of flow from underneath Anatolia lead to the largest reduction in trenchward motion of the overriding plate (up to 9 mmyr −1 in the Aegean) and more northward motion of Anatolia (Fig. 33). LLNL G3D also has the highest sensitivity measures of all experiments in set 3 (Fig. 7). D.7. GAP P4. Aegean crustal motion is reduced and more northward when the GAP P4 tomographic model is used to generate mantle structure (Fig. 34). Anatolian motion is not so much reduced in magnitude, but directed more northward. Strong mantle flow into the area is seen via the Aegean and Cyprus slab window and from underneath eastern Anatolia and Arabia (Fig. 18). Flow underneath Nubia is northward and turns NNW west of -3 • E (in our coordinate system). Interaction of this flow with the slab and overriding plate leads to smaller velocities (by up to 6 mmyr −1 smaller) towards the trench and more northward directed motion of the whole Aegean-Anatolian region.
D.8. SEMUCB. The SEMUCB tomography model results in very smooth temperature variations in the mantle (Fig. 35 bottom left). There is little flow coming from the eastern Anatolian region, although NNE flow from underneath Nubia enters the zone above the slab through the Aegean-Cyprus slab window. Slab motion is reduced and vertical, as seen in the vertical slices of Fig. 35. The combined effect is lessened Aegean motion compared to the reference model (Fig. 35), up to 8 mmyr −1 . Northward motion of Arabia is increased due to N to NNE mantle flow underneath.