Table of Contents

i. Abstract!""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""! #! ii. La déformation de la plaque supérieure dans les zones de subduction en retrait (version française abrégée)! """"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""! $!

Upper plate deformation in retreating subduction zones i. Abstract

The Earth's surface is constantly reshaped by the tectonic plate motion, which is mainly driven by subduction of plates into the deeper mantle. Subduction trenches are also mobile plate boundaries, and are observed to retreat towards the subducting plate or advance towards the upper plate over geological time. Trench retreat has been historically thought to cause extension in the upper plate above the subducting slab.

However, natural subduction systems show several examples of retreating trenches that are associated with upper-plate compression. This thesis explores upper plate (back-arc) deformation in retreating subduction systems. Three techniques are used: large-scale numerical models addressing physical processes; seismic profiles in the Central Aegean addressing basin-scale fault patterns; and field-scale observations clarifying fault kinematics in the Central Aegean. The large-scale thermo-mechanical models deal with viscous deformation of the upper plate, and investigate the relationship between slab pull, slab rollback, trench retreat and upper plate deformation at scales of 100 to 1000 km. They show that asthenosphere flows below the plates (100-200 km depth) can control both trench retreat and upper plate deformation. The type of deformation in the upper plate also depends on the plate's far-field conditions: if the plate is free to move, deformation tends to be compressive, but a fixed upper plate shows extension. The latter is comparable to the Aegean region, an upper plate exhibiting extension above a narrow, retreating subduction zone. Related extensional structures in the central Aegean have been analysed from seismic and field data, revealing co-existing normal, oblique and strike slip faults. These features reflect a combination of rollback-related extension and extrusion-related strike slip activity.

Resulting block rotation and trench retreat re-activate inherited normal faults in obliquenormal slip, while new pure-normal faults are created. We also infer a recent change in stress state possibly related to the slab tear on the western side of the Hellenic slab.

Additionally, accelerated trench retreat and upper plate extension are the cause of the Aegean's high surface heat flow, which makes it potentially suitable for geothermal energy production. As a final perspective on the application of geodynamic research, an assessment of the role of tectonic modelling in predicting geothermal energy potential is presented, using the stretched Aegean upper plate as an example.

ii. La déformation de la plaque supérieure dans les zones de subduction en retrait (version française abrégée)

ii.b. Abstract La surface de la Terre est en permanence remodelée par les mouvements des plaques tectoniques, dont le moteur principal est la subduction, i.e. le plongement de plaques océaniques dans le manteau profond. Les fosses océaniques de subduction constituent également des limites de plaques mobiles, et les observations montrent que, sur des échelles de temps géologiques de plusieurs millions d'années, ces fosses reculent (vers la plaque plongeante) ou avancent (vers la plaque chevauchante/supérieure). Historiquement, le retrait de la fosse a été associé à une extension de la plaque supérieure au-dessus du panneau plongeant. Cependant, les zones de subduction sur Terre montrent plusieurs exemples de fosses en recul associées à des contraintes compressives. Cette thèse étudie la déformation (arrièrearc) de la plaque supérieure pour une subduction en retrait. Trois approches ont été utilisées : des modèles numériques explorant les processus physiques mis en jeu à grande échelle, des profils sismiques en mer Égée centrale permettant d'étudier la répartition des failles à l'échelle du bassin, et des observations de terrain pour caractériser l'évolution temporelle de la déformation de la plaque supérieure en mer Égée centrale. Les modèles thermo-mécaniques à grande échelle reproduisent une déformation visqueuse de la plaque supérieure, et permettent d'analyser les relations entre traction du slab, recul du slab, retrait de la fosse et déformation de la plaque supérieure, à des échelles allant de 100 à 1000 km. Ils montrent que des courants dans le manteau asthénosphérique sous les plaques (vers 100-200 km de profondeur) peuvent contrôler à la fois le mouvement relatif de la fosse et la déformation de la plaque supérieure. Cette dernière dépend également des conditions mécaniques aux limites: si la plaque est libre de bouger, sa déformation sera plutôt compressive ; mais une plaque fixe sera en extension. Ce dernier cas est comparable à la région de la mer Égée, une plaque supérieure montrant de l'extension et associée à une zone de subduction étroite en retrait. Les structures extensives associées ont été analysées grâce à l'observation sur le terrain et à l'étude de profils sismiques, révélant des failles normales, obliques et décrochantes synchrones. Cela est interprété comme résultant de la combinaison de contraintes extensives associées au recul de la fosse et de contraintes décrochantes associées à l'extrusion d'un bloc voisin. La rotation et le recul de la fosse réactivent d'anciennes failles normales dans un mode oblique-extensif, et engendrent des nouvelles failles purement normales. Les données suggèrent également un changement récent de l'état de contrainte mécanique dans la plaque, qui pourrait être dû à une déchirure du panneau plongeant côté Ouest. En sus, l'accélération du recul de la fosse et l'intensification de l'extension de la plaque supérieure expliquent probablement le flux de chaleur élevé en mer Égée, ce qui rend l'énergie géothermique potentiellement exploitable dans cette zone. Une évaluation de l'apport de la modélisation tectonique pour prédire le potentiel géothermique est finalement présentée comme perspective de l'application des recherches en géodynamique, s'appuyant sur l'exemple de la plaque supérieure égéenne amincie.

ii.b. Synthèse

La subduction est le principal moteur des mouvements des plaques tectoniques et créer certaines des caractéristiques les plus importantes sur la surface de la Terre. La fosse, l'interface entre la plaque subductante et la plaque supérieure, est une limite mobile. Cette dernière peut soit avancer vers la plaque supérieure soit reculer vers la plaque subductante (Figure 1). Le mouvement de la fosse est aussi connu pour créer des déformations au niveau de la plaque supérieure, mais ces mouvements sont associés à la fois à des plaques supérieures en compression mais aussi en extension. † Goes, Saskia, Roberto Agrusta, Jeroen van Hunen, and Fanny Garel. Cette thèse a traité le mouvement du retrait de la fosse et la déformation de la plaque supérieure associée en adressant trois questions majeures :

-Est-ce que le tirage de la plaque plongeante contrôle le recul de la fosse ? -Quelle est la relation entre le recul de la fosse et la déformation de la plaque supérieure ?

-La recherche en tectonique peut-elle alimenter directement des applications pouvant être utiles à la société ? Dans le but de répondre à ces questions, une approche multi-échelle a été utilisée pour mieux comprendre les processus physiques dans les zones de subduction en retrait et ce qui contrôle la déformation de la plaque supérieure associée à ces dernières. Nous avons utilisé des modèles numériques à grande échelle pour isoler les effets du tirage ‡ Heuret, Arnauld, and Serge Lallemand. 2005. "Plate Motions, Slab Dynamics and Back-Arc Deformation." Physics of the Earth and Planetary Interiors 149 (1-2 SPEC. ISS.): 31-51. https://doi.org /10.1016/j.pepi.2004.08.022. de la plaque subductante, premièrement, sur le recul de la fosse et deuxièmement, sur la déformation de la plaque supérieure (Figure 4). Nous avons aussi utilisé des données de terrain et de sismiques au niveau de la plaque supérieure Egéenne pour mieux comprendre la disposition des failles fragiles à l'échelle des bassins dans les systèmes de subduction en retrait.

Figure 4: Paramétrisation initiale du modèle numérique

Nos modèles montrent que le recul de la fosse n'est pas seulement contrôlé par le tirage de la plaque plongeante, mais aussi par le flux peu profond du manteau environnant. Les tirages de plaque plongeante les plus forts induisent des subductions plus rapides, ce qui déclenche à son tour un écoulement de manteau plus rapide sous les deux types de plaques. Ces flux de manteau plus rapides se confrontent, et la fosse peut reculer si le retour de flux dans le coin mantellique (sous la plaque supérieure) est plus rapide que celui situé sous la plaque plongeante. Cependant si cette dernière subducte rapidement et entraine le manteau situé juste en dessous de manière à ce qu'il s'écoule plus rapidement que celui dans le coin mantellique, il pourra alors résister au recul de la plaque plongeante et de la fosse même si le tirage est fort et la subduction rapide. Ceci montre que le tirage de la plaque plongeante n'est pas seulement le facteur affectant le mouvement de la fosse mais que le flux mantellique (que cela résulte ou non de la traction de la plaque inférieure) a aussi un fort contrôle sur le mouvement de la fosse.

Nous avons constaté que l'écoulement du manteau est également responsable de la déformation de la plaque supérieure ; ainsi, le mouvement de la fosse ne peut causer indirectement la déformation de la plaque supérieure que s'il déclenche un mouvement plus rapide de l'écoulement du manteau. Nos modèles montrent que si l'écoulement du manteau sous la plaque supérieure est suffisamment fort pour dépasser la résistance de la plaque supérieure, il peut entrainer la plaque supérieure en profondeur et créer de la déformation. Le type de déformation induit dépend alors du fait que la plaque soit entrainée dans son entièreté ou que seul des portions le soient. La plaque peut être entrainée entièrement si elle présente un bord libre, par exemple une dorsale à l'arrière de sa direction d'entrainement. Dans ce cas, la plaque supérieure se déforme en compression. Par ailleurs, si seules des portions de plaque sont entrainées, comme dans le cas présent lorsque la plaque est reliée à l'arrière a un large continent, elle se déforme en extension lorsqu'elle est entrainée au front vers la fosse. Dans les deux cas, la résistance de la plaque supérieure a aussi son importance, puisqu'elle doit être suffisamment faible pour être déformée par le flux mantellique sous-jacent. Une autre façon d'exprimer cette idée est : le flux mantellique sous-jacent doit être assez fort pour dépasser la résistance de la plaque. Nous illustrons ce concept Figure 5.

Figure 5 : Interprétation de processus physiques issus des modèles numériques

Ceci a un autre effet sur le retrait de la fosse. Si la plaque supérieure est fixe, i.e. incapable de translater latéralement, la fosse ne peut reculer que si la plaque supérieure se déforme. Sans déformation, le mouvement de la fosse est bloqué. Ceci fait apparaitre une nouvelle perspective concernant le retrait de la fosse, étant donné qu'il a toujours été considéré comme la cause de la déformation de la plaque supérieure. Notre modèle montre qu'il peut aussi être le résultat de la déformation de la plaque supérieure.

Ceci peut également être pertinent pour le système Égéen, où la fosse se retire dans une étroite zone de subduction, sous une plaque supérieure relativement fixée. Le retrait de la fosse s'est possiblement initié en raison d'une plus forte traction de plaque plongeante, mais ceci a pu engendrer un flux mantellique suffisamment rapide pour qu'il soit présentement le principal moteur. Par conséquent, la déformation dans la mer Égée rend compte de l'évolution du recul de la plaque plongeante, comme le montre le chapitre 3. Ici, nous avons un changement progressif dans le style des failles depuis l'exhumation et la formation localisée de failles à faible pendage, jusqu'à la segmentation et la formation distribuée de failles à fort pendage. D'après les résultats de nos modèles, la segmentation a été rendue possible du fait que la plaque soit préalablement affaiblie par la précédente exhumation, permettant au manteau sousjacent de dépasser la résistance de la plaque et créant une déformation plus pervasive. Nos données sur le domaine Égéen montrent la coexistence de failles normales, obliques et décrochantes (Figure 6). Alors que la présence des failles normales et décrochantes est prévisible avec les processus de recul de la plaque plongeante et de l'extrusion Anatolienne respectivement, la présence de failles obliques ajoute un élément de complexité. Ces dernières suggèrent que les générations précédentes de failles normales ne sont plus alignées avec la direction d'extension, ce qui entraîne un glissement oblique. Notre interprétation de la mer Égée centrale montre que la géométrie des failles lors de la segmentation tourne avec la rotation du recul de la plaque plongeante. Les orientations des failles (toutes les générations) suggèrent qu'il y a eu un changement dans l'état de contrainte en plus (de l'interprétation) de la rotation des blocs. Cela pourrait représenter la rotation des contraintes résultant de la déchirure de la dalle dans la mer Égée occidentale (Figure 7).

Introduction

Preface

The Earth is a complex dynamo with a hot interior and a rigid surface broken up into plates. These plates move and interact with the Earth's interior over geological time, as described by the theory of Plate Tectonics (Le [START_REF] Pichon | Sea-Floor Spreading and Continental Drift[END_REF]. Plates are, thus, defined by mechanical rather than compositional layer designation, and make up the lithosphere, the Earth's outermost mechanical layer, which lies above the asthenosphere [START_REF] Montagner | Earth's Structure, Global[END_REF]. Compositionally, the lithosphere includes the crust and part of the upper mantle. Tectonic plate motion is thought to be largely driven by subduction, the sinking of lithosphere into the asthenosphere (e.g. Forsythe and [START_REF] Forsyth | On the Relative Importance of the Driving Forces of Plate Motion[END_REF]. There are currently tens of active subduction zones on Earth, which are shown in The portion of the plate that has already subducted beneath the surface is called the slab. It is thought that subduction is mostly driven by the slab pulling the rest of the plate since it is heavier than the surrounding mantle. This is called the slab pull force [START_REF] Forsyth | On the Relative Importance of the Driving Forces of Plate Motion[END_REF][START_REF] Davies | Regional Compensation of Subducted Lithosphere: Effects on Geoid, Gravity and Topography from a Preliminary Model[END_REF][START_REF] Conrad | How Mantle Slabs Drive Plate Tectonics[END_REF]. Slab pull occurs because the slab is mainly made of cold lithosphere, which is heavier than the surrounding hot asthenospheric mantle. As the slab subducts through the upper mantle, it starts to warm up, but the increased pressure first metamorphoses the basaltic oceanic crust to eclogite, and then the olivine rich lithospheric mantle to wadsleyite in the mid-upper mantle. This significantly raises the density of the slab, raising slab pull even more (Figure 1.2). resisting forces result from mantle drag, slab bending, upper plate strength and the upper-lower mantle phase transition. (Goes et al., 2017) While slab pull is the largest force in the system, the subduction zone (SZ hereafter) is also affected by several other forces. In addition to slab pull, ridge push acts as a driving force. However, resistance occurs due to lithospheric bending, frictional platecoupling, viscous shear in the mantle, and positive buoyancy due to the ringwooditeto-perovskite/ magnesiowustite phase transition at a depth of 660 km, the bottom of the upper mantle [START_REF] Billen | Modeling the Dynamics of Subducting Slabs[END_REF]. These forces are shown conceptually in Figure 1.2. The variable strength of each of these forces, as well as other variations, for example to farfield plate motions, cause the characteristics of each SZ to vary. For example, SZs are not generally stationary, but move over geological time. The movement is conventionally ascribed to the trench, the interface between the subducting plate and the upper (or overriding) plate (e.g. Chase, 1978). Trenches are usually described to advance, i.e. move towards the upper plate, or retreat, i.e. move towards the subducting plate. Since all plates move relative to each other, a single reference point must be defined in order to describe surface motion. Hot spots are usually used to define different reference frames. (Goes et al., 2017).

Most of the world's trenches move no faster than 5 cm/yr (Heuret and Lallemand., 2005), however, the world's fastest trench, in the northern Tonga SZ, retreats at a rate of 14.5 cm/yr (Lallemand et al., 2008). Figure 1.1 shows the normal component of global trench motion. Statistical data shows that there are approximately as many retreating trenches as there are advancing ones (Heuret and Lallemand, 2005). This challenges our current understanding that slabs ought to spontaneously rollback and retreat because of their negative buoyancy (e.g. [START_REF] Molnar | Interarc Spreading and Cordilleran Tectonics as Alternates Related to the Age of Subducted Oceanic Lithosphere[END_REF][START_REF] Dewey | Episodicity, Sequence, and Style at Convergent Plate Boundaries[END_REF]Garfunkel et al., 1986). The lack of observed correlation between trench retreat and negative buoyancy (i.e. slab pull resulting from plate age, see Figure 1.3), also challenges this notion (e.g. Heuret and Lallemand., 2005;Sdrolias and Müller, 2006;Goes et al., 2017). Contrastingly, trench motion seems to correlate better to upper plate motion, although there are also many exceptions to this statement (e.g. Heuret and Lallemand., 2005;Lallemand et al., 2008). This leads us to question our understanding of the system's driving forces, raising the first question of this thesis: how does slab pull influence trench retreat?

The second question of this thesis relates to an apparent consequence of trench retreat: deformation of the upper plate away from the trench and beyond the volcanic arc, or what is known as the back-arc. Slab rollback and trench retreat have often been assumed to cause extensional stresses in the upper plate leading to back arc spreading (e.g. [START_REF] Molnar | Interarc Spreading and Cordilleran Tectonics as Alternates Related to the Age of Subducted Oceanic Lithosphere[END_REF]. But again, there is no clear correlation between trench retreat and back arc extension, as shown in Figure 1.4 (Heuret and Lallemand., 2005). This observation poses the second question of this thesis: is there a relationship between trench retreat and upper plate deformation? (Heuret and Lallemand, 2005).

In order to answer these two questions, we have approached the problem from two different angles. The first attempts to improve our understanding of the physical processes controlling trench retreat and associated upper plate deformation. To do this, we have constructed thermo-mechanical numerical models on the scale of the upper mantle, dealing with viscous deformation on the time-scale of several million years. The second angle attempts to improve our understanding of the system using observations of brittle deformation from a retreating subduction zone, where the rate of trench retreat and associated upper plate deformation has changed through time. To do this, we use field and seismic data from the Aegean upper plate. Using these approaches allows us to tackle the problem at different scales; the scale of the upper mantle, basin scale, and local scale; while also touching on both the physical processes and natural observations. While we primarily deal with the fundamental scientific concepts relating to upper plate deformation in retreating SZs, we expect results from our multi-scale approach to be more easily disseminated into applied studies, and we pose a third question in this thesis: Can tectonic research (such as performed here) feed back directly into applications that can be useful for society? To answer this third and final question, a brief study was conducted to see if basin-scale tectonic modelling can improve thermo-mechanical predictions of geothermal energy potential, using the Aegean as an example. Sections 1.2 -1.4 below provide the background to each study performed: large scale numerical modelling of upper plate deformation in retreating SZs; case study of brittle extension in the Aegean; and potential use of tectonic modelling in predicting geothermal potential.

Summary of objectives

This thesis attempts to better understand upper plate deformation in retreating subduction zones, by answering the following three questions:

• How does slab pull influence trench retreat?

• Is there a relationship between trench retreat and upper plate deformation?

• Can tectonic research feed back directly into applications that can be useful for society?

We tackle this with a multi-scale approach, using different tools to observe different parts of the retreating subduction system. We, thus, use large-scale numerical models, basin-scale seismic analysis, and local-scale field studies. We then combine basin scale observations with thermomechanical modeling to assess the use of tectonic research in predicting geothermal energy potential.

Thus, rather than focusing on a single specific detail or a particular method, our strategy is to combine multiple angles to gain a better understanding of upper plate deformation in retreating subduction zones, and how the system's main controls relate to each other.

1.2.

Numerical modelling of subduction trench retreat and upper plate deformation

Geodynamic models of subduction

Studying lithosphere and asthenosphere dynamics, such as in subduction, usually falls under Geodynamics, a discipline that combines structural geology, geochemistry and geophysics. This mostly deals with the subsurface, which cannot be directly observed (the deepest borehole drilled to date is 12 km deep, while the mantle extends to 2900 km depth). Therefore, much of the discipline started from mathematical constructs of physical concepts. Recent improvements in data acquisition, such seismic tomography and long-series GPS, has added crucial constraints, improving geodynamic models.

Theoretical models and observations from nature (both direct and indirect) are, thus, strongly complimentary and can produce sophisticated models of natural processes.

Theoretical models are often first inspired by concepts derived from sparse observations. Combining these with physics-based predictions and then comparing the model results to natural observations can largely improve our understanding of the unobserved processes. From the early days of tectonic research, analogue modelling has played an important role in improving our understanding of plate motion and subduction (e.g. [START_REF] Griffiths | A Laboratory Investigation of Effects of Trench Migration on the Descent of Subducted Slabs[END_REF].

For example, analogue models revealed that trench retreat could be responsible for back-arc extension (e.g. Figure 1. 2.1, Faccenna et al., 2001). However, these models generally did not include upper plates, and incorporating upper plates (into numerical models) revealed that back-arc extension might relate to other factors than trench retreat (e.g. Capitanio et a., 2010). While analogue models are useful for studying slab dynamics, retreating subduction with an upper plate remains difficult to model experimentally (but not impossible as shown by Pitard et al., 2018, amongst others). This is more easily handled by numerical models, which have numerous other advantages for tackling upper plate deformation in retreating subduction zones. For example, different scales are more easily captured by numerical models, which is particularly useful for large-scale setups with variable mesh size. Also, model evolution is not limited to kinematics but can also include thermal and compositional evolution. Model setup and parametrisation can also be varied more quickly numerically, allowing us to more efficiently test different scenarios. Numerical model results are also inherently quantified, facilitating analysis of the findings. We have, therefore, chosen to model trench retreat and upper plate deformation numerically.

Numerical modelling of subduction processes

Numerical modelling of subduction processes is usually based on solving conservation equations, often applied to fluid dynamics. The plate and the asthenosphere are, thus, treated as viscous fluids deforming over multi-million year timescales [START_REF] Billen | Modeling the Dynamics of Subducting Slabs[END_REF]. The thickness of the plate is defined by the lithosphere-asthenosphere boundary, or the LAB. The LAB can be defined by different proxies, such as shown in Figure 1.2.2 [START_REF] Eaton | The Elusive Lithosphere-Asthenosphere Boundary (LAB) beneath Cratons[END_REF]. In the models we present in Chapter 2, we define the LAB thermally, since strain rate is calculated from viscosity, which in turn is calculated from temperature.

Thus, in our models, the LAB is defined thermally, and mechanical decoupling between the plates and the lithosphere is maintained; the plates have lower strain rate than the asthenosphere, making them more rigid. There are various types of subduction zone numerical models. The first distinction is between 2-D and 3-D models (e.g. Figure 1.2.3). Trench motion and upper plate deformation are inherently 3-dimensional phenomena. If we also consider mantle flow beneath the plates, the role of lateral influences becomes even more important.

However, 2-D models can be useful to better identify effects that are not related to lateral influences, and these models are less computationally expensive. This makes 2-D models relevant for simplifying a specific problem and gaining a more fundamental understanding of the processes at play. advance [START_REF] Clark | Episodicity in Back-Arc Tectonic Regimes[END_REF].

Regardless of model dimensions, there are several types of numerical models, the choice of which depends on the question addressed. [START_REF] Billen | Modeling the Dynamics of Subducting Slabs[END_REF] summarises the main model types as:

-Instantaneous models: the equations of conservation of mass and momentum are solved, but not the conservation of energy (the thermal field is not advected forward in time). These have the advantage of enabling comparison of results to a large range of observations, e.g. the observed geoid was found to require an increase in viscosity between the upper and lower mantle [START_REF] Hager | Constraints on Mantle Rheology and Flow Residual Geoid : Degree 2-10[END_REF].

-Dynamic models: all three conservation equations are solved (mass, momentum, and energy), updating the thermal and viscosity field (as well as the composition and phase changes where applicable) for each time step. These are useful for exploring slab dynamics, such as rheologic controls on the slab, e.g. [START_REF] Billen | Rheologic Controls on Slab Dynamics[END_REF]. In this category, we can also distinguish models in which slab sinking drives asthenospheric flow (e.g. Garel et al., 2014) from models in which subduction is a consequence of whole-mantle thermal convection (e.g. Crameri and Tackley, 2014).

-A combination of dynamic and kinematic models: the three conservation equations are solved (dynamic), but with part of the model domain having prescribed velocity (kinematic). Models of this type can be purely thermomechanical or include composition, where the latter tend to be more complex. A common use of these models is to combine slab dynamics with prescribed surface kinematics such as observed plate motion (e.g. [START_REF] Tan | Slabs in the Lower Mantle and Their Modulation of Plume Formation[END_REF].

Our work, presented in Chapter 2, uses 2-D, fully dynamic thermo-mechanical models driven by slab sinking. These allow us to identify the effect of slab pull without obscuring plate deformation with prescribed surface velocities.

Previous numerical models of upper plate deformation in retreating subduction systems

Controls on trench motion

Many numerical studies attempt to understand trench motion, as it is not yet clear whether trench motion is a cause or effect of slab dynamics [START_REF] Billen | Modeling the Dynamics of Subducting Slabs[END_REF][START_REF] Gerya | Future Directions in Subduction Modeling[END_REF].

Some studies have found that trench motion correlates with slab dip, where steep slabs are likely to result in trench advance, and gently inclined slabs to trench retreat (e.g. [START_REF] Bellahsen | Dynamics of Subduction and Plate Motion in Laboratory Experiments: Insights into the 'Plate Tectonics' Behavior of the Earth[END_REF][START_REF] Faccenna | Predicting Trench and Plate Motion from the Dynamics of a Strong Slab[END_REF][START_REF] Giuseppe | Slab Stiffness Control of Trench Motion: Insights from Numerical Models[END_REF][START_REF] Stegman | A Regime Diagram for Subduction Styles from 3-D Numerical Models of Free Subduction[END_REF][START_REF] Gerault | The Role of Slabs and Oceanic Plate Geometry in the Net Rotation of the Lithosphere , Trench Motions , and Slab Return Flow To Cite This Version : HAL Id : Insu-00695975 The Role of Slabs and Oceanic Plate Geometry in the Net Rotation of the Lithosphere[END_REF]. This correlation between slab dip and trench motion is also observed in nature (e.g. Lallemand et al., 2005). Other studies link trench motion to mantle dynamics, e.g. slow down of slab rollback and trench retreat due to slab anchoring in the lower mantle (e.g. [START_REF] Billen | Modeling the Dynamics of Subducting Slabs[END_REF]Enns et al., 2005).

Upper plate deformation

Various studies have also looked at upper plate deformation associated with trench retreat. Focusing on retreating systems, UP extension is mechanically expected, according to early studies such as [START_REF] Dvorkin | Narrow Subducting Slabs and the Origin of Backarc Basins[END_REF][START_REF] Faccenna | The Dynamics of Back-Arc Extension: An Experimental Approach to the Opening of the Tyrrhenian Sea[END_REF]Jolivet et al., 1994. Thus, understanding UP compression in retreating systems (e.g. in the South American SZ) has been the focus of many studies. Additionally, upper plate extension has been found to be more complex than being directly caused by trench retreat, as the other forces (Figure 1.2) in the system come into play.

Some studies have focused on the motion of the UP to understand its deformation. For example, Silver et al. (1998) found the South American plate acceleration towards the trench (due to ridge push) caused Andean compression, while [START_REF] Hunen | Various Mechanisms to Induce Present-Day Shallow Flat Subduction and Implications for the Younger Earth: A Numerical Parameter Study[END_REF] also found that this UP motion can cause slab flattening and UP obduction. Sobolev and Babeyko (2005) found that in addition to UP acceleration towards the trench, the crustal structure of the UP, as well as the degree of shear coupling at the plate interface, were also important controls, particularly in the case of the Andes. Other studies attribute UP compressive deformation (again, taking the example of the Andes) to slab buoyancy, e.g. by subduction of ridges (e.g. [START_REF] Gutscher | Geodynamics of Flat Subduction: Seismicity and Tomographic Constraints from the Andean Margin[END_REF][START_REF] Jadamec | Three-Dimensional Numerical Models of Flat Slab Subduction and the Denali Fault Driving Deformation in South-Central Alaska[END_REF]. This leads us to investigate how plate and trench motion could vary in relation to slab pull and mantle flow.

As for retreat-related extension in the UP, several studies have attributed it to the speed of subduction caused by slab buoyancy (e.g. Capitanio et al., 2011). Additionally, [START_REF] Clark | Episodicity in Back-Arc Tectonic Regimes[END_REF] found that in this scenario, back-arc extension tends not to be continuous, but episodic, where the deformation can be pseudo-, quasi-or hyper- For both extensional and compressional deformation in the UP, deformation has also often been attributed to shallow asthenospheric flow induced by subduction (e.g. Wdowinski et al., 1989;[START_REF] Schellart | Cenozoic Tectonics of Western North America Controlled by Evolving Width of Farallon Slab[END_REF]Husson et al., 2012;Schellart and Moresi, 2013;[START_REF] Chen | Does Subduction-Induced Mantle Flow Drive Backarc Extension ?[END_REF]Faccenna et al., 2017). Several studies argue that compression is associated with poloidal flow, where the flows vertically beneath the slab, while extension is associated with toroidal flow, where the mantle flows laterally around the slab (illustrated in Figure 1.2.4). Many of these studies argue that since toroidal flow cannot be captured in 2-D models, 3-D models can more accurately constrain UP extension (e.g. [START_REF] Piromallo | Three-Dimensional Instantaneous Mantle Flow Induced by Subduction[END_REF]Capitanio et al., 2011). These models can also capture more complex slab dynamics, such as slab tearing, which has been show to further enhance trench retreat and mantle flow, dragging the UP from below and introducing extra strain (e.g. Sternai et al., 2014). However, depending on the scenario studied, 2-D models can suffice. For example, by modelling Andean deformation in 3-D, Schellart (2017) showed that poloidal flow produce both extension and compression depending on the size of the convection cell. Here, convection on the scale of the whole mantle scaled with the entire UP, dragging the entire plate and generation compression. However, convection on the scale of the upper mantle (before the slab reaches the lower mantle) is smaller scale and only drags the front of the upper plate, generating back-arc extension. This is an effect that would be captured in 2-D (Figure 1.2.5), as shown by Schellart and Moresi (2013) and Faccenna et al. (2017).

Apart from model dimensions, the choice of boundary conditions is also important. The type of side boundaries can control subduction geometry (e.g. Enns et al., 2005, refer to Chapter 2.4.1 for more). However, the choice of surface boundary conditions can largely influence the observed results when exploring UP deformation. Models benefit from using a free surface as the upper boundary condition to allow a fully dynamic setup where surface deformation and topography can develop freely. Many authors (including our models in Chapter 2) are now using this free-surface for a more accurate analysis (e.g. [START_REF] Bonnardot | Effect of Margin Curvature on Plate Deformation in a 3-D Numerical Model of Subduction Zones[END_REF][START_REF] Schmeling | A Benchmark Comparison of Spontaneous Subduction Models-Towards a Free Surface[END_REF][START_REF] Zhu | Three-dimensional Dynamics of Hydrous Thermalchemical Plumes in Oceanic Subduction Zones[END_REF][START_REF] Zhu | Influences of the Buoyancy of Partially Molten Rock on 3-D Plume Patterns and Melt Productivity above Retreating Slabs[END_REF]Garel et al., 2014;Ficini et al., 2017;[START_REF] Yoshida | Trench Dynamics: Effects of Dynamically Migrating Trench on Subducting Slab Morphology and Characteristics of Subduction Zones Systems[END_REF]Cerpa et al., 2018). Finally, compositional models can also shed light on upper plate deformation [START_REF] Gerya | Future Directions in Subduction Modeling[END_REF], particularly in identifying the role of UP weakening by fluids (e.g. Arcay et al., 2005) and magmatism (e.g., [START_REF] Ueda | Subduction Initiation by Thermal-Chemical Plumes: Numerical Studies[END_REF][START_REF] Sizova | Subduction Styles in the Precambrian: Insight from Numerical Experiments[END_REF]. [START_REF] Gerya | Geodynamic Regimes of Subduction under an Active Margin: Effects of Rheological Weakening by Fluids and Melts[END_REF] found that aqueous fluids mainly affect the forearc and affect plate coupling, while magmatism weakens the lithosphere below the arc, which controls UP extension and compression. The combination of these could result in retreating subduction and back-arc spreading, due to strong rheological weakening combined with weak plate coupling. These compositional models, however, are more complex to build, costly to run, and are not essential when tackling the dynamic issues investigated in this thesis.

Between the large body of work available and the myriad approaches adopted, it is still not clear what is the relationship between slab pull, trench retreat and upper plate deformation. Several studies have invoked various local physical processes influencing UP deformation, as summarised above; however, their relationships (i.e. identifying causality rather than correlation) and their integration into a global mechanical scenario remain elusive. We, therefore, aim to clarify the physical processes that relate slab pull, trench retreat and upper plate deformation.

Choice of model presented in Chapter 2

As we aim to better understand the fundamental relationship between slab pull, trench retreat and UP deformation, our strategy was to produce the a simple model which shows trench retreat and UP deformation occuring during a deviation from a steadystate subduction. We, thus, construct fully dynamic, thermo-mechanical models in 2-D.

We only vary slab pull by introducing a density anomaly in the slab, and observe the resulting changes in trench retreat and UP deformation. We will focus on the Aegean region, where African slab is subducting beneath Eurasia.

The subduction stacked north-dipping nappes, which are now undergoing back-arc extension due to trench retreat [START_REF] Royden | The Tectonic Expression Slab Pull at Continental Convergent Boundaries[END_REF][START_REF] Jolivet | Mediterranean Extension and the Africa-Eurasia Collision[END_REF]Faccenna et al. 2003Faccenna et al. , 2014;;Brun and Faccenna 2008). The trench is currently retreating approximately 3 cm/yr (e.g. Shaw and Jackson, 2010), however, the migration of the volcanic arc, migration of high pressure metamorphism, and slow down of plate convergence (Ring et al., 2010 and references therein) suggest that trench retreat has been accelerating since it first started in the Eocene (Brun et al., 2016). Consequently, the style of deformation in the UP has also evolved through time.

The Aegean shows three main stages of rollback related deformation, as shown by Jolivet and Brun (2010) and [START_REF] Philippon | Deciphering Subduction from Exhumation in the Segmented Cycladic Blueschist Unit (Central Aegean, Greece)[END_REF] The rollback is thought to have caused up to 600 km of extension in the Aegean upper plate to date (Jolivet and Brun 2010;Jolivet et al. 2013). Various studies suggest that this rollback has accelerated over time, either due to the slab flattening on the lower mantle (e.g. Ring et al., 2010) or due to slab tearing (e.g. Wortel and Spakman, 2000;Brun and Sokoutis, 2010;Jolivet et al., 2013). As a result of this acceleration, the style of extensional deformation is thought to have changed from localised exhumation to distributed segmentation (e.g. Figure 1.3.6, Brun et al., 2016). The exhumation stage has already been heavily studied (e.g. Lister, 1992;Jolivet et al., 1996;Avigad et al., 1997;Krohe and Mposkos 2002;Kumerics et al., 2005;Brichau et al., 2007;Brun and Facenna, 2008, etc). The segmentation stage, however, has received little attention so far. We therefore focus on the segmentation stage and how the fault pattern relates to slab rollback and trench retreat.

Recent deformation in the Aegean upper plate

The Aegean upper plate currently shows pervasive high-angle normal faulting. Early efforts by Mascle and Martin (1990) mapped the subsurface structure of the Aegean's offshore faults using seismic reflection profiles (example shown in Figure 1.3.7). They showed that the Aegean is undergoing classical back-arc rifting. However, Anatolia is also extruding westwards [START_REF] Mcclusky | Global Positioning System Constraints on Plate Kinematics and Dynamics in the Eastern Mediterranean and Caucasus[END_REF][START_REF] Hollenstein | Crustal Motion and Deformation in Greece from a Decade of GPS Measurements, 1993-2003[END_REF][START_REF] Müller | Velocity and Deformation Fields in the North Aegean Domain, Greece, and Implications for Fault Kinematics, Derived from GPS Data 1993-2009[END_REF], which overprints the Aegean with a strike-slip regime, as shown in Figure 1.3.8.

Earthquake focal solutions show that there are NE-trending, dextral strike slip faults in the Central and Eastern Aegean (e.g. Taymaz et al., 1991) as shown in Figure 1.3.9. Representative of the entire Aegean Sea. From Mascle and Martin (1990). When this interaction started is not clear, but recent studies have suggested that it is as old as the Mid-Miocene (e.g. [START_REF] Philippon | Tectonophysics The Interaction between Aegean Back-Arc Extension and Anatolia Escape since Middle Miocene[END_REF]. The oldest extension-related sediments are known to be Mid-Miocene (e.g. Mascle and Martin, 1990;Sánchez-Gómez et al., 2002;Beniest et al., 2016), but the age of dispersed strike slip faults is more elusive. The NAF is known to have localised recently, dated around 5

Ma (e.g. [START_REF] Armijo | Westward Propagation of the North Anatolian Fault into the Northern Aegean: Timing and Kinematics[END_REF]Şengör et al., 2005), however, large-scale dextral shear zones in the Aegean may have initiated in the Mid-Miocene [START_REF] Armijo | Quaternary Evolution of the Corinth Rift and Its Implications for the Late Cenozoic Evolutin of the Aegean[END_REF][START_REF] Armijo | Westward Propagation of the North Anatolian Fault into the Northern Aegean: Timing and Kinematics[END_REF]Şengör et al., 2005;[START_REF] Reilinger | Geodetic Constraints on the Tectonic Evolution of the Aegean Region and Strain Accumulation along the Hellenic Subduction Zone[END_REF]. Additionally, there is evidence of NE-trending strike slip faults partly localising magmatism since the Mid-Miocene (Kokkalas and Aydin, 2013). Therefore, there is evidence for strike slip activity as far back as the Mid-Miocene, but whether or not it influenced back arc extension remains unclear.

Focus on the Cyclades and workflow presented in Chapter 3

We attempt to clarify the relationship between slab rollback, extension and strike slip faulting in the central Aegean, and how this relationship has changed through time. To do this, we focus on the Cyclades, where the geology is well known, but the latest segmentation related structures are poorly constrained. The Cyclades block also contains the limit between clockwise and counter-clockwise block rotation (Figure .3.11), indicated from palaeomagnetic studies (Morris and Andersen, 1996;Avigad et al., 1998). This is makes it particularly interesting to understand the role of rotation in a retreating subduction system. Additionally, the Cyclades block is one of the least active areas in the Aegean, making it a suitable candidate for distinguishing older structures.

1
We therefore start by addressing the question: why do the normal faults in the Aegean have variable strike? This was already evident in the early structural maps of Mascle and Martin (1990), and can be seen in Figure 1.3.10. We distinguish multiple normal fault generations and interpret their formation in the context of a rotating block. This is reconciled with slab rollback and extrusion. As we aim to better understand the large-scale pattern of recent high angle normal faults (i.e. segmentation), we will only focus on structures that post-date low-angle detachments. We find that combining field and seismic data to be particularly useful in interpreting faulting pattern and kinematics. We use this to produce a regional interpretation of progression of high angle faults with rotation related to trench retreat and Anatolian extrusion. have not yet conceived of, but today, our main deep energy resources are fossil fuels and geothermal energy. Buried hundreds to thousands of meters beneath the surface, these resources are strongly affected by tectonic activity, and lie where direct observation is not possible. As the petroleum industry matured, it developed various subsurface observation tools on multiple scales, including gravity, magnetic, and seismic surveys, as well various logging tools. Academia has also made large strides in processing seismic tomography data in recent years. These combined with the present availability of satellite data provides us with a suite of tools from which the petrology, temperature and fluid distribution in the lithosphere can be studied. The petroleum industry has made extensive use of this and been able to go as far as exploiting reservoirs in unconventional settings. Consequently, the petroleum industry has produced a wealth of knowledge on geological settings and subsurface processes, particularly related to sedimentary basins, subsurface fluid flow, structural geology and geomechanics. However, as the world's fossil fuel reservoirs deplete and awareness of their environmental impact rises, we must replace energy production from fossil fuels by more sustainable resources. This brings us to our second deep energy resource, geothermal energy.

What is geothermal energy?

Unlike fossil fuels, geothermal energy production is less mature and less predictable, as it requires a better understanding of the heat transport within the lithosphere. It is considered to be a renewable energy source, which exploits the Earth's internal natural heat. This is predominantly primordial heat from the Earth's accreation 4.5 Ga ago, in addition to radiogenic heat generated in the lithosphere and within the Earth's mantle [START_REF] Stacey | Physics of the Earth[END_REF][START_REF] Manzella | General Introduction to Geothermal Energy BT -Geothermal Energy and Society[END_REF]. As the solid Earth is cooling from its surface, the crust is generally coolest and temperature increases with depth. This vertical temperature profile is called the geotherm. In the simplest assumption, the continental lithosphere is assumed to be in thermal equilibrium, i.e. the geotherm is in steady state (Cloetingh et al., 2010). In geologically stable continental regions the ambient heat flux averages around 0.065 W/m 2 [START_REF] Stacey | Physics of the Earth[END_REF] and the geotherm averages around 20-30 °C/km [START_REF] Manzella | General Introduction to Geothermal Energy BT -Geothermal Energy and Society[END_REF]. This heat profile is too cold to be exploited for energy production, which requires a higher geotherm. The latter can result from thermal influences, such as magmatism or tectonic activity changing the thickness of the lithosphere, or inducing metamorphism. The highest geotherms are, thus, observed in volcanic areas and are between 40-80 °C/km (Arndt, 2011). Geothermal energy is already exploited from steam in such volcanic areas, e.g. Iceland, Italy, California and New Zealand [START_REF] Stacey | Physics of the Earth[END_REF]. Tectonic processes such as extension or magmatism also raise local geotherms, such as in the Pannonian Basin, the Aegean and parts of western Turkey (Cloetingh et al., 2010). Interestingly, tectonic processes can also raise local geotherms in amagmatic systems by creating fracture networks.

Such permeable networks allow superheated fluids to migrate from significant depths to the upper crust. An example of this is seen in Taiwan where a geothermal reservoir has formed in fractured metamorphic rocks as a result of oblique collision, creating local hot springs (Gup et al., 2018)

Geothermal energy contribution

On average, global energy use increases by 3% a year [START_REF] Stacey | Physics of the Earth[END_REF]. In 2015, geothermal energy contributed a projected 73,549 GWh to global energy use, a 16% rise from the previous 5 years [START_REF] Bertani | Geothermal Power Generation in the World 2010-2014 Update Report[END_REF]. Figure 1 shows the rise in geothermal energy production from 1950 to a 2020 projection. [START_REF] Bertani | Geothermal Power Generation in the World 2010-2014 Update Report[END_REF] showed that with 16600 GWh, the USA produced the largest amount of geothermal energy 2015, followed by the Philippines, Indonesia, and New Zealand, who produce between 5000-10000 GWh each. Turkey, Kenya, Japan, Costa Rica and El Salvador each produced between 1000-5000 GWh in 2015. The prominence of subduction zones in these values is noteworthy (Philippines, Indonesia, New Zealand, Japan, Costa Rica and El Salvador). projected to 2020 [START_REF] Bertani | Geothermal Power Generation in the World 2010-2014 Update Report[END_REF].

Geothermal energy is used for both heating and producing electricity. In general, a wide range of geothermal reservoir temperatures can be used for heating, but producing electricity usually requires temperatures over 100°C [START_REF] Manzella | General Introduction to Geothermal Energy BT -Geothermal Energy and Society[END_REF]. Presently, global geothermal power production has a capacity of 12.7 GW (IRENA 2018), and annual electricity generation in 2015 reached 80.9 TWh, approximately 0.3% of global electricity generation (IRENA 2017) and 18% of total geothermal energy use [START_REF] Manzella | General Introduction to Geothermal Energy BT -Geothermal Energy and Society[END_REF]. With growing capacity, geothermal energy can provide a more widespread energy source since production is independent of seasonal variations, providing a continuous energy base-load which can respond to energy demand [START_REF] Manzella | General Introduction to Geothermal Energy BT -Geothermal Energy and Society[END_REF].

Geothermal technologies

Geothermal plants use heat directly, or most commonly, extract hot geothermal fluid, circulate it through a heat exchanger, and re-inject it into the subsurface. Geothermal plants must, thus, optimise their locations, manage borehole extraction and injection, and manage the heat exchanger. For heating, temperatures tapped usually range from 70-90°C [START_REF] Manzella | General Introduction to Geothermal Energy BT -Geothermal Energy and Society[END_REF]. Higher temperatures > 100°C are required for electricity production, and the temperature conditions affects the efficiency of electricity generation, as higher temperatures can produce more power for the same flow rates (e.g. [START_REF] Dipippo | Ideal Thermal Efficiency for Geothermal Binary Plants[END_REF], Cloetingh et al., 2010). These usually require deeper drilling (around 5 km for a mean continental heat flow and a thermal gradient of 20°C/km).

Shallower reservoirs are mainly used for heating and cooling. The Enhanced Geothermal System (EGS) EGS exploits deep geothermal reservoirs which are normally too impermeable for commercial production by employing hydraulic and chemical stimulation. EGS involves drilling at least two boreholes (a 'doublet'), one used for production, the other for injection, where both wells are stimulated to artificially enhance the permeability of the reservoir between them, as shown in Figure 1.4.3 (Gérard et al., 2006). The potential for EGS was proven in the Soultz research project in eastern France [START_REF] Gérard | The Deep EGS (Enhanced Geothermal System) Project at Soultz-Sous-Forêts (Alsace, France)[END_REF] and commercially in the Landau project, Palatinate, western Germany (Cloetingh et al., 2010).

Challenges and predictability

The first challenge of geothermal energy use is to accurately predict the geotherm. To reduce drilling costs, relatively high geotherms are targeted. These can be a result of magmatism, e.g. Iceland (Cloetingh et al., 2010); hydrothermal fluids circulating from deeper to shallower levels through open faults and fracture, e.g. Soultz, France [START_REF] Gérard | The Deep EGS (Enhanced Geothermal System) Project at Soultz-Sous-Forêts (Alsace, France)[END_REF]; or recent crustal thinning, e.g. the Pannonian basin, Hungary (Horvath et al., 2015). Predicting geothermal potential, thus, requires an understanding not only of the thermal state of a region, but also its tectonic structures. This is not a trivial task, as the present day thermal and structural conditions of the lithosphere requires inerpretation from surface observations (e.g. geological cross-sections, seismic profiles). Many of these processes also have different effects, for example, basin formation through crustal thinning can raise heat flow, as in the Pannonian Basin, however sedimentation can significantly reduce surface heat flow due to the blanketing effect (e.g. [START_REF] Wangen | The Blanketing Effect in Sedimentary Basins[END_REF]. While this is partially compensated by radiogenic heat production in sediments ( Van Wees et al., 2009), erosion and thinning also reduce heat flow due to the loss of radiogenic heat from crustal material (e.g. [START_REF] Wees | Lithosphere Rheology during Intraplate Basin Extension and Inversion: Inferences from Automated Modeling of Four Basins in Western Europe[END_REF]. The latter two effects will create different geotherms, as sediment blanketing may reduce surface heat flow, but preserve it at depth. In fact, heat flow modelling of continental extensional basins shows that over long time-scales, heat flow is lowered compared to unrifted basins ( Van Wees et al., 2009) and is confirmed by lower heat flow in the deepest centres of certain basins compared with their margins, e.g. in the Black Sea [START_REF] Cloetingh | Thermo-Mechanical Modelling of Black Sea Basin (de) Formation[END_REF]. Thus, taking the example of basins, geotherms and heat flow will vastly differ depending on the age of rifting, the stretching factors and the sedimentation rates. Other geological settings similarly have competing thermal influences. Incorporating these various factors is facilitated by modelling, however, it is a challenge to collectively model both the thermal state, the kinematic evolution and the development of smaller scale influences such as fractures.

Frontier technologies, such as EGS, may be able to exploit a wider range of reservoirs, but also entail further technical challenges. Cloetingh et al. (2010) show that fluid stimulation requires a clear understanding of the reservoir's stress state to optimise production and minimise induced seismicity. This is because inducing fractures in the subsurface preferentially creates open fractures in the direction of maximum stress. If critically stressed fractures exist in this direction, they will also open, enhancing the permeability, however, the presence of critically stressed faults may induce seismicity (e.g. [START_REF] Charléty | Large Earthquakes during Hydraulic Stimulations at the Geothermal Site of Soultz-Sous-Forêts[END_REF]. This is not a new problem, and has been tackled by the petroleum industry for decades, however, the geothermal industry does not yet benefit from the same amount of data for more complete geomechanical analyses.

Interestingly, this is partly addressed by the petroleum industry's data becoming available, but the sub-surface of continental regions remains vastly unexplored. Some factors remain poorly understood, such as the lithosphere layered rheology or stress state. In such situations, tectonic studies could largely aid in geothermal exploration efforts.

The role of tectonic modeling

Tectonic models can address several challenges (see above) and help evaluate a region's potential for geothermal energy production. This was already shown by Cloetingh et al. (2010) who focus on EGS, a geothermal technology that could disseminate geothermal energy to non-volcanic regions (see above). Tectonic modelling could improve conventional geothermal as well as unconventional technologies in a variety of ways, of which a non-exhaustive list is presented below.

Location optimization

This is one the largest controlling factors of a project's success, as it not only determines which temperatures are tapped but also production risk. This is because producing geothermal energy does not only require high temperatures, but also needs particular sub-surface permeabilities and stress conditions. Tectonic modelling, can thus, first aid by identifying regions of high temperature. These are not limited to temperatures elevated by magmatism, but include temperature highs which result from tectonic processes, such as rifting. Many conventional geothermal reservoirs are produceable because they contain hot fluids which are brought to drillable depths via permeable fault and fracture networks (e.g. Fernández et al., 1990;[START_REF] Faulds | Characterizing Structural Controls of Geothermal Fields in the Northwestern Great Basin: A Progress Report[END_REF][START_REF] Sanjuan | Tracer Testing of the Geothermal Heat Exchanger at Soultz-Sous-Forêts (France) between 2000 and 2005[END_REF], e.g. the Upper Rhine Graben reservoir [START_REF] Vidal | Overview of Naturally Permeable Fractured Reservoirs in the Central and Southern Upper Rhine Graben: Insights from Geothermal Wells[END_REF].

Such reservoirs often carry a greater location risk, as the permeable network comprises discrete produceable planes, which if missed, will yield a dry hole. Tectonic modelling here is crucial to better understand the structurally permeable network and its extent.

Reduced risk drilling and production

Understanding regional stresses and lithological failure criteria are crucial to safely drill a borehole. This is the field classically referred to as geomechanics, which uses tectonic information to design borehole drilling plans. These include drilling directions and mud (drilling fluid) weight windows which will keep the borehole stable (the mud must be heavy enough to stop the borehole from caving in, but not so heavy that it fractures the formations). No matter how well the borehole is designed, results can be disastrous if the borehole crosses a critically stressed fault, as this could induce slip on the fault plane, creating seismicity or collapsing the borehole. Similarly, reservoir production -particularly where stimulation is involved-can induce seismicity on critically stressed faults, as was the case in Basel 2006 EGS project, where stimulation caused several seismic events during stimulation, that culminated in a ML 3.4 earthquake after stimulation was stopped (Häring et al., 2008). Thus, understanding the local stress state, particularly in the presence of faults, is critical for safe production.

Reduced cost production

Understanding the stress state can not only improve drilling operations, but also largely reduce the cost of production, especially where stimulation is involved. While stimulation in the vicinity of a critically stressed fault may be dangerous (e.g. [START_REF] Mignan | Including Seismic Risk Mitigation Measures into the Levelized Cost Of Electricity in Enhanced Geothermal Systems for Optimal Siting[END_REF], triggering critically stressed fractures is favourable, as these create a more extensive permeable fracture network. However, if these fractures do not cross the reservoir, stimulation could be useless, since critically stressed fractures form preferential pathways for hydraulic stimulation (Häring et al., 2008). Thus, before hydraulic stimulation, tectonic modelling can be a large aid in predicting the success of stimulation, by first understanding the existing fault and fracture networks, and their stress states. Additionally, a better understanding of the stress state allows a better choice of stimulation direction (e.g. along the maximum stress), which will render stimulation more efficient, reducing its cost.

Prolonged reservoir life

Similarly to reducing stimulation cost, more efficient stimulation will also prolong the lifetime of production, since steadier pressures can be maintained. This also holds true for conventional geothermal technologies, which do not require stimulation. A geothermal reservoir fed by a naturally permeable fracture network will still tend to have higher permeability in a certain direction. Tapping fluid flow in this direction will ensure longer-lived production, as the pressure would quickly drop in other directions.

Additionally, where geothermal fluids flow along open faults, recent structural studies have shown that fluid flow does not necessarily coincide with the positions of maximum slip, but seems to correlate more with fault zone intersections, overlaps or terminations (e.g. [START_REF] Faulds | Structural Controls on Geothermal Systems in Western Turkey: A Preliminary Report[END_REF]. Thus, targeting these fault segments will ensure longer production.

Future perspectives

Conventional geothermal is presently limited to volcanic regions, regions with natural hot springs or permeable deep geothermal reservoirs. However, many parts of the world have sufficiently high geotherms to utilise EGS even in the absence of the aforementioned conditions. For example, EGS is estimated likely to provide 5-10 % of the USA's electricity demand in 2050 (e.g. [START_REF] Tester | Impact of Enhanced Geothermal Systems on US Energy Supply in the Twenty-First Century[END_REF]. Focussing on Europe, Figure 1.4.4 shows temperature gradient values in Europe, (extracted from the international heat flow database http://www.heatflow.und.edu, figure taken from Cloetingh et al., 2010), identifying several regions with high enough geotherms for geothermal energy production, particularly around Europe's mountain belts. These combined with tectonic knowledge (geology, stress state and active seismicity) can be used to identify areas which are critically stressed (Cloetingh et al., 2010). Additionally, onshore first order heat flow appears to relate to Cenozoic tectonics at various depth levels in the lithosphere (Cloetingh et al., 2010). These can be modelled, allowing more accurate prediction of resources. Medium-high enthalpy geothermal is thought to be possible in various areas in Europe, by exploiting active faults that provide fluid pathways from depth (e.g. Fernández et al., 1990;[START_REF] Gartrell | A New Model for Assessing Trap Integrity and Oil Preservation Risks Associated with Postrift Fault Reactivation in the Timor Sea[END_REF]. These are thought to be most promising in extensional regions, where faults and fractures are more likely to be open (e.g. [START_REF] Zoback | Reservoir Geomechanics[END_REF]. This renders the Mediterranean region particularly interesting for geothermal prospecting, as it is marked by active deformation, local magmatism and pronounced lithospheric block movement (e.g. [START_REF] Cloetingh | Lithospheric Memory, State of Stress and Rheology: Neotectonic Controls on Europe's Intraplate Continental Topography[END_REF]. and no electricity is produced from geothermal [START_REF] Mendrinos | Exploring for Geothermal Resources in Greece[END_REF]. This could be addressed by deep conventional geothermal as well as EGS, depending on the reservoir, however, the reservoirs must first be better understood to select the most suitable sites and technology. We show in Chapter 4 that a better understanding of geothermal reservoirs can benefit from incorporating tectonic concepts. While this chapter does not aim to produce realistic geothermal models of the Aegean (that would be an entire PhD in itself), it shows that incorporating tectonic modelling into the thermal modelling workflow can largely change resulting predictions. From this, we conclude that tectonic models and concepts should be systematically included in thermal modelling workflows. This calls for further development of said workflows, but it is evident that this could add large potential to geothermal exploration, reduce its costs, and increase its production. plate deformation partly depends on the upper plate's ability to translate horizontally: an upper plate with a ridge at its trailing edge deforms by shortening, while a fixed upper plate deforms by extension. Finally, in some cases, we observe that upper plate deformation permits trench retreat if the upper plate is weak enough to be deformed by basal drag from underlying asthenospheric flow. Our results provide insights into retreating subduction systems with contrasting upper plate deformation modes, such as the compressive Andes and the extensional Aegean.

Upper plate deformation

Introduction

The slab pull force, generated by the negative buoyancy of the subducting lithosphere relative to the surrounding mantle, is a key component of the force-balance governing tectonic plate motions (e.g. [START_REF] Forsyth | On the Relative Importance of the Driving Forces of Plate Motion[END_REF]. There are tens of subduction zones on Earth, and observations show that most of them have mobile trenches (e.g. Chase, 1978;[START_REF] Heki | The Amurian Plate Motion and Current Plate Kinematics in Eastern Asia[END_REF]. While some trenches, such as the Mariana, advance towards the upper plate, others, including Tonga, retreat into the subducting plate (Heuret and Lallemand, 2005).

In the past, trench migration has been considered to cause deformation in the upper plate (UP hereafter). It is presumed that trench retreat creates extensional stresses in the UP, forming back-arc basins, while trench advance pushes against the UP, leading to UP compression (e.g. [START_REF] Uyeda | Back-Arc Opening and the Mode of Subduction[END_REF]. However, Heuret and Lallemand, 2005 have shown that this correlation is not manifest on Earth. Furthermore, Heuret and Lallemand (2005) and [START_REF] Arcay | Back-Arc Strain in Subduction Zones: Statistical Observations versus Numerical Modeling[END_REF] demonstrated that the bulk upper plate velocity away from the trench (V up ) plays a role in this deformation. For example, trench retreat coupled with back-arc extension is observed in the Mediterranean, New Hebrides, Ryukyu and Scotia; while trench retreat is coupled with UP compression in Chile, Japan, Costa Rica and Manila (Heuret and Lallemand, 2005). This upper plate deformation is thought to result from a complex force balance involving the motions of the sinking plate, surface plates and their respective underlying (asthenospheric) mantle flows (Heuret and Lallemand, 2005).

In a subduction system, the trench and plate motions, as well as the asthenospheric flows, can result from a variety of physical processes. These include buoyancy variation in the subducting plate (e.g. [START_REF] Royden | Subduction Zone Geodynamics[END_REF], deformation through piling and folding of the slab at depth (e.g. [START_REF] Capitanio | Upper Plate Controls on Deep Subduction, Trench Migrations and Deformations at Convergent Margins[END_REF][START_REF] Holt | Trench Migration and Overriding Plate Stress in Dynamic Subduction Models[END_REF], shear drag by asthenospheric flow at the base of the upper plate (e.g. [START_REF] Nakakuki | Dynamics of Slab Rollback and Induced Back-Arc Basin Formation[END_REF], variation in UP thermal and mechanical properties (e.g. [START_REF] Rodríguez-González | The Role of the Overriding Plate Thermal State on Slab Dip Variability and on the Occurrence of Flat Subduction[END_REF], lateral mantle flow, or "mantle wind" (e.g. Ficini et al., 2017), 3-D effects from side plates (e.g. [START_REF] Yamato | Influence of Surrounding Plates on 3D Subduction Dynamics[END_REF] and toroidal mantle flow (e.g. Schellart and Moresi, 2013). In our study, we have designed a simple, 2-D model to explore the physical interactions between some of the dominant controls on trench retreat and upper plate deformation.

Trench retreat is often assumed to be primarily driven by slab pull [START_REF] Forsyth | On the Relative Importance of the Driving Forces of Plate Motion[END_REF], although some recent studies have challenged this notion (e.g. [START_REF] Stotz | Pressure-Driven Poiseuille Flow: A Major Component of the Torque-Balance Governing Pacific Plate Motion[END_REF].

Slab pull exerts a bending torque at the subduction hinge, contributing to a downward folding of the slab. This is thought to generate spontaneous retreat of the trench -or rollback, suggesting that older (and therefore, colder and denser) slabs should have retreating trenches (e.g. Garfunkel et al., 1986;[START_REF] Conrad | Effects of Plate Bending and Fault Strength at Subduction Zones on Plate Dynamics[END_REF]Lallemand et al., 2008). However, in nature, such a relationship is not observed (Heuret and Lallemand, 2005). To better understand this, we model a slab with varying slab pull through time. Varying slab pull and slab sinking rates can be caused by density anomalies within the sinking slab, for example, positive anomalies from eclogitised crust [START_REF] Liu | The Role of Oceanic Plateau Subduction in the Laramide Orogeny[END_REF][START_REF] Arrial | Influence of Geometry and Eclogitization on Oceanic Plateau Subduction[END_REF], or negative anomalies caused by ridges, seamounts or continental crust (e.g. [START_REF] Royden | Subduction Zone Geodynamics[END_REF][START_REF] Magni | How Collision Triggers Backarc Extension: Insight into Mediterranean Style of Extension from 3-d Numerical Models[END_REF].

Our numerical study features a transient increase of slab buoyancy, to explore the effect of changing the slab pull force on subduction dynamics, including trench motion, asthenospheric flow and upper plate deformation. We also explore the effect of UP strength and its role in this relationship. Previous modelling studies on UP deformation have often focused on one particular aspect, such as UP motion (Cerpa et al., 2018), large-scale mantle flow [START_REF] Husson | Trench Migration and Upper Plate Strain over a Convecting Mantle[END_REF], or boundary conditions [START_REF] Capitanio | Upper Plate Controls on Deep Subduction, Trench Migrations and Deformations at Convergent Margins[END_REF].

This paper aims to explore the relationship between some of these aspects, by focusing on 2-D models of retreating subduction zones and observing UP deformation with varying rates of trench retreat and varying UP boundary conditions. This is particularly interesting for systems such as the South American subduction zone, where the trench is retreating, yet the upper plate experiences strong shortening, contributing to the formation of the Andes (e.g. [START_REF] Oncken | Deformation of the Central Andean Upper Plate System ---Facts, Fiction, and Constraints for Plateau Models[END_REF].

Model description

We design 2-D, upper mantle, thermo-mechanical models of subduction, with a freely moving trench and upper plate. We start with quasi-steady-state subduction (e.g [START_REF] Capitanio | Dynamic Models of Downgoing Plate-Buoyancy Driven Subduction: Subduction Motions and Energy Dissipation[END_REF]) and introduce a positive density anomaly into the slab to increase slab pull, allowing us to study the temporal relationships between slab pull, trench retreat, asthenospheric flow, and upper plate deformation without directly forcing plate motion. This setup is illustrated in Figure 2.1.

We use the unstructured, adaptive mesh, Fluidity computational modelling framework (e.g. [START_REF] Davies | Fluidity: A Fully Unstructured Anisotropic Adaptive Mesh Computational Modeling Framework for Geodynamics[END_REF][START_REF] Kramer | An Implicit Free Surface Algorithm for Geodynamical Simulations[END_REF][START_REF] Voci | A Systematic 2-D Investigation into the Mantle Wedge's Transient Flow Regime and Thermal Structure: Complexities Arising from a Hydrated Rheology and Thermal Buoyancy[END_REF] to solve the equations describing the conservation of mass, momentum, and energy for an incompressible Stokes fluid, under the Boussinesq approximation:

! ! ! ! ! ! (1) ! ! ! !" ! !!!! ! (2) !! !! ! ! ! ! ! ! ! ! ! ! ! !(3)
where u and g denote velocity and gravity vectors, respectively, σ ij the stress tensor, T the temperature, κ the thermal diffusivity, and Δ⍴=-⍺⍴ s (T-T s ) the density difference due to temperature, with ⍺ the coefficient of thermal expansion and ⍴ s the nominal density at the surface temperature T s (parameter values are provided in Table 1).

We use a wide box of 6,000 km to minimise potential boundary effects. The domain height is 660 km to simulate the upper mantle, with the first-order approximation that the high-viscosity lower mantle is a barrier to slab penetration. The top boundary is a free surface while the bottom and sides are free slip. Plates are defined thermally using the half-space cooling model, and plate ages increase linearly from 0 Myr at the leftmost and rightmost corners ("ridges") to the limit between the subducting and upper plate ("trench"). The trench is initially located at 3500 km from the left boundary. Note that the simulated UP is free to move laterally towards the trench due to the rightmost ridge, which mimics a partially oceanic upper plate. These ridges also enable asthenopheric mantle material to transform into lithospheric plate material by thermal diffusion from the cold surface. Hence, as in nature, the plates renew and do not have a finite length.

Initial plate thermal ages are 40 Ma for the subducting plate (left) and 20 Ma for the upper plate (right). These ages are chosen to generate a slab deflecting horizontally at the bottom of the upper mantle (see Garel et al. (2014) for a regime diagram of slab morphology as a function of initial plate ages). The subducting plate is prescribed a proto-slab shape that is sufficient to initiate subduction under its negative buoyancy. The numerical simulation evolves for 35 Myr, after which subduction dynamics reach a quasisteady-state, i.e. plate and trench velocities exhibit only small variations through time.

We then restart our simulations from this thermal state of mature subduction, where the slab has already subducted through the upper mantle and has partly flattened on the bottom boundary at 660 km depth (Figure 2.1). Hereafter, the time indicated in the simulations is the time passed since this initial mature state.

Most of the domain material is assigned a "mantle" composite rheology, and only an 8 km thick, weak, decoupling layer at the surface of the subducting plate is assigned a weak "crustal" rheology. The weak layer properties revert to those of normal material below 200 km depth. The "mantle" rheology takes into account 4 mechanisms: pseudobrittle yielding, dislocation creep, diffusion creep, and low-temperature plasticity or "Peierls" creep. Strain-rates for these mechanisms add up, resulting in a composite bulk viscosity:

! ! ! ! !"## ! ! ! !"#$ ! ! ! ! ! ! ! ! !! (4)
with creep viscosities μ diff , μ disl , μ p calculated as 77

! !"## !"#$ ! ! ! !! ! !"# !!!" !"! ! ! !! !!! ! (5)
with A a prefactor, n the stress exponent, E, V the activation and volume energies respectively, P=⍴ s gz the lithostatic pressure, R the gas constant, and the second invariant of the strain-rate tensor. Tr is the temperature obtained by adding an adiabatic gradient of 0.5 K/km to the Boussinesq solution [START_REF] Fowler | The Solid Earth: An Introduction to Global Geophysics[END_REF].

Brittle failure at low lithostatic pressure is approximated through a "yielding" viscosity

! ! ! ! ! ! ! !! !! (6)
where the yield strength τ y is given by

! ! ! !"#! ! ! ! ! ! !! ! !"#$ (7)
with τ 0 the surface yield strength, f c the friction coefficient, and τ ymax the maximum yield strength.

The bulk viscosity calculated in equation ( 4) is limited by lower-and upper-bound values of 10 18 and 10 25 Pa.s, respectively.

Key model parameters are summarised in Table 1. Our geometrical set-up, boundary conditions, rheology parameterization and numerical solution are similar to Garel et al. (2014).

To vary the slab pull force through time, we add a 1000x50 km 2 denser block near the trench in the subducting plate, below the weak decoupling crust (Figure 2.1). This block triggers a change in the overall system dynamics through a transient increase in slab pull, without requiring a prescribed force at the model surface, thus avoiding artificial UP deformation. We vary the block density anomaly from +10 to +100 kg/m 3 in our models to achieve different magnitudes of slab pull and trench retreat rates (Table 2).

Table 2.1: Model parameters

We sometimes modify the upper plate (section 3.5) for a chosen density anomaly (+70 kg/m 3 ) in two increments:

1-To test the influence of the ridge located at the trailing edge of the UP (topright corner in Figure 2.1), we remove the ridge and impose a constant thickness in the UP from the trench to the top-right corner of the model. The thickness used is the average thickness of the ridge-edge models, with an age of ~45 Ma, corresponding to a thickness of ~67 km (calculated via the depth of the 1300 K isotherm derived from the half-space cooling model).

2-To test the influence of upper plate strength in addition to removing the ridge, we reduce the thickness of the upper plate by setting a constant UP age of 15 Ma, which corresponds to a thickness of ~38 km (calculated in the same away as above).

In addition to asthenosphere flow pattern, the main surface velocities analysed are the velocity of the subducting plate (V sp ), upper plate (V up ), and trench (V tr ). The trench is defined as the surface location where the weak, decoupling layer meets the upper plate.

Locating the trench using a change in horizontal velocity direction yields similar, although noisier, results. The trench velocity is then calculated as a smoothed timederivative of the trench location. The range of simulated trench velocities is given in Table 2. As a consequence of the model geometry and closed (free-slip) side boundaries, the subducting plate (hereafter SP) advances towards the right while the upper plate and the trench move towards the left. This creates two asthenosphere flow cells: one beneath the subducting slab (clockwise, hereafter referred to as sub-slab), and another in the mantle wedge between the upper plate and the slab (counter-clockwise). These are largely separated by the slab, except for a small, connecting channel beneath the slab (cf. Figure 2.

Results

Steady state reference simulation

2).

The reference simulation exhibits minimal changes in the asthenosphere flow pattern through time, and is characterised by a slow, steady trench retreat (Figure 2.2a). This is quasi-steady-state subduction, which continues unchanged from the initial condition. This is reflected in the plate and trench surface velocities, which are shown in Figure 2.2.3. Here, the subducting plate moves rightwards with a velocity (V sp ) of ca. 4 cm/yr, while the upper plate (V up ) and trench (V tr ) both move leftwards with a velocity of ca. 1 cm/yr throughout the simulation. kg/m 3 ). Here, the SP moves 5 times faster than the reference at the start of the simulation, and peaks at 17 cm/yr, which is 17 times faster than the reference (Figure 2.3). Therefore, the slab sinks faster (see Appendix A), and the slab tip advances further (c.f. Figure 2.2 time t 3 , when the slab tip in simulation DA70 is 150 km further advanced than the slab tip in the reference simulation). Faster subduction does not alter the geometry asthenospheric flow, but it does increase velocity magnitudes (Figure 2.2). This faster flow results in increased lateral motions of both plates through the ridges on both domain sides (also seen in Faccenna et al., 2017). Thus, the UP and the trench also accelerate in response to the accelerated sinking slab, but unlike the SP, they do not accelerate throughout the simulation. V tr accelerates during the first 3 Myr, peaks around t 2 3.4 Myr at 4.2 cm/yr and subsequently decelerates. V up accelerates during the first 4

Introducing a density anomaly in the slab

Myr, peaks at 5 cm/yr and then decelerates.

The trench decelerates before the upper plate (Figure 2.3 after 3 Myr: V up > V tr ), resulting in a deformed, high strain-rate region under compressive stresses in the UP, as illustrated in Figure 2.4. Deformation is largest when the discrepancy between V up and V tr increases around t 4 6 Myr (Figure 2.3). Figure 2.4 shows the velocity field, strain-rate and horizontal stress of the reference simulation and simulation DA70 (block density anomaly of +70 kg/m 3 ) during trench deceleration at t 3 4.8 Myr. Figure 4.a.ii confirms that introducing a density anomaly only changes the velocity field magnitude, not its pattern (also evident in Figure 2.2). At time t 3 , the maximum velocity (in the asthenosphere near the dipping slab) is raised from 5.7 cm/yr in the reference to 14.2 cm/yr in DA70. The two simulations in Figure 2.4.a exhibit a constant velocity along most of the two plates, i.e. the plates translate rigidly. Figure 2.4.b shows that in both simulations, strainrates are highest in the asthenosphere, particularly at the ridge interfaces due to change in flow direction, and between the plates and the asthenosphere due to viscous shear stresses.

We also observe in both simulations a narrow area in the UP immediately next to the trench with large strain-rates, resulting from the bending of the sinking slab at depth [START_REF] Capitanio | Dynamic Models of Downgoing Plate-Buoyancy Driven Subduction: Subduction Motions and Energy Dissipation[END_REF]) transmitted to the UP through the subduction interface.

In simulation DA70, an intraplate lateral velocity gradient near the trench (Figure 2.4.a.ii) is associated with additional UP deformation (Figure 2.4.b.ii). We track the UP region where the strain-rate at 30 km depth is higher than a threshold of 10 -16 s -1 to quantify deformation width in the UP. Figure 2.4.b.ii shows that the deformed region remains smaller than 225 km wide in the reference simulation, while in DA70 it is 576 km at t 3 4.8 Myr (Table 2).

The mode of deformation associated with these high strain-rates is compressive horizontal stress (τ xx ) in both simulations (Figure 2.4.c). As with strain-rate, simulation DA70 exhibits more intense and widespread horizontal compression in the UP than the reference simulation.

Varying the density anomaly in the slab

To better understand the effect of increased slab pull, we varied the block density anomaly in increments of 10 kg/m 3 between 0 and +100 kg/m 3 . A summary of the results is presented in Table 2, which shows for each density anomaly: the relative increase in slab pull, calculated by multiplying the density anomaly by the 2-D block dimensions and gravity; peak V tr and its associated time; the peak discrepancy between V up and V tr and its associated time; max UP deformation width (where strain rate is higher than 10 -16 s -1 at 30 km depth); the time range during which the UP experiences deformation wider than the reference's (i.e. where the deformed region defined by strain rate higher than 10 -16 s -1 at 30 km depth is wider than 225 km). Appendix B also shows graphs of each density anomaly against its maximum V tr , each density anomaly against its maximum width of upper plate deformation and the maximum width of upper plate deformation against the maximum discrepancy between V up and V tr . Table 2.2: List of simulations and their block density anomalies: additional slab pull (calculated by multiplying density anomaly by 2-D block dimensions and gravity); peak trench retreat (V tr ) and associated time; the peak discrepancy between upper plate and trench retreat velocities (V up -V tr ) and associated time; maximum upper plate deformation width (where strain rate is higher than 10 -16 s -1 at 30 km depth), and the time range during which the UP deformation width is higher than the reference's. The same trench retreat pattern observed in simulation DA70 (acceleration then deceleration) is observed for simulation DA40 (Figure 2.5), and also for all density anomalies higher than +20 kg/m 3 . Note that the velocities for DA70 oscillate after 6 Myr due to slab buckling as the denser block reaches the closed bottom boundary (see Appendix C). We stop our analysis of UP deformation before slab buckling starts. b) Plots of the second invariant of the strain-rate tensor for each simulation at their times of widest upper plate deformation (where strain rate is higher than 10 -16 s -1 at 30 km depth). Magenta lines delimit the deformed region. The plates are outlined with the 1300 K isotherm in black, while the denser block is outlined in white.

Simulation DA10 shows equal V up and V tr , which is 0.5 cm/yr faster than the reference. In DA40, there is a discrepancy between V up and V tr , where V up becomes progressively faster than V tr . For all simulations, trench and UP velocities both increase with increasing density anomaly (Table 2). Peak V tr increases almost linearly with increasing slab pull.

The discrepancy between V up and V tr increases with higher density anomalies, and is associated with stronger and wider UP deformation (see Appendix B and high strain-rate regions in Figure 2.5.b). The plates are outlined in black by the 1300 K contour. The denser block is outlined in white.

Fixing the edge of the upper plate

The deformation observed in our models is strongly dependent on the UP velocity. We now explore the effect of the presence of a ridge at the edge of the UP (top right corner, Figure 2.1) on the asthenosphere return flow and UP deformation. Thus, we design models with an UP of constant thickness and age from the near-trench region to the right side, preventing the rigid translation of the UP by lateral asthenospheric push at the ridge. Note that we preserve the thermal structure of the UP adjacent to the trench (up to ca. 300 km away from the trench) to avoid a change in the corner flow velocity. A block density anomaly of +70 kg/m3 is used to allow direct comparison with simulation DA70.

The constant UP age in simulation OFDA70 (Figure 2.6.b) is 45 Ma (away from the trench), similar to the average of the variable ages along the UP in the "ridge-edge" simulations (e.g. DA70). This simulation shows minimal asthenosphere return flows (e.g. < 3.5 cm/yr at t = 8.8 Myr) below both plates, no trench retreat and no significant UP deformation throughout the simulation.

Changing the strength of the fixed upper plate

In simulation YFD70, the constant UP age is set to 15 Ma to simulate a younger and weaker UP (Figure 2.6.c). With this less resistant UP, the elevated slab pull is able to trigger faster asthenosphere flow, deforming the UP (Figure 2 For both simulations OFDA70 and YFDA70, the UPs exhibit extensional horizontal normal stresses (Figures 6.i.b and 6.i.c). However, in simulation YFDA70, the UP is stretched and thinned by ~10 % due to its lower resistance. Note that the stretched region is localized above the slab tip because of the upward asthenosphere return flow.

Discussion

Limitations of model setup

The geometry of the asthenospheric flow in our model results is partially dictated by the model's 2-dimensionality and its closed free slip side and bottom boundaries. In this section, we provide arguments for the robustness of our interpretations of slab pull variations driving upper deformation. We do, however, acknowledge that since our model does not include a third dimension, some 3-D effects such as lateral slab buoyancy variations or toroidal flow are not captured: these effects can also induce upper plate deformation (e.g. Schellart and Moresi, 2013;[START_REF] Magni | How Collision Triggers Backarc Extension: Insight into Mediterranean Style of Extension from 3-d Numerical Models[END_REF] One of our keys results is that the trench decelerates even though slab pull is relatively high, which we verified not to be imposed by our model geometry. The free slip right side is bound to slow down the slab as the slab tip approaches the boundary. However, this occurs quite late in our simulations, after 6 Myr, whereas the trench starts to decelerate at 3 Myr when the slab is steadily advancing (Appendix E.a). Note the constant acceleration of the slab tip during the elevated slab pull (Appendix E.b), with the slab tip velocity correlated to V sp.

The closed left-side boundary could also hamper rollback since the sub-slab and mantle-wedge flows are isolated from each other, however, this is not the case as the rollback is compensated by slab flattening, elongation of sub-slab channel flow, and to a lesser extent, rising free surface. For example, if we take an average rollback velocity of 3 cm/yr integrated along an inclined slab length of 1000 km, the 2-D volume reduction associated with slab rollback and trench retreat is around 10 -3 m 2 /s. This is mainly accommodated by slab flattening, as trench retreat causes a decrease of slab dip from 50° to about 35° (Appendix F), which creates a sub-slab volume increase of 1.3 x 10 11 m 2 . This corresponds to a 2-D volume rate increase of 0.7 x 10 -3 m 2 /s over 6 Myr.

Additionally, the asthenosphere below the subducting plate moves below the lengthening flat segment of the slab (sub-slab channel at the bottom of the model). If we

Influence of slab buoyancy on subduction dynamics

The force balance governing subduction is not yet entirely clear and has been addressed by numerous studies (e.g. [START_REF] Conrad | Effects of Plate Bending and Fault Strength at Subduction Zones on Plate Dynamics[END_REF][START_REF] Capitanio | Dynamic Models of Downgoing Plate-Buoyancy Driven Subduction: Subduction Motions and Energy Dissipation[END_REF]Husson et al., 2012;Goes et al., 2017). These have attempted to understand how the potential energy of the slab's negative buoyancy is transmitted to slab motion, surface plate motion and asthenospheric motion. They have also questioned how this energy is dissipated through slab bending at the hinge and through drag from the surrounding mantle. Fitting this energetics perspective, we propose a scenario linking the initial slab pull increase to the observed trench motion, slab rollback and upper plate deformation. We conceptually illustrate this in Figure 2.7, and emphasise the importance of competing shallow asthenospheric flows beneath the two plates in controlling slab rollback and trench retreat. This is consistent with the conclusions of [START_REF] Husson | Trench Migration and Upper Plate Strain over a Convecting Mantle[END_REF], who showed that mantle drag is a prominent component of the force-balance governing plate tectonics and trench motion.

Raising slab pull by introducing a dense block in the slab increases the overall speed of subduction (Figure 2.3 and Appendix A) and advance velocity of the slab tip (Appendix E), as seen in other numerical studies (e.g. [START_REF] Holt | Trench Migration and Overriding Plate Stress in Dynamic Subduction Models[END_REF]. The faster slab triggers acceleration in the surrounding asthenosphere on both sides through viscous drag. This has two effects: (i) in the sub-slab region, where it acts as a positive feedback accelerating the subducting plate; and (ii) in the mantle wedge beneath the UP where the negative pressures in the corner flow (see Appendix H) results in a "slab suction" force that drives the UP towards the trench (e.g. [START_REF] Conrad | The Temporal Evolution of Plate Driving Forces: Importance of 'Slab Suction' versus 'Slab Pull' during the Cenozoic[END_REF]. The latter return flow induces shear stresses at the base of the UP, dragging it towards the trench by viscous coupling (as in [START_REF] Husson | Trench Migration and Upper Plate Strain over a Convecting Mantle[END_REF], and Sternai et al., 2014[START_REF] Dal Zilio | The Role of Deep Subduction in Supercontinent Breakup[END_REF] have shown that maximum upper plate deformations occur for maximum lateral gradient of shear stresses. In our free edge models, the UP is also pushed leftwards from its ridge.

Trench retreat initially matches the velocity of the UP (Figure 2.3), suggesting that UP translation controls trench retreat (e.g. Cerpa et al., 2018). However, trench retreat decelerates around halfway through the subduction of the denser block (e.g. after 3.5

Myr for simulation DA70), even though slab pull remains high. This is not caused by changes in slab dip angle, since slab dip adjusts to trench motion with a ~2 Myr delay (see Appendix F). Additionally, in simulation ODA70, the slab progressively steepens as the trench remains stationary, accommodating increased slab pull. Slab dip, therefore, responds to surface trench motion. This is in line with the free subduction models of [START_REF] Capitanio | Dynamic Models of Downgoing Plate-Buoyancy Driven Subduction: Subduction Motions and Energy Dissipation[END_REF] where slab dip adjusts to minimize the energy dissipation associated with bending dissipation, and trench deceleration causes the increase in slab dip angle [START_REF] Royden | Subduction Zone Geodynamics[END_REF]. We suggest that the mantle-wedge return flow beneath the UP promotes slab rollback and trench retreat, while sub-slab flow dragged by the subducting plate opposes this motion. For example, DA70 experiences V tr deceleration at t 3 4.8 Myr. The asthenosphere velocities near the slab at a depth of ~125 km show a sub-slab flow of 8.4 cm/yr, but a mantle wedge flow of 6 cm/yr. Since the sub-slab flow is stronger than the mantlewedge flow, rollback is hampered, causing V tr deceleration. We interpret that this competition modulates trench velocity, as shown conceptually in Figure 2.7. Additionally, the rollback velocities at 200 km depth decelerate before those at 400 km depth, which also suggests that deceleration starts at shallow depths before propagating deeper (Appendix I)

Upper plate deformation during trench retreat

Table 2 shows that the width of the UP deformed region tends to increase with higher slab pull, however, the smallest density anomaly (simulation DA10) does not exhibit extra deformation compared to the reference simulation, despite its faster trench retreat. Thus, for low values of acceleration, the rigid translation of the UP from basal asthenosphere traction is fully accommodated by trench retreat with no discrepancy between V up and V tr (Figure 2.5, simulation DA10). Significant deformation only occurs when trench retreat cannot accommodate a rapid translation of the UP. This effect is seen for higher block density anomalies (simulations DA30 -DA100) due to their faster sub-slab flows, resulting in higher discrepancies between V up and V tr (Table 2, Appendix B). These results show that upper plate deformation in subduction zones can depend on the ability of the trench to retreat in response to a variation in slab sinking rates.

In our free edge models, this results in UP compression (Figure 2.4); however, our fixed edge models show extensional stresses (Figure 2.6). The fixed edge simulations develop a slab "kink" (sharp transition from a vertical to a flat slab) at the bottom of the box (Figure 2.6). This kink might accommodate trench retreat through slab dip variation.

However, we only observe trench retreat in YFDA70 (younger UP) and not in OFDA70 (older UP). In OFDA70, the UP is not free to translate horizontally, and experiences extensional horizontal stresses caused by mantle-wedge flow, but these stressed are insufficient to deform it enough to allow trench retreat (hence, the trench remains stationary since V tr = V up ). Contrastingly, the younger UP in YFDA70 can be deformed by the extensional shear stresses caused by the basal drag, permitting the trench to retreat.

The shear stresses are comparable in both simulations since the elevated slab pull is the same, but the strength of the UP differs, controlling trench retreat in this case. Where the trench retreats in YFDA70, the peculiar slab morphology with mid-mantle dip variation (Figure 2.6) is probably the result of this leftward motion of the trench at the surface.

Therefore, depending on the plate's properties, the same additional slab pull (density anomaly +70 kg/m3) can result in UP plate extension, shortening, or no deformation at all (Figure 2.6). Only a basal drag large enough to deform the UP can cause a discrepancy between UP and trench velocities (e.g. Wdowinski et al., 1998, andSternai et al., 2014). Thus, we consider the UP's motion and deformation to be mainly driven by asthenospheric basal drag, in agreement with the 2-D cylindrical simulations of [START_REF] Gérault | Flat-Slab Subduction, Topography, and Mantle Dynamics in Southwestern Mexico[END_REF]. If the UP is free to move, there is UP translation which trench retreat cannot match, resulting in compression. If the UP is fixed, the basal drag results in UP extension. This far-field influence on UP deformation has also been found by other studies, notably [START_REF] Nakakuki | Dynamics of Slab Rollback and Induced Back-Arc Basin Formation[END_REF], [START_REF] Chen | Overriding Plate Deformation and Variability of Fore-Arc Deformation during Subduction: Insight from Geodynamic Models and Application to the Calabria Subduction Zone[END_REF] and [START_REF] Capitanio | Upper Plate Controls on Deep Subduction, Trench Migrations and Deformations at Convergent Margins[END_REF].

Similar differences are observed in 3-D studies, where trench retreat of wide slabs results in predominantly poloidal flow and UP compression (e.g. Schellart and Moresi, 2013). On the other hand, trench retreat of a narrow slab can create toroidal flow, which results in UP extension (e.g. Schellart and Moresi, 2013;Faccenna et al., 2014). The extension observed in our 2-D fixed-edge simulations (Figure 2.6) in the absence of toroidal flow can be explained by a variable drag exerted by the mantle flow along the UP, as modelled by Schellart et al. (2017). We expect the presence of a 3-D toroidal flow to reinforce the shear drag exerted by mantle wedge flow at the base of the UP, enhancing plate deformation.

Application to natural subduction zones

Our results indicate that trench retreat can be associated with different upper plate deformation styles: no deformation for simulations DA10 and early DA70 (Figures 2 and5), compression at later times in DA70 (Figure 2.4), and extension in YFDA70 (Figure 2.6).

This can partly explain the lack of correlation between trench retreat rates and UP stress regimes observed in nature (Lallemand et al., 2005).

Since our 2-D models neglect the role of toroidal flow, the simulations best represent the dynamics taking place at the centre of a wide slab, where poloidal flow is dominant.

Our fixed-edge simulations reproduce UP extension in 2-D that is also observed in 3-D simulations. We extrapolate two end-member subduction systems exhibiting retreating trenches, summarized in Figure 2.8:

1-Upper plate shortening: requires fast lateral translation of the UP, which could be facilitated by a ridge on the UP's trailing edge (predominantly poloidal flow in the centre of a wide slab). In this case, asthenospheric flow is partitioned into lateral translation of the UP at the ridge, with less energy for basal asthenospheric drag.

A natural example of this case is the South America subduction zone.

2-Upper plate extension: requires asthenospheric basal drag that is sufficient to overcome the strength of the UP, with possible applications to the Mediterranean subduction zones, which have narrow slabs inducing strong 3-D toroidal flows.

The South American subduction zone has a trench that has been retreating since the Cretaceous, a strongly compressive upper plate (the Andes) and the Atlantic ridge in its trailing edge allowing the UP to translate westwards. Compression is thought to have initiated 50 Ma ago, which is synchronous with a discrepancy between V up and V tr , as the trench has been decelerating more than the UP (Faccenna et al., 2017). A hypothesis explaining UP compression is that the slab anchored when it reached the lower mantle, limiting rollback and triggering a large-scale return flow beneath the UP (Faccenna et al., 97 2017). Our results also suggest that compression in the Andes may result from a stronger asthenospheric flow beneath the UP (possibly limited to the upper mantle), triggered by faster subduction related to older slab subduction (e.g. Capitanio et al., 2011). Alternatively, compression could be associated with a larger push from the Mid-Atlantic ridge due to diversion of mantle flow from the African hotspot after Africa's collision with Eurasia (e.g. Silver et al., 1998;Husson et al., 2012). In contrast, the Mediterranean subduction zones (Aegean, Calabrian) show UP extension associated with accelerated trench retreat, which is thought to have resulted from increases in slab density (e.g. [START_REF] Jolivet | The Kinematics of Back-Arc Basins, Examples from the Tyrrhenian, Aegean and Japan Seas[END_REF] or slab tears (Wortel and Spakman, 2000). Our fixed edge models shed some light on these subduction zones, as the fixed UP is analogous to Eurasia, which is much larger than the individual subduction zones and is fixed relative to the Mediterranean slabs. Our results suggest that these fast trench retreats can be a result of UP extension (and the reverse). In the Aegean subduction zone, for example, the Aegean UP was thinned from the Oligocene to the Mid-Miocene (e.g. Brun et al., 2016 and references therein) reducing the crustal thickness by ~100% (Tirel et al., 2004). This means that the UP was weakened, making it more susceptible to deformation by asthenospheric toroidal flow introduced by slab tearing (Faccenna et al., 2014). While our 2D models do not capture these 3D effects, we do observe that a weaker UP can permit trench retreat if the basal flow is strong enough. Thus, if Aegean trench retreat accelerated in the Miocene (Brun et al., 2016), it could be because the UP was weak enough to thin more and allow further trench retreat, creating a positive feedback between UP deformation and trench retreat.

There are many other examples of subduction zones that are affected by dynamic changes through time (e.g. Japan, Banda, Apennines, Scotia, see [START_REF] Clark | Episodicity in Back-Arc Tectonic Regimes[END_REF]). It is hard to disentangle the exact causes for these changes, and our models and other modelling studies show that this can have various causes, including changes in neighbouring plate velocity through time (e.g. Japan, according to [START_REF] Jolivet | The Kinematics of Back-Arc Basins, Examples from the Tyrrhenian, Aegean and Japan Seas[END_REF]), or changes in incoming material (e.g. Banda and Apennines according to [START_REF] Royden | Subduction Zone Geodynamics[END_REF]). The complex relationship between trench motion and UP deformation style also holds for advancing trenches that do not systematically deform in compression, such as the Izu-Bonin or the Kermadec subduction zones (Heuret and Lallemand, 2005).

Our model results show that an increase in slab sinking rates triggers an increase in the spreading rate of ridges located at the end of both the subducting and upper plates. It is well known that slab pull has an effect on the spreading rate of the subducting plate ridge (e.g. [START_REF] Forsyth | On the Relative Importance of the Driving Forces of Plate Motion[END_REF], and our results also suggest a possible influence of subduction on the spreading rate of the oceanic upper plate. For example, the variations of the spreading rate of the North Atlantic ridge may have been influenced by the sinking slab dynamics of the South American slab (Silver et al., 1998).

Conclusion

Our 2-D models of subduction systems shed light on the relationship between slab pull, trench retreat, asthenosphere flow and upper plate deformation. Our results show that depending on the asthenospheric flows, the upper plate strength and its boundary conditions, a retreating trench can be associated with either an upper plate that is undeformed, in compression or extension. We suggest that stresses in the upper plate result from both basal drag by the mantle-wedge flow and from the ability of the trench to retreat, the latter depending on a competition between shallow asthenospheric flows below the subducting and upper plates. For example, if the upper plate is free to move (e.g. if there is a ridge in its trailing edge), deformation can occur in compression if trench retreat is hampered. On the other hand, if the upper plate is fixed, deformation can occur in extension, similar to narrow subduction zones with large toroidal flows. If the velocities of the upper plate and the trench are equal, the trench can retreat without any upper plate deformation. We find that trench retreat also depends on deformation in the upper plate, particularly in fixed edge upper plates, which do not translate rigidly.

Deformation here can permit trench motion. Hence, for the present-day and recent geological past, we emphasize the important role of asthenospheric flow in modulating trench motion and upper plate deformation.

Late, brittle faulting in the

Abstract

Back-arc extension in the Central Aegean has been active since the Late Oligocene and its style of deformation seems to have evolved through time. Extension started with ductile metamorphic rock exhumation in the Late Oligocene, and then high angle normal faulting from the Mid-Miocene, which continues today. The Aegean also experiences strike slip faulting related to the westward extrusion of Anatolia. This study aims to constrain the relationship between strike slip and recent high-angle normal faulting in the Cyclades. We use offshore seismic reflection profiles in the Cyclades, as well as structural data collected on the island of Syros to constrain the pattern and kinematics of faulting. We identify three main sets of faults: NW pure-normal faults, NNW oblique-normal faults, and NNE-NE strike slip faults. Of the two normal fault sets, we suggest that the NW pure-normal set is the youngest, while the NNW obliquenormal set is inherited and reflects block rotation since the Mid-Miocene. We suggest that dextral strike slip faulting has accommodated block rotation, which implies that strike slip faulting was active in the Miocene. We also interpret a more recent change in stress field, possibly due to Hellenic slab tearing since the Pliocene.

** In preparation for submission to Tectonics.

Introduction

The Aegean subduction zone is one of the most active regions in the Mediterranean.

The Hellenic trench is retreating south-westwards at a current rate of approximately 3 cm/yr (e.g. Shaw and Jackson, 2010). Palaeomagnetic studies indicate that the central Aegean has rotated around 23° clockwise north of Naxos, and around 33° counterclockwise south of Naxos since the Mid-Miocene (Morris and Andersen, 1996;Avigad et al., 1998). Hellenic trench retreat is, thus, facilitated by clockwise rotation of the Aegean upper plate about a pole near Scutari-Pec (e.g. van Hinsbergen and Schmid, 2012). The retreat of the subduction zone has created extension in the Aegean upper plate since around 45 Ma (Brun and Sokoutis, 2010). Since this time, extension in the Aegean has been accommodated by a progressive change in deformation style. It started with exhumation of high-pressure, low temperature rocks, followed by the formation of high-temperature metamorphic core complexes (e.g. Lee and Lister, 1992;Jolivet et al., 1996;Avigad et al., 1997;Krohe and Mposkos, 2002;Kumerics et al., 2005;Brichau et al., 2007;Brun and Facenna, 2008). High-pressure, low-temperature exhumation brought the Cycladic Blueschist Unit (CBU) to the surface, which now makes up the Cyclades. Prior to blueschist and eclogite facies metamorphism, the CBU's protolith originated as part of a continental passive margin (Adria) with remnants of ophiolitic melange rocks (Pindos Ocean) [START_REF] Philippon | Deciphering Subduction from Exhumation in the Segmented Cycladic Blueschist Unit (Central Aegean, Greece)[END_REF]. Ductile exhumation in the Aegean ended around the Mid-Miocene, following which the style of deformation changed to high-angle normal faulting (Brun et al., 2016), creating dispersed sedimentary basins (e.g. Mascle and Martin, 1990;Sánchez-Gómez et al., 2002;Beniest et al., 2016). Andersen, 1996;Avigad et al., 1998). Bathymetry from NOAA [START_REF] Amante | ETOPO1 1 Arc-Minute Global Relief Model: Procedures, Data Sources and Analysis[END_REF].

The Aegean is also affected by westward extrusion of Anatolia, mainly accommodated by the North Anatolian Fault since the Pliocene (e.g. Armijo et al., 2004). There is, however, evidence that Anatolian extrusion started earlier in the Mid-Miocene. This is indicated by the initiation of large-scale dextral shear zones around 11-13 Ma in Turkey (Şengör et al., 2005). There are also indications of NE trending strike slip faults in the Central Aegean, partly controlling magma emplacement since the Mid-Miocene (Kokkalas and Aydin, 2013). Therefore, rollback-related extension may have been interacting with extrusion-related strike slip since the Mid-Miocene [START_REF] Philippon | Tectonophysics The Interaction between Aegean Back-Arc Extension and Anatolia Escape since Middle Miocene[END_REF].

We, thus, aim to characterise the style of deformation in the Central Aegean since the Mid-Miocene. We aim to explain how normal faulting has recently evolved with the retreating trench, and how extension interacts with the extrusion regime. We focus our study on the Cyclades block in the Central Aegean. Several studies have been conducted on the neotectonics of continental Greece (e.g. Hatzfel, 1999;Taymaz et l., 1991;Shaw and Jackson, 2010;[START_REF] Sachpazi | Segmented Hellenic Slab Rollback Driving Aegean Deformation and Seismicity[END_REF][START_REF] Vamvakaris | A Detailed Seismic Zonation Model for Shallow Earthquakes in the Broader Aegean Area[END_REF]. Recent offshore faults have also been mapped out in the Aegean Sea (e.g. Mascle and Martin, 1990;Beniest et al., 2016;Sakellario and Tsampouraki-Kraounaki , 2018). However, little work has been done on recent faulting in the Cyclades. We, thus, characterise the present day fault pattern in the Cyclades by identifying multiple fault sets and generation. We reconcile these with the regional tectonic framework, particularly with the plate kinematic reorganisation related to the development of the North Anatolian Fault and the slab tear beneath the Kephalonia Transform Fault since the late Miocene and Pliocene (Royden and Papanikolaou, 2011; Perouse et al., 2012) To do this, we use a combination of offshore seismic data as well as onshore field data.

We conduct our field study on the island of Syros since its geology is well understood and since it has documented late high angle faults (e.g. [START_REF] Keiter | Structural Development of High-Pressure Metamorphic Rocks on Syros Island (Cyclades, Greece)[END_REF]. Offshore data permits the identification of recently active faults and their basin-scale pattern, while onshore data allows us to accurately characterise their kinematics. From this, we can distinguish pure normal faults from oblique slip. With these kinematics clarified fault kinematics, we interpret normal fault evolution in the Cyclades in the context of slab rollback, trench retreat and strike slip faulting. We have performed a structural interpretation, which is mainly based on syn-kinematic sediment packages in faulted blocks. We identify faults as steep features which abut other reflectors. Sediment packages are identified by their closely spaced reflectors (bedding), and their abrupt interruption against steep faults. Syn-kinematic packages additionally have bedding reflectors which do not have continuous thickness, i.e. these packages thicken where accommodation space is created during fault activity. The position of this thickening within a graben indicates which bounding fault is dominant, as the largest amount of accommodation space is created closest to the dominant fault. Unlike sediment packages, the basement (Cycladic Blueschist Unit) does not show many internal reflectors, and the poor data quality limits interpretation inside the basement. We therefore only interpret the top of the basement where possible.

The strike of active faults is interpreted from bathymetry. Where possible, interpreted faults have been correlated from one seismic profile to another using a combination of the bathymetry and the character of the imaged fault. As the seismic digital data files are not available and the lines are not closely spaced, horizon correlation was not possible.

Two main fault families have been interpreted, normal faults and strike slip faults.

Below, we show two representative interpreted profiles, which illustrate the reasoning behind our interpretation. For the complete set of interpreted profiles and their locations, refer to the supplementary materials.

Interpretation of seismic profiles

-Active high angle normal faulting:

We observe normal faulted (half) grabens trending NW-SE to WNW-ESE, based on the bathymetry and correlated from neighbouring seismic lines. -Compressive structures within the grabens:

In addition to normal offset, the imaged half grabens show broad folds within the sediment infill. This is seen in Figure 3.3, where the sediments are slightly folded throughout the imaged depth. The intensity of folding increases with depth in the visible syn-kinematic sediments. This broad folding is seen in all of the imaged profiles except for one (seismic profile 3, Suppl. Mat. We observe evidence for both minor and major strike slip faults. An example of minor strike slip faulting is illustrated in Figure 3.3, where two small negative flower structures can be seen on the eastern edge of the half graben. The subsurface flower, F5, seems to have ceased activity, as it is draped with flat sediments. The easternmost flower, F6, is presently active as it creates a small trough in the sea floor.

Interestingly, the older negative flower, F5, appears to cross cut F7, the eastern binding fault of the graben. The latter is interpreted as a normal fault rather than an eroded surface since the sediments thicken in the middle of the faulted block, not against fault F2. This suggests that the faults were active on both sides of the graben, and that F2 is currently dominant based on sea floor offset. This is supported by the sediments on F7 showing a very high degree of rotation, suggesting it is (or at least has been recently) active but not dominant in the graben. 120 These widely ranging fault orientations of both sets combined with the observation of broad compressive structures in the sediment packages suggests complex fault kinematics. As fault kinematics are not observable from seismic data, we address this by analysing fault exposures onshore. We thus present field data from the island of Syros, in the central Cyclades normal and strike slip faults are observed, and presented below.

Faulting in Syros

Tectonic framework of Syros

Syros is mainly made up of the Cylcadic Blueschist Unit (CBU), a sequence which has undergone high pressure metamorphism during subduction in the Eocene (e.g. [START_REF] Dixon | Glaucophane Schists of Syros, Greece[END_REF]. The sequence experienced retrogression and ductile exhumation in the Oligocene -early Miocene [START_REF] Bröcker | The Geological Significance of 40Ar/39Ar and Rb-Sr White Mica Ages from Syros and Sifnos, Greece: A Record of Continuous (Re)Crystallization during Exhumation?[END_REF] and final exhumation by brittle faulting took place in the mid-late Miocene [START_REF] Ring | Fast Extension but Little Exhumation: The Vari Detachment in the Cyclades, Greece[END_REF][START_REF] Philippon | Tectonics of the Syros Blueschists (Cyclades, Greece): From Subduction to Aegean Extension[END_REF]Soukis and Stöckli, 2013). The Syros CBU is structurally arranged in NE dipping layers. From the structural base to the top, the sequence is made up of albitic micaschists and gneisses, alternating marbles and mica-schists, and metabasites [START_REF] Hecht | Geological Map of Greece 1:50.000, Syros Island[END_REF][START_REF] Keiter | Structural Development of High-Pressure Metamorphic Rocks on Syros Island (Cyclades, Greece)[END_REF][START_REF] Keiter | A New Geological Map of the Island of Syros (Aegean Sea, Greece): Implications for Lithostratigraphy and Structural History of the Cycladic Blueschist Unit[END_REF][START_REF] Philippon | Tectonics of the Syros Blueschists (Cyclades, Greece): From Subduction to Aegean Extension[END_REF].

Since exhumation, Syros has suffered brittle faulting (e.g. [START_REF] Philippon | 3D Geometrical Modelling of Post-Foliation Deformations in Metamorphic Terrains (Syros, Cyclades, Greece)[END_REF][START_REF] Philippon | Tectonics of the Syros Blueschists (Cyclades, Greece): From Subduction to Aegean Extension[END_REF][START_REF] Keiter | A New Geological Map of the Island of Syros (Aegean Sea, Greece): Implications for Lithostratigraphy and Structural History of the Cycladic Blueschist Unit[END_REF] which has received far less attention than its ductile history. In this section, we present results from structural mapping of late stage brittle structures on Syros, summarized in Figure 3 The third set of faults observed in Syros, trending NNW, was already documented by [START_REF] Philippon | 3D Geometrical Modelling of Post-Foliation Deformations in Metamorphic Terrains (Syros, Cyclades, Greece)[END_REF] in the Pyrgos marbles near Papouri (Figure 3.5). At map scale, this fault set is visible in south-western Syros in Galissas (Figure 3.8.a). This fault set is also observed in Palos and correlated to minor faulting observed offshore, shown in

Dextral strike-slip faulting, trending NE-SW

Central Syros is along strike the minor strike slip fault observed on seismic line 3 (F6, Figure 3.3). Following this, we find field evidence for a similarly minor NE-SW trending strike slip fault zone across central Syros. We interpret this to be a minor fault set in Syros, consisting of segmented small-scale structures. The fault planes are poorly exposed here, with the few exposures documented in Figure 3.9. Google Earth satellite imagery shows that this fault zone offsets marble layers with a dextral sense of shear (Figure 3.9.a). One fault plane was found revealing oblique lateral slip, shown in Figure 3.9.a. creating a broad fold in the Pyrgos Marbles [START_REF] Philippon | Tectonics of the Syros Blueschists (Cyclades, Greece): From Subduction to Aegean Extension[END_REF]2015). This fault zone also separates the gently dipping southern half of the island from the steeper northern half, which could be explained by the presence of a strike slip fault. Minor NE-SW strike slip faults have been identified both offshore and onshore in Syros with a dextral sense of shear. One major strike slip zone limits the Cyclades from Attica to the West, with a segment apparently trending NNE-SSW, as already proposed by Sakellariou and Tsampouraki-Kraounaki (2018). The regional strike slip faults in the central Aegean mainly trend NE-SW and have dextral offset (e.g. Taymaz et al., 1991 for earthquake focal mechanisms and Sakellariou and Tsampouraki-Kraounaki, 2018 for offshore faulting patterns). Our interpreted minor strike slip faults on and west of Syros is consistent with this. However, our interpreted strike slip zone in the northern Cyclades has a segment with an apparent trend of NNE-SSW, following the bathymetry. In the north-eastern corner of the Cyclades, the strike slip zone trends NE-SW, but it is not clear from bathymetry alone if these two segments are continuous. To the NE of the Cyclades, the NE-SW trending fault shows a dextral earthquake focal mechanism (Global CMT catalogue, [START_REF] Dziewonski | Determination of Earthquake Source Parameters from Waveform Data for Studies of Global and Regional Seismicity[END_REF][START_REF] Ekström | The Global CMT Project 2004-2010: Centroid-Moment Tensors for 13,017 Earthquakes[END_REF].

Attempting to interpret the NNE-SSW segment as part of the same fault proves problematic, as it strikes at a high angle to the main strike-slip direction, but shows strike slip extension (negative flower) and not compression. It could possibly be a synthetic P structure within the Riedel framework, as interpreted by Sakellariou and Tsampouraki-Kraounaki (2018, but such structures typically have a closer angle to the main strike slip zone. Alternatively, it could be an inherited structure that has mixed kinematics. More data is necessary to reliably interpret this structure. Regardless, this structure does not apparently extend for more than ~65 km before the regional NE-SW dextral strike slip pattern resumes.

The normal faults interpreted from both Syros and the offshore seismic profiles in the Cyclades show two dominant directions: NW-SE and NNW-SSE. Both sets appear to be presently active offshore, as they both offset the sea floor. Slickenlines observed on Syros show pure or near pure vertical slip on the NW-SE fault set (e.g. Figure 3.7), and sinistral oblique slip on the NNW-SSE fault set (e.g. Figure 3.8). Additionally, the pullapart structure in Fabrika (Figure 3.6) shows an older sinistral fault trending NNW-SSE with a later extensional structure trending NW-SE.

We extrapolate this to the offshore faults and interpret similar pure vertical slip on the NW-SE segments, and oblique sinistral slip on the NNW-SSE segments. The sinistral component of these oblique faults could be the cause of the broad compressive structures discussed in section 2.2. Figure 3.10 shows our three interpreted faults in the [START_REF] Gueydan | Sequential Extension as a Record of Corsica Rotation during Apennines Slab Roll-Back[END_REF].

Extrusion and strike slip

In addition to normal faulting, our data shows evidence of dextral strike slip faulting in the Cyclades. These are observed as both major (Figure 3.3) and minor structures (e.g Figure 3.9). The minor structures strike NE-SW, which is consistent with the regional dextral fault pattern. The major structure, however, trends NNE-SSW, at a high angle to the regional pattern. This, combined with the observation of multiple normal fault sets, raises several questions pertaining to the timing and tectonic conditions that control brittle faulting in the Central Aegean. The normal and oblique faults here are a product of back arc extension, which is a function of recent slab rollback (see discussion in section 2.5.2). On the other hand, the presence of strike slip faults (this study and Sakellario and Tsampouraki-Kraounaki , 2018) show that the Aegean is affected by extrusion of Anatolia in addition to back-arc extension (Şengör et al., 2005). It is not clear when strike slip activity started in the Central Aegean, however, there is evidence for regional Mid-Miocene initiation of large-scale dextral shear zones (e.g. Şengör et al., 2005;Kokkalas and Aydin, 2013). If this is the case, these dextral strike slip zones could have facilitated block rotation (represented in is incompatible with NE-SW dextral strike slip faulting. If this stress regime was active in the past and was synchronous with westward-directed extrusion, there must have been strain partitioning between extrusion-related structures and rollback-related structures. Contrastingly, today's extrusion and rollback related stress fields are quite close, which means strain partitioning is not necessary. We have shown that field studies can shed light on offshore fault kinematics in the Aegean. Further such field studies could provide age constraints on the timing of these various fault sets.

Conclusion

We have combined structural field data from Syros with offshore seismic interpretation to characterise the fault pattern and kinematics in the Cyclades. We draw the following conclusions:

• The Cyclades shows three co-existing fault sets that have different kinematics: NW-SE pure normal faults, NNW-SSE oblique (sinistral-normal) faults, and NNE-SWW to NE-SW dextral strike slip faults.

• NNW oblique faults and NW pure normal faults represent multiple generations of normal faulting during Aegean trench retreat.

• The radial retreat of the trench and the dextral strike slip faults (related to Anatolian extrusion) can accommodate block rotation, and subsequently fault rotation.

• This implies that dextral strike slip faulting may have been active since the Mid-Miocene, well before the development of the NAF.

• An additional change in stress pattern is needed to explain the present day stress pattern and resulting NW pure normal faults. This change may be related to recent slab tearing on the western side of the Hellenic slab, which lead to reorganisation of strain in the Pliocene with the development of the NAF and Corinth Gulf.

Acknowledgments

This research has received funding from the European Union's EU Framework Programme for Research and Innovation Horizon 2020 under Grant Agreement No 674899.

Introduction

Plate tectonics is a relatively young discipline, and much of the community's work still lies in the fundamental research stage. However, as plate tectonics shapes the Earth's surface, it influences habitability and resources. While some applications of tectonic studies are fairly advanced, such as mineral resource emplacement, other potential applications are farther from reach, such as earthquake prediction. There are applications which are not yet fully commercial, but are attainable in the near future.

One example of this is geothermal energy production.

Geothermal energy is a renewable and green energy source, which utilises heat in the crust, as introduced in Chapter 1.4. In this chapter, the role of tectonic modelling in geothermal energy production was also demonstrated, particularly for improving geothermal energy predictability in non-volcanic regions. The Aegean is a prime area to explore this (e.g. Cloetingh et al., 2010), as it has widespread high surface heat flow (Figure 4.1), with magmatism and hot springs. Additionally, the complex tectonic history of the Aegean is relatively well understood, rendering it a good candidate to test the ability of tectonic modelling to predict geothermal potential. † † This study was done at TNO Utrecht, in collaboration with Utrecht University. The work was carried out with Kristof Porkalob, Jon Limberger, Jan-Diederik Van Wees and Fred [START_REF] Mendrinos | Exploring for Geothermal Resources in Greece[END_REF]. Continental basins, thus, host significant geothermal resources, most of which are low temperature at the surface (known from hot springs), but higher temperatures at depth (oil exploration wells have proven opposition to high temperature geothermal use. Social awareness needs to be raised before high temperature geothermal resources are accepted again, but low-medium temperature resources are still perceived positively (Karytsas et al., 2019). Therefore, there is an opportunity to develop geothermal production from reservoirs with low surface temperatures, but medium temperatures at depth. As many of these reservoirs lie in continental basins, there is a clear relationship between tectonic activity and geothermal potential in the Aegean, even away from the magmatic arc. In this case, the tectonic activity is in the form of lithospheric extension. This begs the question: can tectonic modelling predict geothermal reservoir potential?

This study aims to address this question using two short thermal modelling workflows, one that is steady state, and another that incorporates extension. We do not attempt to accurately model thermal potential in the Aegean, but only to test whether incorporating extension changes the predicted thermal potential. This will essentially show whether or not it would be valuable to incorporate tectonic modelling into existing thermal modelling workflows. We therefore use the thermal modelling workflow documented in Limberger et al. (2018) and build on the 3D thermotectonic model created by Larede (2018). Larede (2018) has already shown that varying the lithospheric thickness has a first order effect on thermal distribution. In our study, we test the effect of crustal lithological composition in steady state, and the effect of simplified lithospheric stretching over 40 Ma.

Aegean geological setting:

The Aegean subduction zone is one of the most actively deforming regions in the Alpine-Himalayan Belt. Subduction initially created a thickened nappe-stack which later collapsed when the subduction zone started retreating around 55 Ma (Ring et al., 2010 and references therein). During this time, extension started with ductile exhumation and later progressed to brittle normal faulting. The onset of rollback is expected to have largely raised temperatures in the Aegean as the cold slab migrated northwards, and exhumation of high temperature metamorphic core complexes is seen (Ring et al., 2010). The exhumation phase is also thought to have largely thinned the lithosphere (e.g. Brun et al., 2016). The Moho is presently thought to lie around 25 km depth (Tirel et al., 2004), but the thickness of the lithosphere is more poorly constrained. For a more complete overview of the Aegean's geological setting, refer to chapter 1.3. This is the a priori model and starts from a lithological geometrical model. The a priori thermal properties are then populated for each grid cell, to give k (z), a collection of 1D thermal conductivity profiles (W m -1 K -1 ) as a function of depth for the sediments, upper crust, lower crust and lithospheric mantle. The bulk matrix conductivity, k m , (Limberger et al., 2018[START_REF] Sekiguchi | A Method for Determining Terrestrial Heat Flow in Oil Basinal Areas[END_REF] is then used together with lithological properties (Table 1) for the thermal calculations (Hantschel and Kaurauf, 2009;[START_REF] Bär | P3-International PetroPhysical Property Database[END_REF]Limberger et al., 2018). Radiogenic heat production is accounted for and thermal conductivity is also corrected for temperature (Van Wees et al., 2009;Limberger et al., 2018). The a priori model is then used to create a 3D forward model, using a finite-difference approximation, and assuming steady state. This is based on the multi-1D model presented in Van Wees et al. (2009), but with addition of vertical stretching in 3D,

Input data

For the geocellular model input,

we use the crustal model from Larede, 2018, who constructed the model using the software GeoModeller ‡ ‡ , and based the layers on data from Tirel et al., 2004;[START_REF] Hinsbergen | Nappe Stacking Resulting from Subduction of Oceanic and Continental Lithosphere below Greece[END_REF]Brun and Sokoutis, 2010;Jolivet and Brun, 2010;Jolivet et al., 2013;Beniest et al., 2016, and[START_REF] Menant | Kinematic Reconstructions and Magmatic Evolution Illuminating Crustal and Mantle Dynamics of the Eastern Mediterranean Region since the Late Cretaceous[END_REF] This is a simplified lithological model which captures the Hellenides' nappe stack at first order, but does not incorporate a modelled fault network. The main purpose of using this model is capture the different units' thermal properties and their first order distribution across the Aegean. ‡ ‡ www.intrepid-geophysics.com

Generalised workflow

In our study, we first created a present day steady state model, using the Moho (from Tirel et al., 2004) as the bottom boundary condition.

We then created a transient model spanning the last 40 Ma. Here, the lithospheric stretching is forward modelled and employs the principles of non-uniform stretching,

where different stretching factors are used for the crust (⍺) and mantle lithosphere (β) (Royden and Keen, 1980). We consider that stretching in the Aegean started slowly at 40 Ma and has sped up through time (e.g. van Hinsbergen andSchmid, 2012, Brun et al., 2016) The forward model, thus, calculates the evolution of the geotherm with lithospheric stretching, incorporating material properties and thermal velocities ( [START_REF] Wees | Quantitative Modelling of Basin and Rheological Evolution of the Iberian Basin (Central Spain): Implications for Lithospheric Dynamics of Intraplate Extension and Inversion[END_REF] assuming fixed strain-rates.

Synthesis

This thesis has explored trench retreat and upper plate deformation, addressing three main questions:

-How does slab pull influence trench retreat?

-What is the relationship between trench retreat and upper plate deformation?

-Can tectonic research feed back directly into applications that can be useful for society?

In order to answer these questions, a multi-scale approach was used to better understand the physical processes in retreating subduction zones, and what controls upper plate deformation in this setting. We used large-scale (> 100 km) numerical models to isolate the effect of slab pull on trench retreat, and subsequently on upper plate deformation. We used field and seismic data from the Aegean upper plate to better understand basin-scale (> 1 km) brittle fault patterns in retreating subduction systems.

Our models show that trench retreat is not only controlled by slab pull, but also by shallow flow in the mantle below the surface plates and around the subducting slab. Larger slab pull causes faster subduction, which in turn triggers faster mantle flow beneath both the subducting plate and the upper plate. These faster mantle flows compete, and the trench can retreat if the return flow in the mantle-wedge (beneath the upper plate) is faster than the sub-slab flow. However, if the slab subducts faster and drags the sub-slab mantle so that it flows faster than the mantle wedge, it resists slab rollback and trench retreat, even though slab pull is high and subduction is fast. This shows that slab pull is not the only factor affecting the motion of the trench, but that mantle flow (whether or not it is a result of slab pull) is also a strong control on trench motion.

We found that mantle flow can also be responsible for upper plate deformation; thus, trench motion can only indirectly cause upper plate deformation if it triggers faster mantle flow. Our models shows that if mantle flow beneath the upper plate is strong enough to overcome the strength of the upper plate, it can drag the upper plate from below and create deformation. The type of deformation depends on whether the entire plate is dragged or only part of it. The entire plate can be dragged if it has a free edge, for example a ridge in its trailing edge. In this case, the upper plate deforms in compression. On the other hand, if only part of the upper plate is dragged, as is this case if its trailing is fixed because it is part of a large continent, it deforms in extension since its front is dragged away from its trailing edge. In both scenarios, the strength of the upper plate also matters, as it needs to be weak enough to be deformed by the underlying mantle flow. We illustrate this conceptually in This strong mantle flow coupled with a thinned upper plate may be why the Aegean exhibits such elevated heat flow compared to its surroundings. This higher heat flow coupled with hot springs makes the Aegean ripe for geothermal energy exploration, however, understanding the lithospheric structure would largely aid in this effort. The thinned upper plate allows the mantle to flow closer to the surface, however, it also reduces blanketing allow greater diffusion.

Sedimentation, in turn, counteracts this, but the presence of faults which bring up thermal fluids from depth raises the geotherm. Quantifying these effects tectonically can largely improve geothermal predictions, and chapter 4 showed that simply adding lithospheric extension into a thermal modelling workflow significantly changed the predictions.

We often think of geodynamic research as fundamental research far removed from current applications, however, chapter 4 shows that geodynamics has present day, real world applications. The geothermal example is one of a direct application of tectonic research which could introduce geothermal energy to non-volcanic regions. This could increase global renewable energy production. Managing geohazards is the other obvious direct use of geodynamic research. There are many other applications of this kind of research, which are too many to list here, but they are largely contributing to our fast developing technology. Our lives would look very different if we did not benefit from the wealth of resources our planet offers us, and we would not be able to tap into these resources without understanding the planet. So for the sake of our own development, and for the sake of knowledge, the ultimate pursuit, let us always be curious and intrigued by our beautiful Earth.

Figure 1 *

 1 Figure 1 * : Zones de subduction majeures globales, vitesse de subduction de la plaque plongeante en bleu et migration de la fosse en rouge.

Figure 2 †

 2 Figure 2 † : Graphe de l'âge du slab en fonction du mouvement de la fosse qui ne montre aucune corrélation claire.

10 Figure 3 ‡

 103 Figure 3 ‡ : Pas de corrélation claire entre le style de déformation global du domaine arrière-arc et le mouvement absolu de la fosse.

Figure 6 :

 6 Figure 6 : Interprétation des failles de la mer Égée centrale

Figure 7 :

 7 Figure 7 : Interprétation conceptuelle du recul de la fosse, de la rotation de bloc et de la rotation des contraintes en Égée Centrale.

Figure 1

 1 Figure 1.1.

Figure 1 . 1 :

 11 Figure 1.1: Global major subduction zones, trench normal subducting plate velocity in blue, and trench normal trench migration velocity in red. Velocities calculated using the Indo-Atlantic hot spot reference frame (O'Neill et al., 2005). Figure taken from Schellart et al. (2011).

Figure 1 . 2 :

 12 Figure 1.2: Subduction zone simplified force balance: slab pull from negative buoyancy together with the mid-upper mantle phase transition act as positive sinking forces;

Figure 1 . 3 :

 13 Figure 1.3: Plotting slab age against trench motion shows no clear correlation (Goes et

Figure 1 . 4 :

 14 Figure 1.4: No clear correlation between global back-arc deformation style and trench absolute motion(Heuret and Lallemand, 2005).

Figure 1 . 2 . 1 :

 121 Figure 1.2.1: Analogue models of subduction and trench retreat by Faccenna et al., 2001: the slab is simulated using silicone plates and the mantle using glucose syrup mixed with pure white sugar to vary viscosity. The experiments ultimately simulated episodic trench retreat.

Figure 1 . 2 . 2 :

 122 Figure 1.2.2: Different proxies used to define the lithosphere-ashthenosphere boundary (LAB). After Eaton et al. (2009).

Figure 1 . 2 . 3 :

 123 Figure 1.2.3: Examples of 2-D and 3-D numerical models that explore UP dynamics: a) 2D model of flat slab subduction resulting from upper plate acceleration (Van Hunen et al., 2004). b) 3D model of hyper-episodic subduction with alternating trench retreat and

  episodic. Natural systems (e.g. Izu-Bonin Trench, the Mariana Trench, the Japan Trench, the Java-Sunda Trench and the central portion of the Peru-Chile Trench) are thought to be dominantly quasi-episodic, where the back-arc shifts from rifting to spreading to tectonic quiescence. The authors found that this requires fast subducting plate velocities (6-9 cm/yr), as does hyper episodicity, where the trench alternates between advance and retreat. Slower subducting plates, however, are characterised by slower subducting plate motion, such as the in the Mediterranean or Scotia Sea.

Figure 1 . 2 . 4 :

 124 Figure 1.2.4: Illustration of toroidal flow (lateral flow) and poloidal slab (vertical flow) around the slab. Illustration from Schellart and Moresi (2013).

Figure 1 . 2 . 5 :

 125 Figure 1.2.5: The scale of poloidal mantle flow can determine upper plate deformation. a) Upper mantle only, small scale convection can create upper plate extension. b) Whole mantle, large scale convection can create upper plate compression. Illustration from Faccenna et al. (2017).

  zone is part of the complex Alpine-Mediterranean mobile belt, where slabs have been torn, segmented, back-arc basins opened and micro-continents formed (Figure 1.3.1). Only two zones of active subduction remain in the Eastern Mediterranean, the Tyrrhenian and the Aegean. Complex slab dynamics here continue as the slabs tear and rollback further. A simplified illustration of the tectonic evolution of the Eastern Mediterranean is shown in Figure 1.3.2.

Figure 1 . 3 . 1 :

 131 Figure 1.3.1: Major tectonic elements of the Mediterranean system(Royden and Faccenna, 2018)

Figure 1 .

 1 3.3 shows an overview of the main geological units and their simplified tectonic history in cross section.

Figure 1 . 3 . 2 :

 132 Figure 1.3.2: Tectonic evolution of the Eastern Mediterranean since 30 Ma (Royden and Faccenna, 2018). Note the reconstructed trench retreat of the Aegean SZ towards Africa.

Figure 1 . 3 . 3 :

 133 Figure 1.3.3: The main geological units in the Aegean in a) map view, b) cross section, and c) cross sectional tectonic history. After van Hinsbergen et al. (2005).

  : prograde subduction-related deformation; rollback-related exhumation of metamorphic core complexes (Figure 1.3.4); and rollback-related segmentation by high angle normal faults (Figure 1.3.5).

Figure 1 . 3 . 4 :

 134 Figure 1.3.4: Exhumation of (a) high pressure, and (b) high temperature metamorphic core complexes in the Aegean from the Mid-Eocene to the Mid-Miocene (circled numbers indicate age in Ma)(Brun et al., 2016).

Figure 1 . 3 . 5 :

 135 Figure 1.3.5: Segmentation of the Aegean upper plate by high angle normal faults creating widespread sedimentary basins in the Neogene (Brun et al., 2016).

Figure 1 . 3 . 6 :

 136 Figure 1.3.6: Hellenic trench retreat since the Eocene and apparent associated deformation(Brun et al., 2016).

Figure 1 . 3 . 7

 137 Figure 1.3.7 Example of offshore high-angle normal faulting in the Cretan Sea.

Figure 1 . 3 . 8 :

 138 Figure 1.3.8: Major Aegean structures and GPS field in the Aegean, interpreted by Armijo et al., 2004: red and blue shading respectively represent extensional and compressional strain associated with the propagation of the North Anatolian Fault. Rollback related extension is shaded in yellow.

Figure 1 . 3 . 9 :

 139 Figure 1.3.9: Earthquake focal solutions for shallow earthquakes of magnitude >/= 4.0. a) normal faults, b) thrust faults, c) strike-slip faults. From Vamvakaris et al. (2016).

Figure 1 . 3 . 10 :

 1310 Figure 1.3.10: Interpretation of recent the deformation pattern in the Aegean by Sakellariou and Tsampouraki-Kraounaki (2018), which includes normal and strike slip faulting. The insets show the author's interpretation of dextral oblique extension (B) that fits with the Riedel shear framework (A), and creates sinistral pull apart structures in the south Eastern Aegean (C). Volcanic arc marked by red dots, yellow areas cover Plio-Quaternary basins with sediment thickness > 500 m.

Figure 1 . 3 . 11 :

 1311 Figure 1.3.11: Main tectonic elements of the Cyclades with rotation since the Mid-Miocene and an interpreted dextral strike slip accommodating the rotation. After Philippon et al., 2014.

Figure 1 .

 1 Figure 1.3.12: P-T path of the rocks in Syros showing progression from prograde subduction-related metamorphism through retrograde rollback-related extension. From Laurent et al., 2016.

1. 4 . 1 . 4 . 1 .

 4141 The role of tectonic research in predicting geothermal energy Geoscience and energy resourcesGeoscience and energy production have long been closely intertwined, from the early days of coal mining, through the petroleum era and onwards to alternative resources; energy production usually starts from the natural environment in which we live. While several energy sources are derived from the energy Earth's surface, such as solar, wind and hydroelectric power, venturing deeper into the crust has provides a wealth of energy resources. Future generations will undoubtedly uncover further resources we

Figure 1 . 4 . 1 :

 141 Figure 1.4.1: Global geothermal installed capacity and produced energy from 1950

  Conventional deep geothermalDeep geothermal, on the other hand, exploits a high geotherm by drilling into the crust and accessing heated fluids or vapour. There are different types of power plants used in this scenario depending on the fluid temperatures. For temperatures greater than 250°C, dry steam plants are possible. These use geothermal fluids, which are completely vaporised (often due to the drop in a drop which results from drilling) and are piped directly to the power plant. Temperatures greater than 180°C facilitate flash steam systems, where the geothermal fluid is only partially vaporised at the surface (or flashed). For temperatures as low as 110°C, binary cycle technologies can be used, where the geothermal fluid exchanges heat with a working fluid that has a lower boiling point and higher vapour pressure. The working fluid us then vaporised in a heat exchanged.All three types of plants require at least two wells, one for production and the other to inject the geothermal fluids back into the ground, although only the binary cycle type plant is able to re-inject all the of the geothermal fluid, as steam is lost to the atmosphere in the other two plant types. The binary cycle type plant, however, has the disadvantage of lower efficiency, and thus only accounts for 12% of global installed capacity[START_REF] Bertani | Geothermal Power Generation in the World 2010-2014 Update Report[END_REF]. Single flash power plants, on the other hand, account for 41%, while dry steam plants account for 22% of global installed capacity. Other technologies, including hybrid technology plants, account for the remaining 3%. A schematic overview of the geothermal power plants is shown in Figure1.4.2.

Figure 1 . 4 . 2 :

 142 Figure 1.4.2: Simplified overview of geothermal power plants (Manzella, 2019).

Figure 1 . 4 . 3 :

 143 Figure 1.4.3: Conceptual image of the Enhanced Geothermal System (Gerard et al., 2006).

Figure 1 . 4 : 1 . 4 . 8 .

 14148 Figure 1.4: European heat flow map(Cloetingh et al., 2010).

  of Earth Sciences, The Australian National University, Canberra,AustraliaAbstractUpper plate deformation within a subduction zone depends on the complex relationship between surface plate motion, trench motion, slab pull and shallow mantle flow.Previous modelling studies suggest that trench motion rates should be related to slab buoyancy, but this relationship is neither clear nor verified by observations of natural subduction systems. Trench motion is also thought to induce upper plate deformation; however, no clear correlation has been identified between the direction of trench motion and the mode of upper plate deformation. In this study, we construct 2-D thermomechanical models to explore the relationship between slab pull, trench retreat and upper plate deformation, focusing on subduction systems with retreating trenches. We start with quasi-steady-state subduction and introduce a positive density anomaly into the slab to transiently increase slab pull. We vary both the value of the density anomaly and the properties of the upper plate to isolate key controls on trench retreat and upper plate deformation. Our models demonstrate that asthenospheric flow responds to changes in slab pull and influences both trench retreat and upper plate deformation. We find that trench retreat depends on the competition between shallow and opposite asthenospheric flows below the subducting and upper plates, and that a fast sub-slab flow can hamper trench retreat even when slab buoyancy is high. The mode of upper § Submitted to Earth and Planetary Science Letters on 14 Mar 2019. Revised on 02 Jul 2019.

Figure 2 . 1 :

 21 Figure 2.1: Model initial conditions where quasi-steady-state subduction is already achieved. All physical boundaries are free-slip (no normal flow) except for the freesurface top. The initial velocity field is set to null. The beige contour follows the 1300 K isotherm and outlines the surface plates, ridges and slab. The weak decoupling layer is shown in blue, and the 1000x50 km 2 denser block in red. The trench is located at 2940 km from the left side. Subducting and upper plate surface horizontal velocities are retrieved at 1500 km and 4500 km from the left side, respectively.

Figure 2 . 2 :

 22 Figure 2.2: Viscosity plots overlain with velocity glyphs for the same times during reference simulation (no density anomaly, column a) and simulation DA70 (density anomaly +70 kg/m 3 , column b). The glyphs show clockwise motion in the sub-slab flow, and counter-clockwise flow in the mantle wedge beneath the UP. The glyph lengths indicate the velocity magnitude (same scale for all images). Plate material is outlined in black as the 1300 K isotherm. The denser block is outlined in white for simulation DA70. Both simulations have the same starting trench location and slab tip.

Figure 2 . 3 :

 23 Figure 2.3: Surface horizontal velocities for the reference simulation (black curves, no density anomaly) and simulation DA70 (red curves, density anomaly + 70 kg/m 3 ): subducting plate velocity V sp is positive rightwards, upper plate velocity V up and trench velocity V tr are both positive leftwards. The vertical lines indicate the times labelled t 1 to t 4 inFigure 2.

Figure 2 .

 2 Figure 2.2.b and Figure 2.3 show the effect of disturbing the initial equilibrium by introducing a denser block within the subducting plate (DA70, density anomaly of +70

Figure 2 . 4 :

 24 Figure 2.4: Comparison at t 3 4.8 Myr of i) reference simulation with no density anomaly, and ii) simulation DA70 with a block density anomaly of + 70 kg/m 3 . The plates are outlined with the 1300 K isotherm in black, while the denser block is outlined in white. The panels show plots of a) velocity magnitude, b) second invariant of the strain-rate tensor, and c) horizontal normal stresses xx . Panels (a) and (b) are overlain with the velocity pattern. Magenta lines delimit the upper plate deformed region where strain-rate is higher than 10 -16 s -1 .

Figure 2 .

 2 Figure 2.5 illustrates trench and UP surface velocities for simulations with density anomalies: 0 (reference), +10 kg/m 3 (DA10), +40 kg/m 3 (DA40), and +70 kg/m 3 (DA70).

Figure 2 .

 2 Figure 2.5: a) Upper plate and trench retreat velocities for reference simulation with no density anomaly and simulations DA10, DA40 and DA70, which have density anomalies +10 kg/m 3 , +40 kg/m 3 , +70 kg/m 3 respectively.

Figure 2 . 6 :

 26 Figure 2.6: Plots of i. horizontal normal stress τ xx and ii. Second invariant of the strain-rate tensor overlain with velocity glyphs for simulations: a) DA70, block density anomaly of +70 kg/m 3 , unaltered upper plate thermal structure, with a ridge at its trailing edge (same as Figure 2.4). Shown at t 3 4.8 Myr. b) OFDA70, block density anomaly of +70 kg/m 3 , constant upper plate age of 45 Ma (no ridge, "fixed edge"). Shown at t = 8.8 Myr, the time of strongest horizontal stress in the upper plate. c) YFDA70, block density anomaly of +70 kg/m 3 , constant upper plate age of 15 Ma (no ridge, "fixed edge"). Shown at t = 6.2 Myr, the time of strongest horizontal stress in the upper plate.

  .6.c.ii), and the trench is able to retreat. The trench retreats 52 km by 6.2 Myr (Figure 2.6.c) with a peak V tr of 3.3 cm/yr. This differs from the stationary trench in simulation OFD70, which has an old, strong UP (Figure 2.6.b).

Figure 2 . 7 :

 27 Figure 2.7: Cartoon of the main processes influencing trench retreat and upper plate deformation for wide-slab subduction (similar to our 2-D model set-up). Slab pull variations trigger transient asthenospheric flows. The mode of upper plate deformation depends on the upper plate's ability to translate (far-field boundary condition) and to resist shear basal drag by the mantle-wedge flow. Trench motion is controlled by the competition between the opposing asthenospheric flows on each side of the slab and by the deformation of the upper plate.

Figure 2 . 8 :

 28 Figure 2.8: Cartoon of the two end members of upper plate deformation in subduction zones with retreating trenches. a) Upper plate shortening occurs when the sub-slab flow is faster than the mantle wedge flow, hampering trench retreat. This is applicable to some subduction zones with wide slabs. b) Upper plate extension occurs when the mantle wedge flow is larger than the sub-slab flow, enhancing trench retreat. This is applicable to subduction zones with narrow or torn slabs where toroidal flow can contribute to the mantle wedge flow beneath the upper plate.

1

  Geosciences Montpellier, Universite Montpellier, Montpellier, France 2 Hellenic Centre for Marine Research, Athens, Greece

Figure 3 . 1 :

 31 Figure 3.1: Tectonic map of the Aegean region with block rotation (Morris and

1 .

 1 Seismic reflection dataWe use bathymetry data (from NOAA,[START_REF] Amante | ETOPO1 1 Arc-Minute Global Relief Model: Procedures, Data Sources and Analysis[END_REF] combined with 2-D shallow seismic reflection images acquired in the 1980s, and provided by the Hellenic Centre for Marine Research. These lines were used to propose a recent tectonic interpretation of the Aegean (Sakellariou and Tsampouraki-Kraounaki, 2016). The raw and processed digital data are not available; rather, we have used static images of the profiles with no depth/time axes. This means that we are unable to indicate the depth of the interpretation, but the profiles only image the shallow subsurface. The profiles are oriented NE-SW and NW-SE, and their locations are shown on Figure 3.2.

Figure 3 .

 3 3 shows an example of a half graben with eroded basement in the footwall, present day offset of the sea floor and sediment infill in the hanging wall.

Figure 3 .

 3 S.3).

Figure 3 . 2 :-

 32 Figure 3.2: Interpretation of seismic profile 3

Figure 3 . 3 :

 33 Figure 3.3: Interpretation of seismic profile 14

Figure 3 . 3 .

 33 Figure 3.3. Profile 14 trends NW-SE, and thus shares the strike of the grabens described previously. Therefore, it does not capture horizon orientations or sedimentpackages related to extension. Instead, the horizons appear relatively flat (since the profile is acquired along their strike), but to the north, there is zone of complex deformation with a distinct trough in the sea floor.We interpret this to be a major negative flower structure, which creates a bathymetric trough with an overall depth of 300 -350 m, compared to the surrounding footwall which is ~140 m deep to the NW of the fault zone. The major strike-slip zone trends NNE-SSW across the northern Cyclades, and appears to trend more NE-SW to the east and west of the Cyclades (shown as interpreted, dashed lines on Figure3.4).

Figure 3 . 4 :

 34 Figure 3.4: Interpretation of faulting in the Cyclades mainly based on interpreted seismic profiles and bathymetry. Apparently normal and oblique normal faults are shown in red, while strike slip faults are shown in blue. Solid lines used for observed/correlated fault segments from seismic, dashed lines used for fault segments interpreted from bathymetry alone.

. 6 .Figure 3 . 5 :

 635 Figure 3.5: Tectonic map of Syros showing the main units and interpreted fault sets.

Figure 3 . 6 :

 36 Figure 3.6: Normal faulting Komito

3. 3 . 3 .

 33 Figure 3.8: Oblique normal faulting in a) Galissas and b) Palos

Figure 3 .

 3 Figure 3.8.a. The latter is imaged where the seismic profile approaches the Syros coastline (see Supp. Mat Figure 3.S.4 for the complete profile). This part of the coastline is defined by a major fault, visible in the Google Earth image along with tensile joints shown in Figure 3.8.b. This fault plane creates a major cliff (Figure 3.8.e), which shows tensile mineralised veins, and oblique slickenlines (Figure 3.8.c). It also creates a section of fault breccia and gouge in the order of ~20 m (see Figure 3.8.b.iii for an example of small fault breccia). The stereonet plot on Figure 3.8.d. shows that the fault plane has a dominant NNW-SSE trend (maroon) and shows sinistral normal slickenlines. The tensile mineralised veins, plotted in black, also indicate sinistral normal slip. Therefore, this NNW trending fault shows sinistral normal (oblique) slip, similar to the sinistral NNW trending fault observed in Fabrika (Figure 3.7). We therefore distinguish this older oblique NNW fault set from the pure normal NW fault set, suggesting there are at least two generations of normal faulting.

Figure 3 .

 3 9.b shows an exposure of this fault set in eastern Syros near Ermoupouli, where en echelon shear fractures also indicate dextral motion. A single major fault plane is not found; instead, we find small segments of faults and associated fractures, which we plot collectively as associated structures on Figure3.9.b. These have a range of orientations, but are on average NE-SW. We also find a 10m thick NE-SW trending breccia layer within the fault zone in Central Syros, NE of Galissas Bay, which allows us to extend the fault zone drawn on Figure3.5 to Galissas.

Figure 3 . 9 :

 39 Figure 3.9: Dextral strike slip faulting across central Syros

1 .

 1 Fault pattern and kinematicsFrom our field and seismic data, we have identified three main sets of faults, which we show on Figure3.10: NW-SE normal faults (in red), NNW-SSE oblique normal faults (in burgundy), and NNE-SSW to NE-SW strike slip faults (blue).

Figure 3 .

 3 Figure 3.10: a) Interpreted faulting in the Cyclades and stress regime. b) restored block rotation to Mid-Miocene and predicted stress regime

Figure 3 . 11 :

 311 Figure 3.11: Interpretation of fault evolution by trench retreat, extrusion, block rotation and slab tearing.

Figure 3 .

 3 11 by NE-SW trending strike slip faults in blue), while maintaining their orientation since the extrusion stress field does not change. Blocks between dextral strike slip zones rotate clockwise, which is consistent with observed rotation. This could allow new generations of rollbackrelated normal faults to progressively nucleate in the direction of maximum extension, while older generations experience increasingly oblique slip. Oblique strike slip faults may also open during block rotation between two NE-SW strike slip zones, such as our interpreted negative flower structure shown in Figure 3.3, which likely initiated as a sinistral structure, as shown in Figure 3.10.b.However, the older stress regime predicted from the NNE-SSW fault set(Figure 3.10.b)

Figure 3

 3 Figure 3.S.1

4 .

 4 Can tectonic modelling improve modelling geothermal potential? An example from the Aegean. ""

Figure 4 . 1 :

 41 Figure 4.1: European heat flow map, as it appears in Cloetingh et al., 2010.

  This study uses TNO's in-house thermal modeller, Basin 3D, to create a steady state thermal model of the Aegean as well as a transient model which incorporates lithosphere stretching over 40 Ma. The thermal modeller was documented in Limberger et al. (2018): the first model is preprocessed using a 1D steady state thermal solution.

  across the entire modelled region. This updated version of the code was provided by Jan Diederik Van Wees, TNO in early 2019. The code is otherwise similar to that of Limberger et al., 2018, although their full workflow is not followed here. In this study we only create prior and forward models with different input data, as there are no calibration data available for creating posterior models.

Figure 4 . 3 :

 43 Figure 4.3: Map view of model locatoin, model units, and location of cross section A-A'.

Figure 4 . 4 :

 44 Figure 4.4: Grids used for the Moho (after Tirel et al., 2004) and for thre LAB at 40 Ma (interpreted from literature, see main text).

  . We thus use two time-steps where the stretching factors change at 15 Ma,

Figure 4 . 5 :

 45 Figure 4.5: Preliminary model results: a) Steady state model using the Moho as the bottom boundary condition. b) Steady state model using the 40 Ma LAB as the bottom boundary condition. c) Transient stretched model starting from 40 Ma.

Figure 5

 5

Figure 5 . 1 :

 51 Figure 5.1: Shallow asthenospheric flow beneath the plates controls both trench motion and upper plate deformation.

  where there are no seasonal variations. Geothermal energy here is usually produced using a Ground Source Heat Pump (GSHP). This allows both heating and cooling and accounts for 56% of total geothermal energy production.

	Geothermal System (EGS), introduced below. Manzella (2019) provides a concise
	overview of conventional geothermal technologies, summarised here:
	Conventional shallow geothermal
	Apart from direct thermal application, which accounts for 26% of geothermal energy

There are several conventional geothermal technologies which utilise both shallow and deep geothermal reservoirs. A frontier technology is also under development, that is the Enhanced use, geothermal resources are primarily classified by their depth: shallow and deep geothermal. Shallow geothermal exploits the thermally stable layer a few 10 -100 meters deep,

Table 4 .

 4 1: Lithological properties, after[START_REF] Hantschel | Fundamentals of Basin and Petroleum Systems Modeling[END_REF] Kaurauf, 2009, and Limberger et al., 2018.

	IRGranite150	Igneous_Rocks	Granite_150	2.6	1.15	800	2650	3.32	0	0
			_My_old
	MRGneiss	Metamorphic_Rocks	Gneiss	2.7	1.4	800	2740	2.5	0	0
	MRSchist	Metamorphic_Rocks	Schist	2.9	1.35	920	2740	1.44	0	0
	MREclogite	Metamorphic_Rocks	Eclogite	3.55	1.1	750	3400	0.12	0	0
	MRMarble	Metamorphic_Rocks	Marble	2.8	1.02	860	2700	0.34	0	0
	CRDolomiteTy	Carbonate_Rocks_D	typical	4.2	1.06	860	2790	0.29	0.35	0.39
	p	olomite
	MRSerpentinit	Metamorphic_Rocks	Serpentinite	2.6	1.45	785	2900	0.01	0	0
	e
	IRUltramafics	Igneous_Rocks	Ultramafics	3.8	1	900	3310	0	0	0
	IRDiabase	Igneous_Rocks	Diabase	2.6	1	800	2800	0.18	0	0
	CSSandTyp	Clastic_Sediments_S	typical	3.95	1.15	855	2720	0.7	0.41	0.31
		andstone
	CRMarl	Carbonate_Rocks_M	Marl	2	1.45	850	2700	1.18	0.5	0.5
		arl
	CSGypsum	Chemical_Sediments Gypsum	1.5	1.15	1100 2320	0.05	0	0
	CSSandCrich	Clastic_Sediments_S	clay_rich	3.35	1.2	860	2760	1.1	0.4	0.32
		andstone
	CSConglomera	Clastic_Sediments_O	Conglomerat	2.3	1.05	820	2700	0.85	0.3	0.3
	teTyp	ther	e_typical
	CRLimeOoidgr	Carbonate_Rocks_Li	ooid_grainst	3	1.19	835	2740	0.35	0.35	0.01
	ainstone	mestone	one
	CSSalt	Chemical_Sediments Salt	6.5	1.01	860	2740	0.02	0	0
	CRLimeOrgRT	Carbonate_Rocks_Li	org._rich_ty	2	1.95	845	2680	1.4	0.51	0.52
	yp	mestone	pical
	Bouguer	Other_Minerals	Bouguer	1000	1.17	885	2670	0	0	0
	KmLM	IRPeridotite	none	0	0	0	3200	0.02	0	0

Acknowledgements

This study has benefited from several discussions with Jeroen van-Hunen, Serge Lallemand and the other members of the SUBITOP ITN. The paper was largely improved by constructive feedback from Laurent Husson and an anonymous reviewer.

Our work has received funding from the European Union's EU Framework Programme for Research and Innovation Horizon 2020 under Grant Agreement No 674899.

In the results section below, we first present a description of our reference simulation, followed by the results for simulations where we systematically add complexities. These are: (i) introducing a density anomaly in the slab; (ii) varying the value of the density anomaly; (iii) changing the upper plate boundary condition by fixing its edge; and (iv) changing the strength of the fixed upper plate. We then discuss our model results and their implications. take a conservative sub-slab channel thickness of 50 km, which lengthens 600 km over 6 Myr, this yields a volume increase rate of 0.2x10 -3 m 2 /s. Thus, slab flattening and subslab channel flow compensate most of the sub-slab mantle volume reduction caused by rollback, with the free surface elevation above the subducting plate accommodating the rest (Appendix G).

At t 3 4.8 Myr during trench deceleration, dynamic pressure in most of the sub-slab compartment is equivalent to the dynamic pressure in the right compartment near the slab tip (Appendix H). This shows that the sub-slab mantle is not over-pressured and suggests that slab rollback is not limited by the leftward side boundary condition.).

We also tested the subduction dynamics in a set-up with a lower mantle but is otherwise similar to DA70 (Appendix 8). The asthenosphere here flows from the left to the right of the domain throughout the simulation. Since the lower mantle in this test model has a significantly higher viscosity, the velocity magnitudes are all lower than in DA70 (e.g. peak V tr of 2 cm/yr compared to 4 cm/yr in DA70). Despite the different velocity magnitude and slab geometry (the slab anchors in the lower mantle), the trench retreat pattern remains similar to DA70 (trench retreat acceleration and deceleration while slab pull is still relatively high). This shows that V tr deceleration in our upper-mantle simulations is not caused by a closed bottom boundary or the inability of the mantle to flow between the two compartments on each side of the slab.

In our simulations, the slab and asthenosphere cell geometries are forced by the free slip side and bottom boundary conditions. Enns et al. (2005) have shown that the type of side boundary conditions affect the trench velocities, with free-slip (reflective) sides inducing slower rollback velocities than periodic side boundary conditions. Although this suggests that absolute velocity variations inferred from our simulations should be treated with caution, our qualitative interpretation of associated upper plate deformation will likely remain valid.

Supplementary material

Appendix A: Slab sinking rates DA70 surface subducting plate velocity (V sp ) together with the slab sinking velocities at a depth of 250 km: horizontal component (V x ), the vertical component (V y ), and the bulk velocity magnitude (V mag ). The slab sinking velocities exhibit the same temporal pattern as the V sp , and differ from trench retreat velocity (V tr in Figure 2.3). This shows that trench velocity is not solely controlled by slab pull, whereas the motion of the subducting plate is.

Appendix B: Graphical plots of density anomalies and upper plate deformation from Table 2 i) Trench retreat velocity against density anomaly. ii) Maximum width of upper plate deformation against density anomaly. iii) Maximum discrepancy between V up and V tr against maximum width of upper plate deformation. Cyclades, from which we interpret the direction of maximum extension (or sigma 3) to be around NE-SW (e.g. normal to the pure normal faults). The NNW-SSE sinistralnormal fault set and the NE-SW dextral strike slip faults are thus compatible with this stress field. Moreover, this stress field is close to the predicted direction of stretching from geodetic data (e.g. [START_REF] Allmendinger | Strain and Rotation Rate from GPS in Tibet , Anatolia , and the Altiplano[END_REF] show maximum stretching in the central Aegean to be ~NNE-SSW). Since new extensional structures most likely open in the direction of maximum stretching, the presence of oblique extensional faults (NNW-SSE), oriented at a high angle to the direction of stretching, suggests that they act as inherited structures in the present day stress field, and consistent with inferred preexisting structures (e.g. Fabrika sinistral pull-apart, Figure 3.6). These structures require an extension direction closer to WNW-ESE, which differs from today's stress field.

Therefore, either the faults have rotated, the stress field has rotated or both have rotated.

Normal fault rotation since the Mid-Miocene

Palaeomagnetic studies show that the north-western Cycladic area has rotated 22°-23° clockwise since the Miocene (Morris and Andersen, 1996;Avigad et al., 1998). We restore this block rotation in Figure 3 This implies that in addition to block rotation, the stress field must have changed since the Mid-Miocene. We propose this occurred due to the development of the western Hellenic slab tear.

We illustrate this interpretation in Figure 3.11: Mineral Exploration (IGME) of Greece in 1970. Usage, however, is currently limited to direct applications (e.g. for spas or agricultural) from low temperature reservoirs [START_REF] Mendrinos | Exploring for Geothermal Resources in Greece[END_REF]. This is because poorly conducted operations on the Milos Island power plant in the 1970s-1980s led to environmental pollution, and thus, social

Supplementary material

The crustal model consists of four main layers. A map view of the units is shown in The lithological compositions of each of these layers was compiled by Kristof Porkalob (modified from Larede, 2018), and are shown in Table 2. These lithological compositions are used for calculating lithological properties (Table 1) and subsequently thermal conductivity. We found that lithological composition has a first order impact on the thermal model, so we expect a model with increased resolution to give more precise values.

The Moho, the boundary between the crust and the lithospheric mantle, is based on the modelled Moho in Tirel et al., 2004, which The lithosphere-asthenosphere boundary (LAB) is used as the bottom thermal boundary condition. Here, it is considered to be the boundary between mainly conductive and convective mantle, and is represented by the 1200°C isotherm, after [START_REF] Tesauro | A New Thermal and Rheological Model of the European Lithosphere[END_REF].

Since the model is forward modelled, the boundary condition is set for the start of the simulation. Thus, we have created two LAB grids, one for the present day to create the steady state model, and a second for 40 Ma to create the transient, stretched model. This allows us to capture the migrating subduction zone, as the present day offshore Aegean domain is thinned, while at 40 Ma, it would have still been a thickened nappe stack. We interpret the LAB structures from literature reconsructions (e.g. van Hinsbergen and Schmid, 2012;Jolivet and Brun, 2010;[START_REF] Menant | Kinematic Reconstructions and Magmatic Evolution Illuminating Crustal and Mantle Dynamics of the Eastern Mediterranean Region since the Late Cretaceous[END_REF] The relevant results include present day modelled 3D temperature based on the bottom boundary conditions, lithology compositions and geometry, and stretching (for the transient models). In the absence of wellbore data, model validation is limited to comparing modelled and observed surface temperatures.

Preliminary results

We start by a crust-only steady state simulation, using the Moho as the bottom boundary. We then change the bottom boundary condition to the 40 Ma LAB. We then introduce incremental changes in each following simulation. We thus compare the following three models:

Discussion and conclusion

The preliminary models above show that introducing tectonic stretching into thermal model largely affects predicted results. Unfortunately, these results cannot be validated as there are no wellbore data available. These would typically act as additional tie points and could enable further stochastic modelling (as in Limberger et al., 2018). But even without wellbore modelling, there is ample scope to incorporate more tectonic Additionally, the model predictions could be largely improved by incorporating the magmatic arc and its migration through time, an important aspect which is not included in the workflow above. Introducing magmatism is particularly interesting for identifying local areas with elevated temperatures. As the Aegean crust is made up of a nappe stack composed of lithologies with different thermal properties, it would be interesting to see how the migrating arc influences the temperature distribution in the crust.

The preliminary work above shows that there is large potential to improve geothermal model prediction by incorporating large scale tectonic processes. Tectonic modelling could aid thermal prediction in a myriad other ways, particularly when considering different scales. This has already been shown by Cloetingh et al., 2010 who emphasised that tectonic studies can aid in both predicting and producing geothermal energy (refer to chapter 1.4). This is particularly relevant for producing geothermal energy in non-volcanic areas, which is not yet common due to the difficulties in predicting optimum well locations in such areas. Reducing exploration and production risk in non-volcanic areas could largely increase the global potential for geothermal energy, and take it beyond its current status as a resource restricted to volcanic areas.

Acknowledgements

Research presented in this chapter was conducted as part of an industrial secondment